数值实验报告(7)-解初值问题各种方法比较.pdf
数值实验报告数值实验报告七七 一、一、实验名称实验名称 解初值问题各种方法比较 二、二、实验目的实验目的 掌握了解各种解初值问题的方法，体会步长对问题解的影响。 三、三、实验内容实验内容 给定初值问题 2 2 ,12 (0)1 x dy yx ex dxx y 其精确解为 2( ) x yx ee 四、四、实验要求实验要求 分别按 （1）欧拉法，步长0.025,0.1hh； （2）改进的欧拉法，步长0.05,0.1hh； （3）四阶标准龙格-库塔法，步长0.1h ； 编写程序求在节点1 0.1 (1,2,10) k xk k 处的数值解及误差，并比较各方 法的优缺点。 用 MATLAB 中的内部函数 dsolve 求此常微分方程初值问题的解并与 上述结果进行比较。 五、五、实验步骤实验步骤 1 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); 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 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 2、改进的欧拉法、改进的欧拉法 （1）建立主程序 function y = DEModifEuler(f, h,a,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)/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.2032 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 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 = zeros(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); 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 4、使用使用 dsolvedsolve 命令计算命令计算 （1）h=0.1 在 matlab 工作窗口输入 y=dsolve(Dy=2/x*y+x2*exp(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.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.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 1、欧拉法、欧拉法 h=0.1h=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 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.025h=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 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.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.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 2、改进的欧拉法改进的欧拉法 h=0.1h=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 14.2374 14.3231 0.0857 18.5789 18.6831 0.1042 h=0.05h=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 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 3、四阶标准龙格四阶标准龙格- -库塔法库塔法 h=0.1h=0.1 近似解 精确解 绝对误差 0 0 0 0.3459 0.3459 0 0.8666 0.8666 0 1.6072 1.6072 0 2.6203 2.6204 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 的阶为四阶，而改进的欧拉法为二阶，欧拉法为一 阶，可见阶的高低可说明计算法精度的高低。
