华中科技大学计算方法课件-第三章数值微积分

《华中科技大学计算方法课件-第三章数值微积分》由会员分享，可在线阅读，更多相关《华中科技大学计算方法课件-第三章数值微积分（42页珍藏版）》请在装配图网上搜索。

1、第三章数值微积分在现实科学问题中，我们常常涉及到复杂函数的导数及定积分的计算，对其实施精确计算往往存在着异常的困难性，特别就定积分的计算而言，有些积分理论上可证明其原函数存在，但却无法用初等函数明显表出，因而无法获得其精确的积分值。此时，我们需要借助于数值计算。此外，微分方程及微分-积分方程的数值离散问题也必然牵涉到数值微积分。鉴此，本章将探讨导数及定积分的数值计算。3.1 数值微分3.1.1 基本差商逼近公式设函数f 在区间a,b内p+1 阶连续可微，将区间Q沟作均匀剖分：a=XQ Xi X2 t XN=b,这里 Xk=Xk-i+A,h=(b a)/N,k=1,2,则有 Taylor展开式.

2、九）=pE，=o塔皿（士 呼+。（犷+1）,z!=,7V-1,(1)由上式得=/（+2 八加+9，/Lri=+0吟,2h于 （Xk）=+m 2?）+/（跳h）+2）.h2因此，我们有如下差商逼近公式八“回 *1 1,（一阶向前差商）Q）rif J）卜 i）,（一阶向后差商）（3）/CZ一“线+（二阶中心差商）（4）乙/LfMl/（以十1 2器）+f（X m.（二阶中心差商）（5）类似地，我们还可进一步得到高阶差商逼近公式。此外，若二元函数”(巴力)在某区域a,b;c,d上二次连续可微，分划该区域的节点为Xj=a j h,h (b a)/N j=0,1,，N;tk=c-kr,tk+i=c+(k+

3、T=(d c)/M A;=0,1,.,M,则据Taylor展开定理，我们有a”（叼,tk）（叼+i,tQ-u（Xj,tk）dxh+。,dxdu（勺,tQ _dxduXj,udxh（町J+94）一2h+。,叱包十以川),2h+。仇 2）,（6）（8）d2u(xj,tk)_ u(xj+1,tk)-2M叼，4)+(叼 i,4)dx1h2+。优,(1 0)du（勺,tQ udta“（力力力。uJ,k+1)-tk)+9Tj5力k）一乜(IDdtdu(xj,tk)dt&u(叼,tk)dt+O(T),T”(叼，力 k+J (叼，力kJ +2)T_ n(叼,力 卜+i)n(叼)力,一i)+。(7 2)d2ux

4、j,tk)_ n(叼,tk+i)-2矶 叼,4)+(町,4 i)况2T2+。（产）,(1 2)(1 3)(1 4)(1 5)公式（6）-（1 5）将在第九章的偏微分方程数值解中发挥重要作用。3.1.2导数的插值计算公式及外推算法当函数/：Q,“一 R在区间Q阳 内0+1阶连续可微时，由（2.1 6）-（2.1 7）式可知：跳 C Q,“，k=0,1,2,有f（x）=Lp（x）+fx（）,Xi,-,xp,xup+x（x x e a,6,其中“=-叼）（二 二;%+i（”）=HQ -）因此，对于/G Q,4有/=4（/）+、为 幻+1（/）+/比0 也 1,，如 4+1.UJ.JU在上式中取x=x

