常微分方程的数值解法(欧拉法,改进欧拉法,泰勒方法和龙格-库塔法)

1、例1:用欧拉方法与改进的欧拉方法求初值问题dydx2 x37y（o）在区间0,1上取h 0.1的数值解。解欧拉方法的计算公式为yn 1yn2 xnyo1,h 0.1使用excel表格进行运算，相应结果如下例一:欧拉法nxy精确解001110.111.00332220.21.0066671.01315930.31.0198241.02914240.41.0390541.05071850.51.0637541.07721760.61.0932111.10793270.71.1266811.14216580.81.1634431.17927490.91.2028451.2186891011.2443

2、141.259921现用 matlab 编程，程序如下x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0A2);for n=1:9x(n+1)=0.1*(n+1);y( n+1)=y( n)+0.1*2*x( n)/(3*y(门)人2);end ;结果为Columns 1 through 80.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 100.90001.0000Colu mns 1 through 81.00001.00671.01981.03911.06381.09321.

3、12671.1634Colu mns 9 through 101.20281.2443改进的欧拉方法其计算公式为yn(0)1 1ynYnyny。1,hh(2 Xn)h(7P)3 ynh2xn23 Y0.12 Xn 13(W本题的精确解为y(x)即 1 X2，可用来检验数值解的精确度，列出计算结果。使用excel表格进行运算，相应如下例一:改进的欧拉法nx预测y校正y精确解0011110.111.0033331.00332220.21.0099561.013181.01315930.31.0261691.0291711.02914240.41.0480541.0507511.05071850.5

4、1.0749041.0772521.07721760.61.1059761.1079651.10793270.71.1405491.1421941.14216580.81.1779651.1792971.17927490.91.2176461.2187061.2186891011.2591031.259931.259921现用matlab编程，程序如下xO=O;yO=1;x(1)=0.1;ya(1)=y0+0.1*2*x0/(3*y0A2);y(1)=y0+0.05*(2*x0/(3*y0A2)+2*x0/(3*yaA2);for n=1:9x( n+1)=0.1*( n+1);ya( n+1

5、)=ya( n)+0.1*2*x( n)/(3*ya(门)人2);y(n +1)=y( n)+0.05*(2*x( n)/(3*y( n)A2)+2*x( n+1)/(3*ya(门+1)人2);end ;结果为Columns 1 through 80.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 100.9000 1.0000Columns 1 through 81.00001.00991.02611.04791.07481.10591.14071.1783Colu mns 9 through 101.21831

6、.2600例2 :用泰勒方法解分别用二阶、四阶泰勒方法计算点dy2 xdx3y2y（o）1=0.1,0.2, -Xn, 1.0处的数值解，并与精确解进行比较。解：二阶泰勒方法Yn 1 Yn对于本题yxn,yn2x3?y 3y2x3?4x3?4x23?yXn,ynyn 1 yn hynh3yl2Xn h 1 4xn使用excel表格进行运算，相应结果如下nxy精确解001110.11.0033331.00332220.21.0132231.01315930.31.0292911.02914240.41.0509691.05071850.51.0775751.07721760.61.1083881

7、.10793270.71.1427041.14216580.81.1798781.17927490.91.219341.2186891011.2606011.259921现用matlab编程，程序如下xO=O;yO=1;x(1)=0.1;y(1)=y0+0.1/(3*y0A2)*(2*x0+0.1*(1-4*x0A2/(3*y0A3);for n=1:9x(n+1)=0.1*(n+1);y(n +1)=y( n)+0.1/(3*y( n) 2)*(2*x( n)+0.1*(1-4*x( n)A2/(3*y(门)人3);end ;结果为Columns 1 through 90.1000 0.20

8、000.30000.40000.50000.60000.70000.80000.9000Column 101.0000Columns 1 through 91.0033 1.01321.0293 1.05101.07761.10841.14271.17991.2193Column 101.2606四阶泰勒方法x,y2xynf :<n, ynyfXn，yn28x13y29y5ynfXn,yn4Xn2Xn3yn80x；yn4)Xn, yn240 Xn128OX424xy3y23y324x2x3y23y33y228x23y29y516x440x2y9y53y3'C 69y4x4240x

9、2x9y53y39y63y24x80 x33y527 y8480x220x640 x33y59y83y627y9480x220x640 x33y59y83y627y9440x21280x43y53y881 y11x,yx,yx,yy2x3y23y2481 y；1yn 1 yn hynyn hynh2尹h2h3 尹h3h47ryn-yn4)24yb0=2/(3*y0A2)-8*x0A2/(9*y0A5);%二阶导数y。a使用excel表格进行运算，相应结果如下例二：四阶泰勒方法nxyy的一阶导数y的二阶导数y的三阶导数y的四阶导数精确解00100.6666670-1.33333110.11.00

10、33280.0662250.653509-0.13402-1.183061.00332220.21.0131910.1298840.616121-0.2711-0.790411.01315930.31.0292110.1888080.560087-0.40991-0.29471.02914240.41.0508230.2414960.49274-0.543790.1597871.05071850.51.0773460.2871890.421266-0.663420.4828091.07721760.61.1080630.3257850.351405-0.760550.6515451.10793

11、270.71.1422740.3576560.286967-0.830570.6901671.14216580.81.1793390.3834610.229962-0.872960.6414121.17927490.91.2186920.4039830.181038-0.89030.546371.2186891011.259850.4200210.13996-0.886940.4355851.259921现用matlab编程，程序如下xO=O;yO=1;yaO=2*xO/(3*yOA2);% 一阶导数v( L+U)A* L8)/寸 v( l+u)x*08z l-(8v( L+u)A*e)dv(

12、 l+u)x*o 寸+(gv( L+u)A*e)/寸d( l+u)pa=8 v( L+U)A*卜 2)/e v( L+u)x*08-(g v( L+u) A*g( L+u)x*寸d( L+U)OA否 v( L +u)a*6)/2 v( L+U)X*8-(2V( L+u)A*e)dH( L+u)qAm2v( L+u)A*g( L+U)x* 卡(L+u)eAMU)PA* 寸 2/5§0+(u)OA*9/500+(u)qA*2/L00+(u)eA*L0+(u)AH(L+u)Aml+u)*loh(l+u)x6二Hu04M L L v( l)a* L8)EV( L)x*08z L-(8V( L

13、)A*e)dv( L)x*o 寸+(gv( L)A*e)/H( L)PAm8v( l)a* 卜 2)/ev( L)x*08-(gv( L)A*g( l)x*h( l)oaMgv( l)a*6)v( l)x*8-(2v( L)A*g员 L)qA懸吟盒齐齐 m2v( L)A*g( l)x*zh( L)eA6PA* 寸 2/500?OOA*9/500+oqA*2/LOO+oeA*LO+OAH(L)Aloh(l)xM L L VOA* L8)/寸 vox*§z L-(8VOA*e)/2vox*o?(gvoA*e)/HOPAM8V0A* 卜 2)/ev0x*08-(gv0A*e)/0x*寸-h

