整理常微分方程数值解

1、第八章常微分方程数值解姓名学号班级习题主要考察点：欧拉方法的构造，单步法的收敛性和稳定性的讨论，线性多步法中亚当姆 斯方法的构造和讨论。1用改进的欧拉公式，求以下微分方程” , 2xy =y -一“ y xo,iMO） =1（改进的尤拉公式的应用），有史 - 2z = -2x dx的数值解（取步长 h = 0.2），并与精确解作比较。2解：原方程可转化为yyy2 -2x，令z =2解此一阶线性微分方程，可得y = 2x 1。利用以下公式”2人yp = yi +0.2 W )yi2x仏=yi +0.2 5 4 (i =0,1,2, 3, 4)ypyw =扣卩+yj 求在节点x0.2 i ( i

2、= 1, 2,3, 4, 5 )处的数值解yi，其中，初值为x0 = 0 , y0 = 1。MATLAB程序如下：x(1)=0;%初值节点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);%预报值yc=y(i)+0.2*(yp-2*x(i)/yp);%校正值y(i+1)=(yp+yc)/2;%改进值x(i+1)=x(i)+0.2;%节点值yy(i+1)=sqrt(2*x(i+1)+1);%精确解fprintf(x(%d)=%f,y(

3、%d)=%f,yy(%d)=%fn,i+1,x(i+1),i+1,y(i+1),i+1,yy(i+1);end程序运行的结果如下：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.814224, yy(5)=1.612452x(6)=1.000000,

4、y(6)=2.027550, yy(6)=1.7320512用四阶龙格一库塔法求解初值问题丿y y ，取h=0.2，求x=0.2,0.4时的数值解7(00要求写出由h, Xn, yn直接计算yn 1的迭代公式，计算过程保留3位小数。(龙格一库塔方法的应用)解：四阶龙格-库塔经典公式为hyn 1 = Yn(k1 2k22ks k4)6k1 =f(Xn,yn)k21 1= f(Xn 2h,yn 2hk1)k3f(Xn1h, yn 2hk2)2 2k4 = f (Xn h, yn hk3)由于f（x,y） =1 y，在各点的斜率预报值分别为: k1 = 1 - ynk2 =1 (yn ；k1)=1

5、- yn一 yn) =(1 - yn)(1 -；)2 2 2hhhhhk3=1-(ynk2)=1-yn(1-丫.)(1)=(1-丫.)1(1)22222hhhhk4 =1 - (yn hk3)=1 - yn -h(1 -yn)1-;(1 -;) =(1 -yn)1 -；)2222四阶经典公式可改写成以下直接的形式：yn1 二 yn h(1 -yn)(6-3h h26在x =论二0.2处，有y1=0 竺仆 -0)(6-3 0.2(0.2)26(0.2)34)=0.1813在x = x2 = 0.4处，有3y2 =0.1813 一(1 -0.1813)(6-3 0.2(0.2)2 -=0.3297

6、64注：这两个近似值与精确解 y =1 - e在这两点的精确值十分接近。3用梯形方法解初值问题y+ y =01(0) =1证明其近似解为并证明当h )0时，它收敛于原初值问题的准确解y=e。解：显然，y二e是原初值问题的准确解。求解一般微分方程初值问题的梯形公式的形式为hyn 1 =yn - f(Xn,yn) f &n 1, Yn 1)对于该初值问题，其梯形公式的具体形式为hyn 1 二 yn 2(ynhhn1)，(1 尹n(1-2)yn，/ 2 h yn(尸)yn是：*2-hY.0 h严12卄一一I y-l2 + h 丿-hh2 - h、 yn1 5)yn亦即：yn二注意到：xn = 0 n

