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

例1用欧拉方法与改进的欧拉方法求初值问题在区间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.2443141.259921现用matlab编程，程序如下x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y02);for n=1:9 x(n+1)=0.1*(n+1); y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)2); end; x y 结果为x = Columns 1 through 8 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 10 0.9000 1.0000y = Columns 1 through 8 1.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634 Columns 9 through 10 1.2028 1.2443改进的欧拉方法其计算公式为本题的精确解为，可用来检验数值解的精确度，列出计算结果。使用excel表格进行运算，相应如下例一:改进的欧拉法nx预测y校正y精确解0011110.111.0033331.00332220.21.0099561.013181.01315930.31.0261691.0291711.02914240.41.0480541.0507511.05071850.51.0749041.0772521.07721760.61.1059761.1079651.10793270.71.1405491.1421941.14216580.81.1779651.1792971.17927490.91.2176461.2187061.2186891011.2591031.259931.259921现用matlab编程，程序如下x0=0;y0=1;x(1)=0.1;ya(1)=y0+0.1*2*x0/(3*y02);y(1)=y0+0.05*(2*x0/(3*y02)+2*x0/(3*ya2);for n=1:9 x(n+1)=0.1*(n+1); ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)2); y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)2)+2*x(n+1)/(3*ya(n+1)2); end; x y 结果为x = Columns 1 through 8 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 10 0.9000 1.0000y = Columns 1 through 8 1.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783 Columns 9 through 101.2183 1.2600例2用泰勒方法解分别用二阶、四阶泰勒方法计算点0.1, 0.2, , 1.0处的数值解，并与精确解进行比较。解：二阶泰勒方法对于本题故使用excel表格进行运算，相应结果如下nxy精确解001110.11.0033331.00332220.21.0132231.01315930.31.0292911.02914240.41.0509691.05071850.51.0775751.07721760.61.1083881.10793270.71.1427041.14216580.81.1798781.17927490.91.219341.2186891011.2606011.259921现用matlab编程，程序如下x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1/(3*y02)*(2*x0+0.1*(1-4*x02/(3*y03);for n=1:9 x(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)2/(3*y(n)3); end; x y结果为x = Columns 1 through 9 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 Column 10 1.0000y = Columns 1 through 9 1.0033 1.0132 1.0293 1.0510 1.0776 1.1084 1.1427 1.1799 1.2193 Column 10 1.2606四阶泰勒方法使用excel表格进行运算，相应结果如下例二：四阶泰勒方法nxyy的一阶导数y的二阶导数y的三阶导数y的四阶导数精确解00100.666667 0-1.33333110.11.0033280.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.10793270.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编程，程序如下x0=0;y0=1;ya0=2*x0/(3*y02);%一阶导数yb0=2/(3*y02)-8*x02/(9*y05);%二阶导数yc0=-4*x0/(3*y05)-80*x03/(27*y08);%三阶导数yd0=-4/(3*y05)+40*x02/(3*y08)-1280*x04/(81*y011);%四阶导数x(1)=0.1;y(1)=y0+0.1*ya0+0.01/2*yb0+0.001/6*yc0+0.0001/24*yd0;ya(1)=2*x(1)/(3*y(1)2);%一阶导数yb(1)=2/(3*y(1)2)-8*x(1)2/(9*y(1)5);%二阶导数yc(1)=-4*x(1)/(3*y(1)5)-80*x(1)3/(27*y(1)8);%三阶导数yd(1)=-4/(3*y(1)5)+40*x(1)2/(3*y(1)8)-1280*x(1)4/(81*y(1)11);%四阶导数for n=1:9 x(n+1)=0.1*(n+1); y(n+1)=y(n)+0.1*ya(n)+0.01/2*yb(n)+0.001/6*yc(n)+0.0001/24*yd(n); ya(n+1)=2*x(n+1)/(3*y(n+1)2);%一阶导数 yb(n+1)=2/(3*y(n+1)2)-8*x(n+1)2/(9*y(n+1)5);%二阶导数 yc(n+1)=-4*x(n+1)/(3*y(n+1)5)-80*x(n+1)3/(27*y(n+1)8);%三阶导数 yd(n+1)=-4/(3*y(n+1)5)+40*x(n+1)2/(3*y(n+1)8)-1280*x(n+1)4/(81*y(n+1)11);%四阶导数 end; xY结果为x = Columns 1 through 8 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 10 0.9000 1.0000y = Columns 1 through 8 1.0033 1.0132 1.0292 1.0508 1.0773 1.1081 1.1423 1.1793 Columns 9 through 101.2187 1.2598例3用标准四阶RK方法求在区间0, 1上，的数值解以及在区间1, 10上，的数值解，并与精确解进行比较。解：对于本题使用excel表格进行运算，相应结果如下例三：标准四阶RK方法nxyk1k2k3k4精确解00100.0033330.0033220.006623110.11.0033220.0066230.0098690.0098370.0129891.00332220.21.0131590.0129890.016030.0159830.0188831.01315930.31.0291430.0188830.0216320.0215750.0241541.02914240.41.0507180.0241540.026560.02650.0287261.05071850.51.0772170.0287260.0307720.0307150.0325861.07721760.61.1079320.0325860.0342860.0342340.0357721.10793270.71.1421650.0357720.0371550.0371110.038351.14216580.81.1792740.038350.0394540.0394170.0403981.17927490.91.2186890.0403990.0412640.0412350.0419971.2186891011.2599210.0419970.0426630.0426410.0432221.259921现用matlab编程，程序如下x0=0;y0=1;k10=2*0.1*x0/(3*y02);k20=2*0.1*(x0+0.05)/(3*(y0+k10/2)2);k30=2*0.1*(x0+0.05)/(3*(y0+k20/2)2);k40=2*0.1*(x0+0.1)/(3*(y0+k30)2);x(1)=0.1;y(1)=y0+(k10+2*k20+2*k30+k40)/6;k1(1)=2*0.1*x(1)/(3*y(1)2);k2(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k1(1)/2)2);k3(1)=2*0.1*(x(1)+0.05)/(3*(y(1)+k2(1)/2)2);k4(1)=2*0.1*(x(1)+0.1)/(3*(y(1)+k3(1)2);for n=1:9 x(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(n+1)2); k2(n+1)=2*0.1*(x(n+1)+0.05)/(3*(y(n+1)+k1(n+1)/2)2); k3(n+1)=2*0.1*(x(n+1)+0.05)/(3*(y(n+1)+k2(n+1)/2)2); k4(n+1)=2*0.1*(x(n+1)+0.1)/(3*(y(n+1)+k3(n+1)2); end; x y结果为x = Columns 1 through 8 0.1000 0.2000 0.3000 0