数学建模预备知识——微分方程

微分方程建模

微分方程是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。 把实际问题化为求解微分方程的定解问题，大体上可以分为以下步骤

• 根据实际要求确定要研究的量（自变量，未知函数，必要参数）并确定坐标系
• 找出这些量所满足的基本规律
• 运用这些规律列出方程和定解条件

列方程常见的方法有以下几种

1. 按规律直接列方程。
2. 微元分析法与任意区域上取积分的方法
3. 模拟近似法

发射卫星为什么要用三级火箭

采用运载火箭把人造卫星发射到高空轨道上运行，为什么不能采用一级火箭而必须用多级火箭系统？

火箭是一个复杂的系统，为了使问题简单明了，只从动力系统和整体结构上分析，并假设引擎是足够强大的。

为什么不能用以及火箭发射人造卫星

下面用三个数学模型回答

1. 卫星进入 $600km$ 高空轨道的时候，火箭必需的最低速度

   首先将问题理想化，并假设：

   （1）卫星轨道是以地球中心为圆心的某个平面上的圆周，卫星在此轨道上以地球引力作为向心力围绕地球做平面匀速圆周运动。

   （2）地球是固定于空间中的一个均匀球体，其质量集中于球心

   （3）其他星球对卫星的引力忽略不计

   建模与求解：

   设地球半径为 $R$ ，质量为 $M$ ，卫星轨道半径为 $r$ ，卫星质量为 $m$ 。

   根据假设（2）和（3），卫星只受地球引力，由牛顿万有引力定律可知，其引力大小为
   $$F=\frac{GMm}{r^2}\text{\hspace{1em}(1)}$$
   式中：$G$ 为引力常数

   为了消去常数 $G$ ，可以利用 $mg=\frac{GMm}{r^2}$

   经过化简后得到
   $$F=mg(\frac{R}{r}{)}^2\text{\hspace{1em}(2)}$$
   式中： $g=9.81(m/s^2)$ 为重力加速度

   根据假设（1），若卫星围绕地球做匀速圆周运动的速度为 $v$ ，则其向心力为 $m{\omega }^2/r$ ，因为卫星所受的地球引力就是它做匀速运动的向心力，因此
   $$mg=(\frac{R}{r}{)}^2=\frac{m{v}^2}{r}$$
   由此推得卫星据地面距离为 $(r-R)km$ ，必需得最低速度的数学模型为
   $$v=R\sqrt{\frac{g}{r}}\text{\hspace{1em}(3)}$$
   代入数据：$R=6400km,r-R=600km$ 得到
   $$v\approx 7.6km/s$$
   即要把卫星送入离地面 $600km$ 的轨道，火箭的末速度最低为 $7.6km/s$

2. 火箭推进力及升空速度

   火箭的简单模型是由一台发动机和一个燃料舱组成的。燃料燃烧产生大量的气体从火箭末端喷出，给火箭一个向前的推进力。火箭飞行要受到地球引力、空气阻力、地球自转与公转等影响，使火箭升空后做曲线运动。

   为使问题简化，假设：

   （1）火箭在喷气推动下做直线运动，火箭所受的重力与空气阻力忽略不计

   （2）在$t$ 时刻火箭的质量为 $m(t)$ ，速度为 $v(t)$，且均为时间 $t$ 的连续可微函数

   （3）从火箭末端喷出的气体的速度（相对于火箭本身）为常数 $u$

   建模与分析：

   由于火箭在运动的过程中不断喷出气体，使其质量不断减小，在 $(t,t+\Delta t)$ 内的减少量可由泰勒展开式表示为
   $$m(t+\Delta t)-m(t)=\frac{dm}{dt}\Delta t+o(\Delta t)\text{\hspace{1em}(4)}$$
   因为喷出的气体相对于地球的速度为 $v(t)-u$ ，则由动量守恒定律，有
   $$m(t)v(t)=m(t+\Delta t)v(t+\Delta t)-[\frac{dm}{dt}\Delta t+o(t)][v(t)-u]\text{\hspace{1em}(5)}$$
   从式（4）和（5）可以得到火箭推进力的数学模型为
   $$m\frac{dv}{dt}=-u\frac{dm}{dt}\text{\hspace{1em}(6)}$$
   令 $t=0$ 时，$v(0)=v_0,m(0)=m_0$，求式（6），得到火箭升空速度模型为
   $$v(t)=v_0+u\ln \frac{m_0}{m(t)}\text{\hspace{1em}(7)}$$
   式（6）表明火箭所受推力等于燃料的消耗速度与喷气速度的成绩。式（7）表明，在 $v_0,m_0$ 一定的条件下，升空速度 $v(t)$ 由喷气速度 $u$ 以及质量比 $\frac{m_0}{m(t)}$ 决定。这为提高火箭速度找到了正确的途径：从燃料上设法提高 $u$ 的值；从结构上设法减少 $m(t)$。

