[整理]常微分方程數(shù)值解

1、第八章常微分方程數(shù)值解姓名 學(xué)號(hào) 班級(jí)習(xí)題主要考察點(diǎn)：歐拉方法的構(gòu)造，單步法的收斂性和穩(wěn)定性的討論，線性多步法中亞當(dāng)姆 斯方法的構(gòu)造和討論。1用改進(jìn)的歐拉公式，求以下微分方程,2xy =y -一，yxwo,ij（0） =1（改進(jìn)的尤拉公式的應(yīng)用）.dz,有一- 2z = 2x dx的數(shù)值解（取步長(zhǎng) h =0.2）,并與精確解作比較。2解：原方程可轉(zhuǎn)化為yy'= y2 _2x,令z =y-2解此一階線性微分方程，可得y = J2x +1。利用以下公式'2為、yp = y +0.2 (yi ) y2xUc =yi +0.2 (yp ) (i =0,1,2, 3, 4) ypyf =

2、；(yp +丫。)求在節(jié)點(diǎn)xi = 0.2，i (i =1,2,3,4,5)處的數(shù)值解yi,其中，初值為x0=0,y0=1。MATLAB程序如下：x(1)=0;%初值節(jié)點(diǎn)y(1)=1;涮值fprintf( 'x(%d)=%f,y(%d)=%f,yy(%d)=%fn',1,x(1),1,y(1),1,y(1);for i=1:5yp=y(i)+0.2*(y(i)-2*x(i)/y(i);%預(yù)報(bào)值yc=y(i)+0.2*(yp-2*x(i)/yp);% 正值y(i+1)=(yp+yc)/2;%進(jìn)值x(i+1)=x(i)+0.2;初點(diǎn)值yy(i+1)=sqrt(2*x(i+1)+1)

3、;研確解fprintf( 'x(%d)=%f,y(%d)=%f,yy(%d)=%fn',i+1,x(i+1),i+1,y(i+1),i+1,yy(i+1);end程序運(yùn)行的結(jié)果如下：x(1)=0.000000, y(1)=1.000000, yy(1)=1.000000x(2)=0.200000, y(2)=1.220000, yy(2)=1.183216x(3)=0.400000, y(3)=1.420452, yy(3)=1.341641x(4)=0.600000, y(4)=1.615113, yy(4)=1.483240x(5)=0.800000, y(5)=1.814

4、224, yy(5)=1.612452x(6)=1.000000, y(6)=2.027550, yy(6)=1.732051'yy 12用四階龍格庫(kù)塔法求解初值問(wèn)題y y ，取h=0.2，求x = 0.2,0.4時(shí)的數(shù)值解7(0)=0要求寫(xiě)出由h,xn,yn直接計(jì)算yn中的迭代公式，計(jì)算過(guò)程保留3位小數(shù)。(龍格庫(kù)塔方法的應(yīng)用)解：四階龍格-庫(kù)塔經(jīng)典公式為hyn 1 = yn(ki2k2 2k3 kJ6k2k3一1.1、f (xn h, yn -hki) 22一1 .1、f (x- h, yn - hk2)22k4f (xn h, yn hk3)由于f (x, y) =1 - y ,在

5、各點(diǎn)的斜率預(yù)報(bào)值分別為: k1 =1 - ynk2 =1 - (yn hk1)=1 - yn -h(1 - yn) =(1 - yn )(1-7)222hhhhhk3=1-(ynk2)=1-yn 。7口)(1)=(1-丫口)1(1)22222hhhhk4 =1 -(ynhk3) =1 -yn -h(1 -yn)1-(1 -) =(1-yn)1-h(1-(1 -)2222四階經(jīng)典公式可改寫(xiě)成以下直接的形式：h2 hyn 1 = yn(1 - yn)(6 -3h h )64在x = x1 =0.2處，有3y1 =0 02(1 -0)(6 -3 0.2 (0.2)2 -(-) =0.181364在x

6、 =x2 = 0.4處，有30.22 (0.2)y2 =0.1813 (1 -0.1813)(6-3 0.2 (0.2) -(-) -0.329764注：這兩個(gè)近似值與精確解 y =1 -e在這兩點(diǎn)的精確值十分接近。3用梯形方法解初值問(wèn)題y+ y =0：y(0) =1證明其近似解為并證明當(dāng)hT 0時(shí)，它收斂于原初值問(wèn)題的準(zhǔn)確解y = e。解：顯然，y =e”是原初值問(wèn)題的準(zhǔn)確解。求解一般微分方程初值問(wèn)題的梯形公式的形式為 hyn 1 =yn - f (Xn, yn) f (Xn 1, 丫口 .)對(duì)于該初值問(wèn)題，其梯形公式的具體形式為hhh2 - hyn4 =yn +二(yn yn4) , (1

7、 +；)yn 書(shū)=(1 二)yn， yn 書(shū)=(；TZ7)yn2222 h)2 -h2 h)yn但二h仔2 - h儼<2 + h J 九，西,亦即:、. xc2h 111.注息、至 U: xn = 0+nh=nh, n =，令 t = -, = - Wh2 h h t 2xn2 2h V能e,Xn2Xnyn= 1-=(1+t) t 2 =(1+t) t (1+t) 2、2 +h JXnXn從而 lhimQyn= |imQ(1 t) t lim (1 t) 2 = e%即：當(dāng)hT 0時(shí)，yn收斂于原初值問(wèn)題的準(zhǔn)確解y(Xn)=ei。4對(duì)于初值問(wèn)題；y'= -10yJ(0)=1,證明

8、當(dāng)h <0.2時(shí)，歐拉公式絕對(duì)穩(wěn)定。(顯式和隱式歐拉公式的穩(wěn)定性討論)證明：顯式的歐拉公式為yn 1 = yn，hf (Xn, yn) = (1 - 10h)yn隱式的歐拉公式為從而 en4=(110h)en,由于 0<h<0.2, -1 <1-10h <1, en| <|en 因此，顯式歐拉公式絕對(duì)穩(wěn)定。yn 1 =ynhf (Xn 1, Yn 1) = Yn -1。卜丫口 1一ynyn + )1 10henen 1 二1 - 10h因此，隱式的歐拉公式也是絕對(duì)穩(wěn)定的。h .5證明：梯形公式y(tǒng)n =yn+f(xn,yn)+ f(xnMynQ無(wú)條件穩(wěn)定。(梯

9、形公式的穩(wěn)定性討論) 解：對(duì)于微分方程初值問(wèn)題r Vr.y = 一九 y/7 (人 A 0 )、y(0) =1其隱式的梯形公式的具體形式可表示為、，2-h、)%' yn +=(二:F)yn2 hhh、yn +=yn +二一九yn 九yn+，(1 + )yn+ =(1 222 - h、從而 en 1- =()en2 h由h>0,九>0可知,c /2 +九 hJc ° 6書(shū)|<(2)6=異,故隱式的梯形公式無(wú)條件穩(wěn)定。6設(shè)有常微分方程的初值問(wèn)題；y'= f (x,y)=丫幅)二v。，試用泰勒展開(kāi)法，構(gòu)造線性兩步法數(shù)值計(jì)算公式y(tǒng)n由=a(yn + y。)+

10、h(P0fn + P f。)，使其具有二階精度，并推導(dǎo)其局部截?cái)嗾`差主項(xiàng)。(局部截?cái)嗾`差和主項(xiàng)的計(jì)算) 解：假設(shè)yn =y(Xn), yn4=y(Xn)，利用泰勒展式，有y (Xn), 2 y (Xn) .3 .yn"y(Xn)= y(Xn) -丫(Xn)h h - 一二一 h 26fn = f(Xn,yn) = f (Xn, y(Xn) = y (Xn) fnL W yn,)二 f (XnJ, y(XnJ 二 /(心二 y(Xn) - y (Xn)h 號(hào) h2 -yn 1, =2：y(Xn) ( F J -1)y (Xn)h - (- - 、)y (Xn)h2 - (- 一)y (

