数学实验第二次作业——常微分方程数值求解

实验4常微分方程数值解

实验目的：

1.练习数值积分的计算；

2.掌握用MATLAB软件求微分方程初值问题数值解的方法；

3.通过实例学习用微分方程模型解决简化的实际问题；

4.了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。

实验内容：

3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。

解答如下：

这是一个典型的牛顿第二定律问题，分析火箭受力情况；

先规定向上受力为正数

建立数学模型：

A燃料未燃尽前，在任意时刻（t<60s）

火箭受到向上的-F=32000N，

向下的重力G=mg，g=9.8，

向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；

此时火箭收到的合力为F1=（F-mg-f）；

火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;

此刻火箭质量为m=1400-18*t

根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=

(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))

由此可利用龙格-库塔方法来实现，程序实现如下

Function [dx]=rocket[t,x] %建立名为rocket的方程

m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值

dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];

%以向量的形式建立方程

[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式

End;

ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点

x0=[0,0]; %输入x的初始值

[t,x]=ode15s(@rocket,ts,x0); %调用ode15s计算

[t,x];

h=x(:,1);

v=x(:,2);

plot(t,x(:,1)),grid; %绘出火箭高度与时间的关系曲线

title('h-t');

xlabel('t/s');ylabel('h/m'),pause;

plot(t,x(:,2)),grid ; %绘出火箭速度与时间的曲线关系

title('v-t');

xlabel('t/s');ylabel('v/m/s'),pause;

a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t); plot(t,a),grid; %绘出火箭加速度与时间的曲线关系

title('a-t');

xlabel('t/s'),ylabel('a/m^2/s'),pause

火箭高度随时间变化的曲线

火箭速度随时间变化的曲线

火箭加速度随时间变化的曲线

数据过多，故截取部分如下

第一列为时间，第二列为火箭高度，第三列为火箭速度

由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2

B燃料燃尽之后，与A 类似，分析受力如下

火箭受到向上的F=0

向下的重力G=mg，g=9.8，

向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；

此时火箭收到的合力为F2=（-mg-f）；

火箭的初始质量为320kg，恒定

根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320;

程序实现如下

function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数

dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程

ts=60:120; %给出初始时刻，估计终点时刻

x0=[12190,267.26]; %给出x初始值

[t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算

[t,x]

plot(t,x(:,1)),grid; %绘出火箭高度随时间变化的曲线

title('h-t');

xlabel('t/s'),ylabel('h/m'),pause;

plot(t,x(:,2)),grid; %绘出火箭速度随时间的变化曲线

title('v-t');

xlabel('t/s'),ylabel('v/m/s'),pause;

v=x(:,2);

a=-9.8-0.4*v.^2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出火箭加速度随时间变化的曲线

title('a-t');

xlabel('t/s'),ylabel('a/m^2/s'),pause

得到的曲线图形如下

火箭高度随时间的变化曲线

从图中可以大致看出，最高点在13km左右，

火箭速度随时间的变化曲线

加速度随时间变化曲线如下

数据表格大致如下

从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为

-9.8m/m/s^2;

本题总结：

这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a

解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定，先用ode15s 速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长

5．一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。已知河

水流速为v1 与船在静水中的速度v2之比为k。

（1）建立描述小船航线的数学模型，求其解析解；

（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任

意时刻小船的位置及航行曲线，作图，并与解析解比较。

（3）若流速v1=0，0.5，1.5，2（m/s）, 情况又如何