3. 一级火箭末速度上限

   火箭-卫星系统的质量可以分为三部分：$m_p$（有效负载），$m_F$（燃料质量），$m_s$（结构质量）。

   一级火箭末速度上限主要是受目前技术条件的限制，假设：

   （1）目前技术条件为：相对火箭的喷气速度为 $u=3km/s$ 及
   $$\frac{m_s}{m_F+m_s}\ge \frac{1}{9}$$
   （2）初速度 $v_0$ 忽略不计，即 $v_0=0$。

   建模与求解：

   因为升空火箭的最终质量为 $m_p+m_s$，所以由（7）以及假设（2）得到末速度为
   $$v=u\ln \frac{m_0}{m_p+m_s}\text{\hspace{1em}(8)}$$
   令 $m_s=\lambda (m_F+m_s)=\lambda (m_0-m_p)$，带入式（8）可以得到
   $$v=u\ln \frac{m_0}{\lambda m_0+(1-\lambda )m_p}\text{\hspace{1em}(9)}$$
   于是，当卫星脱离火箭的时候，即 $m_p=0$ 时，便得到火箭末速度上限的数学模型为
   $$v^0=u\ln \frac{1}{\lambda }$$
   由假设（1），取 $u=3km,\lambda =\frac{1}{9}$，便得到火箭的速度上限为
   $$v^0=3\ln 9\approx 6.6km/s$$

理想火箭模型

从前面对问题的假设和分析可以看出，火箭推进力自始至终在加速整个火箭，然而随着燃料的不断消耗，所出现的无用结构质量也随之不断加速，做无用功，因而效益低，浪费大。

所谓理想火箭模型，就是能够随着燃料的不断燃烧抛弃火箭的无用结构。下面建立其数学模型

假设在 $(t,t+\Delta t)$ 时段丢弃的结构质量与烧掉的燃料质量以 $\alpha$ 与 $1-\alpha$ 的比例同时进行。

建模与分析：

由动量守恒定律，有
$$m(t)v(t)=m(t+\Delta t)v(t+\Delta t)-\alpha \frac{dm}{dt}\Delta t\cdot v(t)-(1-\alpha )\frac{dm}{dt}\Delta t\cdot [v(t)-u]+o(\Delta t)$$
由上式可以得到理想火箭的数学模型为
$$-m(t)\frac{dv(t)}{dt}=(1-\alpha )\frac{dm}{dt}\cdot u\text{\hspace{1em}(10)}$$
代入初始条件 $v(0)=v_0,m(0)=m_0$

解得
$$v(t)=(1-\alpha )u\ln \frac{m_0}{m(t)}\text{\hspace{1em}(11)}$$
由上式可知，当燃料耗尽，结构质量抛弃完时，便只剩下卫星的质量 $m_p$ ，从而最终速度的数学模型为
$$v(t)=(1-\alpha )u\ln \frac{m_0}{m_p}\text{\hspace{1em}(12)}$$

多级火箭卫星模型

理想火箭时假设把无用结构质量连续抛弃以达到最佳的升空速度，虽然在目前的技术条件下办不到，但是他的确为发展火箭技术指明了奋斗目标。目前已经商业化的多级火箭卫星系统便是朝着这种目标迈进的第一步。多级火箭时自末级开始，逐级燃烧，当第 $i$ 级燃料烧尽的时候，第 $i+1$ 级火箭立即自动点火，并抛弃已经无用的第 $i$ 级。用 $m_i$ 表示第 $i$ 级火箭质量， $m_p$ 表示有效负载。为了简单起见，先做以下假设：

（1）设各级火箭具有相同的 $\lambda$ ，$\lambda m_i$ 表示第 $i$ 级的结构质量，$(1-\lambda )m_i$ 表示第 $i$ 级的燃料质量。

（2）喷气相对火箭的速度 $u$ 相同，燃烧级的初始质量与其负载质量之比保持不变，该比值记为 $k$

先考虑二级火箭。由式（7），当第一级火箭燃烧完时，其速度为
$$v_1=u\ln \frac{m_1+m_2+m_p}{\lambda m_1+m_2+m_p}=u\ln \frac{k+1}{\lambda k+1}$$
在第二级火箭燃烧完时，其速度为
$$v_2=v_1+u\ln \frac{m_2+m_p}{\lambda m_2+m_p}=2u\ln \frac{k+1}{\lambda k+1}\text{\hspace{1em}(13)}$$
仍取 $u=3km/s,\lambda =0.1$ 考虑到阻力等因素，为了达到第一宇宙速度 $7.9km/s$ ，对于耳机火箭，欲使 $V_2=10.5km/s$ ，由式（13）得
$$6\ln \frac{k+1}{0.1k+1}=10.5$$
解得
$$k=11.2$$
这时
$$\frac{m_0}{m_p}=\frac{m_1+m_2+m_p}{m_p}=(k+1{)}^2\approx 149$$
同理可以推出三级火箭
$$v_3=3u\ln \frac{k+1}{\lambda k+1}$$
欲使 $v_3=10.5km/s$ ，应该 $k\approx 3.25$，从而有 $\frac{m_0}{m_p}\approx 77$。

与二级火箭相比，在达到相同效果的情况下，三级火箭的质量几乎节省了一半。

现记 $n$ 级火箭的总质量（包含有效负载 $m_p$）为 $m_0$ ，在相同的假设下（$u=3km/s$,$v_{末}=10.5km/s$，$\lambda =0.1$），可以算出相应的 $\frac{m_0}{m_p}$ 的值，将计算结果列于下表