7、h 二 nh2h2 hyn =Xn彳2h E1 -.2 hXnXnXn(1 t) tXn=(1 -1) t (1 t) 2Xn从而”叫ynjimf t) t lim (1 t)2en即：当h- 0时，yn收敛于原初值问题的准确解y(Xn) =e-Xn。4对于初值问题丿y10y，证明当h.h2十j h由h0,九0可知，en4f .hy = f(x y)6设有常微分方程的初值问题丿，试用泰勒展开法，构造线性两步法数值计算、y(x。)= y。公式yn： (yn ynj) hf-0f -1 fnj)，使其具有二阶精度， 并推导其局部截断误差 主项。(局部截断误差和主项的计算) 解：假设yn =y(Xn

8、)， yn 4 = y(xn j)，利用泰勒展式，有y (Xn), 2 y (Xn) , 3.Yn* =y(XnG =y(Xn) -丫(Xn)hhh2 6fn = f (Xn”n)二 f(Xn, y(Xn)二 y (Xn)n4 = f (Xn4, YnG = f (Xn4, 丫&.4) = y(Xn4)= y (Xn) - y(Xn)hh2yn 4. =2：y(Xn) Co + - ：)y (Xn)h - ( I)y (Xn)h2 (6 (Xn)h3 -2 6 21 2 1 3又 y(Xnl) =y(Xn) y (Xn)h Jy (Xn)h2 y (Xn)h32 6欲使其具有尽可能高的局部截断

9、误差，必须2：二 1,飞1_ ：二1 ,i 二 12 2从而-,飞=7 ,2441 71于是数值计算公式为yn 1(yn ynd) h fnfn)。2 44该数值计算公式的局部截断误差的主项为1s353y(Xn1)-yn1 =(；5 (Xn)h3=三 y (Xn)h3662247已知初值问题y 二 2xy(o) =oy(0.1) -0.01取步长h =0.1，利用阿当姆斯公式yn “ = yn (3fn - fn)，求此微分方程在0，102上的数值解，求此公式的局部截断误差的首项。(阿当姆斯公式的应用)解：假设yn =y(Xn)， yn4二y(Xn)，利用泰勒展开，有yn =y(Xn)， fn

10、 =y (Xn)，仁/ 二 y (x. 4) = y ( X. ) - y (X. )hh 2 -2yn 1 二 y(Xn) y (Xn)h y (Xn)h2 -h3241讯21和3而 y(Xn 1) = y(Xn) y(Xn)h 一讨(Xn)hy(Xn)h2 61 15y(Xn1)-yn1 =(- -)y (Xn)h3(Xn)h35该阿当姆斯两步公式具有 2阶精度，其局部截断误差的主项为y (Xn)h3。12取步长h = 0.1，节点Xn =0.1 n( n = 0,1 , 2，,100)，注意到f (x, y) = 2x，其计算公式可改写为0 1yn 1 = yn(6Xn -2Xn4)=

11、Yn 0.02 n 0.012仅需取一个初值y0 =0，可实现这一公式的实际计算。其 MATLAB 下的程序如下：x0=0; %初值节点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=x1;y0=y1;end运行结果如下：x( 1)=0.10000000,y( 1)=0.01000000x( 2)=0.20000000,y( 2)=0.04000000x( 3)=0.30000000,y( 3)=0.09000000x( 4)=0.

12、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( 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

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)=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)=

14、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)=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.0

15、0000000,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)=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)

16、=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( 43)=18.49000000 x( 44)=4.40000000,y( 44)=19.36000000 x( 45)=4.50000000,y( 45)=20.25000000 x( 46)=4.60000000,y( 46)=21.16000

17、000 x( 47)=4.70000000,y( 47)=22.09000000 x( 48)=4.80000000,y( 48)=23.04000000 x( 49)=4.90000000,y( 49)=24.01000000 x( 50)=5.00000000,y( 50)=25.00000000 x( 51)=5.10000000,y( 51)=26.01000000 x( 52)=5.20000000,y( 52)=27.04000000 x( 53)=5.30000000,y( 53)=28.09000000 x( 54)=5.40000000,y( 54)=29.16000000

18、x( 55)=5.50000000,y( 55)=30.25000000 x( 56)=5.60000000,y( 56)=31.36000000 x( 57)=5.70000000,y( 57)=32.49000000 x( 58)=5.80000000,y( 58)=33.64000000 x( 59)=5.90000000,y( 59)=34.81000000 x( 60)=6.00000000,y( 60)=36.00000000 x( 61)=6.10000000,y( 61)=37.21000000 x( 62)=6.20000000,y( 62)=38.44000000 x( 6

19、3)=6.30000000,y( 63)=39.69000000 x( 64)=6.40000000,y( 64)=40.96000000 x( 65)=6.50000000,y( 65)=42.25000000 x( 66)=6.60000000,y( 66)=43.56000000 x( 67)=6.70000000,y( 67)=44.89000000 x( 68)=6.80000000,y( 68)=46.24000000 x( 69)=6.90000000,y( 69)=47.61000000 x( 70)=7.00000000,y( 70)=49.00000000 x( 71)=7

20、.10000000,y( 71)=50.41000000 x( 72)=7.20000000,y( 72)=51.84000000 x( 73)=7.30000000,y( 73)=53.29000000 x( 74)=7.40000000,y( 74)=54.76000000 x( 75)=7.50000000,y( 75)=56.25000000 x( 76)=7.60000000,y( 76)=57.76000000 x( 77)=7.70000000,y( 77)=59.29000000 x( 78)=7.80000000,y( 78)=60.84000000 x( 79)=7.90000000,y( 79)=62.41000000 x( 80)=8.00000000,y( 80)=64.00000000 x( 81)=8.10000000,y( 81)=65.61000000 x( 82)=8.200000