高等工程數(shù)學(xué)第五講,常微分方程數(shù)值解法

1、1 主要內(nèi)容主要內(nèi)容10.1 Euler方法方法10.2 改進(jìn)改進(jìn)Euler法法10.3 收斂性與穩(wěn)定性收斂性與穩(wěn)定性10.4龍格龍格庫塔庫塔(Runge-Kutta)法法第五講第五講 常微分方程數(shù)值解法常微分方程數(shù)值解法, 許多實(shí)際問題的數(shù)學(xué)模型是微分方程或微分方程組的定解問題,如物體運(yùn)動,電路振蕩,化學(xué)反應(yīng)及生物種群的變化等. 能用解析方法求出精確解的微分方程為數(shù)不多,有的甚至方程即使有解析解,也可能由于解得表達(dá)式非常復(fù)雜而不易計算.因此十分有必要研究微分方程數(shù)值解法.2 目的：目的：用數(shù)值解逼近解析解用數(shù)值解逼近解析解我們考慮一階常微分方程初值問題0( ,),()dyfx yaxbdxy

2、 ay（10.1）( ,)( ,),(10.2)f x yf x yL yy設(shè)f(x,y)連續(xù)，且關(guān)于y滿足Lipschitz條件，存在常數(shù)L,使得由常微分方程理論知，初值問題（10.1）有解且唯一。 3012Naxxxxb 處得近似值 的方法， 稱為問題（9-1）的數(shù)值解。 為步長。 12,Ny yy12,Nyyy1nnnhxx常用幾種方法常用幾種方法:用差商近似導(dǎo)數(shù)；用數(shù)值積用差商近似導(dǎo)數(shù)；用數(shù)值積分方法；用分方法；用Taylor多項(xiàng)式近似。多項(xiàng)式近似。數(shù)值解法-求問題（9-1）的解y(x)在若干點(diǎn)現(xiàn)將微分方程離散化現(xiàn)將微分方程離散化-建立數(shù)值解法建立數(shù)值解法4(1) 用差商近似導(dǎo)數(shù)用差商

