• 工程数学完全征服 第1篇：常微分方程（ODE）
  • 1. 微分方程基础
    • 什么是微分方程？
    • 分类：ODE vs PDE
    • 阶数（Order）与次数（Degree）
    • 线性 vs 非线性
    • 通解、特解、奇异解
    • ODE 在工程中的应用
  • 2. 一阶常微分方程
    • 2.1 变量分离法（Separation of Variables）
    • 2.2 积分因子法（Integrating Factor Method）
    • 2.3 恰当微分方程（Exact ODE）
    • 2.4 伯努利方程（Bernoulli Equation）
  • 3. 二阶线性常微分方程
    • 3.1 齐次方程（Homogeneous）
    • 3.2 RLC 电路完整分析
    • 3.3 非齐次方程（Non-homogeneous）
      • 待定系数法（Method of Undetermined Coefficients）
      • 参数变易法（Variation of Parameters）
  • 4. 微分方程组
    • 4.1 矩阵表示
    • 4.2 特征值/特征向量解法
    • 4.3 电路网络分析应用
  • 5. 拉普拉斯变换（Laplace Transform）
    • 5.1 定义与收敛
    • 5.2 主要变换对
    • 5.3 主要性质
    • 5.4 拉普拉斯逆变换
    • 5.5 用拉普拉斯变换求解 ODE
  • 6. 用 Python 数值求解 ODE
    • 6.1 使用 scipy.integrate.solve_ivp
    • 6.2 用 SymPy 求解析解
    • 6.3 比较不同积分器
    • 6.4 用 Python 自动化拉普拉斯变换
  • 7. 工程应用综合例题
    • 7.1 带电源的 RLC 电路暂态分析
    • 7.2 弹簧-质量-阻尼器共振分析
  • 总结与下一步
  • 参考资料

工程数学完全征服 第1篇：常微分方程（ODE）

对工科学生来说，工程数学是一门绕不开的课程。特别是常微分方程（ODE，Ordinary Differential Equation），是电子电路分析、控制系统设计、机械振动解析等几乎所有工程领域的数学基础。本文将从 ODE 的基础概念到拉普拉斯变换，以工程应用例题为中心系统讲解。

1. 微分方程基础

什么是微分方程？

微分方程是表示未知函数与其导数之间关系的方程。例如：

$\frac{dy}{dx} = 2x, \quad y'' + 3y' + 2y = 0, \quad \frac{\partial u}{\partial t} = k\frac{\partial^2 u}{\partial x^2}$

微分方程是对自然现象和工程系统进行数学建模的核心工具。

分类：ODE vs PDE

常微分方程（ODE）：只包含对一个自变量的导数。

$\frac{d^2x}{dt^2} + 2\frac{dx}{dt} + x = 0$

偏微分方程（PDE）：包含对两个及以上自变量的偏导数。

$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}$

本文聚焦于 ODE。

阶数（Order）与次数（Degree）

阶数是方程中出现的最高阶导数的阶。

• 一阶 ODE：$y' + y = x$（最高阶导数为 $y'$）
• 二阶 ODE：$y'' + 2y' + y = 0$（最高阶导数为 $y''$）

次数（degree）是最高阶导数的指数。在大多数工程问题中，次数为 1。

线性 vs 非线性

线性 ODE：未知函数及其导数只以一次项出现。

$a_n(x)y^{(n)} + a_{n-1}(x)y^{(n-1)} + \cdots + a_1(x)y' + a_0(x)y = f(x)$

非线性 ODE：包含 $y^2$、$y \cdot y'$、$\sin(y)$ 等非线性项。

$y'' + \sin(y) = 0 \quad \text{（单摆方程）}$

通解、特解、奇异解

• 通解：包含任意常数的完整解
• 特解：代入初始条件、确定常数之后得到的解
• 奇异解：无法从通解中得到的特别解

ODE 在工程中的应用

工程系统大多用 ODE 来描述。

RC 充电电路： $RC\frac{dv}{dt} + v = V_s$

弹簧-质量-阻尼器系统： $m\ddot{x} + c\dot{x} + kx = F(t)$

人口增长模型： $\frac{dP}{dt} = kP$

RL 电路： $L\frac{di}{dt} + Ri = E(t)$

以上全部都是 ODE，掌握其求解方法，就能够完整分析工程系统的动态行为。