$n$/级数	1	2	3	4	5	$\cdots$	$\infty$
$\frac{m_0}{m_p}$	x	149	77	65	60	$\cdots$	50

人口模型

Malthus模型

该模型于1789年，由英国神父Malthus提出

模型假设：

（1）设 $x(t)$ 表示 $t$ 时刻的人口数，且 $x(t)$ 连续可微

（2）人口的增长率 $r$ 是常数（增长率=出生率-死亡率）

（3）人口数量的变化是封闭的，即人口数量的增加与减少值取决于人口中个体的生育和死亡，且每一个体都具有同样的生育能力与死亡率

建模与求解：

由假设，$t$ 时刻到 $t+\Delta t$ 时刻人口的增量为
$$x(t+\Delta t)-x(t)=rx(t)\Delta t$$
于是有
$$\{\begin{array}{l}\frac{dx}{dt}=rx\\ x(0)=x_0\end{array}\text{\hspace{1em}(14)}$$
其解为
$$x(t)=x_0{e}^{rt}\text{\hspace{1em}(15)}$$
模型评价：

考虑近百年来的人口增长的实际情况，1961年世界人口总数为 $3.06\times 10^9$，在1961年到1970年这段时间内，每年平均的人口自然增长率为 $2$，则式（15）可以写为
$$x(t)=3.06\times 10^9\cdot e^{0.02(t-1961)}\text{\hspace{1em}(16)}$$
根据1700—1961年间世界人口统计数据，发现这些数据与式（16）的计算结果基本吻合。

但是随着时间的增加，人口数量已知在增加，会达到一个不符合实际的数量，因此，对 $r$ 是常数的假设提出质疑

阻滞增长模型（Logistic模型）

如何对增长率 $r$ 进行修正？在人口比较少的时候，可以把增长率 $r$ 看作一个常数，当人口增加到一定数量之后，就应当视 $r$ 为一个随着人口数量增加而减少的量，即增长率 $r$ 表示为人口 $x(t)$ 的函数，且 $r(x)$ 为 $x$ 的减函数。

模型假设：

1. 设 $r(x)$ 为 $x$ 的线性函数，$r(x)=r-sx$
2. 自然资源与环境条件所能容纳的最大人口数为 $x_m$ ，即当 $x=x_m$ 时，增长率 $r(x_m)=0$

建模与求解：

由假设可知 $r(x)=r(1-\frac{x}{x_m})$ ，于是有
$$\{\begin{array}{l}\frac{dx}{dt}=r(1-\frac{x}{x_m})x\\ x(t_0)=x_0\end{array}\text{\hspace{1em}(17)}$$
此微分方程的解为
$$x(t)=\frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r(t-t_0)}}\text{\hspace{1em}(18)}$$
模型检验：

由式（17）可知
$$\frac{d^{2}x}{d{t}^2}=r^2(1-\frac{x}{x_m})(1-\frac{2x}{x_m})x\text{\hspace{1em}(19)}$$
所以人口总数 $x(t)$ 由如下规律：

1. $\underset{t\to +\infty }{\lim }x(t)=x_m$ ，即无论人口初值 $x_0$ 如何，人口总数都以 $x_m$ 为极限。
2. 当 $0<x_0<x_m$ 时，$\frac{dx}{dt}=r(1-\frac{x}{x_m})x>0$ ，说明 $x(t)$ 是单调增加的。由又式（19）可知，当 $x<\frac{x_m}{2}$ 时，$\frac{d^{2}x}{d{t}^2}>0，x=x(t)$ 为凹函数；当 $x>\frac{x_m}{2}$ 时，$\frac{d^{2}x}{d{t}^2}<0，x=x(t)$ 为凸函数。
3. 人口变化率 $\frac{dx}{dt}$ 在 $x=\frac{x_m}{2}$ 时取得最大值，即人口数达到极限值的一半以前时加速增长状态，超过一半之后，是减速增长状态。

模型推广

