福建农林大学常微分课程论文.doc

福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓 名：XXX系：应用数学专 业：数学与应用数学年 级：2010学 号：102260002XXX指导教师：陈永雪职 称：讲师2011 年 12 月 1 日福建农林大学计算机与信息学院数学类课程实习报告结果评定评语：成绩：指导教师签字：评定日期：19目 录1.实习的目的和任务12.实习要求13.实习地点14.主要仪器设备15.实习内容15.1数值求解Lorenz方程，模拟混沌现象15.1.1数值求解方程25.1.2数值解对初值的敏感性35.1.3模拟混沌现象55.2用不同格式对同一个初值问题的数值求解及其分析75.2.1数值求解75.2.2 用欧拉法求解.95.2.3 用改进欧拉法求解.105.2.4用4阶龙格库塔求解.115.2.5问题讨论与分析.135.3一个算法不同步长求初值问题及其分析.155.3.1取步长.155.3.2分析和结论.196.结束语19参考文献19常微分方程课程实习1.实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。2.实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。3.实习地点学生宿舍4.主要仪器设备计算机、Microsoft Windows XPMatlab 7.65.实习内容 5.1 数值求解Lorenz方程，模拟混沌现象Lorenz方程：，其中a=10,b=8/3,c=28. 5.1.1 数值求解方程 程序代码： 建立函数文件 function xdot = lorenzeq(t,x) a=10;r=28;b=8/3; xdot=a*(x(2)-x(1); r*x(1)-x(2)-x(1)*x(2); x(1)*x(2)-b*x(3);调用函数 global a b c a=10;b=8/3;c=28;%定义参数lorenz=(t,x)a*(x(2)-x(1);c*x(1)-x(1)*x(3)-x( 2);x(1)*x(2)-b*x(3); % 定义函数 T,X=ode45(lorenz,0,20,0;1;2);%数值法解微分方程 Data=T,X plot3(X(:,1),X(:,2),X(:,3),m)%绘图 view(-20,60);%设置视角 xlabel(x);ylabel(y);zlabel(z);%标记坐标轴Data = 0 0 1.0000 2.0000 0.0001 0.0012 0.9999 1.9994 0.0002 0.0025 0.9998 1.9987 0.0004 0.0037 0.9996 1.9980 0.0005 0.0050 0.9995 1.9974. 19.9735 -10.0007 0.2594 38.1095 19.9801 -9.3215 0.8813 37.4048 19.9867 -8.6489 1.4130 36.6808 19.9934 -7.9881 1.8611 35.947920.0000 -7.3435 2.2331 35.2142运行则得到参数a=10，b=8/3，c=28自变量t取值为0,20的Lorenz方程图像（图1）和数据，图中出现一个吸引子，整个图像由螺线型轨道构成。图1 5.1.2 数值解对初值的敏感性将初值x(0)=0稍微改变下，取x(0)=0.001和x(0)=0进行比较程序代码：global a b ca=10;b=8/3;c=28;%定义参数lorenz=(t,x)a*(x(2)-x(1);c*x(1)-x(1)*x(3)-x(2);x(1)*x(2)-b*x(3);%定义函数 subplot(2,2,1);T,X=ode45(lorenz,0,100,0;1;2);%数值法解微分方程plot(T,X(:,1),m-)%绘图xlabel(t);ylabel(x);%标记坐标轴hold onT,X=ode45(lorenz,0,100,0.001;1;2);plot(T,X(:,1),b)xlabel(t);ylabel(x);hold offsubplot(2,2,2);T,X=ode45(lorenz,0,100,0;1;2);plot(T,X(:,2),m)xlabel(t);ylabel(y);hold onT,X=ode45(lorenz,0,100,0.001;1;2);plot(T,X(:,2),b)xlabel(t);ylabel(y);hold offsubplot(2,2,3);T,X=ode45(lorenz,0,100,0;1;2);plot(T,X(:,3),m)xlabel(t);ylabel(z);hold onT,X=ode45(lorenz,0,100,0.001;1;2);plot(T,X(:,3),b)xlabel(t);ylabel(z);hold off050100-20-1001020txX数值解的波动 050100-40-2002040tyy 数值解的波动 0501000204060tzz 数值解的波动 x0=0x0=0.001图2 由图2可知：当初始条件x(0)=0稍微改变为x(0)=0.001其他两个初始条件不变时，Lorenz方程的三个自变量的数值解变化很大，同理可得另外两个初始条件稍微变化后的结果，这里不详细讨论，从而可以简单说明Lorenz方程的数值解对初值很敏感。 5.1.3 模拟混沌现象 将t取不同值分别绘出不同值的函数图像 程序代码： global a b c a=10;b=8/3;c=28;: lorenz=(t,x)a*(x(2)-x(1);c*x(1)-x(1)*x(3)-x(2);x(1)*x(2)-b*x(3); subplot(2,2,1); T,X=ode45(lorenz,0,10,2;3;4); plot3(X(:,1),X(:,2),X(:,3),m) xlabel(x);ylabel(y);zlabel(z); hold on subplot(2,2,2); T,X=ode45(lorenz,0,20,2;3;4); plot3(X(:,1),X(:,2),X(:,3),m) xlabel(x);ylabel(y);zlabel(z); hold on subplot(2,2,3); T,X=ode45(lorenz,0,100,2;3;4); plot3(X(:,1),X(:,2),X(:,3),m) xlabel(x);ylabel(y);zlabel(z); hold on subplot(2,2,4); T,X=ode45(lorenz,0,500,2;3;4); plot3(X(:,1),X(:,2),X(:,3),m) xlabel(x);ylabel(y);zlabel(z); hold off图3 由图3可知：方程有两个平衡点，当初值和参数不变时，随着计算次数的增加，相空间曲线既不趋近于平衡点1，也不趋近于平衡点2，始终在两个平衡点一定范围内无限反复旋转，此时方程进入混沌状态，不管t值取多大，方程依然处于混沌状态，根本无法预测绕平衡点的旋转次数，从而最终得出Lorenz方程对初值非常敏感。 5.2 用欧拉方法，改进欧拉方法，4阶龙格库塔方法分别求下面微分方程的 初 值dy/dx=(1+y2)*sin(x) y(0)=1 x-1.2，1.2 5.2.1求精确解 首先可以求得其精确解为： y= cot(cos(x)-1+1/4*pi) 由matlab可得函数y的图像（图4），由图知该函数在区间上有间断点，但可以取其中一个连续区间-1.2,1.2进行研究。图4程序代码：x=-1.2:0.1:1.2; y= cot(cos(x)-1+1/4*pi); plot(x,y,b*-); Data=x,y Data = -1.2000 6.7186 -1.1000 4.1042 -1.0000 2.9610 -0.9000 2.3198 -0.8000 1.9110 -0.7000 1.6302 -0.6000 1.4285 -0.5000 1.2806 -0.4000 1.1718 -0.3000 1.0936 -0.2000 1.0407 -0.1000 1.0100 0 1.0000 0.1000 1.0100 0.2000 1.0407 0.3000 1.0936 0.4000 1.1718 0.5000 1.2806 0.6000 1.4285 0.7000 1.6302 0.8000 1.9110 0.9000 2.3198 1.0000 2.9610 1.1000 4.1042 1.2000 6.7186 图5 5.2.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输入以下程序： clear all fun=inline(1+y2)*sin(x) ); x,y=cwfa1(fun,-1.2,1.2,7,0.1); x,y plot(x,y,r*-) 图6 5.2.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输入以下程序： clear all fun=inline( (1+y2)*sin(x) ); x,y=cwfa2(fun,-1.2,1.2,7,0.1); x,y plot(x,y,k+-)图7 5.2.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输入以下程序： clear all; fun=inline( (1+y2)*sin(x); x,y=cwfa3(fun,-1.2,1.2,7,0.1); x,y plot(x,y, g+-) 图8 5.2.5 问题讨论与分析 由以上数值分析结果绘制表格：精确解欧拉方法改进的欧拉方法四阶龙格-库塔方法xiyiyi误差yi误差yi误差-1.26.71867-0.28147-0.28147-0.2814-1.14.10422.33981.76444.3814-0.27724.2039-0.0997-12.9611.76281.19823.159-0.1983.015-0.054-0.92.31981.41720.90262.4622-0.14242.3549-0.0351-0.81.9111.18150.72952.019-0.1081.9365-0.0255-0.71.63021.00960.62061.7165-0.08631.6502-0.02-0.61.42850.87950.5491.5008-0.07231.4452-0.0167-0.51.28060.77940.50121.3434-0.06281.295-0.0144-0.41.17180.70230.46951.2283-0.05651.1848-0.013-0.31.09360.64420.44941.1457-0.05211.1056-0.012-0.21.04070.60240.43831.0901-0.04941.0521-0.0114-0.11.010.57530.43471.0579-0.04791.0211-0.0111010.5620.4381.0473-0.04731.0109-0.01090.11.010.5620.4481.0578-0.04781.0211-0.01110.21.04070.57510.46561.0899-0.04921.0521-0.01140.31.09360.60160.4921.1454-0.05181.1056-0.0120.41.17180.64180.531.2277-0.05591.1848-0.0130.51.28060.69680.58381.3426-0.0621.295-0.01440.61.42850.7680.66051.4996-0.07111.4452-0.01670.71.63020.85780.77241.7147-0.08451.6502-0.020.81.9110.96960.94142.0165-0.10551.9365-0.02550.92.31981.10881.2112.4592-0.13942.3548-0.03512.9611.28341.67763.1589-0.19793.0148-0.05381.14.10421.50622.5984.4081-0.30394.2027-0.09851.26.71861.79754.92117.1733-0.45476.9689-0.2503将三种方法图像与精确解作比较程序代码：x=-1.2:0.1:1.2;y=cot(cos(x)-1+1/4*pi);plot(x,y,b*-);hold onfun=inline( (1+y2)*sin(x) ); x,y=cwfa2(fun,-1.2,1.2,7,0.1);plot(x,y,k+-)hold onfun=inline(1+y2)*sin(x) );x,y=cwfa1(fun,-1.2,1.2,7,0.1);plot(x,y,rx-)hold onfun=inline( (1+y2)*sin(x);x,y=cwfa3(fun,-1.2,1.2,7,0.1); plot(x,y, g.-)hold off图9 分析和结论 由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格库塔方法误差相对较小，并且龙格库塔方法误差最小且几乎都跟精确值相同。由欧拉图与精确图相比可清晰看到，随着x由0向两边增加，函数值与精确值的偏差越来越大。5.3 用改进欧拉方法取不同步长分别求下面微分方程的初值dy/dx=(1+y2)*sin(x)， y(0)=1 x-1.2，1.2，并对结果进行分析说明，给出你的结论。 5.3.1 取步长h分别取0.12、0.14、0.16、0.18、0.2程序代码：x=-1.2:0.1:1.2;y=cot(cos(x)-1+1/4*pi);plot(x,y,c-);hold on;fun=inline( (1+y2)*sin(x) ); x,y=cwfa2(fun,-1.2,1.2,7,0.12);x,yplot(x,y,g*-)hold onx,y=cwfa2(fun,-1.2,1.2,7,0.14);x,yplot(x,y,r+:)hold onx,y=cwfa2(fun,-1.2,1.2,7,0.16);x,yplot(x,y,b*-)hold onx,y=cwfa2(fun,-1.2,1.2,7,0.18);x,yplot(x,y,k+:)hold onx,y=cwfa2(fun,-1.2,1.2,7,0.2);x,yplot(x,y,m*-)hold off-1.5-1-0.500.511.5012345678xy改进欧拉法取不同步长与精确解的比较 h=0.14h=精确解h=0.12h=0.16h=0.18h=0.2图10-0.15-0.1-0.0500.050.10.150.50.60.70.80.911.11.21.31.4xy h=0.14精确解h=0.12h=0.16h=0.18h=0.2图110.40.450.50.550.60.650.70.7511.11.21.31.41.51.61.71.81.9xy h=0.14精确解h=0.12h=0.16h=0.18h=0.2图12-1.1-1-0.9-0.8-0.7-0.6-0.5-0.41.41.61.822.22.42.62.833.