11、Xn)h3 2621 .213又 y(Xn 1) = y(Xn) y (Xn)hy (Xn)hy (Xn)h2 6欲使其具有盡可能高的局部截?cái)嗾`差，必須2豆=1,久 +P _口 =1 , - ?1 =- 22從而 口 =1, P。=7, P1 =-1 244一，171于是數(shù)值計(jì)算公式為L(zhǎng)書(shū)=(yn +yn)+h( % 二3)。244該數(shù)值計(jì)算公式的局部截?cái)嗾`差的主項(xiàng)為15y(Xn1)-Yn1 =(-v)y (Xn)h3:丫6 62247已知初值問(wèn)題y y = 2xy(0) =oy(0.1) =0.01取步長(zhǎng)h =0.1 ,利用阿當(dāng)姆斯公式y(tǒng)nd1 = yn +-(3fn - fn)，求此微分方

12、程在0 ,102上的數(shù)值解，求此公式的局部截?cái)嗾`差的首項(xiàng)。(阿當(dāng)姆斯公式的應(yīng)用)解：假設(shè)yn =y(xn) , yn=y(xn)，利用泰勒展開(kāi)，有 yn = y(xn), fn = yH(xn), fn=yxn)=yxn) -y"(xn)h + y 卜)h2 -2yn 1 = y(xn) y (Xn)h 1y (Xn)h2 - y (Xn)h3241 213而 y(xn 1) = y(xn) y (xn)h y (xn)hy (xn)h2 611. Q5. Qy(xn 1) - yn 1 =(- -)y (Xn)hy (xn)h5c該阿當(dāng)姆斯兩步公式具有 2階精度，其局部截?cái)嗾`差的主