5、k得该点导数的插值计算公式P（跳）=4（跳）+/M 孙 ，叼，跳 口（跳 一&）.（1 6）i=0,i壬k该公式将一般函数/（/）的求导问题转化为多项式的求导问题，从而大大地简化了导数的计算。以下，我们介绍导数的Richarson外推算法。若函数/在力处充分可微，则由Taylor展开式有00A岔+无/2）=2=0中 仇/次2.H f（力f（x 九/2）=q2=0 由上述二式有OO/。=庐，2=1其中g=/(印/2)7(1/2),且系数出与九无关。因此，据题1.8给出的结论，导数fx)可近似地由Richarson外推公式计算：(9 o W =g(h),9 m W =4.m_1(/z/2)-。加_

6、1(回/(*-1),m=,(17)其逼近序列 gm(h)的误差为OO2=1这 里 诸 系 数 与 九 无 关。3.2机械求积公式3.2.1 基本数值积分方法定积分求值的困难性往往源于被积函数的复杂性。因此，将复杂被积函数用简单函数近似替代是构造数值积分算法的基本思想。众所周知，从几何观点来看定积分,广/(力)心即为由曲线？/=/，直线/=。、岔=b及力轴所围下窗图形面积的代数和。因此,若用直线段g=/%+(l。河(。(0,1)近似代替曲线段4=(a x 6),则可得矩形积分公式 f(x)dx (6-a)f0a+(1 0)b,0 e 0,1.(18)J a特别，当取e=o,*,i时，我们分别称之

7、为右矩公式，中矩公式及左矩公式。若以过点4(Q,/(Q)、6(bJ(b)的直线段y /(Q)+a),x e a,6 b a近似代替曲线段 =/(/),则得梯形公式pb 1/fx)dx+/(6).(19)考虑过点4(Q J(Q)、。(空J(空)、右在/)的抛物线段y=px+Q X+r,x G a,6,其中p,q,r由方程组pa2+q a+r =/(a),+q(竽)+7 =/(%，pl?qb+r=/(6)确定，用该妈物线段近似代替曲线y=f(x)(a x 1 1/(1)+A2/(2)+4 3 A 3),Jo并由该公式计算定积分:=若 野 心，指出其绝对误差。解 据(22),要公式具3次代数精度，则

8、必有4 +4 +4 +4 =3,Ai+2 A2+34 =I)4 +4 4 +9 4 =9,4 +842+2 7 4 =斗.解之得4 =看4 =A2=1，A3=I，由此即得求积公式r3 3/f d x x-/(0)+3/(1)+3/(2)+/(3),JO 0且当将/(0=/4代入上式时，其不能精确成立，故所得公式具3次代数精度。应用该公式得8-1-F4 9 36e2=4.02404510389840.又其积分的精确值为x exp(jj)(x+I)2e3-1431+x o则上述数值积分的绝对误差为|1-7|=0.00266087310148.3.2.3插值求积法插值求积法是利用被积函数/（力）的插

9、值多项式计算定积分的方法，其根据被积函数在某些节点处的函数值构造一个插值多项式PN,然后用PNQ）近似代替/1（I），而获得积分逼近公式f(x)dx I Pxdx.a这样获得的求积公式称为插值型求积公式。对 于 积 分/(外近，在 区 间a/上 任 取N+1个互异点10,3,xN,构造/(i)的带余项的Lagrange插值公式(23)其中注/(N+1)思G N+1 0)=/-),RN(/,6)=(N+)!=N+I(琉 S e (a,b).将(23)代入积分/(乃冽中得rbNI fxdx=Anfxn)+RN,a n=Q(24)其中G N+1 0)a(1 一 力72屹+(期)dx,RN(/)=x)

10、dx.(25)a在(24)略去余项AN(/)即得插值型求积公式fb N/f(x)d x x 4/出).(26)Ja 2 0若max|/(N+I)I=M v+i，则其余项R v(/)的有如下估计式X (Z,6|RV(/)|(；：+；/QN+13血.(27)例 3.3 取节点g=例4,n=0,1,2,3,4,试利用4 次插值型求积公式计算定积分 I =sin x2)d x,Jo并估计其误差。解 由(25)可计算出求积公式的系数4)=7/90,Ai=1 6/45,A2=2/1 5,A3=1 6/45,A4=7/90.因此，利用4次插值型求积公式有,1 4/sin(x2)dx x V An s in(

11、)=0.3 1 0 261 423 65 3 5 3 74.J。n=0又颂：=m a xo；0,ld(s in z2)dx5m a xx 0,l|3 2x cos(x2)+1 60 x3s in(a?2)1 20 x cos(x2)=3 2r r5 cos (a:2)+1 60 r r3 s in (a?2)1 20 a;cos (a：2)t=1 8.70 89e+0 0 1,则据（27）,其误差估计为四（切/(a:n+i)-/(4)(9n),en e xn,xn+1.求和得b haN l N l I,5 N lE/3%)+2 /升/-/仇).:/ZooUn=0 n=l n=0故有复合Simp

12、son公式Lr b f 3 d x 6hN-l N-l-/(Q)+4 /(7 n+p +2/(g)+/(b),(3 8)n=0 n=l _其余项用2(，)=一 丽 与/仇 人 一 菽/小=一-丽 l/(3)(bH(3)(0).n=0 J a(39)出于计算机编程方面的考虑，我们记逐九=工ni 1 2九 1 =xn_i n=)2)一,)N)而将式(3 8)写成rb h/fx)dx x -N Q)-/的+2 2/(通九 _1)+/(遏九)卜(40)71=1 )其中彳 b-dh=kxn=a+nh,n =1,2,2 N.例3.4取步长/i=I，分别用公式(36)及(40)计算积分/=./o1 smx-

