第8讲(2014-11-19).doc

第4章 数值积分与数值微分上次课要点：4 龙贝格求积公式固定复化求积公式中的步长会导致：（1）步长取得太大，精度难以保证；（2）步长太小，则计算量太大，有时不必要。解决思路：采用变步长方法。一、梯形法的递推化将区间进行等分，设节点为，则梯形法求积公式为再将每一个子区间二等分，记相应的中点为，则此时梯形法求积公式为于是有递推公式。二、龙贝格算法的思路上述算法过于依赖梯形公式，效率较低。问题：能否利用和的某种组合，得到精度更高的求积公式？理查德外推法：用精度较低的近似公式，组合成精度较高的近似公式。龙贝格将它应用到数值积分上。定理4 梯形法二分后的两个积分值的线性组合等于辛普森法积分值，即.用同样方法，可得三、龙贝格算法步1 初始化：计算，置（记录区间的二分次数），；步2（二分） 计算积分值：；步3（加速） 求加速值：计算，；步4 精度检验：对指定的精度，若，则终止计算，并取作为所求的结果；否则置，转步2。5 自适应积分方法问题：如果被积函数在积分区间上，有的部分变化剧烈，有的部分变化平缓，则将整个积分区间一起细分是非理性的。解决方法：根据各个子区间满足精度的情况，采用不同的步长这就是自适应积分方法。8 数值微分研究内容：用数值方法求函数在离散点处的导数值。解决方法：用某些点处函数值的线性组合，近似指定点处的导数值。问题分为两类：（1）求一点处的导数值；（2）求等距分布的多个点处的导数值。一、泰勒展开法1、显式求导公式（用于求一个点处的一阶、二阶导数）（1）一阶向前差商公式（用于求区间左端点的近似导数）：（2）一阶向后差商公式（用于求区间右端点的近似导数）：（3）一阶中心差商公式（用于求区间内点的导数值）：（4）二阶中心差商公式（用于求区间内点的二阶导数值）：注：不能太小，也不能太大。2、隐式求导公式节点：函数值：边界条件：和求：。（1）用表示的近似值，得这是由个未知数、个方程组成的线性方程组。上述方程组主对角占优，解存在且唯一。用追赶法求解是数值稳定的。本次课继续。再看二阶导数。由在处的泰勒展开式，可得而，于是用表示的近似值，得近似等式即再加上两个边值条件，并令便得（大约）具有精度的二阶导数隐格式方程组上述方程组同样是主对角占优的，解存在且唯一，用追赶法求解数值稳定。二、插值型求导公式过已知点，建立插值多项式，用的导数作为的导数的近似值，即。（一）用Lagrange插值多项式求近似导数常用的有三点求导公式和五点求导公式。1、三点公式已知三个点，则的导函数因此，各点处的一阶导数近似公式：2、五点公式已知五个点用三点类似的方法可得（二）用三次样条插值求导（求多个点处的导数值）（略）三、利用数值积分求导设，则上式右边的积分可以采用不同的求积公式，从而得到不同的求导公式。1、中矩形法（用于求一点处的导数）对用中矩形公式则得即。2、辛普森法（用于求多点处的导数）对用辛普森公式则得用表示的近似值，得这是由个未知数、个方程组成的线性方程组。若已知，则可得到前面相同的一阶导数隐格式。四、数值微分的外推算法外推法是一种逐步提高公式精度的方法，使截断误差逐步减小。由泰勒展开式 （1） （2）得 （3）其中只与有关，与无关。令，则（4） （5）由得（6）可见，用精度更高，余项为。令，可以继续外推。一般地，令，可得，误差：。当较大时，截断误差很小。考虑到舍入误差，一般不能太大。五、复数法求近似导数William Squire , George Trapp. Using Complex Variables to Estimate Derivatives of Real Functions. SIAM Review, Vol.40, No.1, PP110-112, March 1998.算法要点：将按复数展开，表示虚数单位。则虚部令 。该式计算的截断误差为。没有两个相近数相减产生的有效数字位数的损失。上述算法与中心差商法进行数值比较 （1） （2）例1 ，公式（1）公式（2）0.1D-010.1D-020.1D-030.1D-040.1D-050.1D-060.1D-070.1D-080.1D-090.1D-0100.1D-0110.1D-0120.1D-0130.1D-0140.1D-0150.18602018344501897D+020.18600824790342307D+020.18600812854818738D+020.18600812735480865 D+020.18600812734248517 D+020.18600812735081185 D+020.18600812723423843 D+020.18600812723423843 D+020.18600814222224926 D+020.18600815332447951 D+020.18600898599174798 D+020.18601231666082185 D+020.18585133432225120 D+020.18596235662471372 D+020.20539125955565396 D+020.18599607128036329D+020.18600800678177631 D+020.18600812613698936 D+020.18600812733054151 D+020.18600812734247702 D+020.18600812734259637 D+020.18600812734259757 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+020.18600812734259759 D+02可以看出，公式（1）在时比较准确。例2 ，公式（1）公式（2）0.100000E-10.100000E-20.100000E-30.100000E-40.100000E-50.100000E-60.362298E+010.362229 E+010.362158 E+010.360012 E+010.357628 E+010.476837 E+010.362109 E+010.362202 E+010.362203 E+010.362203 E+010.362203 E+010.362203 E+01可以看出公式（1）在时比较准确。埃尔米特求积公式的实现中，两个端点的导数用计算即可。第5章 解线性方程组的直接法1 引言本章研究求解线性方程组 （1）的数值方法。这里， 并且是非奇异矩阵，是非零向量。线性方程组的数值解法分类：（1）直接法适用于中小规模的方程组（2）迭代法适用于大型稀疏矩阵方程组。直接法：如果计算过程中没有舍入误差，经过有限步运算，可得到问题的精确解。只有代数问题才可能构造直接法。求解线性方程组的迭代法：用构造近似解序列，使在一定条件下，为方程组的解。2 高斯直接消去法记方程组（1）的增广矩阵： （2）消去法基本思想：（一）消去过程：（二）回代过程：求出上三角矩阵对应的方程组的解。 消去过程第一步消去：在增广矩阵中，假设，以为主元素将增广矩阵的第1列元素，除外全部化为零。记 ，用乘矩阵（2）的第1行，再分别与第行相加（），得到 (3)其中；。第2步消去：再假设，以为主元素将第2列的元素，除最上面的两个元素外全部化为零。记，用乘矩阵(3)的第2行，再与第行相加（），得到 (4)其中；。如此进行下去一般地，第步消去：假设经过步消去后，矩阵中，以为主元素将第列除上面的个元素外全部化为零。记，用乘的第行，再与第行相加，得到(5)其中；。经过步消去后，增广矩阵（2）变为 (6)其中系数矩阵是一个上三角形矩阵。回代过程增广矩阵(6)对应的上三角形方程组 (7)与原方程组(1)是等价的。方程组(7)的解可按下述公式递推得到： 问题：如何保证？答案：系数矩阵的顺序主子式都不等于零。由上述讨论过程可得定理1 若方程组的系数矩阵的顺序主子式不等于零，则可通过高斯消去法求出方程组的唯一解。高斯消去法两个过程的计算量:除法运算次数：乘法运算次数：总计算量：约3 高斯列主元素消去法对于方程组（1），只要系数矩阵是非奇异的，方程组的解就存在且唯一。问题：（一）系数矩阵非奇异，不能保证其顺序主子式都不等于零，从而不能保证高斯消去法中的。（二）即使，但分母与分子（）相比是较小的数，则可能将消去公式中的误差放大倍，导致数值不稳定。 解决措施：（一）全主元高斯消去法：在以为左上角元素的子矩阵中取绝对值最大的元素（当系数矩阵非奇异时，该元素必定不等于零），通过行、列交换到的位置作为主元素。缺点：耗时多。（二）列主元消去法：假设已经进行了次消去过程，得到在增广矩阵的第列，从及其下面的元素中选取绝对值最大的元素作为主元素然后交换第行和第行，并实施消去过程，得到增广矩阵。高斯消去法：要求系数矩阵的顺序主子式不等于零；列主元消去法：只要求方程组的系数矩阵非奇异，并较好地避免了误差的传播。记算法（列主元素高斯消去法）：消元过程：对，执行（1）选主元：；（2）判别：若，停机（为奇异矩阵）；（3）换行：若，则；（4）消元：对，执行（5）判断：如果，停机（为奇异矩阵）。回代过程：对，执行输出方程组的解：例：用高斯主元素消去法解方程组采取8位十进制浮点运算。解：回代求解，得该方程组具有5位有效数字的解（视为精确解）可见，列主元素消去法得到的解精度较好。4 矩阵的三角分解法本节仅考虑与矩阵的直接三角分解法；进一步，还可以构造选主元素的矩阵三角分解法。高斯直接消去法：针对系数矩阵的顺序主子式全不等于零的情况，分两个步骤:消元过程:将方程组对应的增广矩阵化成上三角矩阵;回代过程:通过求解上三角矩阵对应的上三角方程组得到原方程组的解。矩阵的三角分解法与高斯直接消去法等价：第一步： 其中为单位下三角阵， 为上三角阵。第二步：解方程组 和.矩阵LU分解的作用：1）用于求解系列线性方程组。这些方程组的系数矩阵相同，只是右端向量不同（例如反幂法求矩阵按模最小的特征值的情况）2）用于求解特殊类型的方程组。系数矩阵为对称阵；三对角阵一、矩阵的直接三角分解法假设方程组（1）的系数矩阵的顺序主子式都不为零。（一）一般矩阵的三角分解法用将的第1列除外的其余元素化为零的初等变换对应的矩阵为，其中，相应的初等变换可表示为。一般地，用（当系数矩阵的顺序主子式不等于零时，必不等于零）将第列的第个到第个元素化为零的初等变换对应的矩阵为，其中，相应的初等变换为。重复这个过程，得，。从而。注意到可得为单位下三角矩阵。定理2 （矩阵的LU分解）设为阶矩阵，如果的顺序主子式（），则可分解为一个单位下三角矩阵和一个上三角矩阵的乘积，即 且这种分解是唯一的。 上述三角分解叫做Doolittle分解。若为单位上三角矩阵，为下三角矩阵，称为Crout分解。下面给出三角分解法解线性方程组的步骤。首先回顾高斯消去法步骤：消去过程第1步消去：保留增广矩阵的第1行，假设，记 ，令；。得到 (3)据此，可得到的第一行、的第一列。第2步消去：假设，记，令；。得到 (4)据此，可得到的第二行、的第二列。一般地，第步消去：假设矩阵中，记，令；。得到(5)据此，可得到的第行、的第列。经过步消去后，增广矩阵（2）变为 (6)其中系数矩阵是一个上三角形矩阵。回代过程方程组的解可按下述公式递推得到：设非奇异矩阵的顺序主子式都不等于零，将分解为其中，由高斯消去法的实施过程可知，和的元素可通过个循环求得，方法如下：（1）计算的第1行和的第1列元素：，；（2）对，计算的第行，的第列元素：，。进行了上述三角分解后，求解方程组的问题等价于求解两个三角形方程组，先通过第一个方程组求出，然后代入第二个方程组求出。（二） 对称矩阵的三角分解法当矩阵为对称矩阵且顺序主子式不等于零时，将进一步分解为即。又为对称阵，故，再由分解的唯一性，。因此，对于对称矩阵，分解式可表示为。求解方程组，等价于求解三个简单的方程组。特别，当矩阵为对称正定阵时，的对角线元素，令则有，其中为下三角阵。从而有下述定理：定理（对称正定矩阵的三角分解定理）设为对称正定阵，则存在一个非奇异下三角阵，使当限定的对角线元素为正时，这种分解是唯一的。矩阵的上述分解方式叫做Cholesky分解，的元素可按下述方式产生：对于，用Cholesky分解法求解方程组，即等价于求解步骤如下：（1）解方程组，公式为；（2）解方程组，公式为。对于系数矩阵为对称正定矩阵的线性方程组，首先对系数矩阵进行Cholesky分解，然后通过解两个方程组进行求解，叫做解方程组的平方根法，其计算量大约是高斯消去法的一半。（三） 三对角矩阵的