第四章数值积分

1、第四章数值积分定积分的产生是有它重要的应用背景。例如要计算由数据点(xyj (i =0,1,2川丨，n)所围成的平面图形的面积；计算极限n2lim a，这些问题都与定积分有关。在数学分析或 n厂7 n高等数学中已讲过计算定积分的一些方法， 这些方法其最主要的理论基础就是被积函数的原 函数存在。但在实际应用和科学计算过程中， 有些定积分的被积函数的原函数不存在或原函 数比较复杂或不易求出，这时牛顿 -莱布尼茨公式就不好用了。例如定积分1 si n x10 x dx，01 cosx2 dx等其被积函数的原函数不存在。再例如由数据点(Xi，yj (i =0,1,2,|1丨，n)所围成的平面图形的面积

2、不能精确的表示成定积分，但可以近似的表示为数据点(xi, yi) (i =0,1,2,|,n)对应的某个函数的定积分。对这类问题可以用数值积分的方法来讨论和解决。数值积分的应用是较广泛的，尤其在一些实际问题的研究和解决中数值积分法起到了重要的作用，见文献17，20。4.1数值积分初步所谓数值积分就是用函数值的线性组合近似函数的积分值。就是说，如果函数f (x)在n区间a,b上的函数值f(x)(i =0,1,2,山，n)已知，则构造一个数值公式Af(xJ，以i=0b此来近似f(x)dx，即y an(4.1)f (x)dx 八 A f (x)i =0构造数值公式(4.1)的主要方法是利用插值法，即

3、对f(X)构造一个插值多项式p(x)，b用该插值多项式 p(x)的积分近似f(x)dx，即L abbf (x)dx ： p(x)dx( 4.2)aa1梯形公式若函数f (x)在区间a ,b 上的函数值f (a), f (b)已知，那么可以做出过点(a, f (a) , (b, f(b)的线性插值x bx aPi("訂伽市f(b)在区间a,b上用pi(x)代替f (x)得bbb x bxaa f (x)dx：aPi(x)dx二a(rf (a)-f (b)dxaaa a -bb -ab a= -a(f(a) f(b)(4.3 )2b a公式(4.3 )称为梯形公式，记为T(f(a) f

4、(b)。公式(4.3 )的几何意义就是用2梯形面积近似由f(x)所围成的曲边梯形的面积，见图 4.1。2抛物线公式a + ba + b设函数f(x)在区间a,b上的函数值f(a) , f (), f (b)已知，记c，那么2 2可以做出过点(a, f (a), (c, f(c) , (b, f(b)的抛物线，即有二次插值多项式p2(x)(x-b)(x-c) f (a).(x-a)(x-b)f (x-a)(x-c) f (b)(a -b)(a -c)(c_a)(c_b)(b_a)(b_c)b a记h，在区间a,b上用p2(x)代替f (x)得2bbb - aa f(x)dx ：P2(x)dx(f

5、(a) 4f(c) f (b)(4.4 )aa6占f(a) 4f(c)f(b)3公式(4.4)称为Simpson公式，记为(4.5)b aS= = (f(a) 4f(c)f(b)或 S(f(a) 4f (c)f(b)。3公式(4.4 )的几何意义就是用抛物线所围成的曲边梯形面积近似由f (x)所围成的曲边梯形的面积，所以公式(4.4 )也称为抛物线公式，见图4.2。3牛顿-科茨公式如果函数f(x)在区间a,b上的函数值f(x) (i =0,1,2,川,n)已知，则对f(x)可以做出n次Lagrange插值多项式nPn(X)八 h(X)f (Xji =0其中li(x)是 n次Lagrange基插

6、值函数，现用bbPn(x)dx来近似.f (x)dx，所以有aabf (x)dx ：aa Pn(X)dX j Af (Xi)(4.6)其中bA = f li(x)dxL a(4.7)此时把公式(4.6 )称为插值型数值积分公式。现设a = x0 ：为：11( : xn = b，且把区间a,b分成n个相等的小区间Xi,1，记每b a个小区间的长度为hx-x，即h，所以有Xj = a ih (i = 0,1,|1(, n)。令nx = a th,则有0三t乞n，这时有bA 二 li(x)dx-anhnt(t-1)|l|(t-i 1)(t-i -1)川(t-n)0(-1)nhni!(n -i)!(T

7、) h n祐沁讪(t-i Wfzmdt(4.8)引进记号c(n)(-旷n i!(ni)!n0t(t -1)|l|(t -i 1)(t-i -1)川(t-n)dt(4.9)则A =(b -a)c：n)，这时公式(4.6 )可写成bf(x)dx：anPn(x)dx =(b -a)' c(n)f (Xi)i=0(4.10)公式(4.10 )称为牛顿-科茨(Newton-Cotes)公式，ci(n)称为牛顿-科茨系数。利用公式(4.9)可以计算出常用的牛顿-科茨系数c(n)( n <6),计算结果见表4.1。例1试用梯形求积公式、抛物线求积公式、牛顿-科茨求积公式计算定积分解(1)利用梯

8、形求积公式有3 23 _ 1221 In2 xdx (In21 In2 3) = 0 1.206949 = 1.206949a +b(2) 用抛物线公式计算，因c2，所以有23 23 - 1222In2 xdx (In21 4In 2 2 In2 3)1 61 (0 4 0.480453 1.206949) =1.04292(3) 若用n=4的牛顿-科茨求积公式，因Xj=a，ih，= 0.5，并由4 4表4.1得3 2721622216272In xdx (31)( In 1 In 1.5 In 2 In 2.5 In 3)190451545907162167=2(00.1644020.480

9、4530.8395891.206949)9045154590= 1.029817而定积分的具有7位有效数字的准确值为3* In2 xdx = (xln2 x -2x1 n x 2x) 3 =1.029173表4.1nc(n)111221412666133138888716216749045154590192525252519528896144144962884199349941684035280105280358404.2复化数值积分公式从以上例1的计算结果可以看到，一般情况下，抛物线求积公式和牛顿-科茨求积公式要比梯形求积公式好，其主要原因是推导抛物线求积公式和牛顿-科茨求积公式时把积分区间

10、等分的个数（分别为 2个和4个）比梯形公式的区间个数（1个）多。但对某些函数抛物线求积公式并不一定要比梯形求积公式好，如图4.3的情况。但是无论是哪一类的函数，只要被积函数在积分区间上有定义，那么不断增加积分区间的等分的个数时，求积公式得到的结果会逐步的逼近原定积分的精确值。所以考虑到数值求积公式的准确性，先把积分区间a,b分成n个相等的小区间Xi,Xi訂，记每个区间长度为h二人-人，在小区间人，人訂上x 1.-li利用梯形公式得I ' f(x)dx 眉(f(x)+ fg) i2所以有f f (x)dx = f fi+f (x)dx 賂F2 ( f (Xi) + f (Xi J)a7

11、xy2hp(f(X。)2f(Xj 山fX)二(以1)该公式称为复化梯形公式，公式(4.11 )的几何意义是用若干个小梯形面积之和来近似由f(x)所围成的曲边梯形的面积，见图4.4。与推导复化梯形公式 的方法相同，若把区间a,b分成2n个相等的小区间人必訂，b a每个区间长度仍记为 h = x 4 -x，在小区间x2i,x2i 2上利用抛物线(Simpson)2n公式得x2i 2hx f(x)dx(f(X2i) 4f(X21)f(X22)x2i3所以有bn 4 x2i 2n 4 hf f(x)dx=5： J f(x)dxZ ：(f(X2i)+4f(X2i+f(X2iQ)ai 八3hnn 4= (

12、f(X°) 4、f(X2i4)2、f(X2i) f(X2n)=Sn(4.12)3i =1i d该公式称为复化抛物线(Simps on )公式。4.3数值积分公式的误差估计定义对一个一般的数值求积公式(4.13 )nf (x)dx 八 A f (人)7其中A是不依赖于函数 f (x)的常数。若公式(4.13 )中的f(x)为任意一个次数不高于 m次多项式时，其等号成立，即有bna f(X)dX 八 Af (Xji=0而对f (x)是m 1次多项式时，公式(4.3 )不能精确成立，则说数值求积公式(4.13 )具有m次代数精度。根据代数精度定义和数值求积公式的构造过程有以下结果。定理1设

13、f(x)ECa,b , f(f(x)在a,b上存在，则牛顿-科茨公式(4.6 )的代数精度至少为n，当n为偶数时，牛顿-科茨公式的代数精度为 n 1 °证明 设Pn(x)是f (x)的n次Lagrange插值多项式，因f (x)满足以上条件，所以有f(X)二 Pn(x)(n 1)()(n 1)!<b其中(X) =(X -Xo)(X -x!)|l(x -Xn)。当f (x)为n次多项式时，定理结果是明显的。设bf(n 1()a (n 1)! (x)dxf(x)为n1次多项式，最高次项的系数为bm，则得f (n1)(x)二bn.1( n 1)!，由此得bbf f (x)dx _ f

14、 Pn (x)dx =aa0bn .(n 1)!t(t-1)(t-2)|H(t-n)dt(n 1)! 0n二 bn1hn 2 .0t(t -1)(t -2)川(t-n)dt令n = 2k , k为整数，并再次做变量替换u = t - k，则有n0t(t -1)(t -2)|l|(t - n)dt2k=(t(t _1)(t _2)川(t_k)(t _k_1)IH(t_2k+1)(t_2k)dtkk(u k)(u k T)丨l|u(u1)川(u-k 1)(u-k)du令 H (u) = (u k)(u k -1)|u(u -1)lH(u - k 1)(u -k)，则H (_u) -(_u '

15、; k)(-u ' k -1)Hl(-u)(-u -1)IH(-u -k 1)(-u k) = (-1)2k 1H(u -H(u)故H（u）是奇函数，因此有n0t(t-1)(t-2)|(t- n)dtk(u k)(u k -1)|"u(u -1)|(u -k 1)(u -k)du = 0 _k所以得bb上 na f(X)dX 一 乜 Pn(X)dX =bn N 2。t(t - 1)(t - 2川(t - n)dt = 0即当n为偶数时，牛顿-科茨公式的代数精度为 n 1。定理2若f (x) C:口，则对梯形公式(4.3 )有误差估计bb _aR(f,T) = f(x)dx (

16、f(a)f(b)a 2(4.14) Lf ( ). (a,b)12证明因对线性插值p'x）有余项公式所以R(x) = f (x) - P1 (x)=P(x a)(xb) , a 乞乞 bbbb f '牡)R(f ,T) = a f (x)dx- a m(x)dx = a (x -a)(x-b)dx因(xa)(xb)在区间a,b上不变号，且f (x) Cb，由积分中值定理可得bR(f,T) = Jafp(x-a)(x-b)dx = T2 2 (b-a)312b(x _ a)(x _ b)dxa定理3若f（X）E C:,b，则对抛物线公式（4.4 ）有误差估计bb aR(f,S)a

17、f(x)dX-(f(a) 4f(c)f(b)寫宀）5证明根据定理1的结果，抛物线公式（4.4 ）的代数精度是3，所以公式（何不高于3次的代数多项式 p（x）都有bf P(x)dx =ab - a-(p(a) 4p(c)p(b)6（4.15）4.4 ）对任（4.16）a + b其中，若构造满足插值条件P3(a) =f(a), ps(c) = f(c), P3(b)二 f(b), P3(c) = f (c)的3次插值代数多项式 p3（x），则由（4.16 ）得bP3(x)dx 二a故b a= (f( a) 4f(c)f(b)6bR(f,SH ab _ abbf (x)dx(f(a) 4f(c) f

18、(b) f(x)dx-p3(x)dx6、a、a二：P(xa)(x-c)2(x b)dxa 4!因f（x）C：,b，且（x a）（xc）2（x b）在区间a,b上不变号，由积分中值定理可得b f (4) ( ' )2f (4)()(x -a)(x -c)2(x -b)dx =a 4!4!二f()(_ ® £ 一色辺f()120 2880b(x-a)(x-c)2(x-b)dxa4!由此容易得到定理结果。把公式（4.14 ）、（ 4.15 ）分别称为数值积分公式4.3 ）、（4.4 ）的局部截断误差。考虑到误差的大小，把区间a,b分成n个相等的小区间X" 訂，得

19、到了复化梯形公式 （4.11），此时在小区间x,Xj 1上利用梯形公式的截断误差公式（4.14 ）得x 1hXi f(x)dxH(f(x) f(Xi1 沪(Xi 1 -为)123-f ( i) i (xX 1)由此利用微分中值定理得bR(f,Tn)二f(x)dx -Tn力V *二_(x 1 -X)3i=0122 x xf(x)dx：T(f(Xif(Xi1)f (i) 一 Jf()(a,b)即""（艸一“ 一冒公式（4.17 ）称为复化梯形公式（4.11 ）的整体截断误差公式。利用类同方法可得复化h2f ( )(a,b)(4.17)Simpson公式（4.12 ）的整体截断误

20、差 为bR(f,SJ=J f(x)dx& =a益W()(a,b)(4.18)b ab a注意：公式(4.11 )、(4.17 )中 h =，而公式(4.12 )、(4.18 )中 h =。n2n例 2 设 f (0) =1, f (0.5) =2 , f(1) =4 , f (1.5) =6 , f (2) =2，则用复化 Simpson2公式计算jQf (x)dx，若有f兰M，则估计复化Simpson公式的整体截断误差。20 1解因h，则利用复化 Simps on公式有422 h0f(x)dx ： ff(O) 4(f(0.5)f(1.5) 2f(1) f(2)=43石因复化Simps

21、on公式的整体截断误差为R(f ,S4) = f f (x)dx S4 = h4f ( r e (a,b)所以有R(fS)二2 014(4)yr"1802bM =14404.4逐步梯形方法与龙贝格公式复化梯形公式、复化Simpson公式的整体截断误差是比较小，一般情况下都可以满足精度要求。若不满足精度要求， 那么可以把积分区间的等分个数再变大, 的精度。逐步梯形方法就是收敛速度较快的一个提高精度的方法。提高数值积分公式当把区间a,b分成n个相等的小区间Xj,x*时有复化梯形公式bhagX石(伽2畑川2g) g"(4.19)若把区间a,b分成2n个相等的小区间x,x时对应的复

22、化梯形公式为(4.20)bhf(x)d-(f (x0) +2f X) +山 + 2 f (X；nJ + f (X；n) =T2n a2b ab ab a -其中h，因为h22h , x = x2i，所以T2n与的关系为2nn 2nT2n =h2(f(x0)2f(X1)川 2f(X2n" f(X2J)1 h= 22(f(X0)2f(X1)2f(X1)2f(X3)2f(X2)川 2f(Xnj) 2f(X2nT f(Xn)1 h匚Tn 2(f(X1)f(X3)川 f(X2nJ中宁®T)穿(4.21)所以利用公式（4.21 ）可构造出序列Ti,T2,T4,T8j|,T2k，H|其构

23、造序列的公式（4.21 ）可改写为k 11hahaT?k =汀2*丄 亍' f -1）丁），k=1,2川（4.22）则公式（4.22 ）就称为逐步梯形公式。另外，因为有S4T2n - TnSn :4 -1所以利用Tn与T2n可以构造出复化 Simpson公式Sn，即可用代数精度为1的公式推导出代 数精度为3的公式sn，同理再利用Sn与S2n可以构造出复化 Newt on-Cotes公式Cn，即由代数精度为3的公式推导出代数精度为5的Cn，以此类推，可以得到以下序列m (k 1)4 Tm Jm A4 -1k =0,1,2,|l(; m =1,2,川(4.23)由公式（4.23 ）得到的序

24、列就称为龙贝格（Romberg ）序列，（4.23 ）就称为 龙贝格（Romberg ）公式，龙贝格公式的截断误差阶是O（h2m 2）。4.5高斯（Gauss）型求积公式以上求积公式的特点是对给定的数据点（人，） （i =0,1,2jl|, n）作对应的插值多项式bp（x），用插值多项式p（x）的积分近似.f（x）dx，其对应的代数精度一般不超过n 1,逐La步梯形方法与龙贝格公式也不列外。高斯（Gauss）型求积公式的思想是：在节点数目固定为n的条件下，在区间a,b内适当地选择节点x和待定系数A，使求积公式nf (x)dx ：A f (x)具有最高的代数精度。定义 在节点数目固定为n的条件下

25、，在区间a,b内适当地选择节点 xi和待定系数 A，使求积公式f (x)dx 壯瓦 A f (x)a7(4.24 )具有2n -1次的代数精度，此时公式（4.24 ）就称为高斯型求积公式，节点Xi称为区间a,b上的高斯点。首先分析一下对于 n个节点，公式（4.24 ）可以达到的最大代数精度是多少? 假设式（4.24 ）对任意m次多项式mmJ , LPm(x)二amX- amjXa。是准确成立的，于是有am axmdx am.ab m dbbx dx 11 ( ad xdx a0 dxaaan工為 A(amXam必m，V aiXi a。)i A(4.25)b i记i二xidx （i =0,1,|

26、（, m），并重新组合（4.25 ）的右端项得*aamm am Jm71(印叫-& %nnnn=a Axm 飞皿八 AXim71 - a AjXa。'，Ai 1i 1i Ai A(4.26 )从系数ao,ai Jl|,am的任意性，（4.26 ）式成立的充分必要条件是人十人2十川+片=#0A xi + A2X2 + 川 + An Xn =气“ AX2 +A2X； + 川 + AnX：=卩2（427）III HI IIIAxA2Xm III AnXm = Jm在方程组（4.27 ）中有2n个待定常数，最多能给出2n个独立条件，且有 m 1个等式，所以m值最大为2n -1。即对于n

27、个节点的求积公式（4.24 ），可能达到的最大代数精度是 2n -1。并可证明，方程组（4.27 ）当m = 2n-1时是可解的。问题是利用条件（4.27 ）如何选取节点Xi及系数A。下面对n =2的情形介绍方程组（4.27 ）的解法。不失一般性,把积分区间取成-1,1，这是因为利用变换a +b丄b a丄xt2 2总可以将区间a,b变成-1,1，而积分变为b a fgdxp1 t)dt =2b -a""2-对n = 2的情形，问题就变成了如何选择节点x1, x2及系数A , A2使(4.28)1f (x)dx ： a f(Xi) A2 f(X2)32对任何三次多项式 f （

28、x）二a3X - a?x - aix - a°都能精确成立。由式（4.27 ）只要解非线 性方程组(4.29)X + A = 2 An +A?x2 =0A< +A?xf =23A x3 十 A2x； = 0求出节点x1，x2及系数A1，A即可。但这种方法当n较大时就比较困难。此时可用正交多项式的特性来求节点 xi。把三次多项式 f（x） =a3x3 - a2x2 a1x - a0总可以表示为f(X)=( rX ：0)(x -Xi)(X -X2) (bx bo)两边积分得1 1 - - 1f (x)dx 二(,x ：0)(x -xj(x -X2)dx(dx b°)dx(

29、 430 ) 因f(xj =b1x1 b0， f (x2) = dx2 b0，且求积公式(4.28 )对任意一次多项式都精确成 立，所以得1,4(b1x b°)dx = AWN b°) A2(biX2 g) =Af (xj(x?)因此若对任意多项式 -1x恒有1 _ _J（ X + P0）（X - Xj）（X - x2 ）dx = 0（ 4.31）那么由（4.30 ）式就有14 f （x）d = A）f （x1） A2f（x2）（4.28）所以，当节点的选取满足条件 （4.31）时，对任意三次多项式，公式（4.28）是精确成立的。由于（4.31 ）式对任何'-0 ,

30、 '-1都成立，所以必须有1 14（x _片）（x _x2 ）dx 二 0，x（x _xj（xx2）dx 二 0由此得22%x2 =03 1x1 x2 = 01解得Xj - -x2 ： ，再利用求积公式(4.28 )对f(x) =1, f(x)二X准确成立，得到.3A + 人=2A% +A2x2 =0解得Al二A2 = 1，这就得到求积公式把条件(4.31 )称为正交条件。再对一般情形讨论高斯型求积公式，考虑积分bI - a (X) f (X)dX其中co(x)工0是权函数。问题是在区间a,b内如何适当选择节点 X1,X2|,Xn使求积公式(4.32)b/ (X)f(x)dx当f (X

31、)是不高于2n -1次的多项式时精确成立。与前面 n = 2的情形一样，把2n -1次的代数多项式f(x)总可以表示为f(X)二q(x) n(x) r(x)其中''n(x) =(x -X!)(x - x2)IH(x -xn)， q(x) , r(x)都是不超过n -1次的多项式，于是有bbb(x) f (x)dx =，(x)q(x) n(x)dx 亠 i 心(x)r(x)dxaaa如果对任何不超过 n-1次的多项式q(x)都有(4.33)b(x)q(x) n(x)dx = 0a则因求积公式(4.32 )对任何一个不超过 n-1次的多项式精确成立，即有(x)r(x)dx 二、Ar

32、(xjai d而有f(Xi)=r(Xi)，所以得W(x) f(x)dx =迟 A f(Xi)ai4也就是说，只要选取节点满足条件(4.33 ),则求积公式(4.32 )的代数精度就能达到 2n - 1。所以有以下结果。定理对于插值型求积公式(4.24)nf (x)dx 八 A f (x)i=1其节点xi是高斯点的充分必要条件是n(X)=(X -Xi)(X-X2)| H(X -Xn )与任意次数不超过n -1次的多项式q(x)均正交。条件(4.33 )表明，n(x) =(x - Xi)(x - X2)|(x - Xn)和 q(x)是在区间a,b上关于权函数-(x)正交。从正交条件(4.33 )可

33、以解出高斯点 X1, x2|,xn，而由正交多项式的 性质可得一个n次正交多项式的n个零点是实数且不重复，并分布在区间 (a, b)内。所以对 给定的权函数 (x)总可以构造出正交多项式.n(x)，而后再用方程组(4.27 )或用公式b dx(4.34 )a(X Xk)n(Xk)求解待定系数 A ,高斯型求积公式(4.24 )就被确定。1例3试构造高斯型求积公式f (x)dx Af (xi) A2f(x2) AJ(X3)解 先求满足条件(4.33 )的正交多项式 .n(x)-(x-xJd-XzXx-x3)的零点Xi,X2, X3。因条件(4.33 )对任意2次多项式q(x)都成立，所以有1 1

34、J(x xj(x x2)(x x3)dx =0，斗x(x X! )(x x2)(x x3)dx =0，1 2* x (x - 片)(x - x2)(x - x3 )dx = 013(X1+X2*"=011由此得(X/2 %X3 - X2X3)=3 511：(X1+X2+X3)+：XX2X3=0 .53解得X1 X2 X3 =0, X1 X2X3=0，所以至少有一个X=0，不妨设X2=0 ,所以得由于求积公式4.32 )对f(X)=1,X, X2精确成立，即得A + A2 + 人=2< A% + Ax2 + A5X3 =02222Axi +A2X2 +A3X3 =一l3把x2 =

35、 0 , Xi = -X3 = - 一 3代入上式得A + A + 人=2= A A3= 0AA105由此求解得A沁匚8，所以得具有5次代数精度的高斯型求积公式为A +人=一 .9例如用Legendre多项式Pn(X2£，-0，i，2，1"(4.35)3)+9f(0)+if(£应该强调的是：当高斯点的个数较多时用以上例子的方法求高斯点和系数时计算量较大 且有一定的难度。由于Xi,X2H,Xn是区间a,b上的高斯点的充分必要条件是多项式n(x) =(X Xi)(X X2)川(x Xn)是a,b上n i次正交多项式。所以说，求高斯点就是求正交多项式的零点。因此当n较大时最好用已知的正交多项式的零点来构造对应的高斯型求积公式。的零点来构造高斯型求积公式得到Gauss-Legendre求积公式。因Legendre多项式是区间-i,i上关