微分方程数值解课程设计.docx
《微分方程数值解课程设计.docx》由会员分享，可在线阅读，更多相关《微分方程数值解课程设计.docx（14页珍藏版）》请在冰豆网上搜索。
微分方程数值解课程设计
资料范本
本资料为word版本，可以直接编辑和打印，感谢您的下载
微分方程数值解课程设计
地点：
__________________
时间：
__________________
说明：
本资料适用于约定双方经过谈判，协商而共同承认，共同遵守的责任与义务，仅供参考，文档可直接下载或修改，不需要的部分可直接删除，使用时请详细阅读内容
微分方程数值方法课程设计
BasicTheoryofOrdinaryDifferentialEquationsExperimentReport
教务处
2014年7月
课程设计说明书
题目：
常微分方程数值解法课程设计
学院（系）：
理学院
年级专业：
计算科学11-1
学生姓名：
指导教师：
教师职称：
教授
燕山大学课程设计（论文）任务书
院（系）：
理学院教学单位：
说明：
此表一式四份，学生、指导教师、基层教学单位、系部各一份。
年月日
燕山大学课程设计评审意见表
摘要
本文对常微分方程初值问题现有的数值解法进行了综述研究。
主要讨论了几种常用的数值解法：
即欧拉法，改进欧拉法，龙格库塔方法，阿达姆斯PECE,PMECME格式等。
文章最后结合常见数值解法，对较为典型的微分方程模型进行数值求解，探讨了上述数值算法在实际建模问题中的应用。
论文阐述的是常微分方程数值解法的几个问题，通过对以下问题的求解
一.比较Adams四阶PECE模式和PMECME模式。
二.求解贝塞尔方程并与精确解比较。
三.小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
来加强对用数值解法解常微分方程实际问题的能力。
关键词：
常微分方程数值解MATLAB
TOC\o"1-3"\h\zHYPERLINK\l"_Toc293560829"摘要PAGEREF_Toc293560829\hi
HYPERLINK\l"_Toc293560830"绪论PAGEREF_Toc293560830\h1
问题解答…...………………………………………………………………2
HYPERLINK\l"_Toc293560842"总结……………………………………………………………….……...17
HYPERLINK\l"_Toc293560843"参考文献资料……………………………………..…………………..…17
绪论
很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。
建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。
如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。
因此，研究微分方程求解的数值方法是非常有意义的。
问题解答
问题一：
比较Adams四阶PECE模式和PMECME模式。
问题引出：
将阿达姆斯方法显式与隐式方法作一对比，以说明预测——校正格式的必要性。
这些方法的阶及误差常数列表如下：
由于阿达姆斯内插公式是隐式公式，故用它计算时也需用迭代法。
通常把阿达姆斯外插公式与内插公式结合起来使用，先由前者提供初值，再由后者进行修正，所以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;
fork=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）;
end
fork=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）;
end
Un=U;
fork=1:
11
zz（k）=sqrt（1+2*T（k））;
end
e=U-zz;
end
function[Un,e]=pmecme（）
%[Un,e]=pmecme（）
symstu
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;
fork=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）;
end
fork=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）;
end
Un=U;
fork=1:
11
zz（k）=sqrt（1+2*T（k））;
end
e=U-zz;
end
运行结果及图形显示：
运行得：
>>[Un,e]=pece（）
Un=
Columns1through8
1.00001.09541.18321.26491.34161.41421.48321.5492
Columns9through11
1.61251.67331.7321
e=
1.0e-05*
Columns1through8
00.04170.07890.11640.05710.02710.01270.0042
Columns9through11
-0.0013-0.0054-0.0088
pmecme
ans=
Columns1through8
1.00001.09541.18321.26491.33981.41041.47731.5409
Columns9through11
1.60161.65951.7147
编写a.m用于画图显示计算结果，
t0=0:
0.001:
1;
t1=0:
0.1:
1;
U0=sqrt（1+2*t0）;
U2=pece（）;
U1=pmecme（）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'--r','Marker','x'）
plot（t1,U2,'--g','Marker','x'）
legend（'U=sqrt（1+2*t）','U1',’U2’'Location','NorthWest'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下：
放大如下图
可知PMECME模式比PECE模式更为精确
问题二：
求解贝塞尔方程并与精确解比较。
求解x2y’’+xy’+（x2-n2）y=0
y=2,y’=-（贝赛尔方程，令n=0.5），精确解y=sinx
解：
首先将高阶方程装降阶化简为一阶常微分方程组
令y1=y’，令y2=y1’=y’’
将n=0.5代入，则原方程转化为：
y1’=y2
y1=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为步长
编程如下：
clearall
x=[pi/2:
0.1:
pi/2+100-0.05];
y1=2;
y2=-2/pi;
fori=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））;
end
plot（x,y1）,grid
但如果步长选得不一样（横坐标都从开始，到100左右结束，但中间选的点也不太一样），即使采用同样的迭代公式，绘出的曲线却有很大不同：
程序：
clearall
x=[pi/2:
0.01*pi:
30*pi];
y1=2;
y2=-2/pi;
fori=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））;
end
plot（x,y1）,grid
我们认为，可能是因为欧拉公式每算一次都会产生误差，如果取的点正好位置不太合适，会导致误差一步步增大，累加起来，就像蝴蝶效应一样，会产生和真实状况截然不同的结果。
这也是用数值方法求解方程最怕出现的问题，也是应该解决的问题。
这说明向前欧拉公式并不是一种很好的计算方法，误差较大（只有1阶精度），而且容易失真。
这一点在下面还要说明。
（2）龙格—库塔方法
y1’=y2
y1=2,y2=-
编写如下M文件：
functiondx=fangchengzu（t,x）%建立名为fangchengzu的函数M文件
dx=[x
（2）;-x
（2）/t-（1-0.25/t^2）*x
（1）];%表示方程组
用4级5阶龙格—库塔公式进行计算。
编程如下：
clearall
ts=[pi/2:
0.1:
100];
%使用龙格—库塔方法
x0=[2,-2/pi];
[t,x]=ode45（@fangchengzu,ts,x0）;
plot（t,x（:
1））,grid,title（'龙格-库塔方法'）,pause
%绘出精确解y=sinx图像
t=[pi/2:
0.1:
100];
y=sin（t）.*sqrt（2*pi./t）;
plot（t,y）,grid,title（'精确解'）
绘出曲线如下：
比较分析：
从曲线上可以看出，在数值求解贝赛尔方程时，龙格—库塔方法（4级5阶龙格—库塔公式）的结果与精确解y=sinx非常接近，产生震荡，并且随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-kv^2-M0+mg）/（M0+m）QUOTEv=F-kv2-M0+mgM0+m
m=-aQUOTEm=-a
初值QUOTEv=0,m=1080。
由以上各式可以求出t=QUOTET1T1时火箭的速度。
再求解第二阶段：
V=（-kv^2-M0g）/M0QUOTEv=-kv2-M0gM0
m=0QUOTEm=0
可以求出火箭速度降为0的时刻。
将整个过程中的时间向量以及速度向量联合起来，利用第三章所学插值与数值积分的方法可以求得任意时刻火箭的近似高度。
2、方法
求解常微分方程时，我分别采用了自己编写的欧拉公式、改进欧拉公式、4级4阶龙格-库塔公式，以及MATLAB自带的龙格-库塔方法，
结果与分析
1、各种公式的对比
首先，我作出了各种不同公式计算得到的火箭速度随时间变化的图像，图如下：
从图中可以看出，各种公式计算得到的结果基本一致，为确定其区别，将图像放大，放大约2000倍后，得到下图：
分析：
从图中可以看出，自编欧拉公式距离MATLAB自带龙格-库塔公式最远，精度最差；自编的改进欧拉公式和自编的龙格-库塔公式结果基本一致，两者中自编龙格-库塔公式距MATLAB自带龙格-库塔公式的结果稍近。
与之前的分析基本一致。
然而产生自编龙格-库塔公式与MATLAB自带龙格-库塔公式之间的差距的原因还未知。
由于MATLAB自带龙格-库塔公式精度较高，因此以下不在计算其它几项公式的结果。
2、第一阶段
火箭关闭瞬间的速度QUOTEv=267.261240773261m/sv=267.261240773261m/s
关闭前瞬间的加速度a=（F-kv^2-M0g）/M0=0.914286475421320m/s^2QUOTEv=-kv2-M0gM0
此时火箭的高度h=12189.7791507305m
3、第二阶段
由初始速度以及常微分方程可以求得火箭达到最高点的时间约为71.3s；
此时火箭的高度QUOTEh=13115.7148739887mh=13115.748739887m
加速度a=-9.80000406052111m/s^2
4、整体过程
下图为火箭加速度与时间的关系。
由图可以看出，火箭一开始的加速度很大，随着时间的推移，火箭的燃料有所减少，与此同时速度有所上升。
两者中前者使火箭的加速度增大，后者使其减小，综合作用，最终体现为加速度先略有上升，然后慢慢减小。
当火箭中燃料燃尽时，火箭丧失推动力，因而加速度急剧减小为负值。
此后火箭速度不断减小，导致火箭所受阻力逐渐减小，因而加速度的绝对值有所减小，直到最终火箭速度降为零，火箭不受阻力，仅受重力，加速度为重力加速度。
下图为火箭速度与时间的关系：
此图可以看作是由第一幅图对时间积分所得结果，本图的斜率对应第一幅图的值。
第一阶段火箭加速度为正，因此速度不断增加，只是增加的速度不断减慢。
燃料燃尽后，加速度变为负值，因此速度开始急剧下降，与此同时下降的速率不断减小。
最终火箭速度降为0。
下图为火箭高度与时间的关系：
此图可以看作是由第二幅图对时间积分所得结果，本图的斜率对应第一幅图的值。
一开始或减速度较小，因此高度缓慢增加，之后增加的速度不断提升，直到火箭的燃料耗尽，此后速度不断减小，但仍为正值，因此火箭继续向上飞行，只是高度提升的速度逐渐变慢。
知道最终火箭速度降为零。
程序清单
1、第一阶段的微分方程组
functiondx=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、第二阶段微分方程组
functiondx=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、自编欧拉公式
functiony=Euler（fun,ts,x0）
NumOfLine=numel（ts）;
NumOfRow=numel（x0）;
h=ts
（2）-ts
（1）;
y=zeros（NumOfLine,NumOfRow）;
y（1,:
）=x0;
fori=1:
NumOfLine-1
dy=feval（fun,ts（i）,y（i,:
））;
y（i+1,:
）=y（i,:
）+rot90（dy）*h;
end
4、自编改进欧拉公式
functiony=ImprovedEuler（fun,ts,x0）
NumOfLine=numel（ts）;
NumOfRow=numel（x0）;
h=ts
（2）-ts
（1）;
y=zeros（NumOfLine,NumOfRow）;
y（1,:
）=x0;
fori=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））;
end
5、自编龙格-库塔公式
functiony=myRungeKutta（fun,ts,x0）
NumOfLine=numel（ts）;
NumOfRow=numel（x0）;
h=ts
（2）-ts
（1）;
y=zeros（NumOfLine,NumOfRow）;
y（1,:
）=x0;
fori=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）;
end
6、“自启动”脚本1
e=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;
end
y=[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）,y1（:
1）,'y',ts（1:
j）,y2（:
1）,'r'）;
gridon;
legend（'MATLAB自带龙格-库塔公式','自编的欧拉公式','自编的改进欧拉公式','自编的龙格-库塔公式'）;
7、自启动脚本2
e=0.01;
ts1=0:
e:
60;
y0=[0,1080];
[t1,y1]=ode45（@fly,ts1,y0）;
ts2=60:
e:
80;
y0=[y1（（60/e+1）,1）,0];
[t2,y2]=ode45（@fly1,ts2,y0）;
ts=[ts1,ts2];
t=[t1;t2];
i=1;
while（y2（i,1）>0）
i=i+1;
end
y=[y1;y2（1:
i-1,:
）];
y（1,2）=（32000-（320+1080）*9.8）/（320+1080）;
forj=2:
60/e+i
h（j）=mySimpson（t（1:
j）,y（1:
j,1））;
if（j<60/e+1）
y（j,2）=（-0.4*y（j,1）^2+32000-（320+y（j,2））*9.8）/（320+y（j,2））;
else
y（j,2）=（-0.4*y（j,1）^2-320*9.8）/320;
end
end
h=h（1:
j）;
t=t（1:
j）;
plot（t,h,'k'）;
gridon;
总结
本文对常微分方程初值问题现有的数值解法进行了综述研究。
主要讨论了几种常用的数值解法：
即欧拉法，改进欧拉法，龙格库塔方法，阿达
升级会员 