数学实验报告利用MALTAB计算常微分方程数值解Word文档下载推荐.docx
- 文档编号:13250890
- 上传时间:2022-10-08
- 格式:DOCX
- 页数:43
- 大小:411.96KB
数学实验报告利用MALTAB计算常微分方程数值解Word文档下载推荐.docx
《数学实验报告利用MALTAB计算常微分方程数值解Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《数学实验报告利用MALTAB计算常微分方程数值解Word文档下载推荐.docx(43页珍藏版)》请在冰豆网上搜索。
将整个过程中的时间向量以及速度向量联合起来,利用第三章所学插值与数值积分的方法可以求得任意时刻火箭的近似高度。
2、方法
求解常微分方程时,我分别采用了自己编写的欧拉公式、改进欧拉公式、4级4阶龙格-库塔公式,以及MATLAB自带的龙格-库塔方法,求解数值积分时采用辛普森公式。
由于Matlab自带的Simpson公式是自适应的,因此需要使用自己在上一次实验时所编的Simpson公式。
㈢结果与分析
1、各种公式的对比
首先,我作出了各种不同公式计算得到的火箭速度随时间变化的图像,图如下:
从图中可以看出,各种公式计算得到的结果基本一致,为确定其区别,将图像放大,放大约2000倍后,得到下图:
分析:
从图中可以看出,自编欧拉公式距离MATLAB自带龙格-库塔公式最远,精度最差;
自编的改进欧拉公式和自编的龙格-库塔公式结果基本一致,两者中自编龙格-库塔公式距MATLAB自带龙格-库塔公式的结果稍近。
与之前的分析基本一致。
然而产生自编龙格-库塔公式与MATLAB自带龙格-库塔公式之间的差距的原因还未知。
由于MATLAB自带龙格-库塔公式精度较高,因此以下各项实验不再计算其它几项公式的结果。
2、第一阶段
火箭关闭瞬间的速度
关闭前瞬间的加速度
此时火箭的高度
3、第二阶段
由初始速度以及常微分方程可以求得火箭达到最高点的时间约为71.3s;
加速度
4、整体过程
下图为火箭加速度与时间的关系。
由图可以看出,火箭一开始的加速度很大,随着时间的推移,火箭的燃料有所减少,与此同时速度有所上升。
两者中前者使火箭的加速度增大,后者使其减小,综合作用,最终体现为加速度先略有上升,然后慢慢减小。
当火箭中燃料燃尽时,火箭丧失推动力,因而加速度急剧减小为负值。
此后火箭速度不断减小,导致火箭所受阻力逐渐减小,因而加速度的绝对值有所减小,直到最终火箭速度降为零,火箭不受阻力,仅受重力,加速度为重力加速度。
下图为火箭速度与时间的关系:
此图可以看作是由第一幅图对时间积分所得结果,本图的斜率对应第一幅图的值。
第一阶段火箭加速度为正,因此速度不断增加,只是增加的速度不断减慢。
燃料燃尽后,加速度变为负值,因此速度开始急剧下降,与此同时下降的速率不断减小。
最终火箭速度降为0。
下图为火箭高度与实践的关系:
此图可以看作是由第二幅图对时间积分所得结果,本图的斜率对应第一幅图的值。
一开始或减速度较小,因此高度缓慢增加,之后增加的速度不断提升,直到火箭的燃料耗尽,此后速度不断减小,但仍为正值,因此火箭继续向上飞行,只是高度提升的速度逐渐变慢。
知道最终火箭速度降为零。
3、总结
这是自己接触数学实验以来第一次解决相对复杂的实际问题,总体上还是有些不适应。
虽然最终得以成功解决,但是编程的逻辑以及脚本之间的层次没有考虑清楚,需要在之后的学习中不断加强相关能力。
当然也有不小收获:
(1)了解了MATLAB解决常微分方程组的原理与步骤;
(2)进一步熟悉了MATLAB的相关操作,包括画图、调试、编写脚本等;
(3)加强了对之前插值与数值积分的理解,学习将其应用到实际问题的解决中。
㈣程序清单
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)
F=0;
a=0;
-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)
y1=y(i,:
dy1=feval(fun,ts(i)+h,y1);
)+h/2*(rot90(dy+dy1));
5、自编龙格-库塔公式
functiony=myRungeKutta(fun,ts,x0)
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);
)+h/6*rot90(k1+2*k2+2*k3+k4);
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:
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;
y=[y1;
y2(1:
i-1,:
)];
y0=[y01;
y02(1:
y1=[y11;
y12(1:
y2=[y21;
y22(1:
j=numel(y)/numel(y00);
plot(t(1:
j),y(:
1),'
k'
ts(1:
j),y0(:
b'
j),y1(:
y'
j),y2(:
r'
);
gridon;
legend('
MATLAB自带龙格-库塔公式'
'
自编的欧拉公式'
自编的改进欧拉公式'
自编的龙格-库塔公式'
7、自启动脚本2
e=0.01;
y0=[0,1080];
[t1,y1]=ode45(@fly,ts1,y0);
y0=[y1((60/e+1),1),0];
[t2,y2]=ode45(@fly1,ts2,y0);
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
h=h(1:
j);
t=t(1:
plot(t,h,'
二、小船过河
一只小船渡过宽度为d的河流,目标是起点A正对着的另一岸B点。
已知水流速
与船在静水中的速度
之比为k。
建立小船航线的数学模型,求解析解以及数值解。
1、数学模型与解析解
由题图,以B为原点,BA为y轴负半轴建立直角坐标系,设任意时刻小船坐标为(x,y),则易得下列公式
此微分方程组可用于作数值计算。
为求解析解,采用极坐标系分析,设小船到原点的距离为r,与x轴正半轴的夹角为θ
两式相除,分离变量,积分,得
当r=d时,
,代入得
此式即为解析解,经过坐标变换可以得到x,y的关系。
2、数值解
利用上一小节得到的微分方程组,通过龙革-库塔方法即可得到方程组的数值解,此外还可以与解析解所得数值进行比较。
3、各种其他情况
当流速为0,0.5,1.5,2时,可以现根据解析解进行判断,然后再进行数值计算。
1、d=100m,
(1)使用(x
(1)+0.00001计算
函数M文件如下:
functiondz=Cross(t,x)
v1=1;
v2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 实验 报告 利用 MALTAB 计算 微分方程 数值