常微分方程（四）：拉普拉斯变换

拉普拉斯变换把微积分变成代数。 无需繁琐积分、试凑特解，也无需在最后额外处理初值条件——你只需将整个常微分方程（包括方程本身、外力项和初始数据）一次性变换为关于复变量 $$s$$ 的代数方程。解这个方程就像做中学代数题一样简单，再通过反变换即可得到原问题的解。在此过程中，解的形态转化为几何图像：复平面左半部分的极点对应衰减行为，右半部分导致发散，而虚轴上的极点则引发持续振荡。本章将从基本原理出发构建这一图像，并将其与工程中广泛应用的工具——传递函数、Bode 图、PID 控制——联系起来，正是这些工具使拉普拉斯变换成为动力学领域的通用语言。

总结#

定义及对 $e^{-st}$ 作为“衰减探针”的直观理解
微分性质——将 ODE 转化为代数运算的关键
一张精简的变换表，足以应对本章几乎所有反变换问题
部分分式分解：涵盖互异实极点、重极点与共轭复极点的情形
传递函数、极零点图与基于几何位置的稳定性判据
阶跃响应与冲激响应：二者完整刻画系统动态特性
Bode 图：揭示时域与频域之间的对偶关系
PID 控制：P、I、D 三项如何互补短板

前置知识#

第一至三章内容：一阶与二阶 ODE、叠加原理
代数/微积分中的部分分式分解技巧
复数基础知识： $$a + bi$$ 形式、模长与相位
熟悉积分 $\int_0^\infty e^{-st}\,dt$

为什么需要拉普拉斯变换#

RC\,V_c'(t) + V_c(t) = V_s(t), \qquad V_c(0) = V_0.

若 $$V_s$$ 为常数，使用积分因子法仅需两步即可求解。但一旦 $$V_s$$ 变为阶跃信号、短时冲激、斜坡输入，或来自实际信号发生器的分段波形，初等方法便显得笨拙不堪——你会陷入大量积分常数、分段匹配条件和情形讨论的泥潭。

拉普拉斯变换提供了一套统一且高效的流程，适用于所有上述情况：

变换：对 ODE 两边同时进行变换，导数变为 $$s$$ 的多项式，初始条件自动嵌入；
求解：将所得代数方程解出 $$Y(s)$$ ；
反变换：借助一张小表格将 $$Y(s)$$ 变回 $$y(t)$$ 。

繁琐的中间步骤消失不见，取而代之的是清晰的结构划分：谁驱动系统（输入的变换）、系统如何响应（传递函数），以及系统从何处开始（初始条件已自然融入）。

拉普拉斯变换本质上是将信号 $$f(t)$$ 与探针 $e^{-st}$ 进行加权积分。阶跃与冲激是贯穿本章的两类标准输入。

定义与一张实用的变换表#

正变换#

F(s) = \mathcal{L}\{f(t)\} = \int_0^\infty f(t)\,e^{-st}\,dt.

直观理解： 将 $e^{-st}$ 视为一种探针：对每个复频率 $$s$$ ，它询问：“若用衰减速率为 $\operatorname{Re}(s)$ 的指数对 $$f(t)$$ 加权，还能保留多少？”当 $$s$$ 实部很大且为正时，探针主要感知 $$t=0$$ 附近的行为；当 $\operatorname{Re}(s)$ 较小时，则更关注信号的长期尾部。完整的 $$F(s)$$ 即是对所有 $$s$$ 的回答集合，构成了 $$f(t)$$ 的“指纹”。

一张够用的变换表#

$$f(t)$$	$F(s) = \mathcal{L}\{f(t)\}$	收敛域
$$1$$ （或 $$u(t)$$ ）	$$1/s$$	$\operatorname{Re}(s) > 0$
$$t^n$$	$n!/s^{n+1}$	$\operatorname{Re}(s) > 0$
$e^{at}$	$$1/(s-a)$$	$\operatorname{Re}(s) > a$
$\sin\omega t$	$\omega/(s^2+\omega^2)$	$\operatorname{Re}(s) > 0$
$\cos\omega t$	$s/(s^2+\omega^2)$	$\operatorname{Re}(s) > 0$
$e^{at}\sin\omega t$	$\omega/((s-a)^2+\omega^2)$	$\operatorname{Re}(s) > a$
$e^{at}\cos\omega t$	$(s-a)/((s-a)^2+\omega^2)$	$\operatorname{Re}(s) > a$
$\delta(t)$ （冲激）	$$1$$	所有 $$s$$
$$u(t-a)$$ （延迟阶跃）	$e^{-as}/s$	$\operatorname{Re}(s) > 0$

