微分方程数值解课程设计

微分方程数值方法课程设计Basic Theory of Ordinary Differential Equations Experiment Report教 务 处2014年 7 月课 程 设 计 说 明 书题目： 常微分方程数值解法课程设计 学院（系）： 理学院 年级专业： 计算科学11-1 学生姓名： 指导教师： 教师职称： 教授 燕山大学课程设计（论文）任务书院（系）：理学院 教学单位： 学 号XXX学生姓名XXX专业（班级）11计算数学设计题目常微分方程数值解法课程设计设计内容一. 比较Adams四阶PECE模式和PMECME模式。二. 求解贝塞尔方程并与精确解比较。三. 小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。设计要求1、要通过课题背景的介绍，阐明选题的意义；2、运用所学知识，建立问题的数学模型；3、利用Matlab软件进行计算，并给出计算机程序框图及计算机程序；4、运用所学知识，对计算结果进行分析；工作量1、程设计要在两周内完成；2、根据已知题目找出解题方法，然后利用Matlab进行计算；3、正文在3000字左右。工作计划第一天：收集资料，阅读参考文献；第二天：构建数学模型；第三天：设计计算机流程图，编制计算机程序，上机运算，调试程序；第四天：计算结果分析第五天：撰写报告等参考资料1胡建伟汤怀民。微分方程数值方法。科学出版社。1999年01月2 陈渝周璐钱方等。数值方法（MATLAB版）。电子工业出版社。2002.63尹泽明，丁春丽等。精通MATLAB 6。清华大学出版社。2006.64Cleve B.Moler。MATLAB数值计算。机械工业出版社2006.6指导教师签字基层教学单位主任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。年 月 日 燕山大学课程设计评审意见表指导教师评语：成绩： 指导教师： 年 月 日答辩小组评语：成绩： 评阅人： 年 月 日课程设计总成绩：答辩小组成员签字：年 月 日摘 要本文对常微分方程初值问题现有的数值解法进行了综述研究。主要讨论了几种常用的数值解法：即欧拉法，改进欧拉法，龙格库塔方法，阿达姆斯PECE,PMECME格式等。文章最后结合常见数值解法，对较为典型的微分方程模型进行数值求解，探讨了上述数值算法在实际建模问题中的应用。论文阐述的是常微分方程数值解法的几个问题，通过对以下问题的求解一.比较Adams四阶PECE模式和PMECME模式。二.求解贝塞尔方程并与精确解比较。三.小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。来加强对用数值解法解常微分方程实际问题的能力。关键词：常微分方程数值解 MATLAB17目 录摘 要i绪 论1问题解答.2总结.17参考文献资料.171 绪论绪 论很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。因此，研究微分方程求解的数值方法是非常有意义的。燕 山 大 学 课 程 设 计 说 明 书问题解答问题一：比较Adams四阶PECE模式和PMECME模式。问题引出：将阿达姆斯方法显式与隐式方法作一对比，以说明预测校正格式的必要性。这些方法的阶及误差常数列表如下：步数阶误差常数显式隐式显式隐式112223334445由于阿达姆斯内插公式是隐式公式，故用它计算时也需用迭代法。通常把阿达姆斯外插公式与内插公式结合起来使用，先由前者提供初值，再由后者进行修正，所以Adams预测-校正格式既利用了隐式方法较好的稳定性及精确性，有利用了显式方法的简易性，正是把两者结合起来，做到取长补短。Adams四阶预估-校正(PECE)公式： 而PMECME模式公式为对于初值问题u=u-2t/u,u(0)=1,在单区间【0,1】上，用Adams四阶预估-校正算法的PECE模式及PMECME模式求其数值解，取步长h=0.1利用计算结果估计数值解的局部误差主项。（真解u（t）=sqrt(1+2t)）编写matlab程序：function Un,e = pece()% Un,e = pece()f0=u-2*t/u;v=t,u;U=zeros(1,11);T=0:0.1:1;f=zeros(1,11);h=0.1;U(1)=1;f(1)=1;for k=1:3 K1=subs(f0,v,(k-1)*h,U(k); K2=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K1); K3=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K2); K4=subs(f0,v,(k-1)*h+h,U(k)+h*K3); U(k+1)=U(k)+1/6*h*(K1+2*K2+2*K3+K4); f(k+1)=U(k+1)-2*T(k+1)/U(k+1);endfor k=5:11 U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4); f(k)=U(k)-2*T(k)/U(k); U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3); f(k)=U(k)-2*T(k)/U(k);endUn=U;for k=1:11 zz(k)=sqrt(1+2*T(k);ende=U-zz;end function Un,e = pmecme()% Un,e = pmecme()syms t uf0=u-2*t/u;v=t,u;U=zeros(1,11);T=0:0.1:1;f=zeros(1,11);h=0.1;U(1)=1;f(1)=1;for k=1:3 K1=subs(f0,v,(k-1)*h,U(k); K2=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K1); K3=subs(f0,v,(k-1)*h+1/2*h,U(k)+1/2*h*K2); K4=subs(f0,v,(k-1)*h+h,U(k)+h*K3); U(k+1)=U(k)+1/6*h*(K1+2*K2+2*K3+K4); f(k+1)=U(k+1)-2*T(k+1)/U(k+1);endfor k=5:11 U(k)=U(k-1)+h/24*(55*f(k-1)-59*f(k-2)+37*f(k-3)-9*f(k-4); U(k)=U(k)+251/270*(U(k)-U(k-1); f(k)=U(k)-2*T(k)/U(k); U(k)=U(k-1)+h/24*(9*f(k)+19*f(k-1)-5*f(k-2)+f(k-3); U(k)=U(k)-19/270*(U(k)-U(k-1); f(k)=U(k)-2*T(k)/U(k);endUn=U;for k=1:11 zz(k)=sqrt(1+2*T(k);ende=U-zz;end 运行结果及图形显示：运行得： Un,e=pece()Un = Columns 1 through 8 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 Columns 9 through 11 1.6125 1.6733 1.7321e = 1.0e-05 * Columns 1 through 8 0 0.0417 0.0789 0.1164 0.0571 0.0271 0.0127 0.0042 Columns 9 through 11 -0.0013 -0.0054 -0.0088pmecmeans = Columns 1 through 8 1.0000 1.0954 1.1832 1.2649 1.3398 1.4104 1.4773 1.5409 Columns 9 through 11 1.6016 1.6595 1.7147编写a.m用于画图显示计算结果，t0=0:0.001:1;t1=0:0.1:1;U0=sqrt(1+2*t0);U2=pece();U1=pmecme();hold onplot(t0,U0,k)plot(t1,U1,-r,Marker,x)plot(t1,U2,-g,Marker,x)legend(U=sqrt(1+2*t),U1,U2Location,NorthWest)grid onxlable(t)ylable(Un)hold off运行结果如下：放大如下图可知PMECME模式比PECE模式更为精确问题二：求解贝塞尔方程并与精确解比较。求解 x2y+xy+(x2-n2)y=0 y=2,y=- (贝赛尔方程，令n=0.5)，精确解y=sin x解：首先将高阶方程装降阶化简为一阶常微分方程组令y1=y，令y2=y1=y将n=0.5代入，则原方程转化为：y1=y2y1=2, y2=-（1）用向前欧拉公式：y1(n+1)=y1(n)+h*y2(n)y2(n+1)=y2(n)+h*y1(0) =2，y2(0)= -，x(n+1)= +n*h，h为步长编程如下：clear allx=pi/2:0.1:pi/2+100-0.05;y1=2;y2=-2/pi;for i=1:999 y1(i+1)=y1(i)+0.1*y2(i); y2(i+1)=y2(i)-0.1*(y2(i)/x(i)+(1-0.25/x(i)2)*y1(i);endplot(x,y1),grid但如果步长选得不一样(横坐标都从开始，到100左右结束，但中间选的点也不太一样)，即使采用同样的迭代公式，绘出的曲线却有很大不同：程序：clear allx=pi/2:0.01*pi:30*pi;y1=2;y2=-2/pi;for i=1:2950 y1(i+1)=y1(i)+0.01*pi*y2(i); y2(i+1)=y2(i)-0.01*pi*(y2(i)/x(i)+(1-0.25/x(i)2)*y1(i);endplot(x,y1),grid我们认为，可能是因为欧拉公式每算一次都会产生误差，如果取的点正好位置不太合适，会导致误差一步步增大，累加起来，就像蝴蝶效应一样，会产生和真实状况截然不同的结果。这也是用数值方法求解方程最怕出现的问题，也是应该解决的问题。这说明向前欧拉公式并不是一种很好的计算方法，误差较大（只有1阶精度），而且容易失真。这一点在下面还要说明。（2）龙格库塔方法y1=y2y1=2, y2=-编写如下M文件：function dx=fangchengzu(t,x) %建立名为fangchengzu的函数M文件dx=x(2);-x(2)/t-(1-0.25/t2)*x(1); %表示方程组用4级5阶龙格库塔公式进行计算。编程如下：clear allts=pi/2:0.1:100;%使用龙格库塔方法x0=2,-2/pi;t,x=ode45(fangchengzu,ts,x0);plot(t,x(:,1),grid,title(龙格-库塔方法),pause%绘出精确解y=sin x图像t=pi/2:0.1:100;y=sin(t).*sqrt(2*pi./t);plot(t,y),grid,title(精确解)绘出曲线如下：比较分析：从曲线上可以看出，在数值求解贝赛尔方程时，龙格库塔方法（4级5阶龙格库塔公式）的结果与精确解y=sin x非常接近，产生震荡，并且随x的增大振幅逐渐减小；而用欧拉方法（向前欧拉公式）的计算结果却与精确解有较大差异，虽然也产生了震荡，但在x稍大的情况下，振幅随x的增大而增大。这是因为向前欧拉公式只有1阶的精度，而4级5阶龙格库塔公式有5阶的精度，一般情况下，肯定用后者计算更接近精确解。而且从本例可以看出，向前欧拉公式并不是一种很好的计算方法，误差较大，甚至导致结果失真，所以还是应该用龙格库塔方法比较好。问题三：小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭，设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度及火箭达到最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。1、简要分析本题的求解需要用到常微分方程，而整个过程又被分为两个阶段：火箭加速上升阶段和燃料燃尽后减速的阶段。由题目易知第一个阶段持续时间T1=60s列出第一阶段的方程组：设M0为火箭本身质量，m为燃料质量，g为重力加速度 = 9.8m/s2,燃料燃烧率为a，空气阻力的比例系数为k，F为推进力。M0 = 1400-1080 = 320kg; V=(F-kv2-M0+mg)/(M0+m) m=-a初值v=0,m=1080。由以上各式可以求出t=T1时火箭的速度。再求解第二阶段： V=(-kv2-M0g)/M0 m=0可以求出火箭速度降为0的时刻。将整个过程中的时间向量以及速度向量联合起来，利用第三章所学插值与数值积分的方法可以求得任意时刻火箭的近似高度。2、方法求解常微分方程时，我分别采用了自己编写的欧拉公式、改进欧拉公式、4级4阶龙格-库塔公式，以及MATLAB自带的龙格-库塔方法，结果与分析1、各种公式的对比首先，我作出了各种不同公式计算得到的火箭速度随时间变化的图像，图如下：从图中可以看出，各种公式计算得到的结果基本一致，为确定其区别，将图像放大，放大约2000 倍后，得到下图：分析：从图中可以看出，自编欧拉公式距离MATLAB自带龙格-库塔公式最远，精度最差；自编的改进欧拉公式和自编的龙格-库塔公式结果基本一致，两者中自编龙格-库塔公式距MATLAB自带龙格-库塔公式的结果稍近。与之前的分析基本一致。然而产生自编龙格-库塔公式与MATLAB自带龙格-库塔公式之间的差距的原因还未知。由于MATLAB自带龙格-库塔公式精度较高，因此以下不在计算其它几项公式的结果。2、第一阶段火箭关闭瞬间的速度v=267.1m/s关闭前瞬间的加速度a=(F -kv2-M0g)/M0=0.1320m/s2此时火箭的高度h=12189.m3、第二阶段由初始速度以及常微分方程可以求得火箭达到最高点的时间约为71.3s；此时火箭的高度h=13115.m加速度a=-9.111m/s24、整体过程下图为火箭加速度与时间的关系。由图可以看出，火箭一开始的加速度很大，随着时间的推移，火箭的燃料有所减少，与此同时速度有所上升。两者中前者使火箭的加速度增大，后者使其减小，综合作用，最终体现为加速度先略有上升，然后慢慢减小。当火箭中燃料燃尽时，火箭丧失推动力，因而加速度急剧减小为负值。此后火箭速度不断减小，导致火箭所受阻力逐渐减小，因而加速度的绝对值有所减小，直到最终火箭速度降为零，火箭不受阻力，仅受重力，加速度为重力加速度。 下图为火箭速度与时间的关系： 此图可以看作是由第一幅图对时间积分所得结果，本图的斜率对应第一幅图的值。第一阶段火箭加速度为正，因此速度不断增加，只是增加的速度不断减慢。燃料燃尽后，加速度变为负值，因此速度开始急剧下降，与此同时下降的速率不断减小。最终火箭速度降为0。 下图为火箭高度与时间的关系：此图可以看作是由第二幅图对时间积分所得结果，本图的斜率对应第一幅图的值。一开始或减速度较小，因此高度缓慢增加，之后增加的速度不断提升，直到火箭的燃料耗尽，此后速度不断减小，但仍为正值，因此火箭继续向上飞行，只是高度提升的速度逐渐变慢。知道最终火箭速度降为零。程序清单1、第一阶段的微分方程组function dx = fly(t,y)M0 = 320;a = 18; k = 0.4;g = 9.8;F = 32000;dx = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; - a*sign(y(2);2、第二阶段微分方程组function dx = fly1(t,y)M0 = 320; k = 0.4;g = 9.8;F = 0;a = 0;dx = (-k*y(1)2 + F - (M0 + y(2)*g)/(M0+y(2) ; -a;3、自编欧拉公式function y = Euler(fun,ts,x0)NumOfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = x0;for i = 1:NumOfLine-1 dy = feval(fun,ts(i),y(i,:); y(i+1,:)=y(i,:)+rot90(dy)*h;end4、自编改进欧拉公式function y = ImprovedEuler(fun,ts,x0)NumOfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = x0;for i = 1:NumOfLine-1 dy = feval(fun,ts(i),y(i,:); y1 = y(i,:)+rot90(dy)*h; dy1 = feval(fun,ts(i)+h,y1); y(i+1,:)=y(i,:)+h/2*(rot90(dy+dy1);end5、自编龙格-库塔公式function y = myRungeKutta(fun,ts,x0)NumOfLine = numel(ts);NumOfRow = numel(x0);h = ts(2)-ts(1);y = zeros(NumOfLine,NumOfRow);y(1,:) = x0;for i = 1:NumOfLine-1 k1 = feval(fun,ts(i),y(i,:); dy1 = y(i,:)+h/2*rot90(k1); k2 = feval(fun,ts(i)+h/2,dy1); dy2 = y(i,:)+h/2*rot90(k2); k3 = feval(fun,ts(i)+h/2,dy2); dy3 = y(i,:)+h*rot90(k3); k4 = feval(fun,ts(i)+h,dy3); y(i+1,:)=y(i,:)+h/6*rot90(k1+2*k2+2*k3+k4);end6、“自启动”脚本1e =0.1;ts1 = 0 : e : 60;y00 = 0,1080;y01 = Euler(fly,ts1,y00);y11 = ImprovedEuler(fly,ts1,y00);y21 = myRungeKutta(fly,ts1,y00);t1,y1 = ode45(fly,ts1,y00);ts2 = 60 : e :80;y00 = y1(60/e+1),1),0;y000 = y01(60/e+1),1),0;y001 = y11(60/e+1),1),0;y002 = y21(60/e+1),1),0;y02 = Euler(fly1,ts2,y000);y12 = ImprovedEuler(fly1,ts1,y001);y22 = myRungeKutta(fly1,ts1,y002);t2,y2 = ode45(fly1,ts2,y00);ts = ts1,ts2; t = t1;t2;i=1;while(y2(i,1)0) i = i+1;endy = y1;y2(1:i-1,:);y0 = y01;y02(1:i-1,:);y1 = y11;y12(1:i-1,:);y2 = y21;y22(1:i-1,:);j = numel(y)/numel(y00);plot(t(1:j),y(:,1),k,ts(1:j),y0(:,1),b,ts(1:j),y