常微分方程数值求解问题的实习报告.doc
计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓 名：系：信息与计算科学专 业：信息与计算科学年 级：学 号：指导教师：职 称：讲师2010 年 12 月 1 日目 录1.实习的目的和任务12.实习要求13.实习地点14.主要仪器设备15.实习内容15.1 用不同格式对同一个初值问题的数值求解及其分析.15.1.1求精确解15.1.2用欧拉法求解45.1.3用改进欧拉法求解65.1.4用4级4阶龙格库塔法求解8 5.1.5 问题讨论与分析11 5.2 一个算法不同不长求解同一个初值问题及其分析. 146.结束语15参考文献15常微分方程课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。3. 实习地点数学实验室、学生宿舍4. 主要仪器设备计算机、Microsoft Windows XPMatlab 6.55. 实习内容5.1 用欧拉方法，改进欧拉方法，4阶龙格库塔方法分别求下面微分方程的初值dy/dx= y(0)=1 x0，65.1 .1求精确解首先可以求得其精确解为：y=sinx+1/(x+1)程序代码：x=0:0.2:6; y= sinx+1./(x+1) plot(x,y,'b*-');Data=x',y'>> aay = Columns 1 through 9 1.0000 1.0320 1.1037 1.1896 1.2729 1.3415 1.3866 1.4021 1.3842 Columns 10 through 18 1.3310 1.2426 1.1210 0.9696 0.7933 0.5981 0.3911 0.1797 -0.0283 Columns 19 through 27 -0.2251 -0.4035 -0.5568 -0.6793 -0.7664 -0.8151 -0.8238 -0.7923 -0.7222 Columns 28 through 31 -0.6165 -0.4798 -0.3175 -0.1366Data = 0 1.0000 0.2000 1.0320 0.4000 1.1037 0.6000 1.1896 0.8000 1.2729 1.0000 1.3415 1.2000 1.3866 1.4000 1.4021 1.6000 1.3842 1.8000 1.3310 2.0000 1.2426 2.2000 1.1210 2.4000 0.9696 2.6000 0.7933 2.8000 0.5981 3.0000 0.3911 3.2000 0.1797 3.4000 -0.0283 3.6000 -0.2251 3.8000 -0.4035 4.0000 -0.5568 4.2000 -0.6793 4.4000 -0.7664 4.6000 -0.8151 4.8000 -0.8238 5.0000 -0.7923 5.2000 -0.7222 5.4000 -0.6165 5.6000 -0.4798 5.8000 -0.31756.0000 -0.1366 5.1.2 用欧拉法求解程序如下：建立函数文件cwfa1.mfunction x,y=cwfa1(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx=x'y=y'在MATLAB输入以下程序：fun=inline('-y2+2*y*sin(x)+cos(x)-sin(x)2');x,y=cwfa1(fun,0,6,1,0.2); x,yplot(x,y,'r*-')结果及其图象ans = 0 1.0000 0.2000 1.0000 0.4000 1.0676 0.6000 1.1598 0.8000 1.2540 1.0000 1.3358 1.2000 1.3950 1.4000 1.4246 1.6000 1.4200 1.8000 1.3788 2.0000 1.3006 2.2000 1.1867 2.4000 1.0404 2.6000 0.8663 2.8000 0.6703 3.0000 0.4594 3.2000 0.2411 3.4000 0.0235 3.6000 -0.1854 3.8000 -0.3780 4.0000 -0.5471 4.2000 -0.6866 4.4000 -0.7915 4.6000 -0.8581 4.8000 -0.8842 5.0000 -0.8692 5.2000 -0.8141 5.4000 -0.7214 5.6000 -0.5950 5.8000 -0.4401 6.0000 -0.2631 5.1.3用改进欧拉法求解：程序如下：建立函数文件cwfa2.mfunction x,y=cwfa2(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(fun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2;endx=x'y=y'在MATLAB输入以下程序： fun=inline('-y2+2*y*sin(x)+cos(x)-sin(x)2'); x,y=cwfa2(fun,0,6,1,0.2); x,y plot(x,y,'r*-')结果及其图象：ans = 0 1.0000 0.2000 1.0338 0.4000 1.1050 0.6000 1.1898 0.8000 1.2716 1.0000 1.3388 1.2000 1.3827 1.4000 1.3972 1.6000 1.3785 1.8000 1.3248 2.0000 1.2362 2.2000 1.1146 2.4000 0.9635 2.6000 0.7877 2.8000 0.5933 3.0000 0.3871 3.2000 0.1767 3.4000 -0.0302 3.6000 -0.2260 3.8000 -0.4034 4.0000 -0.5557 4.2000 -0.6772 4.4000 -0.7636 4.6000 -0.8117 4.8000 -0.8199 5.0000 -0.7882 5.2000 -0.7180 5.4000 -0.6125 5.6000 -0.4760 5.8000 -0.31426.0000 -0.1338 5.1.4 用4阶龙格库塔求解程序如下：建立函数文件cwfa3.mfunction x,y=cwfa3(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n); k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h/2,y(n)+h/2*k2); k4=feval(fun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x'y=y'在MATLAB输入以下程序：fun=inline('-y2+2*y*sin(x)+cos(x)-sin(x)2'); x,y=cwfa3(fun,0,6,1,0.2); x,y plot(x,y, 'r*-')结果及其图象： 0 1.0000 0.2000 1.0320 0.4000 1.1037 0.6000 1.1896 0.8000 1.2729 1.0000 1.3415 1.2000 1.3866 1.4000 1.4021 1.6000 1.3842 1.8000 1.3310 2.0000 1.2426 2.2000 1.1210 2.4000 0.9696 2.6000 0.7933 2.8000 0.5981 3.0000 0.3911 3.2000 0.1797 3.4000 -0.0283 3.6000 -0.2251 3.8000 -0.4035 4.0000 -0.5568 4.2000 -0.6793 4.4000 -0.7664 4.6000 -0.8151 4.8000 -0.8238 5.0000 -0.7923 5.2000 -0.7222 5.4000 -0.6165 5.6000 -0.4798 5.8000 -0.31766.0000 -0.1366 5.1.5 问题讨论与分析由以上数值分析结果绘制表格：精确解欧拉方法改进的欧拉方法四阶龙格-库塔方法xiyiyi误差yi误差yi误差011010100.21.03210.0321.03380.00181.03200.41.10371.06760.03611.1050.00131.103700.61.18961.15980.02981.18980.00021.189600.81.27291.2540.01891.27160.00131.2729011.34151.33580.00571.33880.00271.341501.21.38661.3950.00841.38270.00391.386601.41.40211.42460.02251.39720.00491.402101.61.38421.420.03581.37850.00571.384201.81.3311.37880.04781.32480.00621.331021.24261.30060.0581.23620.00641.242602.21.1211.18670.06571.11460.00641.12102.40.969581.04040.070820.963490.006090.969571E-052.60.793280.86630.073020.787710.005570.793271E-052.80.598150.670310.072160.593280.004870.598141E-0530.391120.459370.068250.387120.0040.391111E-053.20.179720.241120.06140.17670.003020.179711E-053.4-0.028270.0235210.051789-0.030240.001976-0.028288E-063.6-0.22513-0.185410.03972-0.226030.0009-0.225141E-053.8-0.40352-0.377990.02553-0.403380.00014-0.403531E-054-0.5568-0.547120.00968-0.555670.00113-0.556811E-054.2-0.67927-0.686640.00737-0.677240.00203-0.679281E-054.4-0.76642-0.791530.02511-0.763620.0028-0.766431E-054.6-0.81512-0.858120.043-0.811710.00341-0.815131E-054.8-0.82375-0.884230.06048-0.81990.00385-0.823772E-055-0.79226-0.869240.07698-0.788150.00411-0.792282E-055.2-0.72216-0.814110.09195-0.717990.00417-0.722193E-055.4-0.61651-0.721370.10486-0.612450.00406-0.616543E-055.6-0.47975-0.594960.11521-0.475970.00378-0.479772E-055.8-0.31754-0.440110.12257-0.314210.00333-0.317573E-056-0.13656-0.263130.12657-0.13380.00276-0.136582E-05x=0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.6,1.8,2,2.2,2.4,2.6,2.8,3,3.2,3.4,3.6,3.8,4,4.2,4.4,4.6,4.8,5,5.2,5.4,5.6,5.8,6;y1=1,1.032,1.1037,1.1896,1.2729,1.3415,1.3866,1.4021,1.3842,1.331,1.2426,1.121,0.96958,0.79328,0.59815,0.39112,0.17972,-0.028268,-0.22513,-0.40352,-0.5568,-0.67927,-0.76642,-0.81512,-0.82375,-0.79226,-0.72216,-0.61651,-0.47975,-0.31754,-0.13656;y2=1,1,1.0676,1.1598,1.254,1.3358,1.395,1.4246,1.42,1.3788,1.3006,1.1867,1.0404,0.8663,0.67031,0.45937,0.24112,0.023521,-0.18541,-0.37799,-0.54712,-0.68664,-0.79153,-0.85812,-0.88423,-0.86924,-0.81411,-0.72137,-0.59496,-0.44011,-0.26313;y3=1,1.0338,1.105,1.1898,1.2716,1.3388,1.3827,1.3972,1.3785,1.3248,1.2362,1.1146,0.96349,0.78771,0.59328,0.38712,0.1767,-0.030244,-0.22603,-0.40338,-0.55567,-0.67724,-0.76362,-0.81171,-0.8199,-0.78815,-0.71799,-0.61245,-0.47597,-0.31421,-0.1338;y4=1,1.032,1.1037,1.1896,1.2729,1.3415,1.3866,1.4021,1.3842,1.331,1.2426,1.121,0.96957,0.79327,0.59814,0.39111,0.17971,-0.028276,-0.22514,-0.40353,-0.55681,-0.67928,-0.76643,-0.81513,-0.82377,-0.79228,-0.72219,-0.61654,-0.47977,-0.31757,-0.13658;plot(x,y1,'r+-')hold on,plot(x,y2,'b-') title('精确解与欧拉方法比较')plot(x,y1,'r+ -') hold on,plot(x,y3,'b-') title('精确解与改进欧拉方法比较')plot(x,y1,'r+-')hold on,plot(x,y4,'b-')title('精确解与4阶龙格库塔法比较')由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格库塔方法误差相对较小，并且龙格库塔方法误差最小且大部分值都跟精确值相同。由欧拉图与精确图相比可清晰看到，随着x的增加，函数值与精确值的偏差越来越大。5.2，用改进欧拉方法取不同步长分别求下面微分方程的初值dy/dx= y(0)=1 x0，6，并对结果进行分析说明，得出结论。fun=inline('-y2+2*y*sin(x)+cos(x)-sin(x)2');x1,y1=cwfa2(fun,0,6,1,0.2);x2,y2=cwfa2(fun,0,6,1,0.5);x3,y3=cwfa2(fun,0,6,1,1);plot(x1,y1,'g*-',x2,y2,'b-',x3,y3,'m*-')gtext('h=0.2')gtext('h=0.5')gtext('h=1')由图可以看错出：随着h的取值越小，改进欧拉方程求解趋于光滑曲线，且步长不能太大，否则可能大大偏离精确解，将不能说明问题。结束语在常微分方程中，求近似解的的方法有很多种，但求出的结果也不尽相同，但总有一个是最接近精确解的。计算近似解的过程是很麻烦的，但我们运用计算机帮我们计算，就能既快速又准确的计算出结果。matlab就是这么一个软件，正好我们这学期也学了。我们应更近一步的了解matlab，利用它解决更多的数学问题。进而解决生活中的问题。当然这首先要具备良好的数学基础。参考文献1 张圣勤.MATLAB7.0实用教程. 北京: 机械工业出版社，20042 姜健飞、胡良剑、唐俭.数值分析及其MATLAB实验上海：科学出版社，20043 王高雄，周之铭等编.常微分方程第三版.高等教育出版社.20064 刘卫国主编.MATLAB程序设计与应用第二版.高等教育出版社2006
