第七章: 常微分方程

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

y' = \frac{dy}{dx} = 2xy

一维弹簧振子 假设一个质量为 m 的物体连接在弹簧的一端，弹簧另一端固定. 当物体受到弹簧的恢复力作用时，它会在平衡位置附近做简谐振动. 根据胡克定律，弹簧的恢复力与位移 x(t) 成正比，方向与位移相反，其表达式为：

F_{\text{弹}} = -kx(t)

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

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

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

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

整理得：

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

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

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

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

方程的解

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

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

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

其中，A 和 B 是由初始条件决定的常数.

定解条件
如果初始时刻 t = 0 时，物体的位移为 x(0) = x_0，速度为 v(0) = v_0，我们可以通过代入初始条件求解 A 和 B：

1. 初始位移条件： x(0) = A = x_0

2. 初始速度条件：

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

代入初始条件得到：

x(t) = x_0 \cos(\omega t) + \frac{v_0}{\omega} \sin(\omega t)

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

单摆的描述 单摆由长度为 L 的细绳和质量为 m 的小球组成，小球在重力作用下绕固定点摆动. 假设摆角为 \theta(t)，即小球偏离竖直方向的角度，我们将推导其运动方程.

受力分析 小球所受的力主要包括： 1. 重力 mg，方向竖直向下； 2. 绳子的张力 T，方向沿绳子向内.

将重力分解为两个方向： - 沿绳方向的分力 mg\cos\theta，被张力抵消； - 垂直于绳方向的分力 mg\sin\theta，产生恢复力矩.

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

\tau = -mgL\sin\theta

根据牛顿第二定律的角度形式，力矩等于转动惯量 I 乘角加速度\alpha：

\tau = I \alpha

对于单摆，转动惯量 \displaystyle I = mL^2，角加速度\displaystyle \alpha = \frac{d^2\theta}{dt^2}，因此：

-mgL\sin\theta = mL^2 \frac{d^2\theta}{dt^2}

消去 m 和 L 后，得到：

\frac{d^2\theta}{dt^2} + \frac{g}{L}\sin\theta = 0

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

\frac{d^2\theta}{dt^2} + \frac{g}{L}\theta = 0

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

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

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

\begin{equation} \frac{dy}{dt} = k y \end{equation}

其中 k 为常数. 方程 (4) 是一类非常重要的微分方程, 方程左边是 y 随时间的变化率, 整个方程则表示 y 的瞬时变化率与 y 在当前时刻的值成正比. 这类方程在实际问题中十分常见, 下面我们看两个实例.

情景1: 传染病

考虑一个传染病爆发的场景, 我们用 y(t) 表示 t 时刻感染者的数目. 由于疾病的传播依赖于人和人的接触, 所以当前的感染者 y(t) 约多, 新增的感染人数就会越大. 我们用一个系数 k 来表示这个病的传播能力, 可以理解为在单位时间内一个感染者 (平均) 能够传染 k 个人. 因此在过了 \Delta t 的时间后, t_0+\Delta t 时刻的感染者人数等于 t_0 时刻的感染者人数加上 \Delta t 这段时间内新增的感染者人数, 即

\begin{equation} y(t_0+\Delta t) \approx y(t_0) + k y(t_0) \Delta t \end{equation}

注意这里时约等于, 因为在 \Delta 这段时间内, 感染者的人数并总是不严格等于 y(t_0) 的, 但是如果 \Delta t 取得很小, 那么 y 的变化也非常小 (y(t) 的连续性), 从而可以近似看成常数 y(t_0).

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

\frac{y(t_0+\Delta t) - y(t_0)} {\Delta t} \approx k y(t_0)

根据上面的分析, 当我们让 \Delta t 区域 0 时, 左边趋近于 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+\Delta t) \approx y(t_0) - k y(t_0) \Delta t

注意这里的负号表示原子数目在减小, 而约等于号表示在 \Delta t 时间内，衰变速率近似保持恒定，因为当 \Delta t 很小时，y(t) 变化非常小，y(t_0) 可以被近似看作常数.

对上式整理并移项得到：

\frac{y(t_0+\Delta t) - y(t_0)}{\Delta t} \approx -k y(t_0)

当我们令 \Delta t \to 0 时，左边趋近于 y'(t_0)，约等号也趋近于等号，因此得到：

y'(t_0) = -k y(t_0)

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

下面我们来”蒙”方程(4)的解. 我们知道, \mathrm{e}^t 的导数等于它自己, 因此 \mathrm{e}^t 是满足微分方程 y' = y 的一个解. 对 \mathrm{e}^x 进行适当边形, 可以构造出 (4) 的解为

y(t) = C\mathrm{e}^{kt}

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

在情景1中, 假设初始时刻 t=0 时感染者的人数 y(0) = 100, 可以得到 C = 100. 于是 t 时刻的感染者等于

y(t) = 100\mathrm{e}^{kt}

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

在情景2中，假设初始时刻 t=0 时放射性原子的数目 y(0) = N_0，因此微分方程

y' = -ky

的解为:

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，从而了解放射性物质的衰变特性.

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

\frac{dy}{y} = k dx

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

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

从而

y = \mathrm{e}^{kx +\tilde{C}} = C\mathrm{e}^{kx}

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

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

step 1 将方程的非齐次项 b(x) 去掉，得到齐次方程：

y' = k(x)y

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

\frac{dy}{y} = k(x)dx

积分得到：

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

进一步简化：

y_h = Ce^{\int k(x) dx}

这是齐次方程的一般解.

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

y_p = u(x)e^{\int k(x) dx}

其中 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) dx} + u(x)k(x)e^{\int k(x) dx}

将 y_p 和 y_p' 代入原方程：

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

化简后得：

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

两边同时除以 \displaystyle e^{\int k(x) dx}：

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

积分求解 u(x)：

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

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

y = y_h + y_p = Ce^{\int k(x) dx} + \left(\int b(x)e^{-\int k(x) dx} dx\right)e^{\int k(x) dx}

最后得到方程的解为:

y = Ce^{\int k(x) dx} + e^{\int k(x) dx} \int b(x)e^{-\int k(x) dx} dx

其中 C 是任意常数.

考虑一个质量为 $ m $ 的物体，固定在弹簧上，弹簧常数为 $ k $，并且系统中存在阻力（或阻尼），阻力与速度成正比，比例常数为 $c $. 此外，系统受到一个外力 $ F(t) $ 的作用. 我们的目标是建立该系统的运动方程.

系统的受力分析为:

• 弹簧力：根据胡克定律，弹簧的恢复力是与物体的位移 $ x(t) $ 成正比的，方向与位移相反，表达式为： F_{\text{spring}} = -k \cdot x(t)

• 阻尼力：阻尼力与物体的速度 $ $ 成正比，方向与速度相反，表达式为： F_{\text{damping}} = -c \cdot \frac{dx}{dt}

• 外力：外力 F(t) 是已知的外部作用力.

根据牛顿第二定律，物体的加速度 $ $ 与合力之比为物体的质量：

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

将各个力代入:

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

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

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

这就是有外力的阻尼震荡方程. - 当 F(t) = 0 时，方程为齐次方程，描述的是自由阻尼振动： \frac{d^2x}{dt^2} + \frac{c}{m} \cdot \frac{dx}{dt} + \frac{k}{m} \cdot x(t) = 0

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

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

\frac{dx}{dt} = \alpha x - \beta xy

\frac{dy}{dt} = \delta xy - \gamma y

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

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

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

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

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

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