可以从另一个角度推导出阻滞增长模型，在 Malthus 模型的基础上增加一个竞争项 $-b{x}^2(b>0)$ ，他的作用是使纯增长率减小。建立方程
$$\{\begin{array}{l}\frac{dx}{dt}=x(a-bx)a,b>0\\ x(t_0)=x_0\end{array}\text{\hspace{1em}(20)}$$
其解为
$$x(t)=\frac{a{x}_0}{b{x}_0+(a-b{x}_0)e^{-a(t-t_0)}}\text{\hspace{1em}(21)}$$
由式（21）可得到
$$\frac{d^{2}x}{d{t}^2}=(a-2bx)(a-bx)x\text{\hspace{1em}(22)}$$
对以上三式进行分析

1. 对任意 $t>t_0$ ，有 $x(t)>0$ ，且 $\underset{t\to +\infty }{\lim }x(t)=\frac{a}{b}$
2. 当 $0<x<\frac{a}{b}$ 时，$\frac{dx}{dt}>0$，$x(t)$ 递增；当 $x>\frac{a}{b}$ 时，$\frac{dx}{dt}<0$，$x(t)$ 递减。
3. 当 $0<x<\frac{a}{2b}$时，$\frac{d^{2}x}{d{t}^2}>0$，$x(t)$ 为凹函数；当 $x>\frac{a}{2b}$时，$\frac{d^{2}x}{d{t}^2}<0$，$x(t)$ 为凸函数。

令式（20）第一个方程右侧等于0，得到 $x_1=0,x_2=\frac{a}{b}$ ，称为微分方程的平衡解。易知 $\underset{t\to +\infty }{\lim }x(t)=\frac{a}{b}$ ，故又称 $\frac{a}{b}$ 为式（20）的稳定平衡解。因此，无论人口初值为多少，经过相当长的时间后，人口总数将维持在 $\frac{a}{b}$ 。

例：美国人口预报模型。利用下表中的数据，建立人口预测模型，并预报2010年的人口数量。

建模与求解

记 $x(t)$ 为第 $t$ 年的人口数量，设人口年增长率 $r(x)$ 为 $x$ 的线性函数，$r(x)=r-sx$ 。自然资源与环境条件所能容纳的最大人口数量为 $x_m$，即当 $x=x_m$ 时，增长率 $r(x_m)=0$，可得 $r(x)=r(1-\frac{x}{x_m})$，建立 Logistixc 人口模型
$$\{\begin{array}{l}\frac{dx}{dt}=r(1-\frac{x}{x_m})x\\ x(t_0)=x_0\end{array}$$
其解为
$$x(t)=\frac{x_m}{1+(\frac{x_m}{x_0}-1)e^{-r(t-t_0)}}\text{\hspace{1em}(23)}$$
参数估计

采用线性最小二乘法估计参数
$$\frac{1}{x}\cdot \frac{dx}{dt}=r-sx,s=\frac{r}{x_m}\text{\hspace{1em}(24)}$$
利用向后差分得到差分方程
$$\frac{1}{x(k)}\cdot \frac{x(k)-x(k-1)}{\Delta t}=r-sx(k)\text{\hspace{1em}(25)}$$
其中，取步长 $\Delta t=10$

下面拟合参数 $r$ 和 $s$ ，具体代码如下

x=[3.9 5.3 7.8 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4]';
t=[1790:10:2000]';
a=[ones(21,1),-x(2:end)];
b=diff(x)./x(2:end)/10;
cs=a\b;
r=cs(1),xm=cs(2);

解得 $x_m=373.5135,r=0.0247$

Matlab求解微分方程的符号解

在Matlab中，符号运算工具箱提供了功能强大的求解常微分方程的运算命令 dsolve

在Matlab中求解常微分方程的符号解，首先要定义符号变量，然后调用命令 dsolve

具体调用格式如下
$$[y1,\cdots ,yN]=dsolve(eqns,conds,Name,Value)$$
其中：eqns 为符号微分方程或者符号微分方程组；conds 为初值条件或者边界条件；Name 和 Value 为可选用的成对参数。

试求解常微分方程的通解

试解常微分方程
$$x^2+y+(x-2y)y^{\prime }=0$$
解

syms y(x)
dsolve(x^2+y+(x-2*y)*diff(y)==0)

求解常微分方程的初值问题

试解常微分方程
$$y^{\prime \prime \prime }-y^{\prime \prime }=x,y(1)8,y(1{)}^{\prime }=7,y(2{)}^{\prime \prime }=4$$
解

syms y(x)
dy=diff(y);
d2y=diff(y,2);
y=dsolve(diff(y,3)-diff(y,2)==x,y(1)==8,dy(1)==7,d2y(4)==4)
y=simplify(y);

求解常微分方程组

试求解常微分方程组
$$\{\begin{array}{l}{f}^{\prime \prime }+3g=\sin x\\ g^{\prime }+f^{\prime }=\cos x\end{array}$$
的通解和在初值边界条件为 $f^{\prime }(2)=0,f(3)=3,g(5)=1$ 的解

解

syms f(x) g(x)
df= diff(f);
[f1,g1]=dsolve(diff(f,2)+3*diff(g)==sin(x),diff(g)+diff(f)==cos(x));
f1=simplify(f1),g1=simplify(g1);
[f2,g2]=dsolve(diff(f,2)+3**diff(g)==sin(x),diff(g)+diff(f)==cos(x),df(2)==0，f(3)==3,g(5)==1);
f2=simplify(f2),g2=simplify(g2);

求解线性常微分方程组

• 一阶齐次线性微分方程组
  $${\mathbf{X}}^{\prime }=\mathbf{A}\mathbf{X},\mathbf{X}=\left[\begin{array}{c}{x}_1\\ ⋮\\ x_n\end{array}\right],\mathbf{A}=\left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1n}\\ ⋮& \ddots & ⋮\\ x_{n1}& \cdots & a_{nn}\end{array}\right]$$
  其解为
  $$\mathbf{X}(\mathbf{t})=e^{\mathbf{A}(t-t_0)}{\mathbf{X}}_0$$

  试解初值问题
  $$X^{\prime }=\left[\begin{array}{ccc}2& 1& 3\\ 0& 2& -1\\ 0& 0& 2\end{array}\right]X,X(0)=\left[\begin{array}{c}1\\ 2\\ 1\end{array}\right]$$
  具体程序如下

  syms x(t),y(t),z(t)
  A=[2 1 3;0 2 -1;0 0 2];
  b=[1;2;1];
  X=[x;y;z];
  [x,y,z]=dsolve(diff(X)==A*X,X(0)==B);
  

• 非齐次线性微分方程组

  有参数变易法可求得初值问题
  $$X^{\prime }=AX+f(t),X(t_0)=X_0$$
  的解为
  $$X(t)=e^{A(t-t_0)}{X}_0+{\int }_{t_0}^t{e}^{A(t-s)}f(s)ds$$

  试解初值问题
  $$X^{\prime }=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& -2\\ 3& 2& 1\end{array}\right]X+\left[\begin{array}{c}0\\ 0\\ e^t\cos 2t\end{array}\right],X(0)=\left[\begin{array}{c}0\\ 1\\ 1\end{array}\right]$$
  具体程序如下

  syms x(t) y(t) z(t)
  A=[1 0 0;2 1 -2;3 2 1];
  B=[0;0;exp(t)*cos(2*t)];
  X=[x;y;z];
  x0=[0;1;1];
  [x,y,z]=dsolve(diff(X)==A*X+B,X(0)==x0);
  

