计算方法-数值积分

1、第5章 数值积分,若函数f(x)在区间a, b上连续且其原函数为F(x) , 则可用牛顿莱布尼兹公式，来求定积分,51,求定积分,复习,函数关系由表格或图形表示,无法求出原函数,定积分计算可能遭遇的三种情况,被积函数的原函数不是初等函数,被积函数f(x)没有具体的解析表达式,被积函数f(x)的原函数F(x)不易找到,第5章 数值积分,从几何上看定积分,定积分是曲边梯形的面积,图 5.1,左矩形,右矩形,52,53,图 5.2 梯形面积,图5.3 抛物求积,54,55,第5章 数值积分,近似值,5.1,5.2,5.4,牛顿 柯特斯 (NewtonCotes) 公式,复合求积公式,龙贝格(Romb

2、erg) 积分方法,5.1 牛顿 柯特斯(NewtonCotes) 公式 建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数(x),用(x)代替被积函数f(x),于是有 现用第四章介绍的插值多项式Pn(x)来代替被积函数f(x),即有,将积分区间a,b n等分，则节点是等距分布的，节点x0 ,x1 ,x2 , xn可表示成xk=x0+kh (k=0,1,n)，其中 x0=a, xn=b,称为步长,Newton-Cotes公式,若Ln (x)为Lagrange插值多项式，则由公式,于是,令,5.5,公式(5.6)称为等距节点内插求积公式,则有,5.6,求Ak,在等距节点前提下，做变换

3、,由,可得,而x-xj=(t-j)h (j=0,1,2,n) ，xk-xj=(k-j)h (j,k=0,1,2,n且jk)。于是(5.5)式即为,记,则,5.9,称为牛顿-柯特斯公式。其中Ck(n) 叫Cotes系数，Cotes系数与被积函数及积分区间无关,计算柯特斯系数,n=1时，有两个Cotes系数,n=2时，有三个Cotes系数,类似可得，n=3时有四个Cotes系数,n=4时，有五个Cotes系数,几个常用的牛顿-柯特斯公式,n=1时,此即(5.3)式，为梯形公式,其中,称为Simpson公式,其中 c,d,e为a,b的四等分点，称为Cotes公式,n=2时,n=4时,表 51 柯特斯

4、系数,柯特斯系数C(n)i仅与n和i有关,与被积函数f(x)无关,且满足,515,柯特斯公式对f(x)=1是准确成立的,柯特斯系数的特点,例1 试分别用梯形公式和辛普森公式计算积分 解:利用梯形公式,利用抛物线公式,原积分的准确值,5.1.2 误差估计 现对牛顿柯特斯求积公式所产生的误差作一个分析。牛顿 柯特斯求积公式的余项为 易知,牛顿柯特斯求积公式对任何不高于n次的多项式是准确成立的。这是因为 f(n+1)()0 故 Rn(f)0,510,代数精度 一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)0),而对于某一次数为m+1的多项式并不准确成立(即Rn(f) 0),

5、则称这一求积公式的代数精度为m。 牛顿 柯特斯求积公式的代数精度至少为n，若n为偶数，则至少具有n+1次代数精度。通常在基点个数相等的情况下,代数精度愈高,求积公式愈精确。 梯形公式、辛普森公式、柯特斯公式分别具有1、3、5次代数精度,例5.1 分别利用梯形公式、Simpson公式、Cotes公式计算 , n=1,2,3,4,5，并与用牛顿-莱布尼兹公式计算的结果进行比较,解 计算结果列于表5-2中,证 由式 知,梯形公式的余项为,x-a)(x-b)在区间(a, b)内不变号,f()是x的函数且在a,b上连续,故根据积分第二中值定理参见有关数学分析教材中“一元函数积分学第二中值定理”。 知,存