14、ooaColumns 1 through 80.10000.20000.30000.40000.50000.60000.70000.8000Columns 9 through 100.90001.0000Columns 1 through 81.00331.01321.02921.05081.07731.10811.14231.1793Columns 9 through 101.21871.2598例3 :用标准四阶R K方法求dy 2 X dxy（o）3y21在区间0, 1上,h0.1的数值解以及在区间1, 10上,h1的数值解，并与精确解进行比较。解：对于本题k1hf(Xn,yn)=绘k2

15、hf (Xnk3hf (Xnhh，ynk4hf (Xnh, ynk3)2hXnh23yn0 222hXnh23yn& 222hXnh32yn k3ynyn珈2k262k3k4)使用excel表格进行运算，相应结果如下例三：标准四阶R K方法nXyk1k2k3k4精确解00100.0033330.0033220.00662311.003320.006620.0098610.12390.0098370.0129891.0033221.013150.0129820.2990.016030.0159830.0188831.0131591.029140.018880.0216330.33320.

16、0215750.0241541.0291421.050710.0241540.4840.026560.02650.0287261.0507181.077210.028720.0307750.57620.0307150.0325861.0772171.107930.032580.0342860.62660.0342340.0357721.1079321.142160.035770.0371570.75250.0371110.038351.1421651.179270.0394580.840.0383540.0394170.0403981.1792741.218680.040390.0412690

17、.99940.0412350.0419971.2186891.259920.041990.042661011730.0426410.0432221.259921现用 matlab 编程，程序如下x0=0;y0=1;k10=2*0.1*x0/(3*y0A2);k20=2*0.1*(x0+0.05)/(3*(y0+k10/2)A2);k30=2*0.1*(x0+0.05)/(3*(y0+k20/2)A2);k40=2*0.1*(x0+0.1)/(3*(y0+k30)A2);x(1)=0.1;y(1)=y0+(k10+2*k20+2*k30+k40)/6;k1(1)=2*0.1*x(1)/(3*y(

18、1)A2);k2(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k1(1)/2)A2);k3(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k2(1)/2)A2);k4(1)=2*0.1*(x(1)+0.1)/(3*(y(1)+k3(1)A2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+(k1(n)+2*k2(n)+2*k3(n)+k4(n)/6;k1( n+1)=2*0.1*x( n+1)/(3*y(门+1)人2);k2( n+1)=2*0.1*(x( n+1)+0.05)/(3*(y( n+1)+k1(门+1)/2)人2);k3( n+1)=2*0.1*(x( n+1)+0.05)/