常微分方程（三）高阶线性微分方程
Chen Kai BOSS
  2019-04-14 10:15:00 2019-04-14 10:15      5.8k 字  29 分钟  

当一阶方程不足以描述系统时，我们需要高阶微分方程。弹簧的振动、桥梁的摇晃、电路的谐振——这些现象都需要二阶或更高阶的 ODE 来建模。本章将系统讲解高阶线性 ODE 的理论和求解方法。

从物理直觉开始：弹簧-质量系统

最简单的振动模型

弹簧振动系统示意图

想象一个质量为 的物体挂在弹簧上。根据胡克定律，弹簧的恢复力与位移成正比：

其中 是弹簧常数，负号表示力的方向与位移相反。

根据牛顿第二定律：

整理得到：

其中 是自然角频率。

这是一个二阶常系数齐次线性 ODE！

直觉理解

问题：这个方程的解是什么？

猜测：什么函数的二阶导数等于它自己乘以一个负数？

答案：正弦和余弦函数！

通解：

或等价地写成：

其中 是振幅， 是初相位。

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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

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

def simple_harmonic():
    omega0 = 2 * np.pi  # 频率 1Hz
    t = np.linspace(0, 3, 500)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 不同初始条件
    A, phi = 1.0, 0
    x = A * np.cos(omega0 * t + phi)
    v = -A * omega0 * np.sin(omega0 * t + phi)
    
    # 位移-时间图
    axes[0, 0].plot(t, x, 'b-', linewidth=2.5)
    axes[0, 0].set_xlabel('Time (s)', fontsize=12)
    axes[0, 0].set_ylabel('Displacement x (m)', fontsize=12)
    axes[0, 0].set_title('Simple Harmonic Motion', fontsize=14, fontweight='bold')
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].axhline(y=0, color='k', linewidth=0.5)
    
    # 速度-时间图
    axes[0, 1].plot(t, v, 'r-', linewidth=2.5)
    axes[0, 1].set_xlabel('Time (s)', fontsize=12)
    axes[0, 1].set_ylabel('Velocity v (m/s)', fontsize=12)
    axes[0, 1].set_title('Velocity vs Time', fontsize=14, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].axhline(y=0, color='k', linewidth=0.5)
    
    # 相空间轨迹
    axes[1, 0].plot(x, v, 'g-', linewidth=2.5)
    axes[1, 0].plot(x[0], v[0], 'ro', markersize=10, label='Start')
    axes[1, 0].set_xlabel('Displacement x (m)', fontsize=12)
    axes[1, 0].set_ylabel('Velocity v (m/s)', fontsize=12)
    axes[1, 0].set_title('Phase Portrait (Ellipse)', fontsize=14, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].set_aspect('equal')
    axes[1, 0].legend()
    
    # 能量守恒
    KE = 0.5 * 1.0 * v**2  # 动能，假设 m=1
    PE = 0.5 * (omega0**2) * x**2  # 势能
    TE = KE + PE
    
    axes[1, 1].plot(t, KE, 'r-', linewidth=2, label='Kinetic Energy')
    axes[1, 1].plot(t, PE, 'b-', linewidth=2, label='Potential Energy')
    axes[1, 1].plot(t, TE, 'k--', linewidth=2, label='Total Energy')
    axes[1, 1].set_xlabel('Time (s)', fontsize=12)
    axes[1, 1].set_ylabel('Energy (J)', fontsize=12)
    axes[1, 1].set_title('Energy Conservation', fontsize=14, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig1_shm.png', dpi=150, bbox_inches='tight')
    plt.close()

simple_harmonic()

高阶线性 ODE 的一般理论

定义与标准形式

欠阻尼、临界阻尼和过阻尼对比

阶线性 ODE：

标准形式（首项系数为 1）：

齐次方程： 非齐次方程：

解的结构

定理 1（叠加原理）：如果 和 是齐次方程的解，则 也是解。

定理 2（通解结构）：阶齐次线性 ODE 的通解由 个线性无关的特解构成：

定理 3（非齐次方程）：非齐次方程的通解 = 齐次方程的通解 + 一个特解

线性无关与 Wronskian

定义：函数 在区间 上线性无关，如果

只有 时成立。

Wronskian 行列式：

定理：如果 在某点，则 线性无关。

例：验证 和 线性无关

所以它们线性无关。

常系数齐次方程：特征方程法

核心思想

强迫振动的共振现象

考虑阶常系数齐次线性 ODE：

关键洞察：尝试 的形式

代入方程：

由于，得到特征方程：

三种情况

情况 1： 个不同实根

如果特征方程有 个不同实根：

例 1：求解 特征方程： ，得 通解：

情况 2：有重根

如果 是重根，对应的解是：

例 2：求解 特征方程： ， 是二重根

通解：

情况 3：复根

如果有复根（成对出现），对应的实值解是：

例 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
60
61
62
63
64
65
66
67
68
def characteristic_equation_cases():
    x = np.linspace(0, 5, 300)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Case 1: Different real roots
    # y'' - 5y' + 6y = 0, y = c1*e^2x + c2*e^3x
    y1 = np.exp(2*x)
    y2 = np.exp(3*x)
    axes[0, 0].plot(x, y1, 'b-', linewidth=2, label='$e^{2x}$')
    axes[0, 0].plot(x, y2, 'r-', linewidth=2, label='$e^{3x}$')
    axes[0, 0].plot(x, y1 - y2, 'g--', linewidth=2, label='$e^{2x} - e^{3x}$')
    axes[0, 0].set_xlabel('x', fontsize=12)
    axes[0, 0].set_ylabel('y', fontsize=12)
    axes[0, 0].set_title("Case 1: Distinct Real Roots\n$y'' - 5y' + 6y = 0$", fontsize=12, fontweight='bold')
    axes[0, 0].legend()
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].set_ylim(-10, 50)
    
    # Case 2: Repeated roots
    # y'' - 4y' + 4y = 0, y = (c1 + c2*x)e^2x
    y1 = np.exp(2*x)
    y2 = x * np.exp(2*x)
    axes[0, 1].plot(x, y1, 'b-', linewidth=2, label='$e^{2x}$')
    axes[0, 1].plot(x, y2, 'r-', linewidth=2, label='$x e^{2x}$')
    axes[0, 1].plot(x, (1 + x) * np.exp(2*x), 'g--', linewidth=2, label='$(1+x)e^{2x}$')
    axes[0, 1].set_xlabel('x', fontsize=12)
    axes[0, 1].set_ylabel('y', fontsize=12)
    axes[0, 1].set_title("Case 2: Repeated Root\n$y'' - 4y' + 4y = 0$", fontsize=12, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].set_ylim(-10, 100)
    
    # Case 3: Complex roots - underdamped
    # y'' + 2y' + 5y = 0, r = -1 ± 2i
    y1 = np.exp(-x) * np.cos(2*x)
    y2 = np.exp(-x) * np.sin(2*x)
    env = np.exp(-x)
    
    axes[1, 0].plot(x, y1, 'b-', linewidth=2, label='$e^{-x}\\cos 2x$')
    axes[1, 0].plot(x, y2, 'r-', linewidth=2, label='$e^{-x}\\sin 2x$')
    axes[1, 0].plot(x, env, 'k--', linewidth=1.5, alpha=0.5, label='Envelope$e^{-x}$')
    axes[1, 0].plot(x, -env, 'k--', linewidth=1.5, alpha=0.5)
    axes[1, 0].set_xlabel('x', fontsize=12)
    axes[1, 0].set_ylabel('y', fontsize=12)
    axes[1, 0].set_title("Case 3: Complex Roots (Damped Oscillation)\n$y'' + 2y' + 5y = 0$", fontsize=12, fontweight='bold')
    axes[1, 0].legend()
    axes[1, 0].grid(True, alpha=0.3)
    
    # Case 3b: Pure imaginary roots - undamped
    # y'' + 4y = 0, r = ± 2i
    y1 = np.cos(2*x)
    y2 = np.sin(2*x)
    
    axes[1, 1].plot(x, y1, 'b-', linewidth=2, label='$\\cos 2x$')
    axes[1, 1].plot(x, y2, 'r-', linewidth=2, label='$\\sin 2x$')
    axes[1, 1].plot(x, y1 + 0.5*y2, 'g--', linewidth=2, label='$\\cos 2x + 0.5\\sin 2x$')
    axes[1, 1].set_xlabel('x', fontsize=12)
    axes[1, 1].set_ylabel('y', fontsize=12)
    axes[1, 1].set_title("Case 3b: Pure Imaginary Roots (Undamped)\n$y'' + 4y = 0$", fontsize=12, fontweight='bold')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig2_characteristic.png', dpi=150, bbox_inches='tight')
    plt.close()

