大学本科常微分实验报告

福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓 名：系：信息与计算科学专 业：信息与计算科学年 级：2010学 号：指导教师：职 称：讲师2011年 12 月 1 日福建农林大学计算机与信息学院数学类课程实习报告结果评定评语：成绩：指导教师签字：评定日期：1.实习的目的和任务 1.2.实习要求 1.3.实习地点 1.4.主要仪器设备 1.5.实习内容 1.5.1用不同格式对同一个初值问题的数值求解及其分析 .15.1.1 求精确解1.用欧拉法求解3.用改进欧拉法求解 5.用4级4阶龙格库塔法求解 7.问题讨论与分析 105.2 一个算法不同不长求解同一个初值问题及其分析 .135.3洛伦茨方程模拟混沌现象 186. 结束语 2.1.参考文献 2.1.常微分方程课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB来计算微分方程(组)的数值解；了解常微分方程数值解。任务：通过具体的问题，利用MATLAB件来计算问题的结果，分析问题的结论。2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用 所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。3. 实习地点南 2#4254. 主要仪器设备计算机Microsoft Win dows 7Matlab R2009a5. 实习内容5.1用欧拉方法，改进欧拉方法，4阶龙格一库塔方法分别求下面微分方程 的初值 dy/dx=y*cos(2*x) y(0)=1 x 0，25.1 .1求精确解变量分离方程情形：形如5 = f(x)* g(y)的方程,这里f (x), g(y)分别是x, y的 dx连续函数.如果g(y)=0,我们可将方程改写成 型f(x)dx,这样,变量就”分g(y)离”开来了 ,两边同时积分即可：dy = f(x)dx c,c为任意常数.g(y)用变量分离法可求出其精确为：y=exp(0.5*si n(2*x)程序代码： x=0:0.1:2; y=exp(0.5*si n(2*x) plot(x,y,rs-); Data=x,y 结果及图像： y =Columns 1 through 61.4314 1.52311.0000 1.1044 1.2150 1.3262Columns 7 through 121.5936 1.6368 1.6484 1.62731.57561.4982Columns 13 through 180.9712 0.88011.4018 1.2940 1.1823 1.0731Columns 19 through 210.80150.7364 0.6850Data =0 1.00000.10001.10440.20001.21500.30001.32620.40001.43140.50001.52310.60001.59360.70001.63680.80001.64840.90001.62731.00001.57561.10001.49821.20001.40181.30001.29401.40001.18231.50001.07311.60000.97121.7000 0.88011.8000 0.80151.9000 0.73642.00000.68502.1100 20 40 60 8112141,61 B 2用欧拉法求解设常微分方程的初始问题= f(x,y) dxy(xJ = yo有唯一解。则由欧拉法求初值问题(Tn*+hf (召，yn) n = 0,1,2iy=y(x)程序如下：建立函数文件cwfl.mfun cti on x,y=cwf1(fu n, x_spa n,y0,h)x=x_spa n(1):h:x_spa n(2);y(i)=yo；for n=1:le ngth(x)-1y( n+1)=y( n)+h*feval(fu n, x( n),y( n)endx=x;y=y;在MATLA输入以下程序： clear all fun=inlin e( y*cos(2*x); x,y=cwf1(fu n,0,2,1,0.1); x,y plot(x,y,r*-)结果及图像：ans =0 1.00000.1000 1.10000.20001.20780.30001.31910.40001.42790.50001.52740.60001.60990.70001.66830.80001.69660.90001.69171.00001.65321.10001.58441.20001.49121.30001.38121.40001.26291.50001.14391.60001.0306(1)(2)1)，(2)的数值解的计算公式为:h = xn 1 - Xn)1.7000 0.92781.8000 0.83811.9000 0.7629用改进欧拉法求解：计算公式为：yn+(0)=yn +hf(Xn，yn)预报格式J1yn 1yn 2 h f (Xn，yn) f (Xn 1, yn V)校正格式当然也可以迭代多次：yn/)= yn +hf (Xn，yn)预报格式jynP廿=yn +；h【f(Xn，yn) + f (Xn卅,n出)校正格式程序如下：建立函数文件cwf2.mfun cti on x,y=cwf2(fu n, x_spa n,y0,h)x=x_spa n(1):h:x_spa n(2);y(1)=y0;for n=1:length(x)-1k1=feval(fu n,x( n),y( n);y( n+1)=y( n)+h*k1;k2=feval(fu n,x( n+1),y( n+1);y(n+1)=y(n)+h*(k1+k2)/2;endx=x;y=y;在MATLA输入以下程序： clear all fun=inline( y*cos(2*x); x,y=cwf2(fun,0,2,1,0.1); x,y plot(x,y,b+-) 结果及图像： ans =0 1.00000.10001.10390.20001.21380.30001.32440.40001.42900.50001.52010.60001.59020.70001.63300.80001.64450.90001.62341.00001.57201.10001.49491.20001.39911.30001.29201.40001.18101.50001.07241.60000.97111.70000.88031.80000.80211.90000.73732.0000 0.6859200.20.40.60.311.21.4161.02用4阶龙格一库塔求解 四级四阶的龙格一库塔的一般算式= yn *(& +C2K2 +C3K3 +C4K4)*K!=f(Xn,yn)K2 = f (Xn +a2h,yn +hb2iKi)K3 = f (Xn +a3h,yn +hb3iKi +hb32K2)3= f (Xn +a4h, yn +hb4iKi +hb42K2+b43hK3)公式的截断误差阶为o(h5)。ia2 二 a3 ：经典四级四阶龙格一库塔格式取定2，则得ynyn -(Ki 2K2 2K3 K4)6K1 = f (Xn, yn)hh心=f (Xn +：, yn 十二 KJhhK2K3 = f (Xn + , yn +)22a = f (Xn +h, yn +hK3)这是最为著名的经典四级四阶龙格一库塔格式。 程序如下： 建立函数文件cwf3.mfun cti on x,y=cwf3(fu n, x_spa n,yO,h) x=x_spa n(1):h:x_spa n(2);y(1)=y0;for n=1:le ngth(x)-1k1=feval(fu n,x( n),y( n);k2=feval(fu n,x( n)+h/2,y( n)+h/2*k1);k3=feval(fu n,x( n)+h/2,y( n)+h/2*k2);k4=feval(fu n, x( n+1),y( n)+h*k3);y( n+1)=y( n)+h*(k1+2*k2+2*k3+k4)/6;end x=x;y=y;在MATLA输入以下程序： clear all; fun=inlin e( y*cos(2*x); x,y=cwf3(fu n,0,2,1,0.1); x,y plot(x,y, b+-)结果及其图象：ans =0 1.00000.1000 1.10440.2000 1.21500.3000 1.32620.4000 1.43140.5000 1.52310.6000 1.59360.7000 1.63680.8000 1.64840.9000 1.62731.0000 1.57561.1000 1.49821.2000 1.40181.3000 1.29401.4000 1.18231.5000 1.07311.6000 0.97121.7000 0.88011.8000 0.80151.9000 0.73642.0000 0.68501 800.20.40.60.811.21.41 61.82问题讨论与分析由以上数值分析结果绘制表格:精确解欧拉方法改进的欧拉方法四阶龙格-库塔方法xiyiyi误差yi误差yi误差011010100.11.10441.10.00441.10390.00051.104400.21.2151.20780.00721.21380.00121.21500.31.32621.31910.00711.32440.00181.326200.41.43141.42790.00351.4290.00241.431400.51.52311.5274-0.00431.52010.0031.523100.61.59361.6099-0.01631.59020.00341.593600.71.63681.6683-0.03151.6330.00381.636800.81.64841.6966-0.04821.64450.00391.648400.91.62731.6917-0.06441.62340.00391.627301.01.57561.6532-0.07761.5720.00361.575601.11.49821.5844-0.08621.49490.00331.498201.21.40181.4912-0.08941.39910.00271.401801.31.2941.3812-0.08721.2920.0021.29401.41.18231.2629-0.08061.1810.00131.182301.51.07311.1439-0.07081.07240.00071.073101.60.97121.0306-0.05940.97110.00010.971201.70.88010.9278-0.04770.8803-0.00020.880101.80.80150.8381-0.03660.8021-0.00060.801501.90.73640.7629-0.02650.7373-0.00090.7364020.6850.7026-0.01760.6859-0.00090.6850在MATLA输入以下程序： x=0:0.1:2; y仁1.0000 1.1044 1.2150 1.3262 1.4314 1.5231 1.5936 1.6368 1.6484 1.6273 1.5756 1.4982 1.4018 1.2940 1.1823 1.0731 0.9712 0.8801 0.8015 0.7364 0.6850; y2=1.0000 1.1000 1.2078 1.3191 1.4279 1.5274 1.6099 1.6683 1.6966 1.6917 1.6532 1.5844 1.4912 1.3812 1.2629 1.1439 1.0306 0.9278 0.8381 0.7629 0.7026; y3=1.0000 1.1039 1.2138 1.3244 1.4290 1.5201 1.5902 1.6330 1.64451.6234 1.5720 1.4949 1.3991 1.2920 1.1810 1.0724 0.9711 0.8803 0.8021 0.7373 0.6859; y4=1.0000 1.1044 1.2150 1.3262 1.4314 1.5231 1.5936 1.6368 1.64841.6273 1.5756 1.4982 1.4018 1.2940 1.1823 1.0731 0.9712 0.8801 0.8015 0.7364 0.6850; plot(x,y1,r+-) hold on ,plot(x,y2,b-) plot(x,y1,r+-) hold on ,plot(x,y3,b-) plot(x,y1,r+-) hold on ,plot(x,y4,b-)结果如下：精确解与欧拉方法的比较精确解与4阶龙格一库塔方法比较结果分析：由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格一库塔方法误差相对较小，并且龙格一库塔方法误差最小且大部分值都跟精确值相同。由欧拉图与精确图相比可清晰看到，随着 x的增加，函数值与精确值的偏差越来越大。5.2选择用欧拉方法取不同步长分别求下面微分方程的初值dy/dx=y*cos(2*x )y(0)=1，x 0，2。程序如下：建立函数文件cwfl.mfun cti on x,y=cwf1(fu n, x_spa n,y0,h)x=x_spa n(1):h:x_spa n(2);y(i)=yo；for n=1:le ngth(x)-1y( n+1)=y( n)+h*feval(fu n,x( n),y( n);endx=x;y=y;当h=0.05时在MATLA输入以下程序： clear all fun=inlin e( y*cos(2*x); x,y=cwf1(fu n,0,2,1,0.05); x,y plot(x,y,r*-) 结果及图像： ans =0 1.00000.0500 1.05000.1000 1.10220.1500 1.15630.2000 1.21150.2500 1.26730.3000 1.32290.3500 1.37750.4000 1.43010.4500 1.48000.5000 1.52600.5500 1.56720.6000 1.60270.6500 1.63180.7000 1.65360.7500 1.66770.8000 1.67350.8500 1.67110.9000 1.66030.9500 1.64151.0000 1.61491.0500 1.58131.1000 1.54141.1500 1.49611.2000 1.44621.2500 1.39291.30001.33711.35001.27981.40001.22201.45001.16441.50001.10791.55001.05301.60001.00041.65000.95051.70000.90361.75000.85991.80000.81961.85000.78291.90000.74971.95000.72002.00000.6939当h=0.1时在MATLA输入以下程序： clear all fun=inline( y*cos(2*x) ); x,y=cwf1(fun,0,2,1,0.1); x,y plot(x,y,r*-) 结果及图像： ans =0 1.00000.10001.10000.20001.20780.30001.31910.40001.42790.50001.52740.60001.60990.70001.66830.80001.69660.90001.69171.00001.65321.10001.58441.20001.49121.30001.38121.40001.26291.50001.14391.60001.03061.70000.92781.80000.83811.90000.76292.0000 0.7026在MATLA输入以下程序： clear all fun=inlin e( y*cos(2*x); x,y=cwf1(fu n,0,2,1,0.5); x,y plot(x,y,r*-)结果及图像：ans =0 1.00000.50001.50001.00001.90521.50001.50882.00000.7619根据以上的结果和图像与精确解的比较可知, 越接近精确解。步长越小，用欧拉方法所求的值就5.3数值求解Lorenz方程，模拟混沌现象洛伦茨方程如下:dx dt10( - X + V)(lydldzXV - 8z/3将x，y，z表示y，y(2)，y(3)，即列向量y中三个分量 程序如下：(1) 建立自定义函数用edit命令建立自定义函数名为lorenz.m的内容为:fun cti on dy=lore nz(t,y) dy=zeros(3,1); % 建立一个三维列向量dy(1)=10*(-y(1)+y(2);dy(2)=28*y(1)-y(2)-y(1)*y(3);dy(3)=y(1)*y(2)-8/3*y(3)(2) 用edit命令建立一个命令文件xian.m的内容为t,y=ode45( lorenz ,0 , 30,12,2,9);%表示在030秒内求解，在零时刻设y(1)为12, y(2)为2, y(3)为9plot(t,y(:,1) % 显示y( 1) pause % 暂停 plot(t,y(:,2) %显示y (2)pauseplot(t,y(:,3) %显示y (3)pauseplot3(y(:,1),y(:,2),y(:,3); % view(20,42);%(3) 求解结果即x与时间的关系图即y与时间的关系图即z与时间的关系图显示三分量的关系图设定一个较好的观察角度在matlab窗口中执行xian.m文件，数据偏多就不显示了，其图像如下:201610-2011111051015202530y( 1)即x与时间的关系图y (2)即y与时间的关系图4511II1叶丿V1-5 11IIiii051015202530y (3)即z与时间的关系图显示三分量的关系图(4) 用edit命令建立一个命令文件 wea.m的内容为t,y=ode45( lorenz ,0 ， 30,3,7,18);%表示在030秒内求解，在零时刻设y(1)为12, y(2)为2, y(3)为9.plot(t,y(:,1) % pause%plot(t,y(:,2) % pause plot(t,y(:,3) %显示y (1)即x与时间的关系图 暂停显示y (2)即y与时间的关系图显示y (3)即z与时间的关系图pauseplot3(y(:,1),y(:,2),y(:,3); % view(20,42);%在matlab窗口中执行wea.m文件,显示三分量的关系图设定一个较好的观察角度 其图像如下：显示三分量的关系图结果分析：当参数特定选定时，初值稍有变化，函数图像的绕法就会发生变 化，也就是说微分系统对初值的敏感度相当的大，这就形成了所谓的“混沌”。6结束语通过这段期间的实习，我基本掌握了 MATLA语言的一些基本操作，特别是绘 图功能及编程能力由于有些知识以前没接触过，通过实习提高了我的临时自学 能力，同时提高了自己抽象思维的能力和从数学模型中理解问题所在。在这次的实习中，让我明白了网络的方便与简洁， 对于不懂的问题可以及时得到解决， 同 时我也知道了学校的图书馆原来有那么多的学习资源，以后要多去图书馆吸收知识，充实自己。参考文献：1 张圣勤.MATLAB7.C实用教程.北京：机械工业出版社，20042 姜健飞、胡良剑、唐俭数值分析及其MATLA实验.上海：科学出版社， 20043 许波，刘征.MATLAB工程数学应用M.北京：清华大学出版社，20004 刘华杰.分形艺术（电子版M.长沙：湖南电子音像出版社，1997