微分方程数值解课程设计Word格式.docx
- 文档编号:14258073
- 上传时间:2022-10-21
- 格式:DOCX
- 页数:14
- 大小:21.26KB
微分方程数值解课程设计Word格式.docx
《微分方程数值解课程设计Word格式.docx》由会员分享,可在线阅读,更多相关《微分方程数值解课程设计Word格式.docx(14页珍藏版)》请在冰豆网上搜索。
年级专业:
计算科学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
_Toc293560842"
总结……………………………………………………………….……...17
_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));
Un=U;
zz(k)=sqrt(1+2*T(k));
e=U-zz;
function[Un,e]=pmecme()
%[Un,e]=pmecme()
symstu
U(k)=U(k)+251/270*(U(k)-U(k-1));
U(k)=U(k)-19/270*(U(k)-U(k-1));
运行结果及图形显示:
运行得:
>
[Un,e]=pece()
Un=
Columns1through8
1.00001.09541.18321.26491.34161.41421.48321.5492
Columns9through11
1.61251.67331.7321
e=
1.0e-05*
00.04170.07890.11640.05710.02710.01270.0042
-0.0013-0.0054-0.0088
pmecme
ans=
1.00001.09541.18321.26491.33981.41041.47731.5409
1.60161.65951.7147
编写a.m用于画图显示计算结果,
t0=0:
0.001:
t1=0:
U0=sqrt(1+2*t0);
U2=pece();
U1=pmecme();
holdon
plot(t0,U0,'
k'
)
plot(t1,U1,'
--r'
'
Marker'
x'
plot(t1,U2,'
--g'
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:
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));
plot(x,y1),grid
但如果步长选得不一样(横坐标都从开始,到100左右结束,但中间选的点也不太一样),即使采用同样的迭代公式,绘出的曲线却有很大不同:
程序:
0.01*pi:
30*pi];
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));
我们认为,可能是因为欧拉公式每算一次都会产生误差,如果取的点正好位置不太合适,会导致误差一步步增大,累加起来,就像蝴蝶效应一样,会产生和真实状况截然不同的结果。
这也是用数值方法求解方程最怕出现的问题,也是应该解决的问题。
这说明向前欧拉公式并不是一种很好的计算方法,误差较大(只有1阶精度),而且容易失真。
这一点在下面还要说明。
(2)龙格—库塔方法
编写如下M文件:
functiondx=fangchengzu(t,x)%建立名为fangchengzu的函数M文件
dx=[x
(2);
-x
(2)/t-(1-0.25/t^2)*x
(1)];
%表示方程组
用4级5阶龙格—库塔公式进行计算。
ts=[pi/2:
0.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 课程设计
![提示](https://static.bdocx.com/images/bang_tan.gif)