characteristic_equation_cases()

阻尼振动：完整分析

物理模型

RLC电路分析

考虑有阻尼的弹簧系统： -：质量 -：阻尼系数 -：弹簧常数

标准化为：

其中： -：自然频率 -：阻尼比（ damping ratio）

特征方程分析

特征方程：

根据 的值，分为三种情况：

欠阻尼（ Underdamped）：

，其中 是阻尼频率

特点：振荡并逐渐衰减

临界阻尼（ Critically Damped）：

（重根）

特点：最快回到平衡位置，无振荡

过阻尼（ Overdamped）：

（两个负实根）

特点：缓慢回到平衡位置，无振荡

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
def damping_comparison():
    omega0 = 2 * np.pi  # 自然频率
    t = np.linspace(0, 4, 500)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 不同阻尼比
    zeta_values = [0.1, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0]
    colors = plt.cm.viridis(np.linspace(0, 1, len(zeta_values)))
    
    for zeta, color in zip(zeta_values, colors):
        if zeta < 1:  # 欠阻尼
            omega_d = omega0 * np.sqrt(1 - zeta**2)
            x = np.exp(-zeta * omega0 * t) * np.cos(omega_d * t)
            label = f'ζ = {zeta} (underdamped)'
        elif zeta == 1:  # 临界阻尼
            x = (1 + omega0 * t) * np.exp(-omega0 * t)
            label = f'ζ = {zeta} (critical)'
        else:  # 过阻尼
            r1 = -zeta * omega0 + omega0 * np.sqrt(zeta**2 - 1)
            r2 = -zeta * omega0 - omega0 * np.sqrt(zeta**2 - 1)
            # 初始条件 x(0)=1, x'(0)=0
            A = r2 / (r2 - r1)
            B = -r1 / (r2 - r1)
            x = A * np.exp(r1 * t) + B * np.exp(r2 * t)
            label = f'ζ = {zeta} (overdamped)'
        
        axes[0].plot(t, x, color=color, linewidth=2, label=label)
    
    axes[0].axhline(y=0, color='k', linewidth=0.5)
    axes[0].set_xlabel('Time (s)', fontsize=12)
    axes[0].set_ylabel('Displacement x', fontsize=12)
    axes[0].set_title('Effect of Damping Ratio ζ', fontsize=14, fontweight='bold')
    axes[0].legend(fontsize=9, loc='upper right')
    axes[0].grid(True, alpha=0.3)
    axes[0].set_xlim(0, 4)
    
    # 相空间对比
    for zeta, color in zip([0.2, 1.0, 2.0], ['blue', 'green', 'red']):
        if zeta < 1:
            omega_d = omega0 * np.sqrt(1 - zeta**2)
            x = np.exp(-zeta * omega0 * t) * np.cos(omega_d * t)
            v = -np.exp(-zeta * omega0 * t) * (zeta * omega0 * np.cos(omega_d * t) + omega_d * np.sin(omega_d * t))
            label = f'ζ = {zeta} (underdamped)'
        elif zeta == 1:
            x = (1 + omega0 * t) * np.exp(-omega0 * t)
            v = -omega0**2 * t * np.exp(-omega0 * t)
            label = f'ζ = {zeta} (critical)'
        else:
            r1 = -zeta * omega0 + omega0 * np.sqrt(zeta**2 - 1)
            r2 = -zeta * omega0 - omega0 * np.sqrt(zeta**2 - 1)
            A = r2 / (r2 - r1)
            B = -r1 / (r2 - r1)
            x = A * np.exp(r1 * t) + B * np.exp(r2 * t)
            v = A * r1 * np.exp(r1 * t) + B * r2 * np.exp(r2 * t)
            label = f'ζ = {zeta} (overdamped)'
        
        axes[1].plot(x, v, color=color, linewidth=2, label=label)
        axes[1].plot(x[0], v[0], 'o', color=color, markersize=8)
    
    axes[1].plot(0, 0, 'k*', markersize=15, label='Equilibrium')
    axes[1].set_xlabel('Displacement x', fontsize=12)
    axes[1].set_ylabel('Velocity v', fontsize=12)
    axes[1].set_title('Phase Portrait', fontsize=14, fontweight='bold')
    axes[1].legend(fontsize=10)
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig3_damping.png', dpi=150, bbox_inches='tight')
    plt.close()

