MATLAB与数值分析第二部分—数值积分

1、第二部分第二部分 数值分析基础理论与算法数值分析基础理论与算法四、数值积分四、数值积分 （第五章，张德丰等编教材）（第五章，张德丰等编教材）数值积分数值积分利用利用插值多项式插值多项式 ，则积分易算。则积分易算。)()(xfxPn 问题的提出：问题的提出： f(x)的原函数没有具体的解析表达式或表达式的原函数没有具体的解析表达式或表达式很复杂不适宜计算，只有很复杂不适宜计算，只有f(x)的离散数据点。的离散数据点。如何求如何求 ？( )baf x dx( )baIf x dx寻找近似计算方法寻找近似计算方法 在在a, b上取上取 节点节点a x0 x1 xn b，做，做 f 的的 n 次插值次

2、插值多项式多项式 ，则积分为：，则积分为：0( )() ( )nnkkkPxf xlx00()()()()nbbkkaaknkkkf x d xf xlx d xA f xAk bakjxxxxkdxAjkj)()(由由节点节点决定决定，与与f(x)无关无关。插值型积分公式插值型积分公式对于次数不高于对于次数不高于n n的多项式的多项式f(xf(x) )，其余项为零，其余项为零，求积公式求积公式至少至少具有具有n n次代数精度次代数精度0(1)0 ()() ()()()()()(1)!nbkkakbbnnaannbxkakR ff x dxA f xf xPxdxRx dxfxxdxn产生的误

3、差，即余项为：产生的误差，即余项为：线性插值求积公式线性插值求积公式（即（即梯形积分公式梯形积分公式）00()()()()nbbkkaaknkkkf x d xf xlx d xA f x2 ( )( )b abaf affxbx d称为梯形称为梯形求积公式求积公式1( )( )( )( )xbxaf xP xf af babbaab抛物型求积公式：220 , ,2( )() ( )bb()(b)()()()()b22( )()( )bbbb2()(b)()(b)()()2222kkkaba bablagrangeP xf x lxaaxxxa xxa xbaf aff baaaaaaaba

4、b将区间先先 二 二等等分分： 再再作作二二次次插插值值二次插值求积公式二次插值求积公式（即（即SimpsonSimpson公式）公式）22bbb ()() ( ) 2()() () ()() ( )222()aaaxx b f ax a x b fx a xf bb a 称为称为也称为也称为2 2阶阶Newton-CotesNewton-Cotes公式公式2( )4( )( ) 62( )bbaabaabf aff x dxPf bx dx用更高阶插值来构造数值积分方法，称为用更高阶插值来构造数值积分方法，称为Newton-Cotes方法方法积分为：积分为：Newton-Cotes求积公式求

5、积公式1020301015(4)01202:( ) ()()( )2122:( ) ()4 ()()( )3903: xxxxTrapezoidalhhf x dxf xf xfxxSimpsonhhf x dxf xf xf xfxxSimpson1阶阶阶梯形（）积分公式： ， 其中，积分公式： ， 其中，30405(4)0123037(6)0123403/8 33( ) ()3 ()3 ()()( )8804:28( )7 ()32 () 12 ()32 ()7 ()( )45945xxxxhhf x dxf xf xf xf xfxxhhf x dxf xf xf xf xf xfx阶积

6、分公式： ， 其中， ， 其中，4xf(x) P4(x)f(x) P3(x)插值型求积公式的插值型求积公式的MATLABMATLAB函数实现函数实现定定 义义若某个求积公式所对应的误差若某个求积公式所对应的误差R f 满足：满足：R Pk =0 对对任任意意 k n 阶阶的多项式均成立，且的多项式均成立，且 R Pn+1 0 对对某个某个 n+1 阶多项阶多项式成立，则称此求积公式的式成立，则称此求积公式的代数精度代数精度为为 n 。代数精度：代数精度：(Degree of accuracy)刻画数值积分方法优劣的刻画数值积分方法优劣的一种度量一种度量即：如果某求积公式对于即：如果某求积公式对

7、于但对于但对于，则称此求积公式的，则称此求积公式的 。注意：注意：n n不是插值不是插值节点数节点数, ,是区间数是区间数怎样验证代数精度：怎样验证代数精度：022211( )()( )1, ,., 1()2.,.1()1bnkkkankkknnnkkf x dxA ff xx xxAbaA xbaA xnbanx数积对一一般般地地，若若要要使使求求公公式式： 只只要要 能能具具有有 次次代代精精确确度度精精确确成成立立。即即banx dx梯形公式梯形公式例例1：对于：对于a, b上线性插值，有上线性插值，有1( )( )( )x bx aa bb aP xf af b)()()(2221bf

8、afdxxfAAabbaab 考察其代数精度。考察其代数精度。f(x)abf(a)f(b)解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 f(x) =1： baabdx111 2 ab=代入代入 f(x) = x ：=代入代入 f(x) = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代数精度代数精度 = 1利用梯形公式利用梯形公式22223( )( )1( )( - )-,1 4 1 1-( )6-( )( )2-,-( )4( )62462babababasimpsonsimpsonP x dxf xf xb ab asimpson

9、b af xbaf xxf xb ab asimpsonab aa bf aff b 积为例例 ：明明方方法法的的代代精精度度明明： 分分公公式式利利用用公公式式算算得得：，利利用用公公式式 得得到到： 证数为证当时，计当时( - )()( )2bab a a bbf x33222222-( )( )3,-(-)()4()( )623bababaf xxf xsimpsonbababaaabbabf x，利利用用公公式式 得得到到：当时 3443223333455444( )-( ),42-( - )()4 ()( )624( )-( ),52-4 ()( )62babababaf xxbab

10、 af xsimpsonhb abab a aa bb ababf xf xxbab af xsimpsonhb abaabf x 当时，利用公式计算得到： 当时，利用公式计算得到： 代数精度代数精度 = 3思考：思考： 代数精度是否是越代数精度是否是越高越好？高越好？显然，显然，由于龙格现象，不由于龙格现象，不稳定性越大稳定性越大, , 积分误差积分误差可能增加可能增加。梯形公式的误差梯形公式的误差,32 ( ) ( )( )-( ) - ( ) 12 ( , )( - )a bbbnaaf xCR ff x dxP x dxa bfb a则：定定理理 : : 若 插值求积公式的误差插值求积

11、公式的误差定理证明313( )( )( )()()( , )2( )()()()( ) ()()()( )12( )( )12babaff xxxa xba bPfxa xb dxabfxabaRxb dxfff 明明： 由 由插插值值公公式式余余：证项 辛普森辛普森(Simpson)(Simpson)求积公式的误差：求积公式的误差：4 , 5544( ),: 462( )( ) 902880a bba()()f xbaabR(f)f(x)dxf(a)f()f(b)h (a,b)C(b a)ff 则 若 若 理理 定定 2bahNewton-Cotes求积公式误差求积公式误差3(0022002

12、)( )1-,( - )/( , )( )2( )( , )( ) ( )( )(1).()(2)!nniiinnnbiiaina f xnNewton Cotesxa xb hb ana bnf x nf xCa bhff x dxa f xtttn dtn定理： 设为点公式， 且， 则存在满足： 当 为偶数，且 阶连续，即，时：2(1)100( )1( )( , )( ) ( )( )(1).()(1)!nnnbniiainnnf x nf xCa bhff x dxa f xt ttn dtn 当 为奇数，且 阶连续，即，时： 结论是什么？结论是什么？ 怎么办？怎么办？结论：结论：应用应

13、用高阶型高阶型插值求积公式计算插值求积公式计算 会出现会出现 数值数值不稳定不稳定，而低阶公式（如梯形、辛普生公式），而低阶公式（如梯形、辛普生公式）又因又因积分区间步长过大积分区间步长过大使得使得离散误差大离散误差大( )baf x dx缩小缩小h h，即把积分区间分成若干子区间，在，即把积分区间分成若干子区间，在每个每个子区间子区间上使用上使用低阶积分公式低阶积分公式，再将结果，再将结果加起加起来来，这种公式称为，这种公式称为复合求积公式复合求积公式。复合求积公式：复合求积公式： 复合梯形复合梯形求积公式：求积公式：,(0,., )kbahxak hknn在每个在每个 上用梯形公式：上用梯

14、形公式：,1kkxx 高次插值有高次插值有Runge 现现象，故采用象，故采用分段低次分段低次插值插值nkxfxfxxdxxfkkxxkkkk,., 1,)()(2)(111 h1111111( ) ()() ()()22( ) 2()( )2nnbkkkkkkakknkkxxhf x dxf xf xf xf xhf af xf b= Tn3321112() ()()()121212()( ),( , )12nknnkkkkkfhhhR fffbanhba fa b 中值定中值定理理xk，xk+1上的误差上的误差分段低次合成分段低次合成的复的复合梯形求积公式合梯形求积公式111( ) ()4

15、 ()()62kkxkkkkxxxhf x dxf xff xkx12kkxx1 kx111100 ( )4()2()( )(62)nnkkkkkbaxxhf aff xxfdf xb= Sn)(2180)4(4 fhabfR 复合复合Simpson公式：公式：利用利用Simpson公式，在公式，在 xk，xk+1上有上有112211;2kkkkkkhxxhxxxx复合复合SimpsonSimpson求积公式的求积公式的MATLABMATLAB实现实现复合求积举例：复合求积举例：10?xnnIdxnsimpsonnSTe计积数算算分分 保 保留留五五位位有有效效字字。用用梯梯形形分分公公式式算

16、算 ， ， 和 和分分公公式式算算试积 计积 计222-4-.(,)( ) 12( )( )(0,1)( )-11(,)( )1012122nnxnb aTR fh fTf xfxxfxeeb aR fh fehT误解解 的的差差因因公式为： 又： ，则 -46100.014968,167.3088nhenh为证够保保足足精精度度，分分大大于于，即即则 积区间数 应同理：同理：(4)44-41( )(0,1)( )11( ,)1028802xnfxxfxeeR feSh以以及及 , ，则 两种方法谁好？两种方法谁好？4(4) ( )1802nbahSR ff 误：由由的的差差公公式式1/4-4

17、1440100.4798,12.08443hennh为证够保保足足精精度度，分分大大于于则 积区间数 应 ，所以，。SimpsonSimpson积分公式分段数远远小于梯形公式。积分公式分段数远远小于梯形公式。4.2 4.2 高斯型高斯型积分积分0( )()nbkkakf x dxA f x即：在如下求积公式中：即：在如下求积公式中： 将节点将节点 x0 xn 以及系数以及系数 A0 An 都作为都作为。令令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式具有代入可求解，得到的公式具有。这样的。这样的节点称为节点称为Gauss 点点，公式称为，公式称为是不等距划分插值

18、型是不等距划分插值型求积公式求积公式通过适当选取节点位通过适当选取节点位置，可提高代数精度置，可提高代数精度若构造的若构造的n+1n+1个节点的插值求积公式个节点的插值求积公式 的代数精度是的代数精度是2n2n1 1，称它为，称它为高斯型高斯型积分积分。0( )()nbkkakf x dxA f x. . . 01n0011n2n+12n+12n+12n+10011n1AAAxAA.A.xAA.Ababnabnadxdxxxxdxxxx将将 f (x) = 1, x, x2, , x2n+1 代入求积公式可得，代入求积公式可得，联合求解，有联合求解，有2n1次代数精度次代数精度构造构造n+1个

19、节点的个节点的Gauss 型求积公式型求积公式0( )()nbkkakf x dxA f x注意注意: 高斯积分是不等距划分插值型求积公式高斯积分是不等距划分插值型求积公式2n+22n+2个方程个方程联合求解联合求解练习题：练习题：Gauss Gauss 型求积公式的型求积公式的MATLABMATLAB实现实现具有具有3阶代数精阶代数精度度，比梯形公式，比梯形公式1阶代数精度高阶代数精度高例例1 1：在在两点数值积分两点数值积分公式中，如果积分点也作为未公式中，如果积分点也作为未知量，则有知量，则有4 4个未知量个未知量, ,可以列出可以列出4 4个方程个方程：（以：（以-1,1-1,1为例）

20、为例）1010,xxaa032021111133113001122112001111001110dxxxaxadxxxaxaxdxxaxadxaa可解出：可解出：31,31, 1, 11010 xxaa两点高斯积分公式两点高斯积分公式1111()(3)3(fff x dx很重要！很重要！1.511.51.5111.51.512-1.763407241879019(1.5-1)1.79214( -1.259992724699278)xxxe dxe dxeeeee dxtx计积顿 莱兹积点积变换例例 ：算算分分解解：方方法法一一、牛牛布布尼尼公公式式： 方 方法法二二、梯梯形形插插值值分分： 方

21、 方法法三三、二二高高斯斯分分： 、 、量量代代：令令 1.51(5)/41-1- 33(5)/4(5)/41(5)/433-11-11421111( )8193076 1.76334(-)()3673438xttf x dxffe dxedxedxee则则带点积（） ： ： 、 、入入二二高高斯斯分分公公式式： 精度比较：精度比较： 梯形插值积分：选择梯形插值积分：选择被被积函数端点积函数端点构造线性函构造线性函数，来近似被积函数。数，来近似被积函数。 高斯积分：选择积分高斯积分：选择积分区区间内的点间内的点构造函数，近构造函数，近似被积函数。似被积函数。示意图示意图节点位置节点位置不固定，

22、不固定，精度更高精度更高节点位置节点位置固定固定112-11234 (-1)2 ()3 () (1) ( )3 (2) ( )( )( )( )( )( )112 1 3 1 23 baff xf xf x dxf x dxw f aw f bw faw fbf x 数积参数数尽解：(1)：例例3：3：确确定定下下列列值值分分公公式式中中的的，使使代代精精度度可可能能高高。112-12221212-11222212 ( ) 123 03 ( ) 1232 332310.8110 233f xxxxxdxf xxxxx dxxxxxx 或 1-0.716-0.41101.11655 x 或 12

23、2212342( )1,( )0 1 ( ),( )1 2 ( ),( )2 babaf xfxwwdxbaf xx fxbawaw bwwxdxf xxfxx (2) 3322212343244332231234 223 ( ),( )3 334bababawaw bw aw bx dxf xxfxxbawaw bw aw bx dx122212343322123444332212342223334wwbabawaw bwwbawaw bw aw bbawaw bw aw b12232422()12()12bawbawbawbaw 利用利用Matlab中符号计算得：中符号计算得：Matlab

24、符号计算代码：符号计算代码：syms a bA = 1, 1, 0, 0; a, b, 1, 1; a2, b2, 2 * a, 2 * b; a3, b3, 3 * a2, 3 * b2;b = b - a; (b2 - a2) / 2; (b3 - a3) / 3; (b4 - a4) / 4;w = inv(A) * b;w = simplify(w) 1-111111222211( )5 (- 0.6)8 (0)5 ( 0.6)591( )1 5 (0.6)8 (0)5 ( 0.6)2191 ( ) 50.650.60912 ( ) 5 (0.6)5 ( 0.6) 93 f x dxf

25、fff xfffdxf xxxdxf xxx dx 证数积数为例例4 明4 明： 值 值分分公公式式代代精精度度解解：1333311444411555516661 ( ) 5 (0.6)5 ( 0.6) 0912 ( ) 5 (0.6)5 ( 0.6) 951 ( ) 5 (0.6)5 ( 0.6) 091 ( ) 5 (0.6)5 ( 0.6) 0.249f xxx dxf xxx dxf xxx dxf xxx 16127dx推广：推广：n加权加权GaussGauss积分公式积分公式0( ) ( )()nbkkakx f x dxA f x权函数权函数例：例：求求 的的 2 点点 Gaus

26、s 公式。公式。dxxfx)(10 解：设解：设 ，应有，应有 (2n+1=) 3次代次代数精度数精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3 31130092211200721100521032xAxAxAxAxAxAAA2776. 03891. 02899. 08212. 01010 AAxx不是线性方程组，不是线性方程组，不易求解。不易求解。寻找简单的高斯节点计算方法寻找简单的高斯节点计算方法00( )()( -)(0,1,., )( )( )( ) ( )0nbkkaknkkkbax knxP xP xx dxf x dxA

27、f xx xn 积这点节点点条数过点项项 为定定理理： 求求插插值值型型求求公公式式 的的是是以以些些的的多多式式 与与任任意意次次是是高高斯斯不不超超的的充充分分必必要要件件零零均均正正交交，即即的的多多式式其它的高斯节点计算其它的高斯节点计算n+1阶多项阶多项式式证明：证明：必要性：必要性： 设设x0 xn 为为 Gauss 点点, 即公式即公式 至少有至少有 2n+1 次代数精度。次代数精度。 bankkkxfAdxxfx0)()()( 对任意次数对任意次数不大于不大于n 的多项式的多项式 Pm(x)， Pm(x) w(x)的次数的次数不大于不大于2n+1，则代入公式应则代入公式应精确成

28、立精确成立，即：，即： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 充分性：充分性： 要证明要证明 x0 xn 为为 Gauss 点，即要证公式对任意次点，即要证公式对任意次数数不大于不大于2n+1 的多项式的多项式 Pm(x) 精确成立，即证明：精确成立，即证明：0( )( )()n nb bmmk kmmk ka ak kx x P Px x d dx xA A P Px x bababamdxxrxdxxqxwxdxxPx)()()()()()()( 0 nkkkxrA0)(0()nkmkkA P x 设设)()()()(xrxqxwxPm Pm(x)为不

29、高于为不高于2n+1次的多项式次的多项式00()()()n nn nk kk kk kk kk kk kk kA A q q x x ww x xA A r r x x n+1个节点的插值求个节点的插值求积公式的代数精度至积公式的代数精度至少是少是nPm(x)与与w(x)正交正交正交多项式族正交多项式族 0 0, , 1 1, , , , n n, , 有性质有性质：kxn 推推： 插插值值型型求求公公式式 的的勒勒德德多多式式的的。是是解解阶让项点论积节nn证阶 让项阶 项让项明明： 勒勒 德德多多 式式是是多多 式式，且且勒勒 德德多多 式式相相互互正正交交。若取若取 w(x) 为其中的为

30、其中的 n+1，则，则勒让德多项式的解可勒让德多项式的解可作为高斯节点作为高斯节点练习题：练习题：Gauss-Gauss-勒让德多项式勒让德多项式求积公式的求积公式的MATLABMATLAB实现实现再解上例：再解上例： 101100)()()(xfAxfAdxxfxStep 1：构造正交多项式构造正交多项式 n设设cbxxxaxxx 2210)(,)(, 1)( 53 a0)(10 dxaxx0),(10 215910 cb 102110200)(53(0),(0)(0),(dxcbxxxxdxcbxxx 2 2 即：即：215910)(22 xxx Step 2：求求 2 = 0 的的 2

31、个根，即为个根，即为 Gauss 点点 x0 ，x1221/20)9/10(9/1021;0 xStep 3：代入代入 f (x) = 1, x 以求解以求解 A0 ，A1解线性方程组，解线性方程组，简单。简单。结果与前一方法相同结果与前一方法相同：2776. 0,3891. 0,2899. 0,8212. 01010 AAxx 利用此公式计算利用此公式计算 的值的值 10dxexx2555. 1 10dxexx2899. 08212. 0102776. 03891. 010eeeAeAxx Matlab 积分函数积分函数函数名函数名功能功能quad采用采用Simpson计算积分。精度高，较常

32、用计算积分。精度高，较常用quad8采用采用8样条样条Newton-Cotes公式计算积分。精度公式计算积分。精度高，最常用高，最常用trapz采用梯形法计算积分。精度差，速度快采用梯形法计算积分。精度差，速度快cumtrapz 采用梯形法求一区间上的积分曲线。精度差，采用梯形法求一区间上的积分曲线。精度差，速度快速度快sum等宽矩形法求定积分。精度很差，速度快，一等宽矩形法求定积分。精度很差，速度快，一般不用般不用cumsum等宽矩形法求一区间上的积分曲线。精度很差，等宽矩形法求一区间上的积分曲线。精度很差，速度快，一般不用速度快，一般不用MatlabMatlab 积分函数积分函数符号积分：

33、符号积分：4.3 4.3 积分积分方程方程的数值求解的数值求解( )( , ) ( )( )baFredholmy tk t s y s dsf t类积第第二二分分方方程程：怎么求解？怎么求解？方程中含有积分，方程中含有积分，而积分中又有未而积分中又有未知函数的方程知函数的方程应用广泛应用广泛求解思路：求解思路：n用数值积分代替下面方程中的积分( )( , ) ( )( )bay xk x s y s dsf x1( , ) ( )( ,) ()bnkkkkak x s y s dsA k x xy x即1( )( ,) ()( )nkkkky xA k x xy xf x1( )( ,) ()( )nkkkky xA k x xy xf x1( )( ,) ()( )nikikkiky xA k x xy xf x未知未知i=1,2,n进一步将进一步将x离散化取值离散化取值，则则方程转化为：方程转化为：用矩阵表示可得：1()()IK yfyIKf000101001011110011(,)(,).(,)( ,)( ,).( ,).(,)(,).(,)nnnnnnnnnA k x sAk x sA k x sA k x sAk x sA k x