数值分析试验报告之常微分方程数值解讲解

1、数学与计算科学学院实验报告实验项目名称 常微分方程数值解所属课程名称 数值方法 B实验类 型 验证实验日 期 2013.11.11班级学号姓名成绩一、实验概述：【实验目的】1. 掌握求解常微分方程的欧拉法；2. 掌握求解常微分方程的预估校正法；3. 掌握求解常微分方程的经典的四阶龙格库塔法；4. 能用 C 语言或 MATLAB 将上述三种算法用程序运行出来；5. 将算法实例化，并得出三种算法的相关关系，如收敛性、精度等；6. 附带书中例题的源程序见附录 1。【实验原理】1欧拉格式（1）显式欧拉格式：yn 1 yn hf (xn, yn)局部截断误差：h2h2y(xn 1) yn 1 h y (

2、 ) h y (xn ) o(h2)22（2）隐式欧拉格式：yn 1 yn hf (xn 1,yn 1)局部截断误差：h22y(xn 1) yn 1y (xn) o(h2)22预估校正法预估：yn 1 yn hf (xn,yn)校正：hyn 1 yn f(xn, yn) f(xn 1,yn 1)2统一格式： yn 1 yn h2 f (xn, yn) f (xn h, yn hf (xn,yn)yp yn hf (xn,yn),平均化格式： yc yn hf (xn 1,yp ),1 yn 1 2(yp yc).3四阶龙格库塔方法的格式（经典格式）yn 1 yn 6h(K1 2K2 2K3 K

3、4), K1 f (xn, yn), hhK2 f (xn h,yn hK1),22 K3 f (xn h2,yn 2hK2), K4 f (xn h,yn hK3).【实验环境】1. 硬件环境：HPMicrosoft76481-640-8834005-23929HP CorporationIntel(R) Core(TM)I5-2400 CPU 3.10GHz3.09GHz,3.16GB的内存2. 软件环境：Microsoft Windows XPProfessional 版本 2002 Service Pack 3、实验内容：实验方案】ODE 问题：方案一： 用欧拉法，预估校正法，经典的四

4、阶龙格库塔方法求解下列例题：在区间【 0 ，1】上以 h=0.1 用欧拉法，预估校正法 ,经典的四阶龙格库塔法求解微分方程 dy/dx=-y+x+1 ，初值 y(0)=1 ；其精确解为 y=x+exp(-x), 且将计 算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。方案二： 用欧拉法，预估校正法 , 经典的四阶龙格库塔方法 求解初值问题1dy/dx= xe x y，初值 y(0)=1； 将计算结果与精确解为 1(x2 2)e x 比较在区 2间0,1 上分别取步长 h=0.1; 0.05 时进行计算。对三个算法的收敛性进行分析比较，【实验过程】(实验步骤、记录、数据、分析)注：以下

5、图形是通过 Excel 表格处理数据得出，并未通过 MATLAB 编程序所得！ dy y x 11、 dxy(0) 1由题可知精确解为： y x e ，当 x=0 时， y(x)=0 。h=0.1表 1 h=0.1 时三个方法与精确值的真值表步长Euler 法预估校正法经典四阶库精确值0.11.0100001.0050001.0048381.2490800.21.0290001.0190251.0187311.0554550.31.0561001.0412181.0408181.0912170.41.0904901.0708021.0703201.1318030.51.1314411.1070

6、761.1065311.1768510.61.1782971.1494041.1488121.2260250.71.2304671.1972111.1965861.2790160.81.2874211.2499751.2493291.3355360.91.3486781.3072281.3065701.3953221.01.4138111.3685411.3678801.458127图 1 h=0.1 时三个方法走势图h=0.05 (此时将源程序中 i 的范围进行扩大，即 for(i=0;i20;i+) )表 2 h=0.05 时三个方法与精确值的真值表步长Euler 法预估校正法经典四阶库精

7、确值0.051.0025001.0012501.0012291.0117210.101.0073751.0048771.0048371.0249080.151.0145061.0107641.0107081.0395040.201.0237811.0188021.0187311.0554550.251.0350921.0288851.0288011.0727100.301.0483371.0409151.0408181.0912170.351.0634211.0547951.0546881.1109310.401.0802501.0704361.0703201.1318010.451.0987

8、371.0877521.0876281.1537910.501.1188001.1066621.1065311.1768510.551.1403601.1270871.1269501.2009420.601.1633421.1489541.1488121.2260250.651.1876751.1721931.1720461.2520620.701.2132911.1967361.1965851.2790160.751.2401271.2225201.2223671.3068520.801.2681211.2494851.2493291.3355360.851.2972151.2775721.

9、2774151.3650370.901.3273541.3067281.3065701.3953220.951.3584861.3369001.3367411.4263621.001.3905621.3680391.3678801.458127图 2 h=0.05 时三个方法走势图52、 ddyx xe x y y(0) 1由题可知精确解为：12(x2 2)ex，当 x=0时，y(x)=0。h=0.1表 3 h=0.1 时三个方法与精确值的真值表步长Euler 法预估校正法经典四阶库精确值0.10.9000000.9096250.9094280.9295330.20.8192490.83592

10、70.8355930.8725640.30.7544330.7760810.7756550.8268220.40.7027260.7276710.7271890.7903480.50.6617260.6886360.6881270.7614570.60.6293960.6572250.6567110.7387090.70.6040180.6319570.6314530.7208740.80.5841470.6115820.6111000.7069080.90.5685750.5950500.5945990.6959271.00.5562970.5814870.5810720.687191图 3

11、 h=0.1 时三个方法走势图h=0.05 (此时将源程序中 i 的范围进行扩大，即 for(i=0;i20;i+) )表 4 h=0.05 时三个方法与精确值的真值表步长Euler 法预估校正法经典四阶库精确值0.050.9500000.9524520.9524270.9629240.100.9049040.9094740.9094280.9295330.150.8642840.8706700.8706060.8995110.200.8277410.8356710.8355920.8725640.250.7949080.8041370.8040470.8484190.300.7654470.

12、7757550.7756550.8268220.350.7390430.7502320.7501250.8075380.400.7154070.7273020.7271890.7903480.450.6942720.7067150.7065990.7750500.500.6753940.6882450.6881260.7614570.550.6585460.6716820.6715610.7493970.600.6435190.6568300.6567100.7387090.650.6301240.6435140.6433950.7292470.700.6181850.6315700.6314

13、530.72087470.750.6075410.6208480.6207330.7134660.800.5980460.6112110.6111000.7069080.850.5895650.6025350.6024260.7010940.900.5819760.5947030.5945990.6959270.950.5751670.5876120.5875120.6913201.000.5690350.5811670.5810710.687191图 4 h=0.05 时三个方法走势图【实验结论】（结果）1.预估校正法的精确度比经典的四阶库法的精确度高，库塔法最低；2.从表中数据可知三个算法

14、所得数据与精确值相比，可得出以下结论（针对方 案二）：（1）Euler 法所得值偏离精确值最大，因此可知其精度相对来说最差；（2）经典的四阶库塔法所得值与精确值距离较近，因此可知对于Euler 来说，该法更加有效；（3）预估校正法的数据时距离精确值最近的，其骤减幅度较小，因此对精度上的考虑而言，预估校正法应属于最佳解法；（4）由数据可知，上述两个方程的解的光滑性都比较差，从而导致四阶龙格库 塔法的精度低于预估校正法的精度。3.由图形可知，三个算法所得数据均呈递减趋势，对于他们的收敛性有以下结 论：（1）用上述三法得到的结果大致趋近于 0.581 ，相对于精确值来说，还是存在 较大的误差；（2）

15、就误差最小，应首选预估校正法对问题进行求解，它与经典的四阶库塔法 所得值比较接近，因此误差不会相差太大。【实验小结】（收获体会）由此次试验，我一方面强化了自己的编程能力，另一方面也对（1 ）库塔法，（2）预估校正法，（3）经典的四阶龙格库塔法有了全新的认识，并能巧妙的将他 们运用到数学建模中，解决一些追求高精度的问题。其次在使用上述三种方法时要充分考虑方程的解的光滑性质， 并对照光滑性的好坏选择相应的求解方法，以达到想要的精度的目的。三、指导教师评语及成绩：评语评语等级优良中及格不及格1.实验报告按时完成 ,字迹清楚 , 文字叙述流畅 ,逻辑性强2.实验方案设计合理3. 实验过程（实验步骤详细

16、 ,记录完整 ,数据合理 ,分析透彻）4 实验结论正确 .成 绩：指导教师签名：批阅日期：附录 1：源 程 序1. 书中例题： 精确值： #include #include main() int i,p;float y0,x0,yi,xi,yp,xp,h; xi=0.0;printf( 请输入步长 h:); scanf(%f,&h); for(p=0;p=10;p+)printf(p=%d ,p);10if(p=0)xp=0.0; yp=1.0; yp=sqrt(1+2*xp); printf(x%d=%f,y%d=%fn,p,xp,p,yp) xp+=h;Euler 格式： #include

17、 main()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);for(i=1;i=10;i+)p=i-1; printf(i=%d ,i);if(p=0)xp=0.0; yp=1.0; yi=yp+h*(yp-2*xp/yp); printf(t=%f ,xp/yp); yp=yi; xi+=h; xp=xi;printf(x%d=%f,y%d=%fn,i,xi,i,yi);预估校正格式： #include main()int i,t;float yi,xi,m,xp,n,yt,xt,h;11xi