damping_comparison()

工程应用：汽车悬挂

目标：设计悬挂系统，使汽车在遇到路面颠簸后能平稳、快速地恢复。

最优选择：略微欠阻尼（）

  • 太小的：振荡太多，乘坐不舒适 -：理论上最快，但对参数变化敏感
  • 太大的：恢复太慢
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
def car_suspension():
    # 模拟汽车遇到凸起后的响应
    t = np.linspace(0, 3, 500)
    omega0 = 10  # rad/s
    
    fig, ax = plt.subplots(figsize=(12, 6))
    
    zeta_values = [0.3, 0.5, 0.7, 1.0, 1.5]
    styles = [':', '--', '-', '-.', ':']
    
    for zeta, style in zip(zeta_values, styles):
        if zeta < 1:
            omega_d = omega0 * np.sqrt(1 - zeta**2)
            x = np.exp(-zeta * omega0 * t) * (np.cos(omega_d * t) + 
                                              zeta/np.sqrt(1-zeta**2) * np.sin(omega_d * t))
        elif zeta == 1:
            x = (1 + omega0 * t) * np.exp(-omega0 * t)
        else:
            r1 = -zeta * omega0 + omega0 * np.sqrt(zeta**2 - 1)
            r2 = -zeta * omega0 - omega0 * np.sqrt(zeta**2 - 1)
            A = (r2 + omega0) / (r2 - r1)
            B = -(r1 + omega0) / (r2 - r1)
            x = A * np.exp(r1 * t) + B * np.exp(r2 * t)
        
        ax.plot(t, x, linestyle=style, linewidth=2.5, label=f'ζ = {zeta}')
    
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.axhline(y=0.02, color='gray', linestyle='--', alpha=0.5)
    ax.axhline(y=-0.02, color='gray', linestyle='--', alpha=0.5)
    ax.fill_between(t, -0.02, 0.02, alpha=0.1, color='green', label='2% Settling Band')
    
    ax.set_xlabel('Time (s)', fontsize=12)
    ax.set_ylabel('Displacement (normalized)', fontsize=12)
    ax.set_title('Car Suspension Response to Bump', fontsize=14, fontweight='bold')
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 3)
    ax.set_ylim(-0.5, 1.2)
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig4_suspension.png', dpi=150, bbox_inches='tight')
    plt.close()