放射性废料的处理

问题的提出

美国原子能委员会以往处理浓缩的放射性废料的方法，一直是把它们装入密封的圆筒例，然后扔到水深为90米的海底。生态学家和科学家担心圆筒下沉到海底时与海底发生碰撞而破裂，从而造成核污染。为此，工程师进行了碰撞试验，发现当圆筒以超过 $12.2m/s$ 的速度与海底相撞的时候，圆筒就可能发生破裂。为了避免圆筒破裂，需要计算一下圆筒沉到海底的速度为多少。已知圆筒的质量 $m=239.46kg$，体积 $V=0.2058m^3$，海水密度 $\rho =1035.71kg/m^3$，若圆筒速度小于 $12.2m/s$ 就说明这种方法时安全可靠的，否则就要禁止使用这种方法。假设水的阻力与速度大小成正比例，其正比例系数为 $k=0.6$。先要求建立合理的数学模型，解决如下问题：

1. 判断这种处理废料的方法时候合理。
2. 一般情况下， $v$ 大，$k$ 也大。当 $v$ 很大的时候，常用 $kv$ 来代替 $k$ ，那么这时的速度与时间的关系如何？并求出当速度不超过 $12.2m/s$ 时，圆筒的运动时间 $t$ 和位移 $s$ 应不超过多少

模型的建立与求解

1. 问题一的模型

   以海平面上的一点为为坐标原点，垂直向下为坐标轴的正向建立坐标系。由于圆桶在运动过程中受到本身的重力以及水的浮力 $H$ 和水的阻力 $f$ 的作用，所以根据牛顿运动定律得到圆桶受到的合力 $F$ 满足
   $$F=G-H-f\text{\hspace{1em}(1)}$$
   又因为 $F=ma=m\frac{dv}{dt}=m\frac{d^{2}v}{d{t}^2},G=mg,H=\rho gV$ 以及 $f=kv=k\frac{ds}{dt}$，所以圆桶的位移和速度分别满足以下微分方程
   $$m\frac{d^{2}s}{d{t}^2}=mg-\rho gV-k\frac{ds}{dt}\text{\hspace{1em}(2)}$$

   $$m\frac{dv}{dt}=mg-\rho gV-kv\text{\hspace{1em}(3)}$$

   由方程（2）加上初始条件 $\frac{ds}{dt}{|}_{t=0}=s{|}_{t=0}=0$，求得位移函数为
   $$s(t)=-171511.0+429.744t+171511.0e^{-0.00250564t}\text{\hspace{1em}(4)}$$
   由方程（3）加上初始条件 $v{|}_{t=0}=0$，求得速度函数为
   $$v(t)=429.744-429.744e^{-0.00250564t}\text{\hspace{1em}(5)}$$
   由 $s(t)=90$，求得圆桶到达水深为 $90m$ 的海底需要时间 $t=12.994s$，再把他带入方程（5）求出圆桶到达海底的速度为 $v=13.7720m/s$ 。

   具体代码如下

   syms m V rho g k s(t) v(t)
   ds=diff(s);
   s=dsolve(m*diff(s,2)==m*g-rho*g*V-k*diff(s),s(0)==0,ds(s)==0);
   s=subs(s,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
   s=simply(s);
   s=vpa(s,6);
   v=dsolve(m*diff(v)==m*g-rho*g*V-k*v,v(0)==0);
   v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
   v=simply(v);
   v=vpa(v,6);
   y=s-90;
   tt=solve(y);
   tt=double(tt);
   vv=subs(v,tt);
   vv=double(vv);
   

2. 问题二的模型

   在假设条件下，圆桶受到的阻力应该改为 $f=k{v}^2=k(\frac{ds}{dt}{)}^2$，类似问题一的模型，可以的到圆桶的速度满足如下微分方程
   $$m\frac{dv}{dt}=mg-\rho gV-k{v}^2\text{\hspace{1em}(6)}$$
   根据方程（6），加上初始条件 $v{|}_{t=0}=0$，求得圆桶的速度为
   $$v(t)=20.7303\tanh (0.0519t)$$
   这时若速度要小于 $12.2m/s$，经计算可得圆桶的运动时间就不能超过 $T=13.0025s$，利用 $s(T)={\int }_0^{T}v(t)dt$，计算得到的位移不能超过 $84.8439m$。

   通过这个模型计算，也可以得到处理核废料的方法是不合理的。

   具体代码如下

   syms m V rho g k v(t)
   v=dsolve(m*diff(v)==m*g-rho*g*V-k*v^2,v(0)==0);
   v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});
   v=simply(v);
   v=vpa(v,6);
   T=solve(v-12.2);
   T=double(T);
   s=int(v,0,T);
   