18、=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(t=0 x0=0.000000,y0=1.000000n) for(i=0;i10;i+)t=i+1; printf(t=%d ,t); if(i=0) xi=0.0; yi=1.0; xt=xi+h;m=yi+h*(yi-2*xi/yi); n=yi+h*(m-2*xt/m); yt=(m+n)/2.0; yi=yt; xi+=h;printf(x%d=%f,y%d=%fn,t,xt,t,yt); 经典的四阶龙格库塔方法的格式： #include main()int i,t;float yi,xi,K1,

19、K2,K3,K4,xp,yt,xt,h; xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(t=0 x0=0.000000,y0=1.000000n) for(i=0;i10;i+)t=i+1; printf(t=%d ,t); if(i=0) xi=0.0; yi=1.0;K1=yi-2*xi/yi; K2=yi+h*K1/2.0-(2*xi+h)/(yi+h*K1/2.0);12K3=yi+h*K2/2.0-(2*xi+h)/(yi+h*K2/2.0); K4=yi+h*K3-2*(xi+h)/(yi+h*K3); yt=yi+h*(K1+2*K2

20、+2*K3+K4)/6.0; xt=xi+h;xi+=h;yi=yt;printf(x%d=%f,y%d=%fn,t,xt,t,yt);2.精确解： y x e x ( 方案一)#include #include #define e 2.1828main()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);for(p=0;p=10;p+)printf(p=%d ,p);if(p=0)xp=0.0;yp=1.0;yp=xp+pow(e,-xp);printf(x%d=%f,y%d=%fn,p,xp,p