car_suspension()

非齐次方程：待定系数法

方法概述

特征根与解的关系

对于非齐次方程：

步骤： 1. 解齐次方程得$y_h根据f(x)的形式猜特解y_p通解y = y_h + y_p$

猜解规则

的形式 猜测的
或
（多项式）

重要修正：如果猜测的 的某一项已经是 的解，则需要乘以（或 等）来修正。

例题详解

例 4：求解 Step 1：解齐次方程 特征方程： Step 2：猜特解，按规则猜 但 已经是 的一部分！需要修正为 代入原方程： 代入： 所以 Step 3：通解

共振现象

例 5：求解 齐次解：（因为）

注意： 已经是 的一部分！

猜 经过计算（略去细节）得：

关键观察：特解含有 因子，意味着振幅随时间线性增长！

这就是共振（ resonance）现象——当外力频率等于自然频率时，振幅会无限增大（在无阻尼情况下）。

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
def resonance_demo():
    t = np.linspace(0, 20, 1000)
    omega0 = 2  # 自然频率
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 非共振情况： omega != omega0
    omega = 2.5  # 驱动频率
    # y'' + 4y = cos(2.5t)
    # 特解： y_p = A*cos(2.5t), A = 1/(4 - 2.5^2) = 1/(4-6.25) = -1/2.25
    A = 1 / (4 - omega**2)
    y_nonres = A * np.cos(omega * t)
    
    axes[0, 0].plot(t, y_nonres, 'b-', linewidth=2)
    axes[0, 0].set_xlabel('Time', fontsize=12)
    axes[0, 0].set_ylabel('Displacement', fontsize=12)
    axes[0, 0].set_title(f'Non-Resonance: ω = {omega} ≠ ω₀ = 2', fontsize=12, fontweight='bold')
    axes[0, 0].grid(True, alpha=0.3)
    
    # 共振情况： omega = omega0
    # y'' + 4y = cos(2t)
    # 特解： y_p = (t/4)*sin(2t)
    y_res = (t / 4) * np.sin(2 * t)
    
    axes[0, 1].plot(t, y_res, 'r-', linewidth=2)
    axes[0, 1].plot(t, t/4, 'k--', linewidth=1.5, alpha=0.5, label='Envelope')
    axes[0, 1].plot(t, -t/4, 'k--', linewidth=1.5, alpha=0.5)
    axes[0, 1].set_xlabel('Time', fontsize=12)
    axes[0, 1].set_ylabel('Displacement', fontsize=12)
    axes[0, 1].set_title('Resonance: ω = ω₀ = 2 (Unbounded Growth!)', fontsize=12, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3)
    axes[0, 1].legend()
    
    # 拍频现象： omega ≈ omega0
    omega = 2.1
    A = 1 / (4 - omega**2)
    # 通解加上齐次部分，初始静止
    # 完整解： y = A(cos(ω t) - cos(2t))
    y_beat = A * (np.cos(omega * t) - np.cos(2 * t))
    # 拍频公式： 2A*sin((ω-2)t/2)*sin((ω+2)t/2)
    envelope = 2 * np.abs(A) * np.abs(np.sin((omega - 2) * t / 2))
    
    axes[1, 0].plot(t, y_beat, 'g-', linewidth=2)
    axes[1, 0].plot(t, envelope, 'k--', linewidth=1.5, alpha=0.5, label='Envelope')
    axes[1, 0].plot(t, -envelope, 'k--', linewidth=1.5, alpha=0.5)
    axes[1, 0].set_xlabel('Time', fontsize=12)
    axes[1, 0].set_ylabel('Displacement', fontsize=12)
    axes[1, 0].set_title(f'Beating: ω = {omega} ≈ ω₀ = 2', fontsize=12, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].legend()
    
    # 频率响应曲线
    omega_range = np.linspace(0.1, 4, 200)
    omega0 = 2
    zeta = 0.1  # 小阻尼
    
    # 稳态振幅：|H(ω)| = 1/sqrt((omega0^2 - omega^2)^2 + (2*zeta*omega0*omega)^2)
    H = 1 / np.sqrt((omega0**2 - omega_range**2)**2 + (2*zeta*omega0*omega_range)**2)
    
    axes[1, 1].plot(omega_range, H, 'b-', linewidth=2.5)
    axes[1, 1].axvline(x=omega0, color='r', linestyle='--', label=f'ω₀ = {omega0}')
    axes[1, 1].set_xlabel('Driving Frequency ω', fontsize=12)
    axes[1, 1].set_ylabel('Amplitude |H(ω)|', fontsize=12)
    axes[1, 1].set_title(f'Frequency Response (ζ = {zeta})', fontsize=12, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig5_resonance.png', dpi=150, bbox_inches='tight')
    plt.close()

resonance_demo()

参数变动法（ Variation of Parameters）

为什么需要？

待定系数法只能处理特定形式的（指数、三角、多项式）。对于更一般的情况，需要参数变动法。

二阶方程的公式

对于 如果已知齐次解，则特解为：

其中 是 Wronskian 。

例题

例 6：求解 齐次解：

RLC 电路：高阶 ODE 的经典应用

电路方程

考虑串联 RLC 电路，接入电压源。

基尔霍夫电压定律：

对时间求导（设 可导）：

这是关于电流 的二阶 ODE！

或者用电容电压：

特征方程分析

齐次方程的特征方程：

定义： -$_0 = ：谐振频率= ：衰减系数Q = $：品质因数

根的判别： -：欠阻尼（振荡） -：临界阻尼 -：过阻尼

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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def rlc_circuit():
    # RLC 电路参数
    L = 0.1   # 电感 100mH
    C = 1e-6  # 电容 1 μ F
    
    omega0 = 1 / np.sqrt(L * C)
    print(f"Resonant frequency: {omega0/(2*np.pi):.1f} Hz")
    
    t = np.linspace(0, 0.02, 1000)  # 20ms
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 不同电阻值
    R_values = [10, 100, 632, 1000, 2000]  # 临界阻尼: R = 2*sqrt(L/C) ≈ 632 Ω
    R_critical = 2 * np.sqrt(L / C)
    print(f"Critical resistance: {R_critical:.1f} Ω")
    
    for R in R_values:
        alpha = R / (2 * L)
        discriminant = alpha**2 - omega0**2
        
        if discriminant < 0:  # 欠阻尼
            omega_d = np.sqrt(-discriminant)
            # 初始条件： Vc(0)=0, i(0)=0, 接入阶跃电压
            # 响应： Vc = V0*(1 - exp(-α t)*(cos(ω d*t) + α/ω d*sin(ω d*t)))
            Vc = 1 - np.exp(-alpha * t) * (np.cos(omega_d * t) + alpha/omega_d * np.sin(omega_d * t))
            label = f'R = {R}Ω (underdamped)'
        elif abs(discriminant) < 1e-6:  # 临界阻尼
            Vc = 1 - (1 + alpha * t) * np.exp(-alpha * t)
            label = f'R = {R}Ω (critical)'
        else:  # 过阻尼
            s1 = -alpha + np.sqrt(discriminant)
            s2 = -alpha - np.sqrt(discriminant)
            A1 = s2 / (s2 - s1)
            A2 = -s1 / (s2 - s1)
            Vc = 1 - A1 * np.exp(s1 * t) - A2 * np.exp(s2 * t)
            label = f'R = {R}Ω (overdamped)'
        
        axes[0, 0].plot(t * 1000, Vc, linewidth=2, label=label)
    
    axes[0, 0].axhline(y=1, color='k', linestyle='--', alpha=0.5)
    axes[0, 0].set_xlabel('Time (ms)', fontsize=12)
    axes[0, 0].set_ylabel('Capacitor Voltage (V)', fontsize=12)
    axes[0, 0].set_title('RLC Step Response: Different Damping', fontsize=12, fontweight='bold')
    axes[0, 0].legend(fontsize=9)
    axes[0, 0].grid(True, alpha=0.3)
    
    # 频率响应
    R = 100  # 选择一个欠阻尼情况
    omega = np.linspace(0.01 * omega0, 3 * omega0, 500)
    
    # 传递函数幅值 |H(j ω)| = 1/|1 - (ω/ω 0)^2 + j*2 ζ*(ω/ω 0)|
    zeta = R / (2 * np.sqrt(L / C))
    omega_ratio = omega / omega0
    H_mag = 1 / np.sqrt((1 - omega_ratio**2)**2 + (2*zeta*omega_ratio)**2)
    H_phase = -np.arctan2(2*zeta*omega_ratio, 1 - omega_ratio**2) * 180 / np.pi
    
    axes[0, 1].semilogy(omega / omega0, H_mag, 'b-', linewidth=2.5)
    axes[0, 1].axvline(x=1, color='r', linestyle='--', label='Resonance')
    axes[0, 1].set_xlabel('ω/ω₀', fontsize=12)
    axes[0, 1].set_ylabel('|H(j ω)|', fontsize=12)
    axes[0, 1].set_title(f'Frequency Response Magnitude (R={R}Ω)', fontsize=12, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3, which='both')
    axes[0, 1].legend()
    
    # 相位响应
    axes[1, 0].plot(omega / omega0, H_phase, 'g-', linewidth=2.5)
    axes[1, 0].axvline(x=1, color='r', linestyle='--', label='Resonance (-90 °)')
    axes[1, 0].axhline(y=-90, color='gray', linestyle=':', alpha=0.5)
    axes[1, 0].set_xlabel('ω/ω₀', fontsize=12)
    axes[1, 0].set_ylabel('Phase (degrees)', fontsize=12)
    axes[1, 0].set_title('Phase Response', fontsize=12, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].legend()
    
    # Q 因子对共振峰的影响
    Q_values = [1, 5, 10, 50]
    omega_ratio = np.linspace(0.5, 1.5, 500)
    
    for Q in Q_values:
        H_mag = 1 / np.sqrt((1 - omega_ratio**2)**2 + (omega_ratio/Q)**2)
        axes[1, 1].plot(omega_ratio, H_mag, linewidth=2, label=f'Q = {Q}')
    
    axes[1, 1].set_xlabel('ω/ω₀', fontsize=12)
    axes[1, 1].set_ylabel('|H(j ω)|', fontsize=12)
    axes[1, 1].set_title('Effect of Quality Factor Q', fontsize=12, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig6_rlc.png', dpi=150, bbox_inches='tight')
    plt.close()

rlc_circuit()

高阶方程的系统化处理

转化为一阶系统

任何阶 ODE 都可以转化为 个一阶 ODE 的系统。

对于 令：

得到：

矩阵形式：

例：二阶方程的系统形式

令

矩阵 的特征值就是原方程的特征根！

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
def system_form():
    from scipy.integrate import solve_ivp
    
    # 二阶方程: y'' + 2y' + 5y = 0
    # 系统形式
    def system(t, X):
        x1, x2 = X
        return [x2, -5*x1 - 2*x2]
    
    # 初始条件
    X0 = [1, 0]  # y(0) = 1, y'(0) = 0
    t_span = [0, 6]
    t_eval = np.linspace(0, 6, 500)
    
    sol = solve_ivp(system, t_span, X0, t_eval=t_eval, method='RK45')
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    # 时间响应
    axes[0].plot(sol.t, sol.y[0], 'b-', linewidth=2, label='y(t)')
    axes[0].plot(sol.t, sol.y[1], 'r--', linewidth=2, label="y'(t)")
    axes[0].set_xlabel('Time', fontsize=12)
    axes[0].set_ylabel('Value', fontsize=12)
    axes[0].set_title('System Response', fontsize=12, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 相空间
    axes[1].plot(sol.y[0], sol.y[1], 'g-', linewidth=2)
    axes[1].plot(sol.y[0][0], sol.y[1][0], 'ro', markersize=10, label='Start')
    axes[1].plot(0, 0, 'k*', markersize=12, label='Equilibrium')
    axes[1].set_xlabel('y', fontsize=12)
    axes[1].set_ylabel("y'", fontsize=12)
    axes[1].set_title('Phase Portrait', fontsize=12, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    # 向量场
    y_range = np.linspace(-1.5, 1.5, 15)
    yp_range = np.linspace(-3, 3, 15)
    Y, YP = np.meshgrid(y_range, yp_range)
    
    DY = YP
    DYP = -5*Y - 2*YP
    
    # 归一化
    M = np.sqrt(DY**2 + DYP**2)
    M[M == 0] = 1
    DY_norm = DY / M
    DYP_norm = DYP / M
    
    axes[2].quiver(Y, YP, DY_norm, DYP_norm, M, cmap='coolwarm', alpha=0.6)
    
    # 几条轨迹
    for y0, yp0 in [(1, 0), (-1, 0), (0, 2), (0, -2), (1, 2), (-1, -2)]:
        sol_i = solve_ivp(system, [0, 10], [y0, yp0], t_eval=np.linspace(0, 10, 300))
        axes[2].plot(sol_i.y[0], sol_i.y[1], 'k-', linewidth=1, alpha=0.7)
    
    axes[2].set_xlabel('y', fontsize=12)
    axes[2].set_ylabel("y'", fontsize=12)
    axes[2].set_title('Vector Field with Trajectories', fontsize=12, fontweight='bold')
    axes[2].set_xlim(-1.5, 1.5)
    axes[2].set_ylim(-3, 3)
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig7_system.png', dpi=150, bbox_inches='tight')
    plt.close()

system_form()

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
26
27
from sympy import symbols, Function, dsolve, Eq, exp, sin, cos, classify_ode

x = symbols('x')
y = Function('y')

# 例 1：常系数齐次
eq1 = Eq(y(x).diff(x, 2) - 3*y(x).diff(x) + 2*y(x), 0)
print("方程 1：", eq1)
print("类型：", classify_ode(eq1))
sol1 = dsolve(eq1)
print("解：", sol1)
print()

# 例 2：非齐次
eq2 = Eq(y(x).diff(x, 2) + y(x), sin(x))
print("方程 2：", eq2)
sol2 = dsolve(eq2)
print("解：", sol2)
print()

# 例 3：带初始条件
eq3 = Eq(y(x).diff(x, 2) + 2*y(x).diff(x) + 5*y(x), 0)
ics = {y(0): 1, y(x).diff(x).subs(x, 0): 0}
sol3 = dsolve(eq3, ics=ics)
print("方程 3（带 IVP）：", eq3)
print("初始条件： y(0)=1, y'(0)=0")
print("解：", sol3)

使用 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
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

def forced_oscillator(t, X, omega, zeta, omega_f, F0):
    """
    受迫振动： y'' + 2 ζω₀ y' + ω₀² y = F ₀ cos(ω_f t)
    """
    y, v = X
    dydt = v
    dvdt = F0 * np.cos(omega_f * t) - 2*zeta*omega*v - omega**2 * y
    return [dydt, dvdt]

# 参数
omega0 = 2 * np.pi  # 自然频率 1 Hz
zeta = 0.1  # 阻尼比
F0 = 1.0  # 外力幅值

# 不同驱动频率
omega_f_values = [0.5*omega0, omega0, 1.5*omega0]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

for ax, omega_f in zip(axes, omega_f_values):
    sol = solve_ivp(
        forced_oscillator,
        [0, 20],
        [0, 0],
        args=(omega0, zeta, omega_f, F0),
        t_eval=np.linspace(0, 20, 1000),
        method='RK45'
    )
    
    ax.plot(sol.t, sol.y[0], 'b-', linewidth=2)
    ax.set_xlabel('Time (s)', fontsize=12)
    ax.set_ylabel('Displacement', fontsize=12)
    ax.set_title(f'ω_f/ω₀ = {omega_f/omega0:.1f}', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('第 3 章-高阶线性微分方程/fig8_forced.png', dpi=150, bbox_inches='tight')
plt.close()

实际应用：桥梁振动与塔科马海峡大桥

历史背景

1940 年 11 月 7 日，华盛顿州的塔科马海峡大桥在中等风速下发生剧烈扭转振动，最终坍塌。这是工程史上最著名的共振灾难之一。

简化模型

桥梁振动可以用四阶 ODE 描述：

对于特定模态，简化为：

当风的脉动频率接近桥梁的自然频率时，发生共振。

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
def tacoma_simulation():
    """简化模型：模拟塔科马大桥式的振动增长"""
    
    # 参数
    omega_n = 0.2 * 2 * np.pi  # 自然频率约 0.2Hz
    zeta = 0.005  # 非常小的阻尼
    omega_f = omega_n * 0.99  # 接近共振的驱动频率
    F0 = 0.1
    
    def bridge(t, X):
        q, v = X
        # 非线性气动力（简化）
        F = F0 * np.sin(omega_f * t) * (1 + 0.1 * q)  # 位移相关的非线性
        return [v, F - 2*zeta*omega_{n*}v - omega_{n*}*2 * q]
    
    t_span = [0, 600]  # 10 分钟
    sol = solve_ivp(bridge, t_span, [0, 0], 
                   t_eval=np.linspace(0, 600, 5000),
                   method='RK45')
    
    fig, axes = plt.subplots(2, 1, figsize=(14, 8))
    
    axes[0].plot(sol.t, sol.y[0], 'b-', linewidth=1)
    axes[0].set_xlabel('Time (s)', fontsize=12)
    axes[0].set_ylabel('Displacement', fontsize=12)
    axes[0].set_title('Tacoma-style Resonance (Simplified Model)', fontsize=14, fontweight='bold')
    axes[0].grid(True, alpha=0.3)
    
    # 包络线
    from scipy.signal import hilbert
    analytic_signal = hilbert(sol.y[0])
    amplitude_envelope = np.abs(analytic_signal)
    axes[0].plot(sol.t, amplitude_envelope, 'r-', linewidth=2, label='Envelope')
    axes[0].legend()
    
    # 频谱
    from scipy.fft import fft, fftfreq
    
    N = len(sol.t)
    dt = sol.t[1] - sol.t[0]
    yf = fft(sol.y[0])
    xf = fftfreq(N, dt)
    
    positive_freq = xf > 0
    axes[1].plot(xf[positive_freq], 2/N * np.abs(yf[positive_freq]), 'g-', linewidth=2)
    axes[1].axvline(x=omega_n/(2*np.pi), color='r', linestyle='--', label=f'Natural freq = {omega_n/(2*np.pi):.3f} Hz')
    axes[1].set_xlabel('Frequency (Hz)', fontsize=12)
    axes[1].set_ylabel('Amplitude', fontsize=12)
    axes[1].set_title('Frequency Spectrum', fontsize=12, fontweight='bold')
    axes[1].set_xlim(0, 0.5)
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 3 章-高阶线性微分方程/fig9_tacoma.png', dpi=150, bbox_inches='tight')
    plt.close()

tacoma_simulation()

总结

本章我们深入学习了高阶线性微分方程的理论和应用：

核心方法

方程类型 求解方法
常系数齐次 特征方程
常系数非齐次 待定系数法 + 特征方程
变系数非齐次 参数变动法
复杂方程 转化为系统 + 数值方法

特征根与解的关系

特征根类型 对应解
不同实根
重实根
复根

物理直觉

  • 阻尼比：决定振荡特性 -：欠阻尼（振荡衰减） -：临界阻尼（最快无振荡恢复） -：过阻尼（缓慢恢复）

  • 共振：驱动频率 ≈ 自然频率时，振幅急剧增大

工程应用

  • 机械振动与减震设计
  • 电路分析（ RLC 电路）
  • 结构工程（桥梁、建筑）
  • 控制系统设计

练习题

基础题

  1. 求解 的通解。

  2. 求解（提示：重根）。

  3. 求解，，。

  4. 求解，并画出解的图像。

  5. 用待定系数法求 的特解。

进阶题

  1. 求解，讨论共振现象。

  2. 用参数变动法求 的通解。

  3. RLC 电路中，，：

    • 计算自然频率和阻尼比
    • 判断是欠阻尼还是过阻尼
    • 画出阶跃响应曲线

  4. 证明：对于，当 时，解穿越零点的周期是。

  5. 弹簧-质量-阻尼系统，求稳态响应的振幅和相位。

编程题

  1. 用 Python 绘制不同阻尼比下的单位阶跃响应曲线。

  2. 实现参数变动法的一般算法，测试几个例子。

  3. 模拟双摆系统（这是一个混沌系统），观察对初始条件的敏感性。

  4. 设计一个简单的 PID 控制器，分析其对二阶系统的稳定性影响。

参考资料

  1. Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations. Wiley.
  2. Kreyszig, E. (2011). Advanced Engineering Mathematics. Wiley.
  3. Nagle, R. K., Saff, E. B., & Snider, A. D. (2017). Fundamentals of Differential Equations. Pearson.
  4. MIT OCW 18.03 - Differential Equations (Video Lectures)

下一章预告：《拉普拉斯变换》 — 将微分方程变成代数方程的魔法

  • 本文标题：常微分方程（三）高阶线性微分方程
  • 本文作者：Chen Kai
  • 创建时间：2019-04-14 10:15:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E4%B8%89%EF%BC%89%E9%AB%98%E9%98%B6%E7%BA%BF%E6%80%A7%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论