1、数值实验报告 七 一、 实验名称 解初值问题各种方法比较 二、 实验目的 掌握了解各种解初值问题的方法，体会步长对问题解的影响。 三、 实验内容 给定初值问题 22 ,1 2(0 ) 1xdy y x e xd x xy 其精确解为 2()xy x e e 四、 实验要求 分别按 （ 1）欧拉法，步长 0.025, 0.1hh； （ 2）改进的欧拉法，步长 0.05, 0.1hh； （ 3）四阶标准龙格 -库塔法，步长 0.1h ； 编写程序求在节点 1 0 .1 ( 1, 2 , 1 0 )kx k k 处的数值解及误差，并比较各方法的优缺点。用 MATLAB 中的内部函数 dsolve求此
2、常微分方程初值问题的解并与上述结果进行比较。 五、 实验步骤 1、 欧拉法 （ 1） 建立主程序 function y = DEEuler(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步 长 %a:自变量取值下限 %b:自变量取值上限 %y0:函数初值 %varvec:常微分方程的变量组 format long; N = (b-a)/h; y = zeros(N+1,1); x = a:h:b; y(1) = y0; for i=2:N+1 y(i)=y(i-1)+h*(2/x(i-1)*y(i-1)+x(i-1).2*exp(x(i-1);
3、end format short; （ 2）在 matlab工作窗口输入： syms u v;z=2/u*v+u2*exp(u);yy=DEEuler(z,0.1,1,2,0,u,v) syms u v;z=2/u*v+u2*exp(u);yy=DEEuler(z,0.025,1,2,0,u,v) 运行后得出结果： h=0.1 yy = 0 0.2718 0.6848 1.2770 2.0935 3.1874 4.6208 6.4664 8.8091 11.7480 15.3982 h=0.025 yy = 0 0.0680 0.1445 0.2301 0.3255 0.4311 0.5478
4、 0.6760 0.8165 0.9701 1.1374 1.3192 1.5164 1.7297 1.9601 2.2085 2.4757 2.7629 3.0709 3.4009 3.7539 4.1311 4.5337 4.9630 5.4201 5.9065 6.4235 6.9725 7.5551 8.1728 8.8272 9.5200 10.2529 11.0277 11.8464 12.7107 13.6228 14.5847 15.5985 16.6667 17.7914 2、改进的欧拉法 （ 1）建立主程序 function y = DEModifEuler(f, h,a,
5、b,y0,varvec) %f:一阶常微分方程的一 般表达式的右端函数 %h:积分步长 %a:自变量取值下限 %b:自变量取值上限 %y0:函数初值 %varvec:常微分方程的变量组 format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 v1=2/x(i-1)*y(i-1)+x(i-1).2*exp(x(i-1); t = y(i-1) + h*v1; v2=2/x(i)*t+x(i).2*exp(x(i); y(i) = y(i-1)+h*(v1+v2)
6、/2; end format short; （ 2）在 matlab工作窗口输入： syms u v;z=2/u*v+u2*exp(u);yy=DEModifEuler(z,0.1,1,2,0,u,v) syms u v;z=2/u*v+u2*exp(u);yy=DEModifEuler(z,0.05,1,2,0,u,v) 运行后输出 h=0.1 yy = 0 0.3424 0.8583 1.5927 2.5983 3.9364 5.6789 7.9092 10.7245 14.2374 18.5789 h=0.05 yy = 0 0.1532 0.3449 0.5801 0.8643 1.2
7、032 1.6031 2.0710 2.6142 3.2406 3.9589 4.7785 5.7092 6.7621 7.9487 9.2816 10.7744 12.4418 14.2993 16.3641 18.6542 3、 四阶标准龙格 -库塔法 （ 1） 建立主程序 function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长 %a:自变量取值下限 %b:自变量取值上限 %y0:函数初值 %varvec:常微分方程的变量组 format long; N = (b-a)/h; y = z
8、eros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 =2/x(i-1)*y(i-1)+x(i-1).2*exp(x(i-1); K2 = 2/(x(i-1)+h/2)*(y(i-1)+h/2*K1)+(x(i-1)+h/2).2*exp(x(i-1)+h/2); K3 =2/(x(i-1)+h/2)*(y(i-1)+h/2*K2)+(x(i-1)+h/2).2*exp(x(i-1)+h/2); K4 =2/(x(i-1)+h)*(y(i-1)+h*K3)+(x(i-1)+h).2*exp(x(i-1)+h);
9、y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end format short; （ 2）在 matlab工作窗口输入： syms x y;z=2/x*y+x2*exp(x);yy=DELGKT4_lungkuta(z,0.1,1,2,0,x,y) 输出后结果： h=0.1 yy = 0 0.3459 0.8666 1.6072 2.6203 3.9676 5.7209 7.9638 10.7935 14.3229 18.6829 4、 使用 dsolve 命令计算 （ 1） h=0.1 在 matlab工作窗口输入 y=dsolve(Dy=2/x*y+x2*exp
10、(x),y(1)=0,x);x=1:0.1:2;yy=subs(y,x) 运行后输出结果 yy = 0 0.3459 0.8666 1.6072 2.6204 3.9677 5.7210 7.9639 10.7936 14.3231 18.6831 （ 2） h=0.05 在 matlab工作窗口输入 y=dsolve(Dy=2/x*y+x2*exp(x),y(1)=0,x);x=1:0.05:2;yy=subs(y,x) 运行后输出结果 yy = 0 0.1537 0.3459 0.5818 0.8666 1.2063 1.6072 2.0761 2.6204 3.2480 3.9677 4
11、.7886 5.7210 6.7755 7.9639 9.2987 10.7936 12.4632 14.3231 16.3903 18.6831 （ 3） h=0.025 在 matlab工作窗口输入 y=dsolve(Dy=2/x*y+x2*exp(x),y(1)=0,x);x=1:0.025:2;yy=subs(y,x) 运行后输出结果 yy = 0 0.0723 0.1537 0.2447 0.3459 0.4581 0.5818 0.7177 0.8666 1.0293 1.2063 1.3987 1.6072 1.8327 2.0761 2.3383 2.6204 2.9232 3
12、.2480 3.5958 3.9677 4.3649 4.7886 5.2402 5.7210 6.2322 6.7755 7.3522 7.9639 8.6122 9.2987 10.0253 10.7936 11.6056 12.4632 13.3683 14.3231 15.3297 16.3903 17.5073 18.6831 六、 实验结果分析 1、欧拉法 h=0.1 近似解 精确解 绝对误差 0 0 0 0.2718 0.3459 0.0741 0.6848 0.8666 0.1818 1.277 1.6072 0.3302 2.0935 2.6204 0.5269 3.1874
13、 3.9677 0.7803 4.6208 5.721 1.1002 6.4664 7.9639 1.4975 8.8091 10.7936 1.9845 11.748 14.3231 2.5751 15.3982 18.6831 3.2849 h=0.025 近似解 精确解 绝对误差 0 0 0 0.068 0.0723 0.0043 0.1445 0.1537 0.0092 0.2301 0.2447 0.0146 0.3255 0.3459 0.0204 0.4311 0.4581 0.027 0.5478 0.5818 0.034 0.676 0.7177 0.0417 0.8165
14、0.8666 0.0501 0.9701 1.0293 0.0592 1.1374 1.2063 0.0689 1.3192 1.3987 0.0795 1.5164 1.6072 0.0908 1.7297 1.8327 0.103 1.9601 2.0761 0.116 2.2085 2.3383 0.1298 2.4757 2.6204 0.1447 2.7629 2.9232 0.1603 3.0709 3.248 0.1771 3.4009 3.5958 0.1949 3.7539 3.9677 0.2138 4.1311 4.3649 0.2338 4.5337 4.7886 0.
15、2549 4.963 5.2402 0.2772 5.4201 5.721 0.3009 5.9065 6.2322 0.3257 6.4235 6.7755 0.352 6.9725 7.3522 0.3797 7.5551 7.9639 0.4088 8.1728 8.6122 0.4394 8.8272 9.2987 0.4715 9.52 10.0253 0.5053 10.2529 10.7936 0.5407 11.0277 11.6056 0.5779 11.8464 12.4632 0.6168 12.7107 13.3683 0.6576 13.6228 14.3231 0.
16、7003 14.5847 15.3297 0.745 15.5985 16.3903 0.7918 16.6667 17.5073 0.8406 17.7914 18.6831 0.8917 2、 改进的欧拉法 h=0.1 近似解 精确解 绝对误差 0 0 0 0.3424 0.3459 0.0035 0.8583 0.8666 0.0083 1.5927 1.6072 0.0145 2.5983 2.6204 0.0221 3.9364 3.9677 0.0313 5.6789 5.721 0.0421 7.9092 7.9639 0.0547 10.7245 10.7936 0.0691
17、14.2374 14.3231 0.0857 18.5789 18.6831 0.1042 h=0.05 近似解 精确解 绝对误差 0 0 0 0.1532 0.1537 0.0005 0.3449 0.3459 0.001 0.5801 0.5818 0.0017 0.8643 0.8666 0.0023 1.2032 1.2063 0.0031 1.6031 1.6072 0.0041 2.071 2.0761 0.0051 2.6142 2.6204 0.0062 3.2406 3.248 0.0074 3.9589 3.9677 0.0088 4.7785 4.7886 0.0101
18、5.7092 5.721 0.0118 6.7621 6.7755 0.0134 7.9487 7.9639 0.0152 9.2816 9.2987 0.0171 10.7744 10.7936 0.0192 12.4418 12.4632 0.0214 14.2993 14.3231 0.0238 16.3641 16.3903 0.0262 18.6542 18.6831 0.0289 3、 四阶标准龙格 -库塔法 h=0.1 近似解 精确解 绝对误差 0 0 0 0.3459 0.3459 0 0.8666 0.8666 0 1.6072 1.6072 0 2.6203 2.6204
19、0.0001 3.9676 3.9677 1E-04 5.7209 5.721 1E-04 7.9638 7.9639 1E-04 10.7935 10.7936 1E-04 14.3229 14.3231 0.0002 18.6829 18.6831 0.0002 分析 比较： 1、 分别比较欧拉法的两个表格和改进的欧拉法的两个表格，可以得出：对于相同的计算方法（欧拉法、改进的欧拉法），取得步长越小，则与真实值的误差越小； 2、 当取得相同的步长时（ h=0.1)，比较欧拉法、改进的欧拉法和 四阶标准龙格 -库塔法 的计算结果，可以发现： 四阶标准龙格 -库塔法 的计算结果具有更好的精确度，其次是改进的欧拉法，最末是欧拉法； 3、 比较欧拉法步长为 0.025和改进的欧拉法步长为 0.05的最后结果，欧拉法的误差为 0.8917，而改进的欧拉法的为 0.0289，在欧拉法取得更小步长的情况下，其计算的最后结果的精确度还没有取更大步长的改进的欧拉法好； 4、 通过以上的分析比较，可以得出： 四阶标准龙格 -库塔法 具有更好的精确度，其次是改进的欧拉法，最后是欧拉法。 产生以上结果的原因是， 四阶标准龙格 -库塔法 的局部截断误差关于 h的阶为四阶，而改进的欧拉法为二阶，欧拉法为一阶，可见阶的高低可说明计算法精度的高低。
