微分方程数值解课程设计.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.