2. 一阶常微分方程

2.1 变量分离法（Separation of Variables）

形式：$\frac{dy}{dx} = f(x)g(y)$

求解过程：

$\frac{dy}{g(y)} = f(x)dx$

两边同时积分：

$\int \frac{dy}{g(y)} = \int f(x)dx + C$

例题：RC 充电电路

电路方程： $\frac{dv}{dt} = \frac{V_s - v}{RC}$

分离变量： $\frac{dv}{V_s - v} = \frac{dt}{RC}$

积分： $-\ln|V_s - v| = \frac{t}{RC} + C_1$

代入初始条件 $v(0) = 0$：

$V_s - v = V_s e^{-t/(RC)}$

$\boxed{v(t) = V_s\left(1 - e^{-t/(RC)}\right)}$

该结果表示 RC 电路的暂态响应（transient response）。电压以时间常数 $\tau = RC$ 决定的速度收敛到 $V_s$。

物理意义：

• 当 $t = \tau$ 时：$v \approx 0.632 V_s$（充至 63.2%）
• 当 $t = 5\tau$ 时：$v \approx 0.993 V_s$（几乎完全充满）

2.2 积分因子法（Integrating Factor Method）

线性一阶 ODE 的标准形： $\frac{dy}{dx} + P(x)y = Q(x)$

积分因子：$\mu(x) = e^{\int P(x)dx}$

推导过程：

两边同乘 $\mu(x)$： $\mu(x)\frac{dy}{dx} + \mu(x)P(x)y = \mu(x)Q(x)$

左边由乘积求导法则可写成： $\frac{d}{dx}[\mu(x)y] = \mu(x)Q(x)$

积分后： $\mu(x)y = \int \mu(x)Q(x)dx + C$

$\boxed{y = \frac{1}{\mu(x)}\left[\int \mu(x)Q(x)dx + C\right]}$

例题：RL 电路

$L\frac{di}{dt} + Ri = E$

化为标准形： $\frac{di}{dt} + \frac{R}{L}i = \frac{E}{L}$

积分因子：$\mu(t) = e^{(R/L)t}$

$\frac{d}{dt}\left[e^{(R/L)t} \cdot i\right] = \frac{E}{L}e^{(R/L)t}$

积分： $e^{(R/L)t} \cdot i = \frac{E}{R}e^{(R/L)t} + C$

代入初始条件 $i(0) = 0$：

$\boxed{i(t) = \frac{E}{R}\left(1 - e^{-(R/L)t}\right)}$

电流以时间常数 $\tau = L/R$ 收敛到稳态值 $E/R$。

2.3 恰当微分方程（Exact ODE）

形式：$M(x,y)dx + N(x,y)dy = 0$

恰当条件： $\frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}$

若满足该条件，则存在使 $F(x,y) = C$ 成立的势函数 $F$。

求解方法：

1. 由 $\frac{\partial F}{\partial x} = M$，对 $x$ 积分求出 $F$
2. 用 $\frac{\partial F}{\partial y} = N$ 这一条件确定积分常数
3. 通解：$F(x,y) = C$

例题：

$(2xy + y^2)dx + (x^2 + 2xy)dy = 0$

验证是否恰当： $\frac{\partial M}{\partial y} = 2x + 2y = \frac{\partial N}{\partial x}$

是恰当方程，因此： $F = \int M\,dx = \int (2xy + y^2)dx = x^2y + xy^2 + g(y)$

由 $\frac{\partial F}{\partial y} = x^2 + 2xy + g'(y) = N = x^2 + 2xy$，得 $g'(y) = 0$，$g(y) = C$

通解：$x^2y + xy^2 = C$

2.4 伯努利方程（Bernoulli Equation）

形式： $\frac{dy}{dx} + P(x)y = Q(x)y^n \quad (n \neq 0, 1)$

换元：令 $v = y^{1-n}$，即可化为线性 ODE。

由 $v' = (1-n)y^{-n}y'$，将原方程两边除以 $y^n$：

$y^{-n}\frac{dy}{dx} + P(x)y^{1-n} = Q(x)$

$\frac{1}{1-n}\frac{dv}{dx} + P(x)v = Q(x)$

$\frac{dv}{dx} + (1-n)P(x)v = (1-n)Q(x)$

