计算方法5-1

第五章 数值积分法,由高等数学的有关知识，我们知道，要求得 ,只要求出 的一个原函数 ,然后利用牛顿-莱布尼兹公式：,即可求得 在区间 上的积分值。,但是，在工程技术领域，在实际使用上述求积分方法时，往往会遇到下面情况：,1. 函数f (x)没有具体的解析表达式，只有一些由实验测试数据形成的表格或图形。,2. f (x)的原函数无法用初等函数表示出来，如：,3. f (x) 的结构复杂，求原函数困难，即不定积分难求。,由于以上种种原因，因此有必要研究积分的数值计算方法，进而建立起上机计算定积分的算法，此外，数值积分也是研究微分方程和积分方程的数值解法的基础。,同样，对函数f (x)求导，也有类似的问题，需要研究数值微分方法。,5.1 数值积分的基本概念,1 构造数值求积公式的基本思想,定积分I=ab f (x)dx在几何上为x=a, x=b, y=0和y=f (x)所围成的曲边梯形的面积。定积分计算之所以困难，就在于这个曲边梯形中有一条边y=f (x)是曲边，而不是规则图形。 由积分中值定理，对连续函数f (x)，在区间a, b 内至少存在一点，使：,也就是说，曲边梯形的面积I 恰好等于底为(b-a)、高为f ()的规则图形矩形的面积(右图），f ()为曲边梯形的平均高度, 然而点的具体位置一般是不知道的，,因此难以准确地求出f ()的值。但是，由此可以得到这样的启发，只要能对平均高度f ()提供一种近似算法，便可以相应地得到一种数值求积公式。,如，用两端点的函数值f (a)与f (b)取算术平均值作为平均高度f ()的近似值，这样可导出求积公式：,更一般地，可以在区间a, b 上适当选取某些点xk (k=0,1, , n)，然后用f (xk) 的加权平均值近似地表示f ()，这样得到一般的求积公式：,其中，点xk 称为求积节点，系数Ak 称为求积系数，Ak 仅仅与节点xk 的选取有关，而不依赖于被积函数f (x)的具体形式，即xk决定了，Ak也就相应的决定了。,回顾定积分的定义，积分值I 是和式的极限：,其中xk是a, b 的每一个分割小区间的长度，它与f (x)无关，去掉极限，由此得到近似计算公式：,因此，式（5-1）可作为一般的求积公式，其特点是将积分问题归结为函数值的计算，从而避开了使用牛顿一莱布尼慈公式需要求原函数的困难，适合于函数给出时计算积分，也非常便于设计算法,便于上机计算。 求积公式（5-1）的截断误差为：,Rn也称为积分余项,2 代数精度,数值积分是一种近似方法，但其中有的公式能对较多的函数准确成立，而有的公式只对较少的函数准确成立。为了反映数值积分公式在这方面的差别，引入代数精度的概念。,定义1: 如果某个求积公式对所有次数不大于m的多项式都精确成立，而至少对一个m +1次多项式不精确成立，则称该公式具有m次代数精度。,一般来说，代数精度越高，求积公式越好。为了便于应用，由定义1容易得到下面定理。,定理1: 一个求积公式具有m次代数精度的充分必要条件是该求积公式对 1,x,x2,xm 精确成立，而对xm+1不精确成立。,同理可证明矩形公式的代数精度也是一次的,上述过程表明，可以从代数精度的角度出发来构造求积公式。例如，对于求积公式（5-1），若事先选定一组求积节点xk (k=0, 1, n,)， xk可以选为等距点，也可以选为非等距点，则可令公式对f(x) = 1,x,xn 精确成立，即得：,这是关于A0、A1、An的线性方程组，系数行列式为范德蒙行列式，其值不等于零，故方程组存在唯一的一组解。求解该方程组即可确定求积系数Ak，因此求积公式(5-1)至少具有n次代数精度。,解:求积公式中含有三个待定参数，可假定上述近似式的代数精度为m =2，则当f (x)=1，x，x2时，近似式应准确成立，即有：,代回去可得：,上述近似式不仅对特殊的次数不高于3次的多项式f (x) = 1,x,x2, x3准确成立，而且对任意次数不高于3次的多项式，a0+a1x+a2x2 + a2x3 （f (x)=1,x,x2, x3的线性组合）也准确成立，事实上，令R( f )表式近似式的截断误差：,检查上式对 m = 3 是否成立，为此，令 f(x)=x3 代入上式，此时左边,由于对任意的常数, 和函数f (x)，g (x) 成立:,这表明，误差对f (x)=1, x, x2, x3准确成立，则对它们的任意线性组合a0 + a1x + a2x2+ a3x3也准确成立，所以通常检查一个求积公式是否具有m次代数精度，只需检查对f(x)=1,x,xm 是否准确成立即可。,上述方法称为待定系数法！,上述方法称为待定系数法，在具有尽可能高的代数精度的要求下，利用它可以得出各种求积公式。,说明1：由待定系数法确定的求积公式没有确切的误差估计式，只能从其所具有的代数精度去判定求积公式的准确程度。,说明2：因此，希望由待定系数法确定的求积公式的代数精度越高越好，通常的方法是要确定n +1个待定系数,可设求积公式具有n次代数精度，去建立n +1个方程求解。,设给定一组节点a x0 x1 xn-1xn b，且已知f (x) 在这些节点上的函数值，则可求 得f (x)的拉格朗日插值多项式：,3 插值型求积公式,其中lk(x) 为插值基函数。取f (x) Ln(x)，则有：,记：,则有：,这种求积系数由 所确定的求积公式称为插值型求积公式。,根据插值余项定理，插值型求积公式的求积余项为：,其中a,b 且与x有关。在插值中，因f (x) 不知道，所以无法估计插值误差。而在这里，f (x)作为被积函数，上式却可以用于估计积分的误差。,证：(充分性) 设求积公式（5-1）至少具有n次代数精度，那么，由于插值基函数 li(x) (i=0,1,n)均是次数为n的多项式，故式（5-1）对li(x)精确成立，即:,(必要性) 设求积公式（5-1）是插值型的，则对所有次数不大于n的多项式f (x)，其求积余项Rn = 0，即公式是精确成立的。由定义1知求积公式至少具有n次代数精度。,定理2说明，当求积公式（5-1）选定求积节点xk后，确定求积系数Ak有两条可供选择的途径：求解线性方程组或者计,算积分。由此得到的求积公式都是插值型的，其代数精度均不小于n次。,此例说明三个节点的求积公式不一定具有二次数精度，其原因是此求积公式不是插值型的。,5.2 梯形公式、辛甫生公式、柯特斯公式,5.2.1 梯形公式,几何意义：用梯形AabB 的面积近似代替曲边梯形AabB 的面积，因此，上述数值积分公式又称为梯形公式。,通过这三个数据点，就可以做一个二次插值多项式：,5.2.2 辛甫生公式,假设在积分区间 上，取 ， ， ，这样就得到三个点：,为了便于计算右端的积分，作变量代换：,几何意义：通过曲线上的三点作一抛物线，用以抛物线为曲边的曲边梯形面积近似代替原来的曲边梯形的面积。,所以，公式 就称为辛甫生公式。,5.2.3 牛顿一柯特斯（Newton-Cotes）公式,设将积分区间a, b 划分为n等分，步长h=(b-a)/n，求积分节点取为xk = a+kh(k=0,1,n)，由此构造插值型求积公式，则其求积系数为:,记：,称之为n阶牛顿一柯特斯（Newton-Cotes）公式简记为N-C公式， 称为柯特斯系数。显然， 柯特斯系数与被积函数f (x) 和积分区间a, b 无关，且为多项式积分，其值可以事先求出备用。表5-1中给了了部分柯特斯系数。,柯特斯系数,经计算或查表得到柯特斯系数后，便可以写出对应的牛顿一柯特斯（Newton-Cotes）公式。,当n =1时，按公式有：,即为梯形公式,当n=4时，所得的公式称作柯特斯公式，它有五个节点，其系数：,所以柯特斯公式是：,柯特斯系数的性质：,1、与积分区间无关：当n 确定后，其系数和都等于1，即：,2、对称性：,解：由定理2,辛卜生公式至少具有二次代数精度,因此只需检查对f (x)=x3成立否。当f (x)=x3时：,所以I = S，表明辛卜生公式对于次数不超过三次的多项式准确成立，用同样的方法可以验证对于f (x)=x4，辛卜生公式不成立，因此辛卜生公式的代数精度可以达到三次。,例5:分别用梯型公式、辛卜生公式和柯特斯公式计算积分：,解：由梯形公式得：,由辛卜生公式得：,由柯特斯公式得：,积分的精确值：,与之相比可以看到，柯特斯公式的结果最好，具有七位有效数字；辛卜生公式的结果次之，具有四位有效数字；而梯形公式的结果最差，只有两位有效数字。,5.3 复化求积公式,在实验计算中常用的梯型公式、辛卜生公式和柯特斯公式这些低阶的N-C公式作为数值公式，但若积分区间比较大，直接使用这些求积公式，则精度难以保证； 若增加节点，就要使用高阶的N-C公式，当n 8时，N-C公式的收敛性和稳定性得不到保证，因此不能采用高阶的公式，事实上，增加节点，从插值的角度出发，必然会提高插值多项式的次数，Runge现象表明，一般不采用高次插值，亦即不用高阶N-C公式。 因此，为提高精度，当增加求积节点时，考虑对被积函数用分段低次多项式近似，由此导出复化求积公式。,5.3.1 复化梯形公式,用分段线性插值函数近似被积函数，等于把积分区间分成若干小区间，在每个小区间上以梯形面积近似曲边梯形面积，即用梯形公式求小区间上积分的近似值。如下图所示，这样求得的近似值显然比用梯形公式计算高。,定积分存在定理表明，只要被积函数连续，当小区间长度趋于零时，小梯形面积之和趋于曲边梯形面积的准确值，即定积分的准确值。,它实际上就是用定积分定义计算积分，经等分区间，在每个小区间上以直线近似替代曲顶（线）然后求和，略掉无限细分区间（求极限）这一步而得到的近似值。,这就是复化梯形公式。,由于在每个小区间 上， 是用一次多项式去近似的，由插入法可知：,其中,记,为每个小区间 上用小梯形面积近似代替小曲边梯形面积，所产生的截断误差。,由积分中值定理可知，如果 在区间 上连续，则在积分区间 内存在一点 ，使得 成立。,所以,所以，在区间 上总的截断误差为：,又因为 在区间 上连续，,所以，在区间 内必存在一点 ，使得：,这就是复化梯形公式的截断误差公式。,假设,则：,5.3.2 复化Simpson公式,如果用分段二次插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就导出复化Simpson公式。,若 在区间 上有四阶连续导数，那么辛甫生公式的截断误差为：,因为f (4)(x) 连续，故存在(a, b)，使得：,上式表明，步长h越小，截断误差越小。与复化梯形公式的分析相类似，可以证明，当n 时，用复化Simpson公式所求得的近似值收敛于积分值，而且算法具有数值稳定性。,将区间 分成4等份,节点 ，这样就得到五个数据点：,通过这五个数据点，就可以做一个四次插值多项式 ，那么类似于辛甫生公式的推导，就可以得到：,5.3.3 复化柯特斯公式,为了提高求积精度将区间a, b分成 n 等分，分点为：,在每个子区间上应用柯特斯公式，然后相加，便得到复化的柯特斯公式：,柯特斯公式的截断误差：,梯形公式：,具有1次代数精度。,辛甫生公式：,具有3次代数精度。,柯特斯公式：,具有5次代数精度。,当然，也可以由三个公式的截断误差，推出三个公式的所具有的代数精度。,辛甫生公式的插值节点只比梯形公式多一个，但其代数精度却比梯形公式高2，它们都是最为常用的数值积分公式，尤其是辛甫生公式，逻辑结构简单，且精度又比较高。,例 6: 根据函数表,解：（1）由复化梯形公式,n=8, h=1/8：,（2）由复化Simpson 公式,n=4,h=1/4：,与准确值I=0.9460831比较，显然用复化Simpson公式计算精度较高。,事实上，由误差公式可知：RT (f )=O(h2), RS (f )=O(h4)，故当h比较小时，用复化Simpson公式计算误差较小。,由误差估计公式不仅可以计算所求近似值的误差，反之，亦可由给定的精度估计应取多大步长。,5.4逐次分半算法(变步长方法),当实际应用计算机进行数值积分计算时，一般先提出一定的精度要求，然后根据误差估计式对于预先给定的精度给出步长h或n, 但由于误差估计式中要估计高阶导数，而这一点往往很困难。 因此实际计算时，常采用变步长方法：逐步缩小步长，每次将步长缩小一半，或者说逐次等分区间,反复利用复化求积公式，直到相邻两次计算结果相差不大为止或者满足给定精度为止。,5.4.1 梯形法的递推公式,因此计算梯形序列T2m可按：,4. 设将区间a , b n 等分，共有n+1个分点，如果将积分区间再等分一次，则分点增为2n+1个，将等分前后两个积分值联系起来加以考察：,将每个子区间上的积分值相加得：,此为复化梯形公式的递推公式,按上述逐次分半算法，并利用递推公式，T2m 的计算较容易，那么，上述算法何时停止？,以此作为停止计算的控制。, f (m-1) 与 f (m) 是二阶导数 f (x) 在a, b上2m-1个点与2m个点的算术平均数（每个小区间上取一个点）， 若f (x) 在a, b 的二阶导数连续，则当m较大时：,5.4.2 Simpson公式的逐次分半法,上例说明Tn收敛慢，求T128 要计算64个新增的函数值，而将T8与T4重新组合可构造S8。,在复化梯形公式逐次分半算法中：,我们将此误差估计加到T2m上构成新的近似值：,而在Simpson 逐次分半算法中：,即由Simpson序列可构造出收敛更快的Cotes序列 。,具体做法都是利用控制结束的误差式，构成新的，收敛更快的序列，而由前面的推导可知，下面这些公式具有如下规律性：,类似地，也可以推导出：,5.5 龙贝格（Romberg）求积公式,问题：梯形公式的优点是算法简单，但要达到同样的精度要求，积分区间等分为小区间的个数N 要很大，这样一来，收敛速度就会放慢，能不能有一种算法，在利用梯形公式算法简单的基础上，提高收敛速度？,我们在利用梯形公式 时，一般 和 是 是分别计算的，如果能找到一种方法，使得 ，这样，计算 时，前面已算过的 就可以利用，当然能加快收敛速度。,从面复化梯形公式的递推式可以看出，计算 时， 仍然有用，只需计算新增加的分点处的函数值，从而节省了计算量。,利用梯形公式与辛甫生公式的关系：,它恰好