第7章常微分方程7.1 什么是微分方程7.2 可分离变量方程7.3 常系数一阶线性微分方程7.4 变系数一阶线性微分方程7.5 二阶线性微分方程和微分方程组

第7章常微分方程

Tip

微分方程是解决实际问题的重要数学工具. 现实世界中的很多问题, 如天体运动, 空气和水的流动, 疾病的传播等现象都可以用微分方程来描述, 本章的任务是介绍几个典型的微分方程及其求解方法.

7.1 什么是微分方程

Important

顾名思义, 微分方程就是含有微分的等式. 例如

y^{'} = \frac{d y}{d x} = 2 x y

Caution

微分方程的阶数

$y^{\prime \prime}, y^{\prime \prime \prime}, \cdots, y^{(n)}$ , 方程中所包含的最高阶的导数为几阶, 那么我们就称该微分方程为几阶的.

Caution

微分方程的解

$y(x)$ $y(x)$ 很容易判断它是不是一个微分方程的解, 但是反过来却不容易: 事实上, 大部分的微分方程我们是解不出来的. 接下来我们要介绍几类比较简单的能够求解的微分方程.

Note

自由落体的小球

牛顿第二定律 $y$ $t$ $y(t)$ 满足

\begin{matrix} (1) & y^{''} (t) = g \end{matrix}

方程的右边是一个常数 (重力加速度), 左边是一个二阶导数, 所以这是一个二阶微分方程.

下面我们来求解这个微分方程:

Step 1.

$v(t)=y^{\prime}(t)$ , 则原方程可化为:

\begin{matrix} (2) & y^{″} (t) = v^{'} (t) = g \end{matrix}

$v(t)$ $t$ 做不定积分, 得到

v (t) = g t + C_{1}

$v(t)$ 确实满足上述微分方程, 所以这样我们就求解了该方程.

$C_1$ $C_1$ $C_1$ 物理条件 $C_1$ 的值, 这样的条件也被称为定解条件.

定解条件1 $2m/s$ $t=0$ $v(0)=2$ .

$C_1=2$ , 于是我们得到

\begin{matrix} (3) & v (t) = g t + 2 \end{matrix}

为满足定解条件1的微分方程 (2) 的解.

Step 2.

接下来我们在 (3) 式的基础上继续求解微分方程 (1). 此时

y^{'} (t) = v (t) = g t + 2

对等式两边积分得到

y (t) = \frac{1}{2} g t^{2} + 2 t + C_{2}

$C_2$ 的具体数值, 我们需要第二个定解条件.

定解条件2 $t=0$ $y(0)=1$

同定解条件1一样, 定解条件2也是来自于物理而非方程本身.

定解条件2 $C_2 = 1$ , 从而微分方程 (1) 的解为

y (t) = \frac{1}{2} g t^{2} + 2 t + 1.

此时, 方程的解中已经不含有未知数, 也就是说它是微分方程 (1) 和定解条件1&2所确定的唯一解.

Warning

可能有同学已经注意到了, 上个例子中我们所谓的解微分方程不就是在做不定积分吗? 这是对的. 后面我们会体会到, 很多时候求解微分方程都要归结到积分运算.

Note

简谐振动

一维弹簧振子 $m$ $x(t)$ 成正比，方向与位移相反，其表达式为：

F_{弹} = - k x (t)

$k$ $x(t)$ 为物体相对于平衡位置的位移. 根据牛顿第二定律，力等于质量乘加速度，因此：

F = m a = m \frac{d^{2} x (t)}{d t^{2}}

将两者代入得到弹簧振子的运动方程：

m \frac{d^{2} x (t)}{d t^{2}} = - k x (t)

整理得：

\frac{d^{2} x (t)}{d t^{2}} + \frac{k}{m} x (t) = 0

$\displaystyle \omega^2 = \frac{k}{m}$ ，方程化为：

\frac{d^{2} x (t)}{d t^{2}} + ω^{2} x (t) = 0

这就是简谐振动的微分方程.

方程的解

接下来我们来"蒙"一个微分方程的解. 你没有听错, 因为微分方程太难解了, 以至于"蒙"也成了一个合理的办法, 总之不管用什么方法, 只要能找到解就行.

该二阶线性齐次微分方程的通解为：

x (t) = A \cos (ω t) + B \sin (ω t)

$A$ $B$ 是由初始条件决定的常数.

定解条件
$t = 0$ $x(0) = x_0$ $v(0) = v_0$ $A$ $B$ ：

初始位移条件：

x (0) = A = x_{0}

初始速度条件：

\begin{aligned} v (0) & = x^{'} (0) = - A ω \sin (0) + B ω \cos (0) \\ = B ω = v_{0} \end{aligned}

代入初始条件得到：

x (t) = x_{0} \cos (ω t) + \frac{v_{0}}{ω} \sin (ω t)

这个解描述了弹簧振子在简谐运动中的位移随时间的变化情况.

Note
单摆

单摆的描述 $L$ $m$ $\theta(t)$ ，即小球偏离竖直方向的角度，我们将推导其运动方程.

受力分析 小球所受的力主要包括：

重力 $mg$ ，方向竖直向下；
绳子的张力 $T$ ，方向沿绳子向内.

将重力分解为两个方向：

$mg\cos\theta$ ，被张力抵消；
$mg\sin\theta$ ，产生恢复力矩.

力矩分析 $\tau$ 关于固定点 O 为：

τ = - m g L \sin θ

$I$ $\alpha$ ：

τ = I α

$\displaystyle I = mL^2$ $\displaystyle \alpha = \frac{d^2\theta}{dt^2}$ ，因此：

- m g L \sin θ = m L^{2} \frac{d^{2} θ}{d t^{2}}

$m$ $L$ 后，得到：

\frac{d^{2} θ}{d t^{2}} + \frac{g}{L} \sin θ = 0

小角近似 $\theta$ $\sin\theta \approx \theta$ （弧度制），方程简化为：

\frac{d^{2} θ}{d t^{2}} + \frac{g}{L} θ = 0

这就是单摆在小角度近似下的简谐运动方程.

方程的解
该方程的解跟简谐振动完全一样.

Note

LC电路
$L$ $C$ $q(t)$ $i(t)$ 也随时间变化. 我们来推导该系统的运动方程.

基本方程

电容的关系 $u_C(t)$ $q(t)$ 的关系为：

u_{C} (t) = \frac{q (t)}{C}

电感的关系 $u_L(t)$ $\displaystyle \frac{di(t)}{dt}$ 的关系为：

u_{L} (t) = L \frac{d i (t)}{d t}

$i(t)$ $q(t)$ 的时间导数：

i (t) = \frac{d q (t)}{d t}

回路方程（基尔霍夫电压定律） 在LC电路中，总电压为零，即：

u_{L} (t) + u_{C} (t) = 0

代入电感和电容的表达式：

L \frac{d^{2} q (t)}{d t^{2}} + \frac{q (t)}{C} = 0

LC震荡方程 将上式整理为标准形式：

\frac{d^{2} q (t)}{d t^{2}} + \frac{1}{L C} q (t) = 0

这就是LC电路的微分方程，描述了电容器电荷随时间的简谐振荡过程.

特征频率 $\omega_0$ 为：

ω_{0} = \sqrt{\frac{1}{L C}}

这表明LC电路的振荡频率只由电感和电容决定，与初始电荷或电流无关.

Note

地心电梯

牛顿壳层定理 $R$ $M$ $m$ 位于球壳内部或外部，球壳对质点的引力为

$r \geq R$ 时，球壳对质点的引力为：

F = G \frac{M m}{r^{2}}

$r < R$ 时，球壳对质点的引力为：

F = 0

这表明在球壳内部，质点所受引力为零，而在外部，引力遵循经典的万有引力公式，等效于将球壳质量集中于球心.

牛顿壳层定理的微积分推导

$\displaystyle dA = R^2 \sin\theta \, d\theta \, d\phi$ $\theta$ $\phi$ $\displaystyle \sigma = \frac{M}{4\pi R^2}$ ，因此面积元的质量：

d m = σ d A = \frac{M}{4 π R^{2}} \cdot R^{2} \sin θ d θ d ϕ

$P$ $dm$ 之间的距离为：

s = \sqrt{r^{2} + R^{2} - 2 r R \cos θ}

$dm$ $P$ 产生的引力：

d F = G \frac{d m \cdot m}{s^{2}}

由于球壳的对称性，垂直于径向的引力分量会相互抵消，仅剩下沿径向的分量：

d F_{径} = d F \cdot \frac{r - R \cos θ}{s}

$dF_{\text{径}}$ 对整个球壳积分：

F = \int_{0}^{2 π} \int_{0}^{π} G \frac{M m}{4 π R^{2}} \cdot \frac{R^{2} \sin θ (r - R \cos θ)}{(r^{2} + R^{2} - 2 r R \cos θ)^{3 / 2}} d θ d ϕ

计算结果表明，球壳外部的引力等效于质量集中于球心：

F = G \frac{M m}{r^{2}}

$r$ $P$ $P$ $P$ $P$ 的引力为：

F_{外} = 0

穿过地心的电梯在重力作用下的运动方程推导

假设一部电梯沿直线穿过地球中心，从地球表面一端移动到另一端，忽略空气阻力和其他阻力，仅考虑重力作用. 我们将推导电梯在这种情况下的运动方程.

$R$ $M$ $\displaystyle \rho = \frac{M}{\frac{4}{3} \pi R^3}$ .

$r$ $r$ 的球体部分，质量为：

M (r) = ρ \cdot \frac{4}{3} π r^{3} = M \cdot \frac{r^{3}}{R^{3}}

$r$ 处的引力为：

F (r) = G \frac{M (r) \cdot m}{r^{2}} = G \frac{M \cdot m \cdot r}{R^{3}}

因此，重力加速度为：

g (r) = G \frac{M}{R^{3}} \cdot r

$-g(r)$ ，即：

\frac{d^{2} r}{d t^{2}} = - G \frac{M}{R^{3}} \cdot r

这就是一个简谐振动方程，其形式为：

\frac{d^{2} r}{d t^{2}} + ω^{2} r = 0

$\omega$ 为：

ω = \sqrt{\frac{G M}{R^{3}}}

7.2 可分离变量方程

Tip

这一节我们介绍一类比较简单的微分方程的解法.

Note

速率方程

在这个例子中我们考虑一阶微分方程

\begin{matrix} (4) & \frac{d y}{d t} = k y \end{matrix}

$k$ $y$ $y$ $y$ 在当前时刻的值成正比. 这类方程在实际问题中十分常见, 下面我们看两个实例.

情景1: 传染病

$y(t)$ $t$ $y(t)$ $k$ $k$ $\Delta t$ $t_0+\Delta t$ $t_0$ $\Delta t$ 这段时间内新增的感染者人数, 即

\begin{matrix} (5) & y (t_{0} + Δ t) \approx y (t_{0}) + k y (t_{0}) Δ t \end{matrix}

$\Delta$ $y(t_0)$ $\Delta t$ $y$ $y(t)$ $y(t_0)$ .

通过对公式 (5) 进行变换可以得到

\frac{y (t_{0} + Δ t) - y (t_{0})}{Δ t} \approx k y (t_{0})

$\Delta t$ $y'(t_0)$ , 而约等号也会"趋近于"等号, 因此我们得到

y^{'} (t_{0}) = k y (t_{0})

$t_0$ $t$ , 就得到微分方程 (4).

情景2: 化学反应

$y(t)$ $t$ $k$ 来表示这种衰变率，称为衰变常数，即单位时间内一个原子发生衰变的概率.

$\Delta t$ $y(t_0+\Delta t)$ $t_0$ 时刻的原子数减去这段时间内衰变的原子数：

y (t_{0} + Δ t) \approx y (t_{0}) - k y (t_{0}) Δ t

$\Delta t$ $\Delta t$ $y(t)$ $y(t_0)$ 可以被近似看作常数.

对上式整理并移项得到：

\frac{y (t_{0} + Δ t) - y (t_{0})}{Δ t} \approx - k y (t_{0})

$\Delta t \to 0$ $y'(t_0)$ ，约等号也趋近于等号，因此得到：

y^{'} (t_{0}) = - k y (t_{0})

$t_0$ $t$ ，同样得到微分方程 (4), 只是系数变成了负数.

$\mathrm{e}^t$ $\mathrm{e}^t$ $y' = y$ $\mathrm{e}^x$ 进行适当边形, 可以构造出 (4) 的解为

y (t) = C e^{k t}

$C$ 为一个参数, 依赖于定解条件.

情景1 $t=0$ $y(0) = 100$ $C = 100$ $t$ 时刻的感染者等于

y (t) = 100 e^{k t}

我们经常听到传染病人数"指数型增长"的可怕说法, 背后对应的就是这个数学公式. 当然, 实际情况下指数增长只是在传染病爆发的初期出现, 随着时间的增长, 未感染的群体数目会减小, 还有一部分感染者将被治愈, 这时微分方程 (4) 就不再成立了, 需要用一个新的微分方程来描述.

情景2 $t=0$ $y(0) = N_0$ ，因此微分方程

y^{'} = - k y

的解为:

y (t) = N_{0} e^{- k t}

这个公式描述了放射性原子数目随时间的指数型衰减. 与传染病的"指数型增长"相反，放射性衰变是一种指数型减少的过程.

$y'(t) = -k y(t)$ 半衰期 $\displaystyle T_{\frac{1}{2}}$ $\displaystyle y(T_{\frac{1}{2}}) = \frac{N_0}{2}$ ，可以得到：

T_{\frac{1}{2}} = \frac{\ln 2}{k}

$k$ ，从而了解放射性物质的衰变特性.

Important

可分离变量方程的求解

方程 (4) 的求解其实可以不通过"蒙", 而是通过一种叫做 "separation of variables" 的技巧求解. 注意到 (4) 可以变形为

\frac{d y}{y} = k d x

$y$ $x$ , 具有这个性质的微分方程就是可分离变量的微分方程. 这类方程的求解方法为分别对等式左右两边积分, 对上式而言, 可以得到

\ln y = k x + \tilde{C}

从而

y = e^{k x + \tilde{C}} = C e^{k x}

这正是我们刚才猜到的解的形式.

Note

例

$\displaystyle y' = \frac{y^2}{x^2}$

解: 方程可转换为变量分离的形式:

\frac{d y}{y^{2}} = \frac{d x}{x^{2}}

对等式两边积分得到

- \frac{1}{y} = - \frac{1}{x} + C

即

y = \frac{x}{1 + k x}

$k$ 为常数.

7.3 常系数一阶线性微分方程

Tip

$y' = ky + b$ 一阶线性微分方程 $b$ $k$ $b$ $k(x)$ $b(x)$ $x$ 变化的情形.

Note

例: 恒定速率的降温过程(牛顿冷却定律)

模型描述
牛顿冷却定律描述了物体与周围环境之间的温度变化过程，其微分方程模型为：

\frac{d T}{d t} = - k (T - T_{env})

其中：

$\displaystyle T(t)$ $t$ 的温度.
$\displaystyle T_{\text{env}}$ ：环境温度（常数）.
$k$ $k > 0$ .

求解过程（分离变量法）

将方程整理为：

\frac{d T}{T - T_{env}} = - k d t

对两边积分：

\int \frac{d T}{T - T_{env}} = \int - k d t

积分结果为：

\ln | T - T_{env} | = - k t + C

$T$ ：

T - T_{env} = C_{1} e^{- k t}

$\displaystyle C_1 = e^C$ 为任意常数.

$T(0) = T_0$ $C_1$ ：

T_{0} - T_{env} = C_{1}

最终解为：

T (t) = T_{env} + (T_{0} - T_{env}) e^{- k t}

解释 $T_{\text{env}}$ $k$ 控制.

Note

例: 简单贷款偿还模型

模型描述 贷款偿还模型描述了借款人以恒定利率和固定还款额偿还贷款的过程，其微分方程为：

\frac{d B}{d t} = r B - P

其中：

$B(t)$ $t$ 时剩余的贷款余额.
$r$ ：年利率（常数）.
$P$ ：每年固定还款额.

求解过程（分离变量法）

将方程整理为：

\frac{d B}{r B - P} = d t

对两边积分：

\int \frac{d B}{r B - P} = \int d t

积分结果为：

\frac{1}{r} \ln | r B - P | = t + C

$B$ ：

r B - P = C_{1} e^{r t}

$C_1 = e^{rC}$ .

$B(0) = B_0$ $C_1$ ：

r B_{0} - P = C_{1}

最终解为：

B (t) = \frac{P}{r} + (B_{0} - \frac{P}{r}) e^{r t}

解释 $P > rB_0$ $B(t)$ $\displaystyle P \leq rB_0$ ，贷款余额将持续增长.

7.4 变系数一阶线性微分方程

Tip

$k(x)$ $b(x)$ $x$ 一阶线性微分方程 $\displaystyle y' = k(x)y + b(x)$ .

Note

例: 城市热岛效应模型讲义

模型描述

城市热岛效应是指由于人类活动和城市结构导致城市温度显著高于周边郊区的现象. 为了描述这种现象，可以使用以下一阶线性微分方程模型：

\frac{d T}{d t} = - α (t) T + Q (t)

其中：

$T(t)$ ：城市的平均温度（随时间变化）.
$\alpha(t)$ $\alpha(t)$ $\alpha(t)$ 较大.
$Q(t)$ ：时间相关的外部热源，包括：太阳辐射：随时间变化呈现周期性. 人类活动：如车辆排放、工业生产、空调使用等。通常在白天达到高峰.

模型扩展(通常用在数学建模中)

加入风速影响： $\alpha(t)$ $\alpha(t, v_w)$ $v_w$ 为风速. 风速较大时，散热效果增强.
$\frac{d T}{d t} = - α (t, v_{w}) T + Q (t)$
考虑湿度和绿化率： 可以通过引入湿度和绿化率来模拟更复杂的热交换过程：
$\frac{d T}{d t} = - (α (t) + β H + γ G) T + Q (t)$
$H$ $G$ $\beta$ $\gamma$ 为对应的影响系数.

Important

常数变异法

求解变系数的一阶线性微分方程可以使用常数变异法. 过程如下:

step 1 $b(x)$ 去掉，得到齐次方程：

y^{'} = k (x) y

这个方程是可分离变量的，解法为：

\frac{d y}{y} = k (x) d x

积分得到：

\ln | y | = \int k (x) d x + C

进一步简化：

y_{h} = C e^{\int k (x) d x}

这是齐次方程的一般解.

step 2 假设非齐次方程的特解形式为：

y_{p} = u (x) e^{\int k (x) d x}

$u(x)$ 是待定函数.

step 3 $\displaystyle y_p = u(x)e^{\int k(x) dx}$ $y' = k(x)y + b(x)$ ，利用链式法则：

y_{p}^{'} = u^{'} (x) e^{\int k (x) d x} + u (x) k (x) e^{\int k (x) d x}

$y_p$ $y_p'$ 代入原方程：

u^{'} (x) e^{\int k (x) d x} + u (x) k (x) e^{\int k (x) d x} = k (x) u (x) e^{\int k (x) d x} + b (x)

化简后得：

u^{'} (x) e^{\int k (x) d x} = b (x)

$\displaystyle e^{\int k(x) dx}$ ：

u^{'} (x) = b (x) e^{- \int k (x) d x}

$u(x)$ ：

u (x) = \int b (x) e^{- \int k (x) d x} d x + C

通解为齐次解与特解的叠加：

y = y_{h} + y_{p} = C e^{\int k (x) d x} + (\int b (x) e^{- \int k (x) d x} d x) e^{\int k (x) d x}

最后得到方程的解为:

y = C e^{\int k (x) d x} + e^{\int k (x) d x} \int b (x) e^{- \int k (x) d x} d x

$C$ 是任意常数.

7.5 二阶线性微分方程和微分方程组

Tip

实际问题中经常出现形如

y^{″} (x) + Q (x) \cdot y^{'} (x) + R (x) \cdot y (x) = f (x)

的微分方程. 这类方程属于二阶线性微分方程.

Note

带有阻尼和外力的振动方程

$m$ $k$ $c$ $F(t)$ 的作用. 我们的目标是建立该系统的运动方程.

系统的受力分析为:

弹簧力 $x(t)$ 成正比的，方向与位移相反，表达式为：
$F_{spring} = - k \cdot x (t)$
阻尼力 $\frac{dx}{dt}$ 成正比，方向与速度相反，表达式为：
$F_{damping} = - c \cdot \frac{d x}{d t}$
外力 $F(t)$ 是已知的外部作用力.

$\frac{d^2x}{dt^2}$ 与合力之比为物体的质量：

m \cdot \frac{d^{2} x}{d t^{2}} = F_{spring} + F_{damping} + F (t)

将各个力代入:

m \cdot \frac{d^{2} x}{d t^{2}} = - k \cdot x (t) - c \cdot \frac{d x}{d t} + F (t)

将上式整理得到标准形式的二阶线性微分方程：

\frac{d^{2} x}{d t^{2}} + \frac{c}{m} \cdot \frac{d x}{d t} + \frac{k}{m} \cdot x (t) = \frac{F (t)}{m}

这就是有外力的阻尼震荡方程.

$F(t) = 0$ 时，方程为齐次方程，描述的是自由阻尼振动：
$\frac{d^{2} x}{d t^{2}} + \frac{c}{m} \cdot \frac{d x}{d t} + \frac{k}{m} \cdot x (t) = 0$
$F(t) \neq 0$ $F(t) = h\sin( \omega t)$ , 方程为非齐次方程，描述的是有外力作用下的阻尼振动.
$\frac{d^{2} x}{d t^{2}} + \frac{c}{m} \cdot \frac{d x}{d t} + \frac{k}{m} \cdot x (t) = \frac{h}{m} \sin (ω t)$

Warning

前面我们考虑的简谐振动方程可以看成是二阶线性微分方程的特例, 但一般的二阶线性微分方程求解是非常困难的, 需要通过引入新变量的方法将高阶(二阶)线性微分方程转化为一阶线性线性微分方程组.

例如对方程

y^{″} + p y^{'} + q y = 0

$v(t) = y'(t)$ , 则

\begin{matrix} {\begin{aligned} y^{'} = u \\ u^{'} = - p u - q y \end{aligned} \end{matrix}

求解这类一阶线性线性微分方程组需要在高维线性空间中进行, 需要一些线性代数的知识, 因此这部分内容我们将在下学期再介绍.

Tip

由多个微分方程联立得到微分方程组, 刚才的例子就是一个微分方程组. 实际问题中微分方程组也十分的常见.

Note

Lotka-Volterra 模型

Lotka-Volterra 模型是一个经典的捕食-猎物模型，用于描述捕食者与猎物之间的相互作用. 该模型由两个微分方程组成，分别描述猎物和捕食者的动态变化. 模型的标准形式为：

\frac{d x}{d t} = α x - β x y

\frac{d y}{d t} = δ x y - γ y

其中：

$x(t)$ 是猎物（兔子）的数量，
$y(t)$ 是捕食者（狼）的数量，
$\alpha$ 是猎物的自然增长率（不受捕食者影响时），
$\beta$ 是捕食者捕食猎物的效率，
$\delta$ 是捕食者的繁殖率（每捕食一只猎物，捕食者的种群数量增加的比例），
$\gamma$ 是捕食者的死亡率（不受食物影响时）.

这个模型是一个非线性方程组，常用于描述生态系统中捕食者和猎物的动态交互关系.

Note

行星运动方程

\frac{d^{2} x}{d t^{2}} = - G \frac{M}{(x^{2} + y^{2} + z^{2})^{3 / 2}} x

\frac{d^{2} y}{d t^{2}} = - G \frac{M}{(x^{2} + y^{2} + z^{2})^{3 / 2}} y

\frac{d^{2} z}{d t^{2}} = - G \frac{M}{(x^{2} + y^{2} + z^{2})^{3 / 2}} z

该方程组描述了行星在三维空间中的运动.

第7章 常微分方程

7.1 什么是微分方程

7.2 可分离变量方程

7.3 常系数一阶线性微分方程

7.4 变系数一阶线性微分方程

7.5 二阶线性微分方程和微分方程组

第7章常微分方程