《数值分析实验报告之常微分方程数值解.doc》由会员分享，可在线阅读，更多相关《数值分析实验报告之常微分方程数值解.doc（17页珍藏版）》请在三一文库上搜索。
1、数学与计算科学学院实 验 报 告实验项目名称 常微分方程数值解 所属课程名称 数值方法B 实 验 类 型 验证 实 验 日 期 2013.11.11 班 级 学 号 姓 名 成 绩 一、实验概述：【实验目的】1.掌握求解常微分方程的欧拉法；2.掌握求解常微分方程的预估校正法；3.掌握求解常微分方程的经典的四阶龙格库塔法；4.能用C语言或MATLAB将上述三种算法用程序运行出来；5.将算法实例化，并得出三种算法的相关关系，如收敛性、精度等；6.附带书中例题的源程序见附录1。【实验原理】1欧拉格式（1）显式欧拉格式：局部截断误差：（2）隐式欧拉格式：局部截断误差：2预估校正法预估：校正：统一格式：
2、平均化格式：3四阶龙格库塔方法的格式（经典格式）【实验环境】1.硬件环境：HPMicrosoft76481-640-8834005-23929HP CorporationIntel(R) Core(TM)I5-2400 CPU 3.10GHz3.09GHz,3.16GB的内存2.软件环境：Microsoft Windows XPProfessional版本 2002Service Pack 3二、实验内容：【实验方案】方案一：用欧拉法，预估校正法，经典的四阶龙格库塔方法求解下列ODE问题：例题：在区间【0，1】上以h=0.1用欧拉法，预估校正法,经典的四阶龙格库塔法求解微分方程 dy/dx=-
3、y+x+1，初值y（0）=1；其精确解为y=x+exp(-x),且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。方案二：用欧拉法，预估校正法, 经典的四阶龙格库塔方法求解初值问题 dy/dx=，初值y（0）=1； 将计算结果与精确解为比较在区间0,1上分别取步长h=0.1; 0.05时进行计算。对三个算法的收敛性进行分析比较，【实验过程】（实验步骤、记录、数据、分析）注：以下图形是通过Excel表格处理数据得出，并未通过MATLAB编程序所得！1、由题可知精确解为：，当x=0时，y(x)=0。h=0.1表1 h=0.1时三个方法与精确值的真值表步长Euler法预估校正法经典四阶
4、库精确值0.11.0100001.0050001.0048381.2490800.21.0290001.0190251.0187311.0554550.31.0561001.0412181.0408181.0912170.41.0904901.0708021.0703201.1318030.51.1314411.1070761.1065311.1768510.61.1782971.1494041.1488121.2260250.71.2304671.1972111.1965861.2790160.81.2874211.2499751.2493291.3355360.91.3486781.307
5、2281.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法预估校正法经典四阶库精确值0.051.0025001.0012501.0012291.0117210.101.0073751.0048771.0048371.0249080.151.0145061.0107641.0107081.0395040.201.0237811.0188021.0187311.0554550.2
6、51.0350921.0288851.0288011.0727100.301.0483371.0409151.0408181.0912170.351.0634211.0547951.0546881.1109310.401.0802501.0704361.0703201.1318010.451.0987371.0877521.0876281.1537910.501.1188001.1066621.1065311.1768510.551.1403601.1270871.1269501.2009420.601.1633421.1489541.1488121.2260250.651.1876751.1
7、721931.1720461.2520620.701.2132911.1967361.1965851.2790160.751.2401271.2225201.2223671.3068520.801.2681211.2494851.2493291.3355360.851.2972151.2775721.2774151.3650370.901.3273541.3067281.3065701.3953220.951.3584861.3369001.3367411.4263621.001.3905621.3680391.3678801.458127图2 h=0.05时三个方法走势图2、由题可知精确解为
8、：，当x=0时，y(x)=0。h=0.1表3 h=0.1时三个方法与精确值的真值表步长Euler法预估校正法经典四阶库精确值0.10.9000000.9096250.9094280.9295330.20.8192490.8359270.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
9、.7208740.80.5841470.6115820.6111000.7069080.90.5685750.5950500.5945990.6959271.00.5562970.5814870.5810720.687191图3 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.870
10、6700.8706060.8995110.200.8277410.8356710.8355920.8725640.250.7949080.8041370.8040470.8484190.300.7654470.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
11、.7493970.600.6435190.6568300.6567100.7387090.650.6301240.6435140.6433950.7292470.700.6181850.6315700.6314530.7208740.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
12、.5690350.5811670.5810710.687191图4 h=0.05时三个方法走势图【实验结论】（结果）1.预估校正法的精确度比经典的四阶库法的精确度高，库塔法最低；2.从表中数据可知三个算法所得数据与精确值相比，可得出以下结论（针对方案二）：（1）Euler法所得值偏离精确值最大，因此可知其精度相对来说最差；（2）经典的四阶库塔法所得值与精确值距离较近，因此可知对于Euler来说，该法更加有效；（3）预估校正法的数据时距离精确值最近的，其骤减幅度较小，因此对精度上的考虑而言，预估校正法应属于最佳解法；（4）由数据可知，上述两个方程的解的光滑性都比较差，从而导致四阶龙格库塔法的精度
13、低于预估校正法的精度。3.由图形可知，三个算法所得数据均呈递减趋势，对于他们的收敛性有以下结论：（1）用上述三法得到的结果大致趋近于0.581，相对于精确值来说，还是存在较大的误差；（2）就误差最小，应首选预估校正法对问题进行求解，它与经典的四阶库塔法所得值比较接近，因此误差不会相差太大。【实验小结】（收获体会）由此次试验，我一方面强化了自己的编程能力，另一方面也对（1）库塔法，（2）预估校正法，（3）经典的四阶龙格库塔法有了全新的认识，并能巧妙的将他们运用到数学建模中，解决一些追求高精度的问题。其次在使用上述三种方法时要充分考虑方程的解的光滑性质，并对照光滑性的好坏选择相应的求解方法，以达到
14、想要的精度的目的。三、指导教师评语及成绩：评 语评语等级优良中及格不及格1.实验报告按时完成,字迹清楚,文字叙述流畅,逻辑性强2.实验方案设计合理3.实验过程（实验步骤详细,记录完整,数据合理,分析透彻）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);if(p=0)xp=0.0;yp=1.0;yp=s
15、qrt(1+2*xp); printf(x%d=%f,y%d=%fn,p,xp,p,yp)；xp+=h;Euler格式：#include 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);预估校正
16、格式：#include main()int i,t;float yi,xi,m,xp,n,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;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);经典的四阶龙格库塔方
17、法的格式：#include main()int i,t;float yi,xi,K1,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);K3=yi+h*K2/2.0-(2*xi+h)/(yi+h*K2/2.0);K4=yi+h*K3-2*(x
18、i+h)/(yi+h*K3);yt=yi+h*(K1+2*K2+2*K3+K4)/6.0;xt=xi+h;xi+=h;yi=yt;printf(x%d=%f,y%d=%fn,t,xt,t,yt);2.精确解：（方案一）#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=
19、%f,y%d=%fn,p,xp,p,yp); xp+=h;Euler格式：#include main()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=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,
20、yt,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;xi+=h;printf(x%d=%f,y%d=%fn,t,xt,t,yt);经典的四阶龙格库塔方法的格式：#include main()int i,t;float yi,xi,K1,K2,K3,K
21、4,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;yi=yt;printf(x%d=%f,y%d
22、=%fn,t,xt,t,yt);3.精确解：（方案二）#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*xp+2)*pow(e,-xp)/2.0;printf(x%d=%f,y%d=%fn,p,xp,p,yp);xp+=h;Euler格式：#include #include #define e 2.1828m
23、ain()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=%fn,i,xi,i,yi);预估校正格式：#include #include #define e 2.1828main()int
24、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;printf(x%d=%f,y%d=%fn,t,xt,t,yt);经典的四阶龙格库塔方法的格式：#inclu
25、de #include #define e 2.1828main()int i,t;float yi,xi,K1,K2,K3,K4,yt,xt,h;xi=0.0;printf(请输入步长h:);scanf(%f,&h);for(i=0;i10;i+)t=i+1;printf(t=%d ,t);if(i=0)xi=0.0;yi=1.0;K1=xi*pow(e,-xi)-yi;K2=(xi+h/2.0)*pow(e,-(xi+h/2.0)-(yi+h*K1/2.0);K3=(xi+h/2.0)*pow(e,-(xi+h/2.0)-(yi+h*K2/2.0);K4=(xi+h)*pow(e,-(xi
26、+h)-(yi+h*K3);yt=yi+h*(K1+2*K2+2*K3+K4)/6.0;xt=xi+h;xi+=h;yi=yt;printf(x%d=%f,y%d=%fn,t,xt,t,yt);附录2：实验报告填写说明 1实验项目名称：要求与实验教学大纲一致。2实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。3实验原理：简要说明本实验项目所涉及的理论知识。4实验环境：实验用的软、硬件环境。5实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。概括整个实验过程。对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设计思路和设计方法，再配以相应的文字说明。对于创新性实验，还应注明其创新点、特色。6实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。7实验结论（结果）：根据实验过程中得到的结果，做出结论。8实验小结：本次实验心得体会、思考和建议。9指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。16