这是关于 $v$ 的线性一阶 ODE，可用积分因子法求解。

应用：逻辑斯蒂增长模型 $\frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right)$ 就是 $n=2$ 的伯努利方程。

3. 二阶线性常微分方程

3.1 齐次方程（Homogeneous）

标准形： $ay'' + by' + cy = 0$

特征方程的推导：代入 $y = e^{rx}$。

$a r^2 e^{rx} + b r e^{rx} + c e^{rx} = 0$

$\boxed{ar^2 + br + c = 0}$

根据判别式 $\Delta = b^2 - 4ac$ 的符号，会出现三种情形。

情形 1：两个不同实根（$\Delta > 0$）

若 $r_1 \neq r_2$： $y = C_1 e^{r_1 x} + C_2 e^{r_2 x}$

情形 2：重根（$\Delta = 0$）

若 $r_1 = r_2 = r = -b/(2a)$： $y = (C_1 + C_2 x)e^{rx}$

之所以要加上 $C_2 x e^{rx}$，是因为需要构造出另一个线性无关的解，使朗斯基行列式（Wronskian）不为零。

情形 3：共轭复根（$\Delta < 0$）

若 $r = \alpha \pm i\beta$（其中 $\alpha = -b/(2a)$，$\beta = \sqrt{4ac - b^2}/(2a)$）： $y = e^{\alpha x}(C_1\cos\beta x + C_2\sin\beta x)$

这是借助欧拉公式 $e^{i\beta x} = \cos\beta x + i\sin\beta x$，转换为实数解之后的形式。

3.2 RLC 电路完整分析

电路方程（以电荷 $q(t)$ 表示）： $L\frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{q}{C} = 0$

化为标准形： $q'' + \frac{R}{L}q' + \frac{1}{LC}q = 0$

特征方程： $r^2 + \frac{R}{L}r + \frac{1}{LC} = 0$

$r = -\frac{R}{2L} \pm \sqrt{\left(\frac{R}{2L}\right)^2 - \frac{1}{LC}}$

定义阻尼系数 $\alpha = R/(2L)$、谐振频率 $\omega_0 = 1/\sqrt{LC}$，则：

$r = -\alpha \pm \sqrt{\alpha^2 - \omega_0^2}$

过阻尼（Overdamped）：$\alpha > \omega_0$（$R > 2\sqrt{L/C}$）

$q(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t} \quad (r_1, r_2 < 0)$

由两个指数项之和构成，没有振荡，单调递减。

临界阻尼（Critically Damped）：$\alpha = \omega_0$（$R = 2\sqrt{L/C}$）

$q(t) = (C_1 + C_2 t)e^{-\alpha t}$

没有振荡，以最快速度回到平衡状态。在控制系统设计中极为重要。

欠阻尼（Underdamped）：$\alpha < \omega_0$（$R < 2\sqrt{L/C}$）

定义衰减振荡频率 $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$，则：

$q(t) = e^{-\alpha t}(C_1\cos\omega_d t + C_2\sin\omega_d t)$

一边振荡一边衰减。在通信电路、振动系统中很常见。

物理意义总结：

3.3 非齐次方程（Non-homogeneous）

形式： $ay'' + by' + cy = f(x)$

完全解：$y = y_h + y_p$

• $y_h$：齐次方程的通解（余函数，complementary function）
• $y_p$：特解（particular solution）

待定系数法（Method of Undetermined Coefficients）

根据 $f(x)$ 的形式，猜测 $y_p$ 的形式。

注意：若 $y_p$ 的猜测形式与 $y_h$ 中的项重合，需要乘以 $x$ 来修正。

例题：$y'' + 4y = 3\cos 2x$

齐次解：$y_h = C_1\cos 2x + C_2\sin 2x$

由于 $y_p$ 的猜测形式与 $y_h$ 重合，需要修正： $y_p = x(A\cos 2x + B\sin 2x)$

代入后计算得：$A = 0$，$B = 3/4$

$y_p = \frac{3}{4}x\sin 2x$

这是共振（resonance）现象，振幅随时间成正比增大。

受迫振动响应：

$m\ddot{x} + c\dot{x} + kx = F_0\cos\omega t$

稳态特解： $x_p(t) = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}\cos(\omega t - \phi)$

其中频率比 $r = \omega/\omega_n$，阻尼比 $\zeta = c/(2\sqrt{mk})$。