6、在某一(a, b)使,定理2 (抛物线公式的误差)设f(x)在a, b上有连续的四阶导数,则抛物线公式的误差为,定理1 (梯形公式的误差)设f(x)在区间a, b上具有连续的二阶导数,则梯形求积公式的误差为,如果在每个子区间上使用梯形公式，就得到复合梯形公式。将积分区间a,bN等分后的节点记为xk，xk=a+kh(k=0,1,2,N )，在每个子区间xk ,xk+1 (k=0,1,2,，N-1)上应用梯形公式,1.复合梯形公式,5.2 复合求积公式,再求和得,1.复合梯形公式,其中xk=a+kh (k=0,1,2,N,1.复合梯形公式,复合梯形公式,在每个 上用梯形公式,Tn,2.复合Simp

7、son公式,如果在每个子区间上使用Simpson公式，就得到复合Simpson公式。将N等分后的每个子区间再对分一次，于是共有2N+1个节点， (k=0,1,2,2N)，在每个N等分的子区间x2k , x2k+2 (k=0,1,2,N-1)上应用Simpson公式,再求和得,2.复合Simpson公式,其中 (k=0,1,2,2N,2.复合Simpson公式,复合 Simpson 公式,Sn,注：为方便编程，可采用另一记法：令 n = 2n 为偶数， 这时 ，有,其中 (k=0,1,2,4N,3.复合Cotes公式,4、复合Simpson公式算法,1) 输入a,b,N,2,3) 当 i=1,2

8、, ,N时 做循环,x=x+h,s=s+4f(x,x=x+h,s=s+2f(x,4,例 5.2：利用数据表,计算积分,这个问题有明显的答案,取n = 8用复合梯形公式,取n=4,用辛普森公式,二、复合求积公式的余项,梯形公式的余项为,对于复合梯形公式则有,若 在a,b上连续，则存在 ，使,1、复合梯形公式的余项,所以,由 在a,b上连续可知， 在a,b上有界，于是存在常数M2，使,1、复合梯形公式的余项,故,2、复合Simpson公式的余项,同理,由 在a,b上连续可知， 在a,b上有界，于是存在常数M4，使,故,3、复合Cotes公式的余项,由 在a,b上连续可知， 在a,b上有界，于是存在

9、常数M6，使,同理,故,当 时， ，于是从这些余项公式可以看出， 当时，复合求积公式TN ,SN , CN都收敛于定积分值I，而且收敛速度一个比一个快,二、复合求积公式的余项,例5.3 用复合梯形公式、复合Simpson公式、复合Cotes公式在取相同节点的情况下，计算定积分 的近似值。设把区间8等分,1) 用复合梯形公式计算，相当于取,2) 用复合Simpson 公式计算，相当于取N=4，把区间0,1N等分，然后在每个子区间上使用Simpson公式,3) 用复合Cotes 公式计算，相当于取N=2，把区间0,1N等分，然后在每个子区间上使用Cotes公式,的准确值为0.9460831,回顾：

10、复合求积公式的余项,1. 复合梯形公式的余项,2. 复合辛普森公式的余项,3. 复合柯特斯公式的余项,一、变步长梯形公式,1.当把区间a,b 等分时，步长,复合梯形公式为,2.当把区间a,b 等分时,步长,复合梯形公式为,5.3 变步长求积公式,改写上式得,复合梯形公式的递推公式,二、变步长梯形公式算法,1. 输入a,b,精度eps,2. h=b-a,3. 做循环,s=s+f(x,4. 则返回 3,5. 输出T,1. 将a,bN等分后复合梯形公式的余项，h,2. 将a,b2N等分后复合梯形公式的余项，h,设 在a,b上变化不大，即有,于是,整理得,同理，由复合辛普森公式的余项可得,同理，由复合

11、柯特斯公式的余项可得,5.4 龙贝格求积公式,一、龙贝格求积公式 由变步长的求积公式可以看出，利用前后两次计算结果进行适当的线性组合，可以构造出精度更高的计算公式，这就是龙贝格求积公式的基本思想,一、龙贝格求积公式,一、龙贝格求积公式,对于复合Simpson公式，设将区间a,b分成 等份，即步长为 ，节点为 (k=0,1,2,2k,一、龙贝格求积公式,即,同理，由复合Simpson公式的前后两次计算结果作线性组合可以得到精度更高的复合Cotes公式,一、龙贝格求积公式,由复合Cotes公式的前后两次计算结果作线性组合，必可得到精度更高的公式,龙贝格(Romberg)求积公式,龙贝格求积过程,T

12、1,T8,T4,T2,S1,R1,S2,C1,C2,S4,T16,S8,C4,R2,龙贝格求积过程：T数表,引入记号Tk, i，其中i 表示外推的次数，k表示区间a,b对分的次数，即把a,b分成2k等份,如果f(x)充分光滑，那么T数表每一列的元素及对角线元素均收敛到所求的积分值，即,i固定,并且后者的收敛速度比前者快,龙贝格求积过程：T数表,复合梯形公式的递推公式,龙贝格求积过程：T数表,k=1,2,外推公式,龙贝格求积过程：T数表,k=0,1,；i=1,2,例5.3用龙贝格积分方法求 的近似值，精度要求为,解：令 , a=2,b=8,1)在2,8 上用梯形公式计算 k=0 h=b-a=6,

13、2)将区间二等分，此时 k=1 h=(b-a)/2=3,计算新增节点处的函数值,3)将区间四等分k=2 , h=(b-a)/4=3/2,计算新增节点处的函数值,例5.3用龙贝格积分方法求 的近似值，精度要求为,4)将区间八等分k=3 , h=(b-a)/8=3/4,计算新增节点处的函数值,例5.4用龙贝格积分方法求 的近似值，精度要求为,达到了精度要求，故取T0,3 作为积分的近似值，即,例5.3用龙贝格积分方法求 的近似值，精度要求为,2,3,4) 当 i=1,2,k时,5,则返回(3)；否则输出T0,k，结束,龙贝格积分算法,1) 输入积分上、下限a、b，精度要求eps,j=k-i,第5章 小结,第6次作业,用龙贝格求积公式求定积分 的近似值， 要求写出每一步计算的公式，精度要求10-5，每一步的计算结果都至少保留6位小数,作业：用龙贝格积分方法求 的近似值，精度要求为,解：令 , a=1,b=2,1)在1,2 上用梯形公式计算 k=0 h=b-a=1,2)将区间二等