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

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

2、(2 )隐式欧拉格式：yn 1 二 yn hf (Xn 1,yn 1)局部截断误差：h22y(Xn 1)- Yn 1 ： -(Xn) =0(h )22 预估校正法预估：yn 1 =yn hf (Xn,yn)校正：丄h丄_yn 1 =yn A f (Xn, Yn )f ( X 1 , Yn 1)2统一格式：yn =yn +刁f (Xn, y.) + f(Xh,n + hf(X., y.)ryp in +hf (Xn，yn),平均化格式：丿yc =yn+hf (Xnyp),1 =2(yp +yc).3 四阶龙格库塔方法的格式（经典格式）hyn 卅二yn+Ki+2K2+2K3+K4),Ki = f

3、(Xn, Yn),hhKf(X-,y-Ki),2 2K3 = f (Xn 十/An +孑 K2),iK4= f(Xn +h,yn +h&).【实验环境】1. 硬件环境：HPMicrosoft76481-640-8834005-23929 HP Corporati onIntel(R) Core(TM)I5-2400 CPU 3.10GHz3.09GHz,3.16GB 的内存2. 软件环境：Microsoft Win dows XPProfessi onal 版本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 -y，初值y (0) =1 ;将计算结果与精确解为(x2 2)e 比较在区2间0,1上分别取步长h=0.1;0.05时进行计算。对三个算法的收敛性进行分析比较，【实验过程】(实验步骤、记录、数据、分析)注：以下图形是通过 Excel表格处理数据得出，并未通过MATL

5、AB编程序所得!业=_y +x +11、 dx.y(0) =1由题可知精确解为：y = x /，当x=0时，y(x)=O。h=0.10.8表1h=0.1时三个方法与精确值的真值表rH x.步长Euler 法预估校正法经典四阶库精确值0.11.0100001.0050001.0048381.2490800.21.0290001.0190251.0187311.0554550.3 1.0561001.0412181.0408181.0912170.41.0904901.0708021.0703201.1318030.51.1314411.1070761.1065311.1768510.61.178

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

7、.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.0987371.087

8、7521.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.2774151

9、.3650370.901.3273541.3067281.3065701.3953220.951.3584861.3369001.3367411.4263621.001.3905621.3680391.3678801.4581270.80.60.40.201,6 Eulrj去T-预估枝正法 止经典四阶库一精确值00.10.2 0304 0.50.607 0.80.91图2h=0.05时三个方法走势图2、d：iyy(0) J由题可知精确解为:；(x2 2)，当 x=0 时，y(x)=。h=0.1表3h=0.1时三个方法与精确值的真值表rH x.步长Euler 法预估校正法经典四阶库精确值f 0.

10、10.9000000.9096250.9094280.9295330.20.8192490.8359270.8355930.8725641 0.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

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

12、.8355920.872564；0.250.7949080.8041370.804047(0.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.7493970.600.6435190.6568300.6567100.7

13、38709：0.65 :0.6301240.6435140.643395:0.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.691320：1.000.5690350.5811670.5810710.687191 E

14、uIe法十预估校正法 T-经典四阶库 精确值图4h=0.05时三个方法走势图【实验结论】（结果）1. 预估校正法的精确度比经典的四阶库法的精确度高，库塔法最低；2. 从表中数据可知三个算法所得数据与精确值相比，可得出以下结论（针对方案二）：（1）Euler法所得值偏离精确值最大，因此可知其精度相对来说最差；（2 ）经典的四阶库塔法所得值与精确值距离较近，因此可知对于Euler来说,该法更加有效；（3）预估校正法的数据时距离精确值最近的，其骤减幅度较小，因此对精度上的考虑而言，预估校正法应属于最佳解法；（4）由数据可知，上述两个方程的解的光滑性都比较差，从而导致四阶龙格库 塔法的精度低于预估校正

15、法的精度。3.由图形可知，三个算法所得数据均呈递减趋势，对于他们的收敛性有以下结 论：（1）用上述三法得到的结果大致趋近于 0.581，相对于精确值来说，还是存在 较大的误差；（2）就误差最小，应首选预估校正法对问题进行求解，它与经典的四阶库塔法 所得值比较接近，因此误差不会相差太大。【实验小结】（收获体会）由此次试验，我一方面强化了自己的编程能力，另一方面也对（1 ）库塔法，（2）预估校正法，（3）经典的四阶龙格库塔法有了全新的认识，并能巧妙的将他 们运用到数学建模中，解决一些追求高精度的问题。其次在使用上述三种方法时要充分考虑方程的解的光滑性质，并对照光滑性的好坏选择相应的求解方法，以达到

16、想要的精度的目的。三、指导教师评语及成绩：评语评语等级优良中及格不及格1.实验报告按时完成，字迹清楚，文字叙述流畅，逻辑性强2.实验方案设计合理3.实验过程（实验步骤详细，记录完整，数据合理，分析透彻）4实验结论正确.成绩：指导教师签名：批阅日期：附录1源程序1.书中例题： 精确值：#i nclude #in elude mai n()int i,p;float yO,xO,yi,xi,yp,xp,h; xi=0.0;printf(请输入步长h:);sca nf(%f,&h); for(p=0;pv=10;p+) prin tf(p=%d ,p);if(p=O)xp=O.O;yp=1.0;yp

17、=sqrt(1+2*xp);prin tf(x%d=%f,y%d=%fn,p,xp,p,yp) xp+=h;Euler格式：#i nclude mai n()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);for(i=1;i=10;i+)p=i-1;prin tf(i=%d ,i);if(p=0)xp=0.0;yp=1.0;yi=yp+h*(yp-2*xp/yp);prin tf(t=%f ,xp/yp);yp=yi;xi+=h;xp=xi;prin tf(x%d=%f,y%d=%fn,i,xi,i,

18、yi); 预估校正格式：#i nclude mai n()int i,t;float yi,xi,m,xp, n, yt,xt,h;xi=O.O;printf( 请输入步长h:); sea nf(%f,&h);prin tf(t=O xO=O.OOOOOO,yO=1.OOOOOOn) for(i=0;i10;i+)t=i+1;prin tf(t=%d ,t);if(i=O)xi=O.O;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; prin tf(x%d=%f,y%d=%fn,t,