在 $r = 1$（共振条件）时振幅达到最大。

参数变易法（Variation of Parameters）

用于待定系数法无法适用的一般 $f(x)$。

已知 $y_h = C_1 y_1 + C_2 y_2$ 时：

$y_p = -y_1 \int \frac{y_2 f}{W}dx + y_2 \int \frac{y_1 f}{W}dx$

朗斯基行列式（Wronskian）：$W = y_1 y_2' - y_1' y_2$

该公式适用于任意连续函数 $f(x)$。

4. 微分方程组

4.1 矩阵表示

由 $n$ 个一阶 ODE 组成的方程组：

$\frac{d\mathbf{X}}{dt} = A\mathbf{X}$

其中 $\mathbf{X} = [x_1, x_2, \ldots, x_n]^T$，$A$ 为 $n \times n$ 系数矩阵。

4.2 特征值/特征向量解法

解的形式：$\mathbf{X} = \mathbf{v}e^{\lambda t}$

代入后得：$\lambda \mathbf{v} = A\mathbf{v}$

也就是说，$\lambda$ 是 $A$ 的特征值，$\mathbf{v}$ 是对应的特征向量。

特征方程：$\det(A - \lambda I) = 0$

例题：两个相连的水箱系统

$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -2 & 1 \\ 1 & -2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$

特征方程：$(\lambda+2)^2 - 1 = 0$

$\lambda_1 = -1, \quad \lambda_2 = -3$

特征向量：$\mathbf{v}_1 = [1, 1]^T$，$\mathbf{v}_2 = [1, -1]^T$

通解： $\mathbf{X} = C_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix} e^{-t} + C_2 \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-3t}$

4.3 电路网络分析应用

带有两个回路的电路（基尔霍夫电压定律）：

$L_1\frac{di_1}{dt} + R_1 i_1 - M\frac{di_2}{dt} = E_1(t)$ $L_2\frac{di_2}{dt} + R_2 i_2 - M\frac{di_1}{dt} = E_2(t)$

其中 $M$ 为互感。将其整理为矩阵形式后，可用特征值法或拉普拉斯变换求解。

5. 拉普拉斯变换（Laplace Transform）

5.1 定义与收敛

定义： $\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty f(t)e^{-st}dt$

$s = \sigma + j\omega$ 是复变量，积分要收敛，需要满足 $\text{Re}(s) > \sigma_c$（收敛域右侧边界）。

为什么要用拉普拉斯变换？

• 微分运算转化为代数运算
• 自动处理初始条件
• 把复杂的 ODE 变成简单的代数方程来求解
• 是传递函数（Transfer Function）概念的基础

5.2 主要变换对

推导示例：$\mathcal{L}\{e^{at}\}$

$\mathcal{L}\{e^{at}\} = \int_0^\infty e^{at} e^{-st} dt = \int_0^\infty e^{-(s-a)t} dt = \left[\frac{e^{-(s-a)t}}{-(s-a)}\right]_0^\infty = \frac{1}{s-a}$

当 $s > a$ 时收敛。

5.3 主要性质

线性性： $\mathcal{L}\{af(t) + bg(t)\} = aF(s) + bG(s)$

时移（第一移位定理）： $\mathcal{L}\{e^{at}f(t)\} = F(s-a)$

频移（第二移位定理）： $\mathcal{L}\{f(t-a)u(t-a)\} = e^{-as}F(s)$