这张表已足够处理本章几乎所有的反变换问题。

两个基于定义的推导#

\int_0^\infty e^{-st}\,dt = \left[-\frac{1}{s}\,e^{-st}\right]_0^\infty = \frac{1}{s}, \qquad \operatorname{Re}(s) > 0.

\int_0^\infty e^{at}\,e^{-st}\,dt = \int_0^\infty e^{-(s-a)t}\,dt = \frac{1}{s-a}, \qquad \operatorname{Re}(s) > a.

这一技巧——将 $e^{at}$ 合并进核函数——使得后续的频移性质变得显而易见。

核心性质：真正干活的工具#

线性#

\mathcal{L}\{a f + b g\} = a F(s) + b G(s).

微分性质——整个方法的基石#

\mathcal{L}\{f'(t)\} = sF(s) - f(0),

\mathcal{L}\{f''(t)\} = s^2 F(s) - s f(0) - f'(0),

\mathcal{L}\{f^{(n)}(t)\} = s^n F(s) - s^{n-1} f(0) - \cdots - f^{(n-1)}(0).

为何关键？ 对时间 $$t$$ 求导转化为对 $$s$$ 的乘法，而初始条件直接作为公式的一部分出现，无需事后匹配。一个 $$n$$ 阶线性 ODE 因此被转化为一个关于 $$s$$ 的 $$n$$ 次代数方程，该方程天然包含了 $y(0), y'(0), \dots$ 等信息。

\int_0^\infty f'(t)\,e^{-st}\,dt = \left[f(t)\,e^{-st}\right]_0^\infty + s\int_0^\infty f(t)\,e^{-st}\,dt = sF(s) - f(0),

前提是当 $t \to \infty$ 时 $f(t)e^{-st} \to 0$ ，这正是收敛域的物理含义。

频移#

\mathcal{L}\{e^{at}f(t)\} = F(s-a).

在时域乘以 $e^{at}$ ，相当于在 $$s$$ 域将整个变换平移 $$a$$ 。这也解释了为何变换表中所有“阻尼”项都只是无阻尼项将 $$s$$ 替换为 $$s-a$$ 的结果。

时移（延迟）#

\mathcal{L}\{f(t-a)\,u(t-a)\} = e^{-as}F(s), \qquad a > 0.

时域的延迟在 $$s$$ 域表现为相位因子 $e^{-as}$ 。注意必须乘以门函数 $$u(t-a)$$ ，以确保变换的是延迟且截断后的信号，而非原始信号。

卷积#

\mathcal{L}\{(f * g)(t)\} = F(s)\,G(s), \qquad (f * g)(t) = \int_0^t f(\tau)\,g(t-\tau)\,d\tau.

在线性时不变（LTI）系统中，时域响应是卷积，而在 $$s$$ 域则简化为乘积。这正是工程框图能够成立的代数基础。

终值定理#

\lim_{t\to\infty} f(t) = \lim_{s\to 0} sF(s),

前提是极限存在，且 $$sF(s)$$ 的所有极点均位于开左半平面。利用该定理，可直接从 $$Y(s)$$ 读取稳态值，无需进行反变换。

求解 ODE：两个示例展示完整流程#

一阶方程 + 指数激励#

求解 $y' + 2y = e^{-t}$ ， $$y(0) = 1$$ 。

sY(s) - 1 + 2Y(s) = \frac{1}{s+1}.

(s+2)Y(s) = 1 + \frac{1}{s+1}, \qquad Y(s) = \frac{1}{s+2} + \frac{1}{(s+1)(s+2)}.

Y(s) = \frac{1}{s+2} + \frac{1}{s+1} - \frac{1}{s+2} = \frac{1}{s+1}.

\boxed{\; y(t) = e^{-t}.\;}

验证： $y' + 2y = -e^{-t} + 2e^{-t} = e^{-t}$ ，且 $$y(0) = 1$$ ，正确。

二阶共振#

求解 $y'' + \omega_0^2\,y = \cos\omega_0 t$ ，初值 $$y(0) = y'(0) = 0$$ 。

s^2 Y + \omega_0^2 Y = \frac{s}{s^2 + \omega_0^2}, \qquad Y(s) = \frac{s}{(s^2 + \omega_0^2)^2}.

y(t) = \frac{t}{2\omega_0}\,\sin\omega_0 t.

关键在于显式的 $$t$$ 因子。代数上，它源于 $s = \pm j\omega_0$ 处的重极点：高重数极点在时域产生“多项式 × 正弦”形式。物理上，这意味着振幅无界增长——即共振现象。

$$Y(s)$$ 中的重极点是共振的代数标志：反变换会引入 $$t$$ 因子，导致响应随时间线性增长。

部分分式：唯一真正需要掌握的技术#

一旦 $$Y(s)$$ 表示为有理函数，反变换的核心工作就是将其拆分为若干可查表的基本项之和。

互异实极点#

\frac{P(s)}{(s-a)(s-b)} = \frac{A}{s-a} + \frac{B}{s-b}.

使用遮盖法： $$A$$ 等于将 $$(s-a)$$ 遮住后，剩余部分在 $$s=a$$ 处的取值（即 $$P(s)/(s-b)$$ 在 $$s=a$$ 的值）。

重极点#

\frac{P(s)}{(s-a)^3} = \frac{A_1}{s-a} + \frac{A_2}{(s-a)^2} + \frac{A_3}{(s-a)^3}.

每个重数为 $$k$$ 的极点，在时域对应一项 $t^{k-1} e^{at}/(k-1)!$ 。

共轭复极点#

\frac{B s + C}{(s-\alpha)^2 + \beta^2}\;\xrightarrow{\;\mathcal{L}^{-1}\;}\; e^{\alpha t}\big(B \cos\beta t + D \sin\beta t\big),

其中 $D = (C + \alpha B)/\beta$ ，通过将分子改写为 $B(s-\alpha) + (C + \alpha B)$ 即可得出。

$Y(s) = \dfrac{3s+5}{(s+1)(s+2)} = \dfrac{2}{s+1} + \dfrac{1}{s+2}$ ，因此 $y(t) = 2 e^{-t} + e^{-2t}$ 。每个极点贡献一个模态，它们线性叠加。

传递函数与稳定性的几何视角#

定义#

H(s) \;=\; \frac{Y(s)}{U(s)}\quad\text{（在零初始条件下）}.

$$H(s)$$ 仅由系统本身决定，与输入或初始状态无关。

H(s) = \frac{1}{RCs + 1} = \frac{1}{\tau s + 1}, \qquad \tau = RC.

该系统在 $s = -1/\tau$ 处有一个实极点。

极点、零点与稳定性#

零点：使 $$H(s) = 0$$ 的 $$s$$ 值（即分子的根）；
极点：使 $H(s) \to \infty$ 的 $$s$$ 值（即分母的根）。

极点分布完全决定了系统的自由响应。每个极点在 $\mathcal{L}^{-1}\{H(s)\}$ 中贡献一个模态：

极点位置	时域模态	行为
实数， $$s = -a$$ ， $$a > 0$$	$e^{-at}$	衰减至零
实数， $$s = a > 0$$	$e^{at}$	无界增长 — 不稳定
共轭复数 $\alpha \pm j\beta$ ， $\alpha < 0$	$e^{\alpha t}(\cos\beta t,\sin\beta t)$	阻尼振荡
纯虚数 $\pm j\beta$	$\cos\beta t,\,\sin\beta t$	等幅振荡
共轭复数， $\alpha > 0$	$e^{\alpha t}(\cos\beta t,\sin\beta t)$	发散振荡 — 不稳定

几何稳定性判据： 系统渐近稳定的充要条件是： $$H(s)$$ 的所有极点严格位于左半平面 $\operatorname{Re}(s) < 0$ 。

极点的实部决定衰减速率，虚部决定振荡频率。稳定性问题归结为：“所有极点是否都在左半平面？”

阶跃响应与冲激响应#

两类标准输入值得单独命名：

冲激响应： $h(t) = \mathcal{L}^{-1}\{H(s)\}$ ，即输入为 Dirac 冲激 $\delta(t)$ 时的输出；
阶跃响应： $s(t) = \mathcal{L}^{-1}\{H(s)/s\}$ ，即输入为单位阶跃 $$u(t)$$ 时的输出；
关系： $$h(t) = s'(t)$$ 。二者包含等价信息，因为阶跃是冲激的积分。

一旦已知 $$h(t)$$ ，系统对任意输入 $$u(t)$$ 的响应即为卷积 $$y(t) = (h * u)(t)$$ ；等价地，在 $$s$$ 域有 $$Y(s) = H(s) U(s)$$ 。

同一系统的两扇窗：时域与频域#

将 $s = j\omega$ 代入 $$H(s)$$ ，得到频率响应 $H(j\omega)$ 。其模 $|H(j\omega)|$ 表示系统对各正弦频率的增益，相位 $\arg H(j\omega)$ 表示相移。在对数–对数坐标下绘制这两者，即为 Bode 图。

H(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2},

阻尼比 $\zeta$ 主导一切行为：当 $\zeta < 1$ 时，极点为复数，阶跃响应出现超调与振铃；当 $\zeta = 1$ 时，极点合并为 $-\omega_n$ 处的二重实极点，响应最快且无超调；当 $\zeta > 1$ 时，极点分离于实轴，响应变得迟缓。

三种视角，同一系统：时域的超调、频域的峰值、极点的几何分布，都是无量纲参数 $\zeta$ 的不同体现。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

sys = signal.TransferFunction([1.0], [1.0, 0.5, 1.0])

t, y = signal.step(sys)
w, mag, phase = signal.bode(sys)

fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
axes[0].plot(t, y);            axes[0].set_title("阶跃响应")
axes[1].semilogx(w, mag);      axes[1].set_title("Bode 幅频 (dB)")
axes[2].semilogx(w, phase);    axes[2].set_title("Bode 相频 (deg)")
for ax in axes: ax.grid(True, alpha=0.3, which="both")
plt.tight_layout()

PID 控制：每项弥补彼此的不足#

u(t) = K_p\,e(t) + K_i\!\int_0^t e(\tau)\,d\tau + K_d\,\frac{de}{dt}.

C(s) = K_p + \frac{K_i}{s} + K_d\,s.

每一项专门解决其他项无法克服的缺陷。

项	功能	优势	缺陷
P（比例）	响应当前误差	快速初步校正	留下稳态误差
I（积分）	累积历史误差	消除稳态误差	降低系统响应速度，可能引发振荡
D（微分）	预测误差变化趋势	抑制超调	放大测量噪声

T(s) = \frac{C(s)\,G(s)}{1 + C(s)\,G(s)}.

调节 $$K_p, K_i, K_d$$ 的过程，本质上是在 $$s$$ 平面上移动闭环极点。目标是将极点尽可能深地推入左半平面（以保证稳定性和快速性），同时避免虚部过大（以防振铃）。

仅用 P 控制会留下稳态偏差；加入 I 项可消除偏差，但会使系统变慢并增加超调；再加入 D 项则能有效抑制超调。完整的 PID 控制可使系统快速进入 $\pm 5\%$ 误差带并保持稳定。

Python 实践：符号与数值计算#

使用 SymPy 进行符号变换#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
from sympy import (symbols, laplace_transform, inverse_laplace_transform,
                   exp, sin, apart)

s, t = symbols("s t")
a, omega = symbols("a omega", positive=True)

print("L{1}        =", laplace_transform(1, t, s))
print("L{exp(-at)} =", laplace_transform(exp(-a*t), t, s))
print("L{sin(wt)}  =", laplace_transform(sin(omega*t), t, s))

F = (3*s + 5) / ((s + 1) * (s + 2))
print("partial fractions:", apart(F, s))
print("L^{-1}{F}        =", inverse_laplace_transform(F, s, t))

使用 SciPy 进行极零点分析#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import numpy as np
from scipy import signal

sys = signal.TransferFunction([1, 2], [1, 3, 2])

poles = np.roots([1, 3, 2])
zeros = np.roots([1, 2])
stable = all(p.real < 0 for p in poles)

print(f"poles  = {poles}")
print(f"zeros  = {zeros}")
print(f"stable = {stable}")

总结#

五步工作流#

变换：对 ODE 两边进行拉普拉斯变换，自动吸收初始条件；
求解：代数求解 $$Y(s)$$ ；
分解：对 $$Y(s)$$ 进行部分分式分解；
反变换：对照表格逐项反变换；
验证：将结果代回原方程检验。

必须牢记的核心性质#

性质	公式	应用场景
微分	$\mathcal{L}\{f'\} = sF(s) - f(0)$	将 ODE 转化为代数方程
频移	$\mathcal{L}\{e^{at}f\} = F(s-a)$	处理阻尼振荡、指数激励
时移	$\mathcal{L}\{f(t-a)u(t-a)\} = e^{-as}F(s)$	处理延迟或开关输入
卷积	$\mathcal{L}\{f * g\} = F(s)\,G(s)$	系统框图、响应计算
终值	$\lim_{t\to\infty} f = \lim_{s\to 0} sF(s)$	直接获取稳态值

一句话总览#

线性时不变系统天然存在于 $$s$$ 平面：极点即其“基因”，左半平面代表稳定，时域运算在 $$s$$ 域化为代数操作。

练习题#

基础

利用频移性质求 $\mathcal{L}\{t^2 e^{-3t}\}$ 和 $\mathcal{L}\{e^{2t}\sin 3t\}$ 。
对 $F(s) = \dfrac{2s+1}{(s+1)(s+3)}$ 进行反变换。
使用拉普拉斯变换求解 $y' + 3y = e^{-2t}$ ，初值 $$y(0) = 1$$ ，并验证结果。
从定义出发，通过分部积分证明 $\mathcal{L}\{f'\} = sF(s) - f(0)$ ，并明确 $t \to \infty$ 时所需的条件。

进阶

求解 $y'' + 4y = \delta(t)$ ，初值 $$y(0) = y'(0) = 0$$ ，并解释为何冲激响应为 $\tfrac{1}{2}\sin 2t$ 。
对 $H(s) = \dfrac{10}{s^2 + 2s + 10}$ ，求其极点、冲激响应与阶跃响应，并识别 $\zeta$ 与 $\omega_n$ 。
利用终值定理求 $\lim_{t\to\infty} y(t)$ ，其中 $Y(s) = \dfrac{5}{s(s^2+3s+2)}$ ，并验证定理前提是否满足。

编程

绘制 $H(s) = \dfrac{100}{s^2 + 10s + 100}$ 的 Bode 幅频与相频图，从图中估计谐振峰位置，并与理论值 $\omega_n\sqrt{1 - 2\zeta^2}$ 比较。
使用 SciPy 的 lsim 函数模拟具有纯时延的系统 $H(s) = \dfrac{e^{-s}}{s+1}$ 的阶跃响应，需先用 Padé 近似展开延迟项。

下一步#

下一章面对一类不老实的 ODE——系数不是常数、解不能写成初等函数。这时候级数法上场。我们会展示 Frobenius 方法和正则奇点的处理，并最终遇到 Bessel、Legendre、Hermite 这些’特殊函数’——它们其实就是某些经典 ODE 的级数解。

参考文献#

Kreyszig, E. Advanced Engineering Mathematics（第 10 版），Wiley (2011) — 第 6 章是拉普拉斯变换的标准教材处理。
Ogata, K. Modern Control Engineering（第 5 版），Pearson (2010) — 深入讲解极零点分析、Bode 图与 PID 整定。
Oppenheim, A. V. & Willsky, A. S. Signals and Systems（第 2 版），Prentice Hall (1997) — 从信号工程视角阐述卷积等概念。
Strang, G. Differential Equations and Linear Algebra，Wellesley-Cambridge (2014) — 更具几何直觉的入门读物。
SciPy: scipy.signal — 提供传递函数、阶跃/冲激响应及 Bode 图的 Python 工具。