3、近似導(dǎo)數(shù)11()()()()()(,(),(0,1,)代入初值問題中的方程nnnnnnny xy xyxhy xy xf xy xnh可由初值y0逐次算出 .稱離散化問題（9-3）為差分方程初值問題。12,yy1()()(,()nnnny xy xhfxy x10(,),(0,1,),(10.3)( )nnnnyyhf xynyy a1(,),(0,1,)nnnnyyhf xyn5Euler方法-用差分方程初值問題的解近似微分方程初值問題（9-1）。10(,),(0,1,),(10.3)( )nnnnyyhf xynyy a可由初值y0逐次算出.12,yy用Euler方法求下面初值問題的數(shù)值解

4、,01(0)0yxyxy 10.1 Euler方法方法例例1、6解：解：Eulern+1nnn公式的具體形式為y=y +h(x -y ),(n=0,1,)1,xyxe易 知 精 確 解 為若取h=0.1,N=10,可由初值y0=0逐次算出.見下表()()00.0000000.0000000.0000000.10.0000000.0048370.0048370.20.0100000.0187310.0087310.30.0290000.0408180.0118180.90.2874200.3065700.0191501.00.3486780.3678790.019201nnnnnxyy xy x

5、y72,01(0)0.0.1.xyh 22例 用Euler方法求解初值問題y =1+x +y取解：22i記f(x,y)=1+x +y ,x =ih=0.1h,0i10計算可得012()0.10000,()0.20200,()0.31008,yhyhyh221002221122322y(0.1)y1+x +yy(0.2)y1+x +yy(0.3)y1+x +y86789()0.90759,()1.13896,()1.43268,()1.81894,yhyhyhyh227662287722988221099y(0.7)y1+x +yy(0.8)y1+x +yy(0.9)y1+x +yy(1.0)y

6、1+x +y345()0.42870,()0.56307,()0.71987,yhyhyh224332254422655y(0.4)y1+x +yy(0.5)y1+x +yy(0.6)y1+x +y9二二,隱式歐拉公式隱式歐拉公式三三,歐拉中點(diǎn)公式歐拉中點(diǎn)公式)5.10(,)()()(1i11111nnnnnnnyxhfyyxyxyxyhf11111i()()()22,(10.6)nnnnnnnfy xy xyxhyyhfxy10若用梯形公式代替微分方程（10.3）右端積分，有11(),()nnnnyy xyy x得計算公式-梯形公式。1111(,)(,)2nnnnnnyyfxyfxy梯形公式

7、的局部截斷誤差31111()( ),12nnnnnhRy xyyxx 梯形公式為二階方法。10.2 改進(jìn)的改進(jìn)的Euler方法方法一一, ,梯形公式梯形公式111( , ( ) (, ()(, ()2nnxnnnnxhf x y x dxf xy xf xy x11 由于梯形公式仍然是隱式方法，一般用迭代法由于梯形公式仍然是隱式方法，一般用迭代法求解，若每一步只代一次，相當(dāng)將求解，若每一步只代一次，相當(dāng)將EulerEuler公式與梯公式與梯形公式結(jié)合使用。形公式結(jié)合使用。1111(,)(10.9) (,)(,)2nnnnnnnnnnyyhf xyhyyf xyf xy預(yù)測校正為上機(jī)編程方便，改

8、寫為為上機(jī)編程方便，改寫為1(,)(,)(10.10)1()2pnnnqnnpnpqyyhfxyyyhfxhyyyy二二, 改進(jìn)改進(jìn)Euler法法(梯形法梯形法)12改進(jìn)Euler法求解下面初值問題的數(shù)值解。（取h=0.1,N=10）01(0)0,yxyxy 用公式（9-7） (,)0.1()0.1(0.1)pnnnnnnqnnpyyhf xyyxyyyxy若取h=0.1,N=10,可由初值y0=0逐次算出.見下表11()2npqyyy例例2、解：解：13()()00.0000000.0000000.0000000.10.0050000.0048370.0001630.20.0190250.0

9、187310.0002940.30.0412180.0408180.0004000.90.3072280.3065700.0006581.00.3685410.3678790.000662nnnnnxyy xy xy 改進(jìn)Euler法優(yōu)于Euler法.(二階方法)1410.3 收斂性與穩(wěn)定性收斂性與穩(wěn)定性11()()(, ()nnnnnRy xy xhf xy xn+1ny(x)=y(x +h)上式截斷誤差右邊的第一項(xiàng)由Taylor展開式得2( )nn1=y(x )+hy (x )+2h y21( ),nnn1y(x )+hf(x ,y(x )+2nnh yxx截斷誤差截斷誤差由Euler公式

10、收斂性收斂性:1521( ),n+1n+1n+11R=y(x)-y=2nnh yxx11() ()(,()nnnnnRy xy xhf xy x22( ),(0,1,2,)22hhyMn212nhRM由Euler公式111()nnney xy（2）截斷整體誤差）截斷整體誤差（不含舍入誤差）16111()nnney xy配項(xiàng)111111() ()(, () (,) () (, ()(,)nnnnnnnnnnnnnnnnnny xyyyRy xhf xy xyhf xyRy xyh f xy xf xy1()(,),nnnny xhf xy記y有111()()(1)()(1)nnnnnnnnnnR

11、y xyhL y xyRhL y xyRhL e172221(1)2(1)(1)22nnhMhL ehhMhLMhL e20211(1)2(1)1(1)12(1) 12nkknnhMhLh MhLhMhLhLL18NEuler方法的整體截斷誤差與h同階(h0,e0).p+1if數(shù)值方法的局部截斷誤差為O(h),稱其是p階的。Euler方法是一階的(精度不高）()112L b anhMeeL1(1)(1)(1)b anhnhbahLhL由結(jié)合前面有23()()()121(1)hLb ab ahLhLhhhLhLehLhLehLe而！3！19穩(wěn)定性穩(wěn)定性()yy nmm用數(shù)值方法求解為復(fù)數(shù) ，當(dāng)h

12、步長一定，由計算節(jié)點(diǎn)值y 時產(chǎn)生誤差 所引起后面節(jié)點(diǎn)值y的誤差的絕對值均不超過，則稱對h和 是絕對穩(wěn)定的。絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域- 使得方法絕對穩(wěn)定的h和 的全體。絕對穩(wěn)定區(qū)域包含左半平面絕對穩(wěn)定區(qū)域包含左半平面-該方法稱為該方法稱為A穩(wěn)定穩(wěn)定。 模型) 1( ,yy2011,計算y 產(chǎn)生 ，y產(chǎn)生nnn（1）Euler方法穩(wěn)定性方法穩(wěn)定性 由Euler方法知，計算公式為1(1)nnyhy1(1)nh 111記yy，yynnnnn111(1)(1)(1) yyyynnnnnhhh2110.3 收斂性與收斂性與穩(wěn)定性穩(wěn)定性()yy nmm用數(shù)值方法求解為復(fù)數(shù) ，當(dāng)h步長一定，由計算節(jié)點(diǎn)值y 時

13、產(chǎn)生誤差 所引起后面節(jié)點(diǎn)值y的誤差的絕對值均不超過，則稱對h和 是絕對穩(wěn)定的。絕對穩(wěn)定區(qū)域絕對穩(wěn)定區(qū)域- 使得方法絕對穩(wěn)定的h和 的全體。絕對穩(wěn)定區(qū)域包含左半平面絕對穩(wěn)定區(qū)域包含左半平面-該方法稱為該方法稱為A穩(wěn)定穩(wěn)定。 收斂性收斂性22則Euler方法的絕對穩(wěn)定區(qū)域?yàn)?1h特點(diǎn)：特點(diǎn)： 在 h平 面 上 以 -1為 中 心 的 單 位 圓 。1）當(dāng)h 在單位圓內(nèi)，Euler方法是穩(wěn)定的；2）當(dāng)h 在單位圓外，Euler方法是不穩(wěn)定的；oRe()hIm()h-1-223（2）梯形公式的穩(wěn)定性）梯形公式的穩(wěn)定性 由梯形公式知 11112()212nnnnnhhyyyyyh類似于前面梯形公式的方法

14、121,12hh241122hh絕對穩(wěn)定區(qū)域?yàn)樽蟀肫矫?。Re()0hoRe()hIm()h253-2 (Runge-Kutta)法的構(gòu)造法的構(gòu)造一般有 11111(,)(,)(2,3, )pnniiinniininijjjyyhc KKf xyKf xa h yhb Kip確定相應(yīng)參數(shù)的原則：確定相應(yīng)參數(shù)的原則：近似公式在（Xn,Yn) 處的泰勒展式與y(x) 在Xn處的泰勒展式的前面的項(xiàng)盡可能重合，以提高精度。2611122122211()(,)(,)nnnnnnyyh c Kc KKf xyKf xa h yb hK當(dāng)P=2時的模型33()0 ()O hhn+1n1nn2n2n21nnn1

15、nn2nn2xnn21ynnnn2n12nn22 xnn21ynnnny= y + hcf(x ,y )+ c f(x + a h,y + hb f(x ,y )= y + h cf(x ,y )+ cf(x ,y )+ a hf(x ,y )+ b hf(x ,y )f(x ,y )= y +(c + c )f(x ,y )h + ca f(x ,y )+ b hf(x ,y )f(x ,y )h上式在（Xn,Yn)處展開272312323()()()()()nnnnnnnxnnynnnn)+y (x )h+2+f(x ,y )h+f(x ,y )2+f(x ,y )h+f(x ,y )+f

16、(x ,y )f(x ,y )2nnnnny xy xy xhO hhyO hhyO h比較系數(shù)有：122 22 21111 ,22ccc ac b3有無窮多組解，其解回代得模型公式，且局部誤差均為O（h ）是二階方法。28 11212111()22(,)(,)nnnnnnyyhKKKf xyKf xh yhK221,1 ;ab121特別?。海?）c =c =2改進(jìn)改進(jìn)Euler公式公式 29（二階）中點(diǎn)公式（二階）中點(diǎn)公式 12121(,)(,)22nnnnnnyyhKKf xyhhKf xyK221(2),.ab121c =0,c =1230（3 3）三階三階(Runge-Kutta) 公

17、式公式 11231213121(4)6(,)11(,)22(,2)nnnnnnnnyyh KKKKf xyKf xh yhKKf xh yhKhK311123412132431(22)6(,)11(,)2211(,)22(,)nnnnnnnnnnyyh KKKKKf xyKf xh yhKKf xh yhKKf xh yhK（4 4）四階經(jīng)典四階經(jīng)典(Runge-Kutta) 公式公式算法算法9-2（見教材（見教材212） 32用四階經(jīng)典(Runge-Kutta)法求解下面初值問題的數(shù)值解。（取h=0.2,N=5）01(0)0,yxyxy 依題意四階經(jīng)典(Runge-Kutta)公式應(yīng)為112

18、3412132431(22)3011()0.91()0.12211()0.91()0.0922()0.818()0.182nnnnnnnnnnnnnnnnyyKKKKKxyKxhyhKxyKxhyhKxyKxhyhKxy解：解：例例3、33由初值y0=0逐次算出下表()()00.0000000.0000000.0000000.20.0187330.0187310.0000020.40.0703240.0703200.0000040.60.1488170.148812.00000050.80.2493350.2493290.0000061.00.3678660.3678790.000007nnn

19、nnxyy xy xy特點(diǎn)：特點(diǎn)：用四階經(jīng)典(Runge-Kutta)法計算的結(jié)果，其精度高于改進(jìn)Euler法.34例例1.2(0)10,1(0.2).xyyyyh用經(jīng)典的四階龍格 庫塔方法求解初值問題在上的數(shù)值解 解解:利用經(jīng)典的四階龍格-庫塔公式有351211322433112342(,)2 (0 .1)0 .10 .12 (0 .1)0 .10 .12 (0 .2 )0 .20 .20 .1(22)3nnnnnnnnnnnnnnnnxkfxyyyxkykykxkykykxkykykyykkkk36計算結(jié)果如下表。nx 0 0.2 0.4 0.6 0.8 1.0ny)(nxy 1 1.18

20、323 1.34167 1.48324 1.61251 1.73205 1 1.18322 1.34164 1.48328 1.61245 1.73205 該結(jié)果具有四位有效數(shù)字，和本章第一節(jié)例1比，可見四階龍格-庫塔方法的精度高于Euler公式和改進(jìn)Euler（預(yù)估-校正）公式。37%a=0;y0=1;h=0.2;n=0;x_rec(1)=a;y_rec(1)=y0;fprintf(x%2.0f=%10.6f y%2.0f=%10.6fn,n,x_rec(1),n,y_rec(1)x=a;y=y0;while x1 n=n+1 k1=y-2*x/y; k2=(y+h*k1/2)-(2*(x+

21、h/2)/(y+h*k1/2); k3=(y+h*k2/2)-(2*(x+h/2)/(y+h*k2/2); k4=(y+h*k3)-(2*(x+h)/(y+h*k3); y=y+h*(k1+2*k2+2*k3+k4)/6; x=a+n*h; y_rec(n+1)=y;x_rec(n+1)=x;fprintf(x%2.0f=%10.6f y%2.0f=%10.6fn,n,x_rec(n+1),n,y_rec(n+1) endplot(x_rec,y_rec)38x 0= 0.000000 y 0= 1.000000n = 1x 1= 0.200000 y 1= 1.183229n = 2x 2=

22、 0.400000 y 2= 1.341667n = 3x 3= 0.600000 y 3= 1.483281n = 4x 4= 0.800000 y 4= 1.612514n = 5x 5= 1.000000 y 5= 1.73214200.511.522.533.544.5500.20.40.60.811.21.41.61.8xyR-K(X,Y)00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8XYR-K:(X,Y)39%a=0;y0=2;h=0.25;n=0;x_rec(1)=a;y_rec(1)=y0;fprintf(x%2.0

23、f=%10.6f y%2.0f=%10.6fn,n,x_rec(1),n,y_rec(1)x=a;y=y0;while x5 n=n+1 k1=-x*y2; k2=-(x+h/2)*(y+h*k1/2)2; k3=-(x+h/2)*(y+h*k2/2)2; k4=-(x+h)*(y+h*k3)2; y=y+h*(k1+2*k2+2*k3+k4)/6; x=a+n*h; y_rec(n+1)=y;x_rec(n+1)=x;fprintf(x%2.0f=%10.6 y%2.0f=%10.6fn,n,x_rec(n+1),n,y_rec(n+1) end例例1.2,(05)(0)20,5(0.25)

24、.yxyxyh 用經(jīng)典的四階龍格庫塔方法求解初值問題在上的數(shù)值解 40plot(x_rec,y_rec)x 0= 0.000000 y 0= 2.000000n = 1x 1= 0.250000 y 1= 1.882308n = 2x 2= 0.500000 y 2= 1.599896n = 3x 3= 0.750000 y 3= 1.279948n = 4x 4= 1.000000 y 4= 1.000027n = 5x 5= 1.250000 y 5= 0.780556n = 6x 6= 1.500000 y 6= 0.615459n = 7x 7= 1.750000 y 7= 0.492

25、374n = 8x 8= 2.000000 y 8= 0.400054n = 9x 9= 2.250000 y 9= 0.329940n = 10 x10= 2.500000 y10= 0.275895n = 11x11= 2.750000 y11= 0.233602n = 12x12= 3.000000 y12= 0.200020n = 13x13= 3.250000 y13= 0.172989n = 14x14= 3.500000 y14= 0.150956n = 15x15= 3.750000 y15= 0.132790n = 16x16= 4.000000 y16= 0.117655n

26、 = 17x17= 4.250000 y17= 0.104924n = 18x18= 4.500000 y18= 0.094123n = 19x19= 4.750000 y19= 0.084885n = 20 x20= 5.000000 y20= 0.0769274100.511.522.533.544.5500.20.40.60.811.21.41.61.82xyR-K方 法(x,y)42主要內(nèi)容主要內(nèi)容10.1 Euler方法方法10.2 改進(jìn)改進(jìn)Euler法法10.3 收斂性與穩(wěn)定性收斂性與穩(wěn)定性10.4龍格龍格庫塔庫塔(Runge-Kutta)法法小結(jié)本講內(nèi)容小結(jié)本講內(nèi)容思考與練習(xí)見教

27、材思考與練習(xí)見教材PP:22043思考與練習(xí)012(0)1xyyy在 ，上以h=0.1為步長，分別用Euler法和改進(jìn)的Euler法求如下初值問題的數(shù)值解。dy其中dx解：（1）該方程是伯努利方程，其解為y= 1+2x2nnxyn+1nnEuler公式的具體形式為y=y +h(y -)44122 ()nnnnkkxkyxhkkyk hn+1n121n2n111y=y +h()22=h(y -)=h(y + h-)改進(jìn)的Euler公式的具體形式為4500（2）取步長h=0.1,x =0,y =1,計算有如下結(jié)果()()0011110.11.11.0959091.09544520.21.1918181.1840971.18321630.31.2774381.2662011.26491190.91.7177791.6781661.6733201011.7847711.7378671.732051nnnnxyEuleryEuler改進(jìn)準(zhǔn)確值近似值與準(zhǔn)確值比較，Euler法