导数变换（最重要的性质）： $\mathcal{L}\{f'(t)\} = sF(s) - f(0)$ $\mathcal{L}\{f''(t)\} = s^2F(s) - sf(0) - f'(0)$ $\mathcal{L}\{f^{(n)}(t)\} = s^nF(s) - s^{n-1}f(0) - s^{n-2}f'(0) - \cdots - f^{(n-1)}(0)$

初始条件会在 $s$ 域中被自动纳入。

积分变换： $\mathcal{L}\left\{\int_0^t f(\tau)d\tau\right\} = \frac{F(s)}{s}$

卷积定理（Convolution Theorem）： $\mathcal{L}\{(f*g)(t)\} = F(s)G(s)$

其中 $(f*g)(t) = \int_0^t f(\tau)g(t-\tau)d\tau$

该定理从数学上证明了：系统的冲激响应与输入信号的卷积，就是系统的输出。

终值定理： $\lim_{t\to\infty} f(t) = \lim_{s\to 0} sF(s)$

不必在 $t$ 域求解，就能直接得到稳态值。

初值定理： $\lim_{t\to 0^+} f(t) = \lim_{s\to\infty} sF(s)$

5.4 拉普拉斯逆变换

部分分式分解（Partial Fraction Decomposition）：

将 $F(s) = P(s)/Q(s)$ 的分母因式分解，再拆解为变换表中各项之和。

情形 1：不同的实数极点

$F(s) = \frac{N(s)}{(s-p_1)(s-p_2)\cdots(s-p_n)} = \frac{A_1}{s-p_1} + \frac{A_2}{s-p_2} + \cdots$

$A_k = \lim_{s\to p_k}(s-p_k)F(s)$

情形 2：重极点（$s = p$ 为 $m$ 重根）

$F(s) = \cdots + \frac{A_m}{(s-p)^m} + \frac{A_{m-1}}{(s-p)^{m-1}} + \cdots + \frac{A_1}{s-p}$

$A_k = \frac{1}{(m-k)!}\left[\frac{d^{m-k}}{ds^{m-k}}(s-p)^m F(s)\right]_{s=p}$

情形 3：共轭复极点

$\frac{As + B}{s^2 + 2\alpha s + (\alpha^2 + \beta^2)} \Rightarrow e^{-\alpha t}(C_1\cos\beta t + C_2\sin\beta t)$

亥维赛展开定理（Heaviside Expansion）：

当分母是不同一次因式的乘积时： $F(s) = \sum_{k=1}^n \frac{N(p_k)}{Q'(p_k)} \cdot \frac{1}{s - p_k}$

5.5 用拉普拉斯变换求解 ODE

例题：RLC 串联电路暂态响应

$L\frac{d^2i}{dt^2} + R\frac{di}{dt} + \frac{i}{C} = E_0\delta(t)$

初始条件：$i(0) = 0$，$i'(0) = 0$（电容上无电荷）

求解：

两边应用拉普拉斯变换：

$I(s) = \frac{E_0/L}{s^2 + (R/L)s + 1/(LC)}$

定义 $\alpha = R/(2L)$、$\omega_0 = 1/\sqrt{LC}$、$\omega_d = \sqrt{\omega_0^2 - \alpha^2}$（欠阻尼情形）：

$I(s) = \frac{E_0/L}{(s+\alpha)^2 + \omega_d^2}$

拉普拉斯逆变换（查表：$\mathcal{L}\{e^{-\alpha t}\sin\omega_d t\} = \omega_d/[(s+\alpha)^2+\omega_d^2]$）：

$i(t) = \frac{E_0}{L\omega_d}e^{-\alpha t}\sin\omega_d t$

这就是 RLC 电路的冲激响应，同时也是它的格林函数。

6. 用 Python 数值求解 ODE

6.1 使用 scipy.integrate.solve_ivp

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# RLC 电路仿真（2阶 ODE -> 1阶 ODE 方程组）
# L * q'' + R * q' + q/C = Vs（阶跃输入）
# 状态变量: y[0] = q（电荷），y[1] = i = dq/dt（电流）

def rlc_circuit(t, y, R, L, C, Vs):
    q, dq_dt = y
    # q'' = (Vs - R*i - q/C) / L
    d2q_dt2 = (Vs - R * dq_dt - q / C) / L
    return [dq_dt, d2q_dt2]

# 电路参数
R = 100      # 电阻 (Ohm)
L = 0.1      # 电感 (H)
C = 1e-6     # 电容 (F)
Vs = 10.0    # 电源电压 (V)

# 计算谐振频率与阻尼比
omega_0 = 1.0 / np.sqrt(L * C)
alpha = R / (2 * L)
zeta = alpha / omega_0
print(f"谐振频率: {omega_0/(2*np.pi):.1f} Hz")
print(f"阻尼比 zeta: {zeta:.3f}")
if zeta > 1:
    print("-> 过阻尼")
elif zeta == 1:
    print("-> 临界阻尼")
else:
    print("-> 欠阻尼")

# 数值求解
y0 = [0.0, 0.0]  # 初始条件: q(0) = 0, i(0) = 0
t_span = (0, 0.01)
t_eval = np.linspace(0, 0.01, 2000)

sol = solve_ivp(
    rlc_circuit, t_span, y0,
    t_eval=t_eval,
    args=(R, L, C, Vs),
    method='RK45',
    rtol=1e-8,
    atol=1e-10
)

# 计算电压: v_C = q/C
v_C = sol.y[0] / C
i = sol.y[1]

# 可视化
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(sol.t * 1e3, v_C, 'b-', linewidth=2, label='电容电压')
axes[0].axhline(y=Vs, color='r', linestyle='--', label=f'稳态 {Vs}V')
axes[0].set_xlabel('时间 (ms)')
axes[0].set_ylabel('电压 (V)')
axes[0].set_title('RLC 串联电路 - 阶跃响应')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(sol.t * 1e3, i * 1e3, 'g-', linewidth=2, label='电流')
axes[1].set_xlabel('时间 (ms)')
axes[1].set_ylabel('电流 (mA)')
axes[1].set_title('RLC 电路电流')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_response.png', dpi=150)
plt.show()

6.2 用 SymPy 求解析解

import sympy as sp

# 定义符号
t, s = sp.symbols('t s')
R_sym, L_sym, C_sym = sp.symbols('R L C', positive=True)
Vs_sym = sp.Symbol('Vs')

# RC 电路解析解（用 sympy 求解微分方程）
v = sp.Function('v')
tau = sp.Symbol('tau', positive=True)

# dv/dt = (Vs - v) / (R*C)
ode = sp.Eq(v(t).diff(t), (Vs_sym - v(t)) / (R_sym * C_sym))

# 初始条件 v(0) = 0
solution = sp.dsolve(ode, v(t), ics={v(0): 0})
print("RC 电路解:")
print(solution)
sp.pprint(solution)

# 计算拉普拉斯变换
f = sp.exp(-t) * sp.sin(2*t)
F = sp.laplace_transform(f, t, s, noconds=True)
print(f"\nL{{e^(-t) sin(2t)}} = {F}")

# 拉普拉斯逆变换
F_inv = sp.inverse_laplace_transform(
    1 / (s**2 + 2*s + 5), s, t
)
print(f"L^(-1){{1/(s^2+2s+5)}} = {F_inv}")

6.3 比较不同积分器

import numpy as np
import matplotlib.pyplot as plt

# 欧拉法 vs RK4 比较
# 例题: y' = -y, y(0) = 1, 精确解 y = e^(-t)

def f(t, y):
    return -y

def euler_method(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros(n+1)
    t[0], y[0] = t0, y0
    for i in range(n):
        y[i+1] = y[i] + h * f(t[i], y[i])
        t[i+1] = t[i] + h
    return t, y

def rk4_method(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros(n+1)
    t[0], y[0] = t0, y0
    for i in range(n):
        k1 = h * f(t[i], y[i])
        k2 = h * f(t[i] + h/2, y[i] + k1/2)
        k3 = h * f(t[i] + h/2, y[i] + k2/2)
        k4 = h * f(t[i] + h, y[i] + k3)
        y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
        t[i+1] = t[i] + h
    return t, y

h = 0.5
n = 10
t0, y0 = 0.0, 1.0

t_euler, y_euler = euler_method(f, t0, y0, h, n)
t_rk4, y_rk4 = rk4_method(f, t0, y0, h, n)
t_exact = np.linspace(0, n*h, 200)
y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))
plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='精确解')
plt.plot(t_euler, y_euler, 'r--o', label=f'欧拉法 (h={h})')
plt.plot(t_rk4, y_rk4, 'b--s', label=f'RK4 (h={h})')
plt.xlabel('t')
plt.ylabel("y(t)")
plt.title("欧拉法 vs RK4 比较: y' = -y")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 误差比较
print(f"t=5时欧拉法误差: {abs(y_euler[-1] - np.exp(-5)):.6f}")
print(f"t=5时RK4误差: {abs(y_rk4[-1] - np.exp(-5)):.6f}")

6.4 用 Python 自动化拉普拉斯变换

import sympy as sp

s, t = sp.symbols('s t', real=True)
a, omega = sp.symbols('a omega', positive=True)

# 自动生成拉普拉斯变换表
functions = {
    '1': sp.Integer(1),
    't': t,
    't^2': t**2,
    't^n': t**3,
    'e^(at)': sp.exp(a*t),
    'sin(wt)': sp.sin(omega*t),
    'cos(wt)': sp.cos(omega*t),
    't*e^(at)': t * sp.exp(a*t),
    'e^(at)*sin(wt)': sp.exp(a*t) * sp.sin(omega*t),
}

print("拉普拉斯变换表:")
print("-" * 50)
for name, func in functions.items():
    F = sp.laplace_transform(func, t, s, noconds=True)
    print(f"L{{{name}}} = {sp.simplify(F)}")

7. 工程应用综合例题

7.1 带电源的 RLC 电路暂态分析

问题：一个 $R = 50\,\Omega$、$L = 0.2\,\text{H}$、$C = 100\,\mu\text{F}$ 的 RLC 串联电路，在 $t=0$ 时接入 $V_s = 100\,\text{V}$ 的直流电源。初始条件为 $v_C(0) = 0$，$i(0) = 0$。

用拉普拉斯变换求解：

基尔霍夫电压定律： $L\frac{di}{dt} + Ri + \frac{1}{C}\int i\,dt = V_s u(t)$

用 $i = C dv_C/dt$ 代换，整理为电荷 $q$ 的方程： $L\ddot{q} + R\dot{q} + \frac{q}{C} = V_s$

拉普拉斯变换： $\left(Ls^2 + Rs + \frac{1}{C}\right)Q(s) = \frac{V_s}{s}$

$Q(s) = \frac{V_s}{s\left(Ls^2 + Rs + \frac{1}{C}\right)} = \frac{V_s/L}{s\left(s^2 + \frac{R}{L}s + \frac{1}{LC}\right)}$

代入数值：$\alpha = R/(2L) = 125\,\text{s}^{-1}$，$\omega_0 = 1/\sqrt{LC} \approx 223.6\,\text{s}^{-1}$

由于 $\omega_0 > \alpha$，属于欠阻尼。

$\omega_d = \sqrt{\omega_0^2 - \alpha^2} \approx 185.0\,\text{s}^{-1}$

经部分分式分解后做拉普拉斯逆变换：

$v_C(t) = V_s\left[1 - e^{-\alpha t}\left(\cos\omega_d t + \frac{\alpha}{\omega_d}\sin\omega_d t\right)\right]$

$= 100\left[1 - e^{-125t}\left(\cos 185t + 0.676\sin 185t\right)\right]\,\text{V}$

7.2 弹簧-质量-阻尼器共振分析

问题：$m = 1\,\text{kg}$，$k = 100\,\mathrm{N/m}$，$c = 4\,\mathrm{N \cdot s/m}$，$F(t) = 10\cos\omega t\,\mathrm{N}$

固有频率：$\omega_n = \sqrt{k/m} = 10\,\text{rad/s}$

阻尼比：$\zeta = c/(2\sqrt{mk}) = 0.2$（欠阻尼）

稳态响应振幅： $X = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}$