13、ax.x解 记/(乃=乎，并取/(o)=L n等=1,则由公式(36)得 x 0/11 -167/(0)+2 磴+/7 1=1由公式(40)得/J，。)E4 -2/(n=l L=0.9456908635827014.2n-1 n=0.9460833108884721.将上述两个近似值与/的准确值0.9460831 一相比较可知：复合梯形公式仅有2位有效数字，而Simpson公式有6位有效数字。故在计算量基本相同的情况下，后者精度高于前者。3.4变步长求积法梯形公式与Simpson公式的复合使其精度获得改善，但两者均属于定步长公式，若要求达到某个计算精度，其步长的选择则成为一件困难的事情。为此，

14、本节介绍一种在计算机上自动选择步长的变步长方法，同时也将涉及其加速收敛技巧。3.4.1 变步长梯形求积法变步长梯形求积法即是在求积区间上通过步长逐次减半使用复合梯形公式计算定积分的方法。当将求积区间Q,4等分k次时，其复合梯形公式的求积值为_ b C L c Xn(Tl(b d)n/7 Tk=+)+/(与n=l )(41)其满足递推关系211b-Q 2 7 2 -1.、。=尹-1+2 1 Q -a)72=1 Lk=1,2,(42)该式称为变步长梯形求积公式。由(3 7)有(43)/小 1 片 徽尸尸,这里T厂 7 b aI=以*d%h=J d据(43)可得其事后误差估计/Tk+i g(Tk+i

15、 TkY(44)鉴于上式是一个近似估计，因此我们可保守地以|n+1 -Tk作为当前步近似值+1 的误差。若预定精度为：|/-Tk+1 s,则在实际运行变步长梯形求积公式时，不等式|+1 -|=tol&k=1000k=k+l;for n=l:2 7 k-l)g(n)=f(a+(b a)*(2*n 1 )/2Ak);ends=sum(g);tO=t 1 ;tl=t0/2+s*(b-a)/2 k;endtlk例3.5应用变步长梯形求积法计算定积分/exp并要求其计算精度满足：Tk 10-8。解 取a=0,6=3,tol=1 0弋/=，5 八）（力+1产运行算法3.1,经 15次迭代后获得满足精度要求

16、：7V 1 I +髡。-喇,2,72=1/7 4 m 7（+1）T）Tm=手Li-I Z =0,1，；m =1,2,.（45）其计算步骤可按如下箭头所示方向依次进行：I 需 一 方 一 球雪I ）一 7一 球 一碑）若要求逼近精确积分值/的计算精度满足田，)-“心 则可利用终止准则I础)-T巴|=to 11=1+1;h=h/2;for n=1 :mg(n)=f(a+h*(2*n 1);ends=sum(g);t(l+l,l)=t(l ,l)/2+s*h;m=2*m;for k=1:1t(1+1 ,k+l)=(4k*t(1+1 ,k)-t(1,k)/(4 k-l);ende rr=a b s(t

17、(l+l,l +l)t(l,1);endt(l+l,l+l)例 3.6 用Romberg算法计算积分f3 x exp(re)7I=/7-Jo 0 +1)2并要求其计算精度满足：|T俨-康J 10一8。解 取,”、xexp(x)Q=0,6=3,tol 10 f (x)-.(x+I)2据 Romberg算 法 3.2,可计算得满足精度：|我)外%|0n=oIT知公式(50)对2N+2次多项式不精确成立，故(50)为Gauss型求积公式，即3,为，v 为 Gauss 点。由定理3.2可知,若能找到满足(48)的由+1次多项式G N+13)，则公式的Gauss点就确定了，从而确定了一个Gauss型求积

18、公式，为解决这一问题，我们引入Legendre多项式及其相关结论。定义3.3 一个仅以区间 1,1 上的Gauss点 出(=0,1,2,，N)为零点的N+1 次多项式称为Legendre多项式。定理3.4 首项系数为1 的 Legendre多项式可唯一地表示为GN+IQ)(N+1)!1)N+1(2N+2)!dxN+1N=0,1,.证 明 考虑2N+2 次多项式/,、N+ixdxdx dx,7M 力其满足=3N+I(C),)(1)=0,i=0,1,2,N.设M l)为任意N 次多项式，则由(51)得(51)/u(力)3N+i(6)d/=/vxuN+1x)dxJ-l J-l=o(/)(N)(7)/

19、vNx)vx)dx=0(N)_/uNxx)vr,(x)dx=(继续逐次分部积分)Nrl=2(-1)%(，T +/ux)v+Yx)dx=0J-lN=(T)。/N T).i=0而据定理3.3 知/v(x)UN+i(x)dx=0,N则 5 2(1)%(N T)=0.2=0由此及。(%)的任意性有u(NT)(i)=o,i=0 1 2 ，N.(52)(51)和(52)表明/=1均为(乃的N+1重零点，故“=C(1尸+。为待定常数.而仞V+1 Q)的首项系数为1，则由(51)的第一式有y (N +1)!,(N +I)!d+1(，1 严 1(2N +2)!A (2N +2)!dxN ,例 3.7 试构造Ga

20、u s s型求积公式/f(x)dx-Aof(xo)+4 (g)+A2f(x2并 由 此 计 算 姬 磊 山，要求精确到I C T4.解通过求三次Legendre多项式小侬)=/3-5的零点获Gauss点其公式的系数由(25)第一式得.5.8.4)=个 人1 =个 人2y y故所求公式为,j 八 I，-冉+?(0)+/(y i(53)据(46)作变换7=2(力 9，且由公式(53)得(1+力 产Viy/x+1 7-axV2 0，则由(54)看归|2 m ax|sn|.0nN因此Gauss公式的误差是可控的，即具较好的稳定性能。此外，若 危)(-1,1),则Gauss求积公式必收敛，且有N1lim 5 Anfxn)=/f(x)dx.Noo z I in=0 J T