13、項(xiàng)為-5y7xn)h3o取步長(zhǎng)h = 0.1 ,節(jié)點(diǎn)Xn =0.1n ( n = 0,1 , 2，100),注意到f (X, y) = 2x ,其計(jì)算公 式可改寫(xiě)為0.1 - 一 、 yn 1 = yn k(6xn -2xn4) = yn 0.02n 0.012僅需取一個(gè)初值y0 = 0 ,可實(shí)現(xiàn)這一公式的實(shí)際計(jì)算。其 MATLAB 下的程序如下：x0=0; %初值節(jié)點(diǎn)y0=0; %初值for n=0:99y1=y0+0.02*n+0.01;x1=x0+0.1;,n+1,x1,n+1,y1);fprintf( 'x(%3d)=%10.8f,y(%3d)=%10.8fn'x0=x

14、1;y0=y1;end運(yùn)行結(jié)果如下：x( 1)=0.10000000,y( 1)=0.01000000x( 2)=0.20000000,y( 2)=0.04000000x( 3)=0.30000000,y( 3)=0.09000000x( 4)=0.40000000,y( 4)=0.16000000x( 5)=0.50000000,y( 5)=0.25000000x( 6)=0.60000000,y( 6)=0.36000000x( 7)=0.70000000,y( 7)=0.49000000x( 8)=0.80000000,y( 8)=0.64000000x( 9)=0.90000000,y

15、( 9)=0.81000000x( 10)=1.00000000,y( 10)=1.00000000x( 11)=1.10000000,y( 11)=1.21000000x( 12)=1.20000000,y( 12)=1.44000000x( 13)=1.30000000,y( 13)=1.69000000x( 14)=1.40000000,y( 14)=1.96000000x( 15)=1.50000000,y( 15)=2.25000000x( 16)=1.60000000,y( 16)=2.56000000x( 17)=1.70000000,y( 17)=2.89000000x( 18

16、)=1.80000000,y( 18)=3.24000000x( 19)=1.90000000,y( 19)=3.61000000x( 20)=2.00000000,y( 20)=4.00000000x( 21)=2.10000000,y( 21)=4.41000000x( 22)=2.20000000,y( 22)=4.84000000x( 23)=2.30000000,y( 23)=5.29000000x( 24)=2.40000000,y( 24)=5.76000000x( 25)=2.50000000,y( 25)=6.25000000x( 26)=2.60000000,y( 26)=

17、6.76000000x( 27)=2.70000000,y( 27)=7.29000000x( 28)=2.80000000,y( 28)=7.84000000x( 29)=2.90000000,y( 29)=8.41000000x( 30)=3.00000000,y( 30)=9.00000000x( 31)=3.10000000,y( 31)=9.61000000x( 32)=3.20000000,y( 32)=10.24000000x( 33)=3.30000000,y( 33)=10.89000000x( 34)=3.40000000,y( 34)=11.56000000x( 35)=

18、3.50000000,y( 35)=12.25000000x( 36)=3.60000000,y( 36)=12.96000000x( 37)=3.70000000,y( 37)=13.69000000x( 38)=3.80000000,y( 38)=14.44000000x( 39)=3.90000000,y( 39)=15.21000000x( 40)=4.00000000,y( 40)=16.00000000x( 41)=4.10000000,y( 41)=16.81000000x( 42)=4.20000000,y( 42)=17.64000000x( 43)=4.30000000,y

19、( 43)=18.49000000x( 44)=4.40000000,y( 44)=19.36000000x( 45)=4.50000000,y( 45)=20.25000000x( 46)=4.60000000,y( 46)=21.16000000x( 47)=4.70000000,y( 47)=22.09000000x( 48)=4.80000000,y( 48)=23.04000000x( 49)=4.90000000,y( 49)=24.01000000x( 50)=5.00000000,y( 50)=25.00000000x( 51)=5.10000000,y( 51)=26.010

20、00000x( 52)=5.20000000,y( 52)=27.04000000x( 53)=5.30000000,y( 53)=28.09000000x( 54)=5.40000000,y( 54)=29.16000000x( 55)=5.50000000,y( 55)=30.25000000x( 56)=5.60000000,y( 56)=31.36000000x( 57)=5.70000000,y( 57)=32.49000000x( 58)=5.80000000,y( 58)=33.64000000x( 59)=5.90000000,y( 59)=34.81000000x( 60)=

21、6.00000000,y( 60)=36.00000000x( 61)=6.10000000,y( 61)=37.21000000x( 62)=6.20000000,y( 62)=38.44000000x( 63)=6.30000000,y( 63)=39.69000000x( 64)=6.40000000,y( 64)=40.96000000x( 65)=6.50000000,y( 65)=42.25000000x( 66)=6.60000000,y( 66)=43.56000000x( 67)=6.70000000,y( 67)=44.89000000x( 68)=6.80000000,y

22、( 68)=46.24000000x( 69)=6.90000000,y( 69)=47.61000000x( 70)=7.00000000,y( 70)=49.00000000x( 71)=7.10000000,y( 71)=50.41000000x( 72)=7.20000000,y( 72)=51.84000000x( 73)=7.30000000,y( 73)=53.29000000x( 74)=7.40000000,y( 74)=54.76000000x( 75)=7.50000000,y( 75)=56.25000000x( 76)=7.60000000,y( 76)=57.76000000x( 77)=7.70000000,y( 77)=59.29000000x( 78)=7.80000000,y( 78)=60.84000000x( 79)=7.90000000,y( 79)=62.41000000x( 80)=8.000