当 $r = \omega/\omega_n = 1$（共振）时： $X_{\mathrm{res}} = \frac{F_0/k}{2\zeta} = \frac{10/100}{2 \times 0.2} = 0.25\,\text{m}$

是无共振时静态挠度 $F_0/k = 0.1\,\text{m}$ 的 2.5 倍。这就是共振的危险之处。

总结与下一步

本文涉及的内容：

1. 微分方程基础：分类、阶数、线性性、解的种类
2. 一阶 ODE 解法：变量分离法、积分因子法、恰当 ODE、伯努利方程
3. 二阶线性 ODE：特征方程、三种情形、RLC 电路完整分析
4. 非齐次方程：待定系数法、参数变易法、受迫振动
5. 微分方程组：矩阵表示、特征值/特征向量解法
6. 拉普拉斯变换：定义、变换对、性质、逆变换、求解 ODE
7. Python 实现：用 scipy、sympy 求数值解与解析解

下一篇将讲解傅里叶级数/变换与偏微分方程（PDE）。这是信号处理与电磁学的核心数学工具。

参考资料

• Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley
• Boyce, W. & DiPrima, R. "Elementary Differential Equations", 11th Edition, Wiley
• Simmons, G. "Differential Equations with Applications and Historical Notes", 3rd Edition
• SciPy ODE 求解器官方文档
• SymPy 拉普拉斯变换文档