21、,yp); xp+=h;Euler 格式：#include main()int i,p;float y0,x0,yi,xi,yp,xp,h;13xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);for(i=1;i=11;i+)p=i-1; printf(p=%d ,p); if(p=0) xp=0.0; yp=1.0;yi=yp+h*(-yp+xp+1); yp=yi;printf(x%d=%f,y%d=%fn,p,xi,p,yi); xi+=h;xp=xi; 预估校正格式： #include main()int i,t;float yi,xi,m,xp,n,yt,

22、xt,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(t=0 x0=0.000000,y0=1.000000n); for(i=0;i20;i+)t=i+1; printf(t=%d ,t);if(i=0)xi=0.0;yi=1.0; xt=xi+h; m=yi+h*(-yi+xi+1); n=yi+h*(-m+xt+1); yt=(m+n)/2.0; yi=yt;14 xi+=h; printf(x%d=%f,y%d=%fn,t,xt,t,yt);经典的四阶龙格库塔方法的格式： #include main()int i,t;float yi,x

23、i,K1,K2,K3,K4,yt,xt,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(t=0 x0=0.000000,y0=1.000000n); for(i=0;i10;i+)t=i+1; printf(t=%d ,t); if(i=0) xi=0.0; yi=1.0;K1=-yi+xi+1;K2=-(yi+h*K1/2.0)+xi+h/2.0+1;K3=-(yi+h*K2/2.0)+xi+h/2.0+1; K4=-(yi+h*K3)+xi+h+1; yt=yi+h*(K1+2*K2+2*K3+K4)/6.0; xt=xi+h;xi+=h;

24、yi=yt; printf(x%d=%f,y%d=%fn,t,xt,t,yt);3. 精确解： 1(x2 2)e x (方案二)2 #include #include #define e 2.1828 main() int i,p;15float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);for(p=0;p=10;p+)printf(p=%d ,p); if(p=0)xp=0.0;yp=1.0; yp=(xp*xp+2)*pow(e,-xp)/2.0; printf(x%d=%f,y%d=%fn,p,xp,p,yp);

25、 xp+=h;Euler 格式： #include #include #define e 2.1828 main() int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(i=0 x0=0.000000,y0=1.000000n); for(i=1;i=10;i+)p=i-1; printf(i=%d ,i);if(p=0)xp=0.0;yp=1.0; yi=yp+h*(xp*pow(e,-xp)-yp); yp=yi; xi+=h; xp=xi; printf(x%d=%f,y%d=%f

26、n,i,xi,i,yi);16预估校正格式： #include #include #define e 2.1828 main()int i,t;float yi,xi,m,xp,n,yt,xt,h;xi=0.0;printf( 请输入步长 h:); scanf(%f,&h);printf(i=0 x0=0.000000,y0=1.000000n); for(i=0;i20;i+)t=i+1; printf(t=%d ,t); if(i=0) xi=0.0; yi=1.0; xt=xi+h; m=yi+h*(xi*pow(e,-xi)-yi); n=yi+h*(xt*pow(e,-xt)-m); yt=(m+n)/2.0; yi=yt; xi+=h; pri