常微分方程（四）拉普拉斯变换
Chen Kai BOSS
  2019-04-19 15:45:00 2019-04-19 15:45      6.3k 字  33 分钟  

拉普拉斯变换是工程师的秘密武器：它能将令人头疼的微分方程转化为简单的代数方程。从电路分析到控制系统，从信号处理到机械振动，拉普拉斯变换无处不在。本章将揭开这个数学工具的神秘面纱。

从直觉开始：为什么需要拉普拉斯变换？

微分方程的痛点

拉普拉斯变换的几何意义

考虑一个简单的 RC 电路问题：

如果 是阶跃函数、脉冲函数或者任意复杂波形，直接求解非常繁琐。

拉普拉斯变换的魔力： 1. 将微分方程变成代数方程 2. 将卷积变成乘法 3. 自动处理初始条件 4. 统一处理各种输入信号

核心思想：从时域到频域

拉普拉斯变换的本质是将时域函数 映射到复频域函数。Extra close brace or missing open brace\mathcal{L}\{f(t)} = F(s) = \int_0^\infty f(t)e^{-st}dt

其中 是复变量。

直觉： 是一个"探测器"，检测 中各个"频率分量"的含量。

拉普拉斯变换的定义与性质

正式定义

常用函数的拉普拉斯变换表

单边拉普拉斯变换：Extra close brace or missing open braceF(s) = \mathcal{L}\{f(t)} = \int_0^\infty f(t)e^{-st}dt

逆变换：Extra close brace or missing open bracef(t) = \mathcal{L}^{-1}\{F(s)} = \frac{1}{2\pi j}\int_{c-j\infty}^{c+j\infty}F(s)e^{st}ds （实际应用中，逆变换通常用查表或部分分式法）

存在条件 存在的充分条件：

1. 在 上分段连续 2. 存在常数 使得（指数阶）

常用函数的拉普拉斯变换

收敛域
(单位阶跃)
$ (s) > 0e^{at}(s) > at(s) > 0t(s) > 0e^{at}t(s) > ae^{at}t(s) > a(t)$(脉冲)
(延迟阶跃)

推导几个基本变换

例 1： Extra close brace or missing open brace\mathcal{L}\{1} = \int_0^\infty 1 \cdot e^{-st}dt = \left[-\frac{1}{s}e^{-st}\right]_0^\infty = 0 - \left(-\frac{1}{s}\right) = \frac{1}{s} （当 时收敛）

例 2：Extra close brace or missing open brace\mathcal{L}\{e^{at}} Extra close brace or missing open brace\mathcal{L}\{e^{at}} = \int_0^\infty e^{at}e^{-st}dt = \int_0^\infty e^{-(s-a)t}dt = \frac{1}{s-a} （当 时收敛）

例 3： 使用欧拉公式： Extra close brace or missing open brace\mathcal{L}\{\sin \omega t} = \frac{1}{2j}\left(\frac{1}{s-j\omega} - \frac{1}{s+j\omega}\right) = \frac{1}{2j} \cdot \frac{2j\omega}{s^2 + \omega^2} = \frac{\omega}{s^2 + \omega^2} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