初值问题的Matlab数值解

Matlab的工具箱提供了几个解常微分方程的功能函数，如ode45、ode23、ode113。其中，ode45 采用四五阶龙格库塔方法（RK方法），是解非刚性微分方程的首选方法；ode23 采用二三阶RK方法；ode113 采用多步法，效率较高。

在化学工程或者自动化控制等领域中，所涉及的常微分方程组初值问题常常是所谓的刚性问题。具体说就是对于一阶微分方程组，有
$$\frac{dy}{dx}=\mathbf{A}y+\Phi (x)\text{\hspace{1em}(1)}$$
式中，$y,\Phi \in {\mathbf{R}}^m$；$\mathbf{A}$ 为 $m$ 阶方阵。

若矩阵 $\mathbf{A}$ 的特征值为 ${\lambda }_i(i=1,2\cdots ,m)$ 满足关系
$$\begin{gathered} Re{\lambda }_i<0,i=1,2,\cdots ,m \\ \underset{1\le i\le m}{\max }|Re{\lambda }_i|\gg \underset{1\le i\le m}{\min }|Re{\lambda }_i| \end{gathered}$$
则称方程式（1）为刚性方程组或者Stiff方程组，称数
$$s=\underset{1\le i\le m}{\max }|Re{\lambda }_i|/\underset{1\le i\le m}{\min }|Re{\lambda }_i|$$
为刚性比。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长 $h$ 不做任何限制。

Matlab工具箱提供了几个求解刚性常微分方程的功能函数，如ode15s、ode23s、ode23t、ode23tb。

ode23、ode45、ode113的使用