19、xt,t,yt);经典的四阶龙格库塔方法的格式：#i nclude mai n()int i,t;float yi,xi,K1,K2,K3,K4,xp,yt,xt,h;xi=0.0;printf( 请输入步长h:); sca nf(%f,&h);prin tf(t=O x0=0.000000,y0=1.000000n) for(i=O;i1O;i+)t=i+1;prin tf(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

20、+h*K2/2.0);K4=yi+h*K3-2*(xi+h)/(yi+h*K3); yt=yi+h*(K1+2*K2+2*K3+K4)/6.0; xt=xi+h;xi+=h;yi=yt;prin tf(x%d=%f,y%d=%fn,t,xt,t,yt);2.精确解：讨e (方案一)#i nclude #in clude #define e 2.1828mai n()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);for(p=0;p=10;p+)prin tf(p=%d ,p);if(p=0)xp=0.

21、0;yp=1.0;yp=xp+pow(e,-xp);prin tf(x%d=%f,y%d=%fn,p,xp,p,yp); xp+=h;Euler格式：#i nclude mai n()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=O.O;printf(请输入步长h:);sea nf(%f,&h);for(i=1;i=11;i+)P=i-1;prin tf(p=%d ,p);if(p=0)xp=0.0;yp=1.0;yi=yp+h*(-yp+xp+1);yp=yi;prin tf(x%d=%f,y%d=%fn,p,xi,p,yi); xi+=h;xp=xi;预估校正格

22、式：#i nclude mai n()int i,t;float yi,xi,m,xp, n, yt,xt,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);prin tf(t=0 x0=0.000000,y0=1.000000n); for(i=0;i20;i+)t=i+1;prin tf(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;prin tf(x%d=%f,y%d=%fn,t,xt,t,yt);经典的四阶龙

23、格库塔方法的格式：#i nclude mai n()int i,t;float yi,xi,K1,K2,K3,K4,yt,xt,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);prin tf(t=0 x0=0.000000,y0=1.000000n);for(i=0;i10;i+)t=i+1;prin tf(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

24、*(K1+2*K2+2*K3+K4)/6.0; xt=xi+h;xi+=h;yi=yt;prin tf(x%d=%f,y%d=%fn,t,xt,t,yt);3.精确解：1(x2 2)e(方案二)2#i nclude #in clude #define e 2.1828mai n()int i,p;float yO,xO,yi,xi,yp,xp,h;xi=0.0;printf(请输入步长h:);sea nf(%f,&h);for(p=0;pv=10;p+)prin tf(p=%d ,p);if(p=0)xp=0.0;yp=1.0;yp=(xp*xp+2)*pow(e,-xp)/2.0;prin

25、tf(x%d=%f,y%d=%fn,p,xp,p,yp); xp+=h;Euler格式：#i nclude #in clude #define e 2.1828mai n()int i,p;float y0,x0,yi,xi,yp,xp,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);prin tf(i=0 x0=0.000000,y0=1.000000n); for(i=1;i=10;i+)p=i-1;prin tf(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=x

26、i;prin tf(x%d=%f,y%d=%fn,i,xi,i,yi);预估校正格式：#i nclude #in elude #define e 2.1828mai n()int i,t;float yi,xi,m,xp, n, yt,xt,h;xi=0.0;printf(请输入步长h:);sca nf(%f,&h);prin tf(i=0 x0=0.000000,y0=1.000000n); for(i=0;i20;i+)t=i+1;prin tf(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