def basic_transforms():
    fig, axes = plt.subplots(2, 3, figsize=(15, 8))
    t = np.linspace(0, 5, 500)
    
    # 1. 阶跃函数
    axes[0, 0].plot(t, np.ones_like(t), 'b-', linewidth=2.5)
    axes[0, 0].set_xlabel('t', fontsize=12)
    axes[0, 0].set_ylabel('f(t)', fontsize=12)
    axes[0, 0].set_title('Unit Step:$f(t) = 1$\n$F(s) = 1/s$', fontsize=12, fontweight='bold')
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].set_ylim(-0.5, 2)
    
    # 2. 指数衰减
    a = -1
    axes[0, 1].plot(t, np.exp(a * t), 'r-', linewidth=2.5)
    axes[0, 1].set_xlabel('t', fontsize=12)
    axes[0, 1].set_ylabel('f(t)', fontsize=12)
    axes[0, 1].set_title('$f(t) = e^{-t}$\n$F(s) = 1/(s+1)$', fontsize=12, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3)
    
    # 3. 正弦函数
    omega = 2 * np.pi
    axes[0, 2].plot(t, np.sin(omega * t), 'g-', linewidth=2.5)
    axes[0, 2].set_xlabel('t', fontsize=12)
    axes[0, 2].set_ylabel('f(t)', fontsize=12)
    axes[0, 2].set_title('$f(t) = \\sin(2\\pi t)$\n$F(s) = 2\\pi/(s^2 + 4\\pi^2)$', fontsize=12, fontweight='bold')
    axes[0, 2].grid(True, alpha=0.3)
    
    # 4. t 函数
    axes[1, 0].plot(t, t, 'm-', linewidth=2.5)
    axes[1, 0].set_xlabel('t', fontsize=12)
    axes[1, 0].set_ylabel('f(t)', fontsize=12)
    axes[1, 0].set_title('$f(t) = t$\n$F(s) = 1/s^2$', fontsize=12, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    
    # 5. 衰减正弦
    axes[1, 1].plot(t, np.exp(-t) * np.sin(omega * t), 'c-', linewidth=2.5)
    axes[1, 1].plot(t, np.exp(-t), 'k--', linewidth=1.5, alpha=0.5, label='Envelope')
    axes[1, 1].plot(t, -np.exp(-t), 'k--', linewidth=1.5, alpha=0.5)
    axes[1, 1].set_xlabel('t', fontsize=12)
    axes[1, 1].set_ylabel('f(t)', fontsize=12)
    axes[1, 1].set_title('$f(t) = e^{-t}\\sin(2\\pi t)$\n$F(s) = 2\\pi/((s+1)^2 + 4\\pi^2)$', fontsize=12, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()
    
    # 6. 脉冲函数（近似）
    t_pulse = np.linspace(-0.5, 5, 500)
    epsilon = 0.05
    pulse = np.where(np.abs(t_pulse) < epsilon, 1/(2*epsilon), 0)
    axes[1, 2].plot(t_pulse, pulse, 'orange', linewidth=2.5)
    axes[1, 2].set_xlabel('t', fontsize=12)
    axes[1, 2].set_ylabel('f(t)', fontsize=12)
    axes[1, 2].set_title('$f(t) = \\delta(t)$(approx)\n$F(s) = 1$', fontsize=12, fontweight='bold')
    axes[1, 2].grid(True, alpha=0.3)
    axes[1, 2].set_xlim(-0.5, 2)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig1_basic.png', dpi=150, bbox_inches='tight')
    plt.close()

basic_transforms()

拉普拉斯变换的重要性质

线性性Extra close brace or missing open brace\mathcal{L}\{af(t) + bg(t)} = aF(s) + bG(s)

卷积定理示意图

微分性质（最重要！）Extra close brace or missing open brace\mathcal{L}\{f'(t)} = sF(s) - f(0)

Extra close brace or missing open brace\mathcal{L}\{f''(t)} = s^2F(s) - sf(0) - f'(0)

一般地：Extra close brace or missing open brace\mathcal{L}\{f^{(n)}(t)} = s^nF(s) - s^{n-1}f(0) - s^{n-2}f'(0) - \cdots - f^{(n-1)}(0)

这就是拉普拉斯变换解微分方程的关键：微分变成乘以，初始条件自动包含！

积分性质

频移性质（ 域平移）Extra close brace or missing open brace\mathcal{L}\{e^{at}f(t)} = F(s-a)

应用：如果已知Extra close brace or missing open brace\mathcal{L}\{\cos \omega t} = \frac{s}{s^2+\omega^2}，则：Extra close brace or missing open brace\mathcal{L}\{e^{at}\cos \omega t} = \frac{s-a}{(s-a)^2+\omega^2}

时移性质（延迟）Extra close brace or missing open brace\mathcal{L}\{f(t-a)u(t-a)} = e^{-as}F(s) \quad (a > 0)

其中 是延迟阶跃函数。

卷积定理Extra close brace or missing open brace\mathcal{L}\{f * g} = F(s) \cdot G(s)

其中卷积定义为：

意义：时域的卷积 = 频域的乘法（大大简化计算！）

终值定理

如果 存在且 的所有极点在左半平面（除了可能在 的单极点），则：

初值定理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def laplace_properties():
    """演示拉普拉斯变换的性质"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    t = np.linspace(0, 10, 1000)
    
    # 1. 微分性质演示
    # f(t) = sin(t), f'(t) = cos(t)
    f = np.sin(t)
    f_prime = np.cos(t)
    
    axes[0, 0].plot(t, f, 'b-', linewidth=2, label='$f(t) = \\sin t$')
    axes[0, 0].plot(t, f_prime, 'r--', linewidth=2, label="$f'(t) = \\cos t$")
    axes[0, 0].set_xlabel('t', fontsize=12)
    axes[0, 0].set_ylabel('f(t)', fontsize=12)
    axes[0, 0].set_title('Differentiation Property\n$\\mathcal{L}\\{f\'\}  = sF(s) - f(0)$', fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].text(5, 0.5, '$F(s) = \\frac{1}{s^2+1}$\n$sF(s) - 0 = \\frac{s}{s^2+1}$', fontsize=11,
                   bbox=dict(boxstyle='round', facecolor='wheat'))
    
    # 2. 频移性质演示
    # f(t) = sin(t), e^(-t)sin(t)
    f1 = np.sin(t)
    f2 = np.exp(-t) * np.sin(t)
    
    axes[0, 1].plot(t, f1, 'b-', linewidth=2, label='$\\sin t$')
    axes[0, 1].plot(t, f2, 'r-', linewidth=2, label='$e^{-t}\\sin t$')
    axes[0, 1].plot(t, np.exp(-t), 'k--', linewidth=1, alpha=0.5)
    axes[0, 1].plot(t, -np.exp(-t), 'k--', linewidth=1, alpha=0.5)
    axes[0, 1].set_xlabel('t', fontsize=12)
    axes[0, 1].set_ylabel('f(t)', fontsize=12)
    axes[0, 1].set_title('Frequency Shift Property\n$\\mathcal{L}\\{e^{at}f(t)\}  = F(s-a)$', fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    
    # 3. 时移性质演示
    a = 2
    f_original = np.sin(t)
    f_shifted = np.where(t >= a, np.sin(t - a), 0)
    
    axes[1, 0].plot(t, f_original, 'b-', linewidth=2, label='$f(t) = \\sin t$')
    axes[1, 0].plot(t, f_shifted, 'r--', linewidth=2, label='$f(t-2)u(t-2)$')
    axes[1, 0].axvline(x=a, color='green', linestyle=':', linewidth=2, label=f'Delay = {a}')
    axes[1, 0].set_xlabel('t', fontsize=12)
    axes[1, 0].set_ylabel('f(t)', fontsize=12)
    axes[1, 0].set_title('Time Shift Property\n$\\mathcal{L}\\{f(t-a)u(t-a)\}  = e^{-as}F(s)$', fontsize=12, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # 4. 卷积演示
    # f(t) = exp(-t), g(t) = sin(t)
    # 卷积结果
    f = np.exp(-t)
    g = np.sin(t)
    
    # 数值卷积
    dt = t[1] - t[0]
    conv_result = np.convolve(f, g, mode='full')[:len(t)] * dt
    
    axes[1, 1].plot(t, f, 'b-', linewidth=2, label='$f(t) = e^{-t}$')
    axes[1, 1].plot(t, g, 'g-', linewidth=2, label='$g(t) = \\sin t$')
    axes[1, 1].plot(t, conv_result, 'r-', linewidth=2.5, label='$(f * g)(t)$')
    axes[1, 1].set_xlabel('t', fontsize=12)
    axes[1, 1].set_ylabel('Value', fontsize=12)
    axes[1, 1].set_title('Convolution Theorem\n$\\mathcal{L}\\{f * g\}  = F(s) \\cdot G(s)$', fontsize=12, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig2_properties.png', dpi=150, bbox_inches='tight')
    plt.close()

laplace_properties()

用拉普拉斯变换解微分方程

基本步骤

传递函数与系统响应
  1. 变换：对方程两边取拉普拉斯变换
  2. 代数求解：解出（变成代数方程！）
  3. 逆变换：将 变回

例 1：一阶方程

求解， Step 1：变换Extra close brace or missing open brace\mathcal{L}\{y'} + 2\mathcal{L}\{y} = \mathcal{L}\{e^{-t}}

Step 2：解出

Step 3：逆变换

验证：✓，✓

例 2：二阶方程

求解，， Step 1：变换

Step 2：解出

Step 3：逆变换（需要技巧）

查表或用部分分式得：

注意：这正是共振情况！解中含有 因子，振幅线性增长。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def solve_ode_laplace():
    """用数值方法验证拉普拉斯变换的解"""
    from scipy.integrate import solve_ivp
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 例 1： y' + 2y = exp(-t), y(0) = 1
    def ode1(t, y):
        return np.exp(-t) - 2*y
    
    t_span = [0, 5]
    t_eval = np.linspace(0, 5, 200)
    sol1 = solve_ivp(ode1, t_span, [1], t_eval=t_eval)
    
    # 解析解
    y_analytical1 = np.exp(-t_eval)
    
    axes[0].plot(sol1.t, sol1.y[0], 'b-', linewidth=2.5, label='Numerical')
    axes[0].plot(t_eval, y_analytical1, 'r--', linewidth=2, label='Analytical:$y = e^{-t}$')
    axes[0].set_xlabel('t', fontsize=12)
    axes[0].set_ylabel('y(t)', fontsize=12)
    axes[0].set_title("Example 1:$y' + 2y = e^{-t}$,$y(0) = 1$", fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 例 2： y'' + 4y = sin(2t), y(0) = 0, y'(0) = 0
    def ode2(t, Y):
        y, v = Y
        return [v, np.sin(2*t) - 4*y]
    
    t_span = [0, 15]
    t_eval = np.linspace(0, 15, 500)
    sol2 = solve_ivp(ode2, t_span, [0, 0], t_eval=t_eval)
    
    # 解析解（共振情况）
    y_analytical2 = (1/4) * (np.sin(2*t_eval) - 2*t_eval * np.cos(2*t_eval))
    
    axes[1].plot(sol2.t, sol2.y[0], 'b-', linewidth=2.5, label='Numerical')
    axes[1].plot(t_eval, y_analytical2, 'r--', linewidth=2, label='Analytical (Resonance!)')
    axes[1].plot(t_eval, t_eval/2, 'k:', linewidth=1.5, alpha=0.5, label='Envelope growth')
    axes[1].plot(t_eval, -t_eval/2, 'k:', linewidth=1.5, alpha=0.5)
    axes[1].set_xlabel('t', fontsize=12)
    axes[1].set_ylabel('y(t)', fontsize=12)
    axes[1].set_title("Example 2:$y'' + 4y = \\sin 2t$(Resonance!)", fontsize=12, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig3_solve_ode.png', dpi=150, bbox_inches='tight')
    plt.close()

solve_ode_laplace()

部分分式展开

为什么需要？

逆变换的留数方法

逆拉普拉斯变换时，通常是复杂的有理函数。部分分式展开将其分解为简单项之和，便于查表求逆变换。

情况 1：不同实极点

留数法（覆盖法）：

例：展开

逆变换：

情况 2：重极点

求系数：

例：展开

逆变换：

情况 3：复共轭极点

技巧：配方后查表

逆变换： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def partial_fractions():
    """部分分式展开的例子"""
    from sympy import symbols, apart, inverse_laplace_transform, exp, sin, cos
    
    s, t = symbols('s t', real=True, positive=True)
    
    # 例子们
    examples = [
        ("Different real poles", (3*s + 5)/((s + 1)*(s + 2))),
        ("Repeated pole", (s + 3)/(s + 1)**2),
        ("Complex poles", (s + 2)/(s**2 + 2*s + 5)),
    ]
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    t_vals = np.linspace(0, 5, 200)
    
    for ax, (title, F_s) in zip(axes, examples):
        # 数值逆变换的近似（用 scipy）
        # 这里我们直接给出解析解
        pass
    
    # 例 1：不同实极点
    f1 = 2*np.exp(-t_vals) + np.exp(-2*t_vals)
    axes[0].plot(t_vals, f1, 'b-', linewidth=2.5)
    axes[0].set_xlabel('t', fontsize=12)
    axes[0].set_ylabel('f(t)', fontsize=12)
    axes[0].set_title('Different Real Poles\n$F(s) = \\frac{3s+5}{(s+1)(s+2)}$\n$f(t) = 2e^{-t} + e^{-2t}$', 
                     fontsize=10, fontweight='bold')
    axes[0].grid(True, alpha=0.3)
    
    # 例 2：重极点
    f2 = (1 + 2*t_vals) * np.exp(-t_vals)
    axes[1].plot(t_vals, f2, 'r-', linewidth=2.5)
    axes[1].set_xlabel('t', fontsize=12)
    axes[1].set_ylabel('f(t)', fontsize=12)
    axes[1].set_title('Repeated Pole\n$F(s) = \\frac{s+3}{(s+1)^2}$\n$f(t) = (1+2t)e^{-t}$', 
                     fontsize=10, fontweight='bold')
    axes[1].grid(True, alpha=0.3)
    
    # 例 3：复共轭极点
    # (s+2)/(s^2+2s+5) = (s+2)/((s+1)^2+4)
    # = (s+1)/((s+1)^2+4) + 1/((s+1)^2+4)
    # = e^{-t}cos(2t) + (1/2)e^{-t}sin(2t)
    f3 = np.exp(-t_vals) * (np.cos(2*t_vals) + 0.5*np.sin(2*t_vals))
    axes[2].plot(t_vals, f3, 'g-', linewidth=2.5)
    axes[2].plot(t_vals, np.exp(-t_vals), 'k--', linewidth=1, alpha=0.5, label='Envelope')
    axes[2].plot(t_vals, -np.exp(-t_vals), 'k--', linewidth=1, alpha=0.5)
    axes[2].set_xlabel('t', fontsize=12)
    axes[2].set_ylabel('f(t)', fontsize=12)
    axes[2].set_title('Complex Conjugate Poles\n$F(s) = \\frac{s+2}{s^2+2s+5}$\n$f(t) = e^{-t}(\\cos 2t + 0.5\\sin 2t)$', 
                     fontsize=10, fontweight='bold')
    axes[2].legend()
    axes[2].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig4_partial_fractions.png', dpi=150, bbox_inches='tight')
    plt.close()

partial_fractions()

传递函数与系统分析

什么是传递函数？

对于线性时不变（ LTI）系统，传递函数定义为：输出的拉普拉斯变换输入的拉普拉斯变换 （假设初始条件为零）

例： RC 电路

取拉普拉斯变换（零初始条件）：

极点与零点

  • 零点：使 的 值
  • 极点：使 的 值

极点决定系统特性： - 极点在左半平面：稳定 - 极点在右半平面：不稳定 - 极点在虚轴：临界稳定（持续振荡）

极点与时域响应的关系

极点类型 时域响应
负实极点 （衰减）
正实极点 （发散）
复极点 （衰减振荡）
纯虚极点 （等幅振荡） 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def pole_zero_analysis():
    """极点零点分析"""
    
    fig = plt.figure(figsize=(16, 10))
    
    # 创建子图
    ax1 = fig.add_subplot(2, 3, 1)
    ax2 = fig.add_subplot(2, 3, 2)
    ax3 = fig.add_subplot(2, 3, 3)
    ax4 = fig.add_subplot(2, 3, 4)
    ax5 = fig.add_subplot(2, 3, 5)
    ax6 = fig.add_subplot(2, 3, 6)
    
    # 例子系统
    systems = [
        ("Stable (2 real poles)", [-1, -2], [], ax1, ax2),
        ("Underdamped", [-1+2j, -1-2j], [], ax3, ax4),
        ("Unstable", [1], [], ax5, ax6),
    ]
    
    t = np.linspace(0, 10, 500)
    
    for name, poles, zeros, ax_pz, ax_resp in systems:
        # 极零点图
        ax_pz.axhline(y=0, color='k', linewidth=0.5)
        ax_pz.axvline(x=0, color='k', linewidth=0.5)
        
        # 画极点
        for p in poles:
            ax_pz.plot(np.real(p), np.imag(p), 'rx', markersize=15, markeredgewidth=3)
        
        # 画零点
        for z in zeros:
            ax_pz.plot(np.real(z), np.imag(z), 'bo', markersize=12, markerfacecolor='none', markeredgewidth=2)
        
        ax_pz.set_xlabel('Re(s)', fontsize=10)
        ax_pz.set_ylabel('Im(s)', fontsize=10)
        ax_pz.set_title(f'{name}\nPole-Zero Map', fontsize=10, fontweight='bold')
        ax_pz.grid(True, alpha=0.3)
        ax_pz.set_xlim(-3, 2)
        ax_pz.set_ylim(-3, 3)
        
        # 填充不稳定区域
        ax_pz.axvspan(0, 2, alpha=0.1, color='red')
        ax_pz.text(0.5, 2.5, 'Unstable', fontsize=9, color='red')
        
        # 冲激响应
        if name == "Stable (2 real poles)":
            h = np.exp(-t) - np.exp(-2*t)
        elif name == "Underdamped":
            h = np.exp(-t) * np.sin(2*t)
        else:
            h = np.exp(t)
            h[h > 10] = np.nan  # 限制显示
        
        ax_resp.plot(t, h, 'b-', linewidth=2)
        ax_resp.axhline(y=0, color='k', linewidth=0.5)
        ax_resp.set_xlabel('t', fontsize=10)
        ax_resp.set_ylabel('h(t)', fontsize=10)
        ax_resp.set_title('Impulse Response', fontsize=10, fontweight='bold')
        ax_resp.grid(True, alpha=0.3)
        
        if name == "Unstable":
            ax_resp.set_ylim(-1, 15)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig5_poles.png', dpi=150, bbox_inches='tight')
    plt.close()

pole_zero_analysis()

频率响应与 Bode 图

频率响应函数

将 代入传递函数得到频率响应： -：幅频特性 -：相频特性

Bode 图

Bode 图是工程中最常用的频率响应表示方法： - 横轴：频率（对数刻度） - 纵轴（幅值）：（ dB） - 纵轴（相位）：（度）

常见环节的 Bode 图

一阶低通：

转角频率： 在 处： - 幅值下降 3dB - 相位为-45 °

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def bode_plot():
    """Bode 图演示"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 频率范围
    omega = np.logspace(-2, 2, 500)
    
    # 一阶低通： H(s) = 1/(1+s)
    tau = 1
    H1 = 1 / (1 + 1j*omega*tau)
    mag1 = 20 * np.log10(np.abs(H1))
    phase1 = np.angle(H1, deg=True)
    
    axes[0, 0].semilogx(omega, mag1, 'b-', linewidth=2.5, label='Exact')
    # 渐近线
    omega_c = 1/tau
    asymp_low = np.zeros_like(omega[omega < omega_c])
    asymp_high = -20 * np.log10(omega[omega >= omega_c] * tau)
    axes[0, 0].semilogx(omega[omega < omega_c], asymp_low, 'r--', linewidth=2, label='Asymptote')
    axes[0, 0].semilogx(omega[omega >= omega_c], asymp_high, 'r--', linewidth=2)
    axes[0, 0].axvline(x=omega_c, color='green', linestyle=':', label=f'Corner freq = {omega_c}')
    axes[0, 0].axhline(y=-3, color='gray', linestyle=':', alpha=0.5)
    axes[0, 0].set_xlabel('Frequency ω (rad/s)', fontsize=12)
    axes[0, 0].set_ylabel('Magnitude (dB)', fontsize=12)
    axes[0, 0].set_title('First-Order Low-Pass: Magnitude\n$H(s) = 1/(1+s)$', fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3, which='both')
    
    axes[0, 1].semilogx(omega, phase1, 'b-', linewidth=2.5)
    axes[0, 1].axhline(y=-45, color='green', linestyle=':', label='Phase at corner')
    axes[0, 1].axhline(y=-90, color='gray', linestyle=':', alpha=0.5)
    axes[0, 1].set_xlabel('Frequency ω (rad/s)', fontsize=12)
    axes[0, 1].set_ylabel('Phase (degrees)', fontsize=12)
    axes[0, 1].set_title('First-Order Low-Pass: Phase', fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3, which='both')
    
    # 二阶系统： H(s) = ω_n^2 / (s^2 + 2 ζω_{n*}s + ω_n^2)
    omega_n = 1
    zetas = [0.1, 0.3, 0.5, 0.7, 1.0, 2.0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(zetas)))
    
    for zeta, color in zip(zetas, colors):
        H2 = omega_{n*}*2 / (omega_{n*}*2 - omega**2 + 2j*zeta*omega_{n*}omega)
        mag2 = 20 * np.log10(np.abs(H2))
        phase2 = np.angle(H2, deg=True)
        
        axes[1, 0].semilogx(omega, mag2, color=color, linewidth=2, label=f'ζ = {zeta}')
        axes[1, 1].semilogx(omega, phase2, color=color, linewidth=2, label=f'ζ = {zeta}')
    
    axes[1, 0].axvline(x=omega_n, color='red', linestyle=':', label='Natural freq')
    axes[1, 0].set_xlabel('Frequency ω (rad/s)', fontsize=12)
    axes[1, 0].set_ylabel('Magnitude (dB)', fontsize=12)
    axes[1, 0].set_title('Second-Order System: Magnitude', fontsize=12, fontweight='bold')
    axes[1, 0].legend(fontsize=9)
    axes[1, 0].grid(True, alpha=0.3, which='both')
    axes[1, 0].set_ylim(-40, 30)
    
    axes[1, 1].axhline(y=-90, color='gray', linestyle=':', alpha=0.5)
    axes[1, 1].set_xlabel('Frequency ω (rad/s)', fontsize=12)
    axes[1, 1].set_ylabel('Phase (degrees)', fontsize=12)
    axes[1, 1].set_title('Second-Order System: Phase', fontsize=12, fontweight='bold')
    axes[1, 1].legend(fontsize=9)
    axes[1, 1].grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig6_bode.png', dpi=150, bbox_inches='tight')
    plt.close()

bode_plot()

处理不连续输入：阶跃函数和脉冲函数

单位阶跃函数

Extra close brace or missing open brace\mathcal{L}\{u(t)} = \frac{1}{s}

单位脉冲函数（ Dirac delta）

性质： - 当 - - Extra close brace or missing open brace\mathcal{L}\{\delta(t)} = 1

脉冲响应与阶跃响应的关系

脉冲响应： 阶跃响应： 关系： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def step_impulse_response():
    """阶跃响应和脉冲响应"""
    from scipy import signal
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 一阶系统： H(s) = 1/(s+1)
    num1 = [1]
    den1 = [1, 1]
    sys1 = signal.TransferFunction(num1, den1)
    
    t1, y1_step = signal.step(sys1)
    t1, y1_impulse = signal.impulse(sys1)
    
    axes[0, 0].plot(t1, y1_step, 'b-', linewidth=2.5, label='Step Response')
    axes[0, 0].axhline(y=1, color='r', linestyle='--', alpha=0.5)
    axes[0, 0].set_xlabel('Time (s)', fontsize=12)
    axes[0, 0].set_ylabel('Response', fontsize=12)
    axes[0, 0].set_title('First-Order System: Step Response\n$H(s) = 1/(s+1)$', fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    
    axes[0, 1].plot(t1, y1_impulse, 'r-', linewidth=2.5, label='Impulse Response')
    axes[0, 1].set_xlabel('Time (s)', fontsize=12)
    axes[0, 1].set_ylabel('Response', fontsize=12)
    axes[0, 1].set_title('First-Order System: Impulse Response', fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    
    # 二阶欠阻尼系统： H(s) = 4/(s^2 + 0.8s + 4)
    omega_n = 2
    zeta = 0.2
    num2 = [omega_{n*}*2]
    den2 = [1, 2*zeta*omega_n, omega_{n*}*2]
    sys2 = signal.TransferFunction(num2, den2)
    
    t2, y2_step = signal.step(sys2, T=np.linspace(0, 15, 500))
    t2, y2_impulse = signal.impulse(sys2, T=np.linspace(0, 15, 500))
    
    axes[1, 0].plot(t2, y2_step, 'b-', linewidth=2.5, label='Step Response')
    axes[1, 0].axhline(y=1, color='r', linestyle='--', alpha=0.5)
    # 过冲标注
    overshoot = (np.max(y2_step) - 1) * 100
    t_peak = t2[np.argmax(y2_step)]
    axes[1, 0].plot(t_peak, np.max(y2_step), 'go', markersize=10)
    axes[1, 0].annotate(f'Overshoot: {overshoot:.1f}%', xy=(t_peak, np.max(y2_step)),
                       xytext=(t_peak+1, np.max(y2_step)+0.1), fontsize=10,
                       arrowprops=dict(arrowstyle='->', color='green'))
    axes[1, 0].set_xlabel('Time (s)', fontsize=12)
    axes[1, 0].set_ylabel('Response', fontsize=12)
    axes[1, 0].set_title(f'Second-Order System (ζ={zeta}): Step Response', fontsize=12, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    axes[1, 1].plot(t2, y2_impulse, 'r-', linewidth=2.5, label='Impulse Response')
    axes[1, 1].set_xlabel('Time (s)', fontsize=12)
    axes[1, 1].set_ylabel('Response', fontsize=12)
    axes[1, 1].set_title('Second-Order System: Impulse Response', fontsize=12, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig7_step_impulse.png', dpi=150, bbox_inches='tight')
    plt.close()

step_impulse_response()

控制系统应用： PID 控制器

PID 控制器简介

PID 控制器是工业中最常用的控制器： -：比例增益（ Proportional） -：积分增益（ Integral） -：微分增益（ Derivative）

拉普拉斯域表示

各部分的作用

项 作用 优点 缺点
P（比例） 与误差成比例 快速响应 有稳态误差
I（积分） 累积误差 消除稳态误差 可能振荡
D（微分） 预测误差趋势 减少过冲 放大噪声 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
def pid_control():
    """PID 控制器演示"""
    from scipy import signal
    from scipy.integrate import solve_ivp
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 被控对象：二阶系统 G(s) = 1/(s^2 + s + 1)
    # 用单位反馈
    
    # 仿真 PID 控制
    def simulate_pid(Kp, Ki, Kd, t_span, setpoint=1):
        """简化的 PID 仿真"""
        def system(t, state):
            y, v, integral = state
            e = setpoint - y
            
            # PID 控制律
            u = Kp * e + Ki * integral + Kd * (-v)  # 用-v 近似 de/dt
            
            # 被控对象： y'' + y' + y = u
            dydt = v
            dvdt = u - v - y
            dintegral = e
            
            return [dydt, dvdt, dintegral]
        
        sol = solve_ivp(system, t_span, [0, 0, 0], 
                       t_eval=np.linspace(t_span[0], t_span[1], 500),
                       method='RK45')
        return sol.t, sol.y[0]
    
    t_span = [0, 20]
    
    # 仅 P 控制
    t1, y1 = simulate_pid(2, 0, 0, t_span)
    axes[0, 0].plot(t1, y1, 'b-', linewidth=2, label='P control (Kp=2)')
    axes[0, 0].axhline(y=1, color='r', linestyle='--', alpha=0.5, label='Setpoint')
    axes[0, 0].set_xlabel('Time (s)', fontsize=12)
    axes[0, 0].set_ylabel('Output', fontsize=12)
    axes[0, 0].set_title('P Control: Steady-State Error', fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    
    # PI 控制
    t2, y2 = simulate_pid(2, 1, 0, t_span)
    axes[0, 1].plot(t2, y2, 'g-', linewidth=2, label='PI control (Kp=2, Ki=1)')
    axes[0, 1].axhline(y=1, color='r', linestyle='--', alpha=0.5, label='Setpoint')
    axes[0, 1].set_xlabel('Time (s)', fontsize=12)
    axes[0, 1].set_ylabel('Output', fontsize=12)
    axes[0, 1].set_title('PI Control: Eliminates Steady-State Error', fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    
    # PID 控制
    t3, y3 = simulate_pid(2, 1, 0.5, t_span)
    axes[1, 0].plot(t3, y3, 'm-', linewidth=2, label='PID control (Kp=2, Ki=1, Kd=0.5)')
    axes[1, 0].axhline(y=1, color='r', linestyle='--', alpha=0.5, label='Setpoint')
    axes[1, 0].set_xlabel('Time (s)', fontsize=12)
    axes[1, 0].set_ylabel('Output', fontsize=12)
    axes[1, 0].set_title('PID Control: Reduced Overshoot', fontsize=12, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # 比较
    axes[1, 1].plot(t1, y1, 'b-', linewidth=2, label='P only')
    axes[1, 1].plot(t2, y2, 'g-', linewidth=2, label='PI')
    axes[1, 1].plot(t3, y3, 'm-', linewidth=2, label='PID')
    axes[1, 1].axhline(y=1, color='r', linestyle='--', alpha=0.5, label='Setpoint')
    axes[1, 1].set_xlabel('Time (s)', fontsize=12)
    axes[1, 1].set_ylabel('Output', fontsize=12)
    axes[1, 1].set_title('Comparison of P, PI, PID Control', fontsize=12, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 4 章-拉普拉斯变换/fig8_pid.png', dpi=150, bbox_inches='tight')
    plt.close()

pid_control()

Python 实战：符号计算与数值仿真

使用 SymPy 进行符号计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
from sympy import symbols, Function, laplace_transform, inverse_laplace_transform
from sympy import exp, sin, cos, Heaviside, DiracDelta, simplify
from sympy import apart, factor

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

# 基本变换
print("=== 基本拉普拉斯变换 ===")
print(f"L{{ 1 }} = {laplace_transform(1, t, s)}")
print(f"L{{ t }} = {laplace_transform(t, t, s)}")
print(f"L{{ exp(-at) }} = {laplace_transform(exp(-a*t), t, s)}")
print(f"L{{ sin(ω t) }} = {laplace_transform(sin(omega*t), t, s)}")

# 逆变换
print("\n=== 逆拉普拉斯变换 ===")
F1 = 1 / (s + 2)
print(f"L^(-1){{ 1/(s+2) }} = {inverse_laplace_transform(F1, s, t)}")

F2 = (s + 3) / (s**2 + 4*s + 5)
print(f"L^(-1){{ (s+3)/(s^2+4s+5) }} = {simplify(inverse_laplace_transform(F2, s, t))}")

# 部分分式
F3 = (3*s + 5) / ((s + 1) * (s + 2))
print(f"\n 部分分式展开: {apart(F3, s)}")

使用 SciPy 进行系统分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# 定义传递函数： H(s) = (s+2) / (s^2 + 3s + 2)
num = [1, 2]  # 分子系数
den = [1, 3, 2]  # 分母系数

sys = signal.TransferFunction(num, den)

# 极点和零点
poles = np.roots(den)
zeros = np.roots(num)
print(f"极点: {poles}")
print(f"零点: {zeros}")

# 阶跃响应
t, y = signal.step(sys)

# Bode 图
w, mag, phase = signal.bode(sys)

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# 极零点图
axes[0, 0].plot(np.real(poles), np.imag(poles), 'rx', markersize=12, markeredgewidth=2, label='Poles')
axes[0, 0].plot(np.real(zeros), np.imag(zeros), 'bo', markersize=10, markerfacecolor='none', markeredgewidth=2, label='Zeros')
axes[0, 0].axhline(y=0, color='k', linewidth=0.5)
axes[0, 0].axvline(x=0, color='k', linewidth=0.5)
axes[0, 0].set_xlabel('Real', fontsize=12)
axes[0, 0].set_ylabel('Imaginary', fontsize=12)
axes[0, 0].set_title('Pole-Zero Map', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_aspect('equal')

# 阶跃响应
axes[0, 1].plot(t, y, 'b-', linewidth=2)
axes[0, 1].set_xlabel('Time (s)', fontsize=12)
axes[0, 1].set_ylabel('Response', fontsize=12)
axes[0, 1].set_title('Step Response', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)

# Bode 幅值
axes[1, 0].semilogx(w, mag, 'b-', linewidth=2)
axes[1, 0].set_xlabel('Frequency (rad/s)', fontsize=12)
axes[1, 0].set_ylabel('Magnitude (dB)', fontsize=12)
axes[1, 0].set_title('Bode Plot: Magnitude', fontsize=12, fontweight='bold')
axes[1, 0].grid(True, alpha=0.3, which='both')

# Bode 相位
axes[1, 1].semilogx(w, phase, 'r-', linewidth=2)
axes[1, 1].set_xlabel('Frequency (rad/s)', fontsize=12)
axes[1, 1].set_ylabel('Phase (degrees)', fontsize=12)
axes[1, 1].set_title('Bode Plot: Phase', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('第 4 章-拉普拉斯变换/fig9_scipy.png', dpi=150, bbox_inches='tight')
plt.close()

总结

核心概念

概念 定义 意义
拉普拉斯变换 时域→频域
逆变换 查表或部分分式 频域→时域
传递函数 系统特性
极点 的分母零点 决定稳定性
零点 的分子零点 影响频率响应

重要性质

性质 公式 应用
微分 Extra close brace or missing open brace\mathcal{L}\{f'} = sF(s) - f(0) 解 ODE
积分 Extra close brace or missing open brace\mathcal{L}\{\int f} = F(s)/s 处理积分方程
频移 Extra close brace or missing open brace\mathcal{L}\{e^{at}f} = F(s-a) 阻尼系统
时移 Extra close brace or missing open brace\mathcal{L}\{f(t-a)u(t-a)} = e^{-as}F(s) 延迟系统
卷积 Extra close brace or missing open brace\mathcal{L}\{f*g} = F \cdot G 系统响应

解题流程

  1. 对方程取拉普拉斯变换
  2. 解代数方程得
  3. 部分分式展开
  4. 查表求逆变换
  5. 验证答案

练习题

基础题

  1. 求以下函数的拉普拉斯变换： - - -$f(t) = (t-1)u(t-1)求以下F(s)$ 的逆变换： - - -$F(s) = 用拉普拉斯变换求解：y' + 3y = e^{-2t}，y(0) = 1验证微分性质：{tf(t)} = -F'(s)$

进阶题

  1. 用拉普拉斯变换求解：

  2. RC 电路，初始：

    • 当 时，求
    • 当 时，求稳态响应

  3. 证明卷积定理：

  4. 求传递函数 的：

    • 极点和零点
    • 冲激响应
    • 阶跃响应

  5. 设计一个 PID 控制器，使二阶系统 的闭环响应满足：

    • 过冲小于 10%
    • 稳态误差为 0

  6. 用终值定理求，其中

编程题

  1. 用 Python 实现拉普拉斯逆变换的数值算法（ Talbot 方法）。

  2. 编写程序绘制任意传递函数的 Bode 图。

  3. 仿真带有时滞的系统： 的阶跃响应。

  4. 用 Python 实现 Ziegler-Nichols 方法自动整定 PID 参数。

参考资料

  1. Kreyszig, E. (2011). Advanced Engineering Mathematics. Wiley.
  2. Ogata, K. (2010). Modern Control Engineering. Pearson.
  3. Oppenheim, A. V., & Willsky, A. S. (1997). Signals and Systems. Prentice Hall.
  4. SciPy Documentation: scipy.signal
  5. MIT OCW 18.03 - Differential Equations

下一章预告：《微分方程组与矩阵指数》 — 从单个方程到耦合系统

  • 本文标题：常微分方程（四）拉普拉斯变换
  • 本文作者：Chen Kai
  • 创建时间：2019-04-19 15:45:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%9B%9B%EF%BC%89%E6%8B%89%E6%99%AE%E6%8B%89%E6%96%AF%E5%8F%98%E6%8D%A2/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论