对于简单的一阶方程的初值问题
$$\{\begin{array}{l}{y}^{\prime }=f(x,y)\\ y(x_0)=y_0\end{array}$$
Matlab的函数形式为
$$[t,y]=solver{(}^{\prime }{F}^{\prime },tspan,y0)$$
这里solver为ode23、ode45、ode113，输入的参数 F 是用.m文件定义的微分方程 $y^{\prime }=f(x,y)$ 右端的函数。$tspan=[t0,tfinal]$ 为求解的区间。y0为初始值列向量。

注意： 函数 $solver{(}^{\prime }{F}^{\prime },tspan,y0)$ 也可以返回一个值，返回一个值的时候，返回的是结构数组，利用Matlab的命令 deval 和返回的结构数组，可以计算任一点的函数值。

例：用RK方法求解
$$y^{\prime }=-2y+2x^2+2x,0\le x\le 0.5,y(0)=1$$
解：

编写Matlab代码如下

yx=@(x,y) -2*x+2*x^2+2*x;
[x,y]=ode45(yx,[0,0.5],1);
sol=ode45(yx,[0,0.5],1);
y2=deval(sol,x);

高阶微分方程 $y^{(n)}=f(t,y,y^{\prime },\cdots ,y^{(n-1)})$ 和一阶微分方程组的解法

高阶常微分方程，必需做变量替换，化为一阶常微分方程组，才能使用Matlab求数值解。

一阶微分方程组的解法和简单的一阶方程解法是一样的。

例1：考虑初值问题
$$y^{\prime \prime \prime }-3y^{\prime \prime }-y^{\prime }y=0,y(0)=0,y^{\prime }(0)=1,y^{\prime \prime }(0)=-1$$
解：

（1）设 $y_1=y,y_2=y^{\prime },y_3=y^{\prime \prime }$，则有
$$\{\begin{array}{l}{y}_1^{\prime }=y2,y_1(0)=0\\ y_2^{\prime }=y3,y_2(0)=1\\ y_3^{\prime }=3y_3+y_1{y}_2,y_3(0)=-1\end{array}$$
初值问题可以写成 $Y^{\prime }=F(t,Y),Y(0)=Y_0$ 的形式，其中 $Y=[y_1,y_2,y_3{]}^T$。

（2）把一阶方程组写成接受两个参数t和y，返回一个列向量的.m文件

function dy=F(t,y);
dy=[y(2);y(3);3*y(3)+y(1)*y(2)];

（3）利用Matlab求解此问题的函数为

[T,Y]=ode45('F',[0,1],[0;1;-1])；

例2：求解 van der Pol 方程
$$\{\begin{array}{l}{y}^{\prime \prime }-1000(1-y^2)y^{\prime }+y=0\\ y(0)=2\\ y^{\prime }(0)=0\end{array}$$
的数值解，这里是刚性微分方程

解：

（1）先化为常微分方程组，令$y_1=y,y_2=y^{\prime }$，则
$$\{\begin{array}{l}{y}_1^{\prime }=y_2,y_1(0)=2\\ y_2^{\prime }=1000(1-y_1^2)y_2-y_1,y_2(0)=0\end{array}$$
（2）求解数值解，并利用图形输出解的结果。

dy=@(t,y) [y(2);1000*(1-y(1)^2)*y(2)-y(1)];
[t,y]=ode15s(dy,[0 3000],[2;0]);
plot(t,y(:,1),'*');

例3：Lorenz方程是一个三阶的非线性系统，它是由描述大气动力系统的Navier-Stoker偏微分方程演化而来的。自由系统如下：
$$\{\begin{array}{l}\dot{x}=\sigma (y-x)\\ \dot{y}=\beta x-y-xz\\ \dot{z}=-\lambda z+xy\end{array}$$
若系统参数 $\sigma ,\beta ,\lambda$ 在一定范围内，系统就出现混沌，如 $\sigma =10,\beta =28,\lambda =8/3$ 时，出现混沌现象。求在初值条件 $[x(0),y(0),z(0){]}^T=[5,13,17{]}^T$ 时，方程组的数值解，并画出解的图像

解：

sigma=10,beta=28,lambda=8/3;
f=@(t,Y)[sigma*(Y(2)-Y(1));beta*Y(1)-Y(2)-Y(1)*Y(3);-lambda*Y(3)+Y(1)*Y(2)];
[t,y]=ode45(f,[0 30],[5;13;17]);
subplot(2,2,1)
plot(t,y(:,1),'*')
subplot(2,2,2)
plot(t,y(:,2),'+')
subplot(2,2,3)
plot(t,y(:,3),'o')
subplot(2,2,4)
plot3(y(:,1),y(:,2),y(:,3))

边值问题的Matlab数值解

Matlab中用 bvp4c 和 bvpinit 命令求解常微分方程的两点边值问题，微分方程的标准形式为
$$y^{\prime }=f(x,y),bc(y(a),y(b))=0$$
或者是
$$y^{\prime }=f(x,y,p),bc(y(a),y(b),p)=0$$
式中：p 为有关的参数； y，f 可以为向量函数；[a,b]为求解的区间 ;``bc`为边界条件

一般来说，边值问题在计算上比初值问题困难得多，特别是由于边值问题的解可能是多值，bvp4c需要提供猜测的初始值。

一般的，bvp4c 的调用格式如下
$$sol=bvp4c(@odefun,@bcfun,solinit,options,p1,p2,\cdots );$$
函数 odefun 的格式如下
$$yprime=odefun(x,y,p1,p2,\cdots );$$
函数 bcfun 的格式如下
$$res=bcfun(ya,yb,p1,p2,\cdots );$$
下面首先给出一个简单的例子。

例1：考察描述在水平面上一个小水滴横截面形状的标量方程
d 2 d x 2 h ( x ) + [ 1 − h ( x ) ] { 1 + [ d d x h ( x ) ] 2 } 3 2 = 0 , h ( − 1 ) = h ( 1 ) = 0 \frac{d^2}{dx^2}h(x)+[1-h(x)]\{1+[\frac{d}{dx}h(x)]^2\}^\frac{3}{2}=0,\quad h(-1)=h(1)=0 dx2d2​h(x)+[1−h(x)]{1+[dxd​h(x)]2}23​=0,h(−1)=h(1)=0
式中：$h(x)$ 为 $x$ 处水滴的高度

解：

设 $y_1(x)=h(x),y_2(x)=h^{\prime }(x)$ ,把上述微分方程写成两个一阶微分方程组，即
$$\{\begin{array}{l}\frac{d}{dx}{y}_1(x)=y_2(x)\\ \frac{d}{dx}{y}_2(x)=[y_1(x)-1][1+y_2(x{)}^2{]}^{\frac{3}{2}}\end{array}$$
上述微分方程组可由如下函数表示

function yprime=drop(x,y);
yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];

边界条件通过残差函数指定，边界条件通过如下函数表示

function res=dropbc(ya,yb);
res=[ya(1);yb(1)];

如，使用 $y_1(x)=\sqrt{1-x^2}$ 和 $y_2(x)=\frac{-x}{0.1+\sqrt{1-x^2}}$ 作为初始猜测解，由如下函数定义

function yinit=dropinit(x);
yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];

利用如下程序求解微分方程的边值问题并画出图像

solinit=bvpinit(linspace(-1,1,20),@ dropinit);
sol=bvp4c(@ drop,@ dropbc,solinit);
fill(sol.x,sol.y(1,:),[0.7,0.7,0.7]);
xlabel('x','Fontsize',12);
ylabel('h','Rotation',0,'Fontsize',12);

这里通过调用函数bvpinit，计算区间$[-1,1]$ 上等间距的20个点的数据，然后调用函数 bvp4c ，得到数值解的结构 sol 。

初始猜测结构 solinit 由两个域，solinit.x 提供初始猜测的 $x$ 值，排列顺序为从左到右排列，其中 solinit.x(1) 和 solinit.x(end) 分别表示 $a$ 和 $b$。对应的，solinit.y(:,i) 给出点 solinit.x(i) 处的初始猜测解。

输出参数 sol 时包含数值解的一个结构，其中 sol.x 给出计算数值解的 $x$ 点，sol.x(i) 处的数值解由 sol.y(:,i) 给出，类似的，sol.x(i) 处数值解的一阶导数由 sol.yp(:,i) 给出。

可以把上面的所有函数放在一个文件中

function sol=example;
solinit=bvpinit(linspace(-1,1,20),@ dropinit);
sol=bvp4c(@ drop,@ dropbc,solinit);
fill(sol.x,sol.y(1,:).[0.7 0.7 0.7]);
xlabel('x','Fontsize',12);
ylabel('h','Rotation',0,'Fontsize',12);

function yprime=drop(x,y);
yprime=[y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];

function res=dropbc(ya,yb);
res=[ya(1);yb(1)];

function yinit=dropinit(x);
yinit=[sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];

也可以使用匿名函数的形式

yprime=@ (x,y) [y(2);(y(1)-1)*(1+y(2)^2)^(3/2)];
res=@ (ya,yb) [ya(1);yb(1)];
yinit=@ (x) [sqrt(1-x.^2);-x./(0.1+sqrt(1-x.^2))];
solinit=bvpinit(linspace(-1,1,20),yinit);
sol=bvp4c(yprime,res,solinit);
fill(sol.x,sol.y(1,:).[0.7 0.7 0.7]);
xlabel('x','Fontsize',12);
ylabel('h','Rotation',0,'Fontsize',12);

例2：描述 $x=0$ 处固定，$x=1$ 处有弹性支持，沿着 $x$ 轴平衡位置以均匀角速度旋转的绳的位移方程
$$\frac{d^2}{d{x}^2}y(x)+\mu y(x)=0$$
具有边界条件
$$y(0)=0,[\frac{d}{dx}y(x)]{|}_{x=0}=1,[y(x)+\frac{d}{dx}y(x)]{|}_{x=1}=0$$
解：这个边值问题是一个特征问题，必需找到参数 $\mu$ 的值使得方程的解存在。如果提供参数 $\mu$ 的猜测值和对应解的猜测值，可以利用函数 bvp4c 求解特征问题。上述微分方程可以写成下面的微分方程组：
$$\{\begin{array}{l}\frac{d}{dx}{y}_1(x)=y_2(x)\\ \frac{d}{dx}{y}_2(x)=-\mu y_1(x)\end{array}$$
使用 $\mu =5,y_1(x)=\sin x,y_2(x)=\cos x$ 作为初始猜测解，编写程序如下

eq=@ (x,y,mu) [y(2);-mu*y(1)];
bd=@ (ya,yb,mu) [ya(1);ya(2)-1;yb(1)+yb(2)]
guess=@ (x) [sin(x);cos(x)];
guess_structure=bvpinit(linspace(0,1,10),guess,5);
sol=bvp4c(eq,bd,guess_structure);
plot(sol.x,sol.y(1,:),'-',sol.x,sol.yp(1,:),'--','LineWidth',2);
xlabel('x','Fontsize',12);
legend('y_1','y_2');

例3：求解微分方程组
$$\{\begin{array}{l}{u}^{\prime }=0.5u(w-u)/v\\ v^{\prime }=-0.5(w-u)\\ w^{\prime }=[0.9-100(w-y)-0.5w(w-u)]/z\\ z^{\prime }=0.5(w-u)\\ y^{\prime }=-100(y-w)\end{array}$$
边界条件为 $u(0)=v(0)=w(0)=1,z(0)=-10,w(1)=y(1)$。

解：

使用如下猜测
$$\{\begin{array}{l}u(x)=1\\ v(x)=1\\ w(x)=-4.5x^2+8.91x+1\\ z(x)=-10\\ y(x)=-4.5x^2+9x+0.91\end{array}$$
编写程序如下

eq=@ (x,y) [0.5*y(1)*(y(3)-y(1))/y(2)
-0.5*(y(3)-y(1))
(0.9-1000*(y(3)-y(5))-0.5*y(3)*(y(3)-y(1)))/y(4)
0.5*(y(3)-y(1))
100*(y(3)-y(5))];
bd=@ (ya,yb) [ya(1)-1;ya(2)-1;ya(3)-1;ya(4)+10;yb(3)-yb(5)];
guess=@ (x) [1;1;-4.5*x^2+8.91*x+1;-10;-4.5*x^2+9*x+0.91];
guess_structure=bvpinit(linspace(0,1,5),guess);
sol=bvp4c(eq,bd,guess_structure);
plot(sol.x,sol.y(1,:),'-*',sol.x,sol.y(2,:),'-D',sol.x,sol.y(3,:),':S',sol.x,sol.y(4,:),'-.O',sol.x,sol.y(5,:),'--P');
legend('u','v','w','z','y','Location','southwest');