常微分方程（六）线性微分方程组
Chen Kai BOSS
  2019-05-01 11:00:00 2019-05-01 11:00      5.8k 字  28 分钟  

当多个变量相互影响时，单个微分方程就不够用了。捕食者与猎物的数量、电路中的电流与电压、化学反应中的各物质浓度——它们之间的关系需要用微分方程组来描述。本章我们将看到，线性代数与微分方程的结合产生了极其优美的理论：矩阵指数、基解矩阵、相空间分析，这些工具让复杂系统的行为变得清晰可见。

从一个生态学问题说起

想象一个简化的生态系统：兔子和狼。兔子吃草繁殖，狼吃兔子繁殖。设 为兔子数量， 为狼的数量。

直觉建模： - 兔子的增长率取决于自身数量（越多繁殖越快）减去被狼吃掉的数量 - 狼的增长率取决于有多少兔子可吃减去自然死亡

一个简化的线性模型： $

这就是一个线性微分方程组。如何求解它？如何理解 和 的长期行为？

向量形式与矩阵记号

向量化表示

将方程组写成向量形式：

简记为：

其中 是状态向量， 是系数矩阵。

一般形式

齐次线性系统：

非齐次线性系统：

本章主要关注常系数情况： 不依赖于。

回顾：一维情况

对于一维方程$ x' = ax，解是 x(t) = x_0 e^{at}$。

自然猜想：对于，解是否为？

是的！但首先需要定义矩阵指数。

矩阵指数

定义

回忆标量指数函数的泰勒展开： $$

e^x = 1 + x + + + = _{k=0}^{} $$

矩阵指数的定义： $$

e^A = I + A + + + = _{k=0}^{} $$

这个级数对任何方阵都收敛！

基本性质

性质 1：（零矩阵的指数是单位矩阵）

性质 2： 这正是我们需要的微分性质！

性质 3：如果（ 和 可交换），则 警告：如果，则一般！这是矩阵与标量的重要区别。

性质 4： 总是可逆的， 性质 5：

线性系统的解

定理：初值问题

的唯一解是：

验证： - ✓ - ✓

如何计算矩阵指数？

直接用定义计算很麻烦。有几种更实用的方法。

方法一：对角化

如果 可对角化：，其中。

则： $$

A^k = PDkP{-1} $

$

而对角矩阵的指数很简单： $$

e^{Dt} = (e^{_1 t}, , e^{_n t}) $$

方法二：特征值与特征向量

设 有特征值 和对应的特征向量。

如果特征向量线性无关，则通解为：

为什么？ 因为，所以 是方程 的解：

例子： 2 × 2 矩阵

考虑开头的生态模型： $$

A =

$$

步骤 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
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
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.linalg import expm

# 定义系统
A = np.array([[2, -1], 
              [1, 0.5]])

# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)
print("特征值:", eigenvalues)
print("特征向量:\n", eigenvectors)

# 定义微分方程
def system(x, t, A):
    return A @ x

# 求解多个初值
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

t = np.linspace(0, 5, 500)

# 左图：相空间轨迹
ax1 = axes[0]
initial_conditions = [
    [1, 0], [0, 1], [-1, 0], [0, -1],
    [1, 1], [-1, 1], [1, -1], [-1, -1],
    [2, 0], [0, 2]
]

for x0 in initial_conditions:
    sol = odeint(system, x0, t, args=(A,))
    ax1.plot(sol[:, 0], sol[:, 1], linewidth=1.5, alpha=0.7)
    ax1.scatter(x0[0], x0[1], s=50, zorder=3)

# 画特征向量方向
for i, (ev, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
    if np.isreal(ev):
        scale = 2
        ax1.arrow(0, 0, scale*vec[0].real, scale*vec[1].real, 
                 head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
                 linewidth=2, label=f'$\\lambda={ev:.2f}$')

ax1.set_xlabel('x (rabbits)', fontsize=12)
ax1.set_ylabel('y (wolves)', fontsize=12)
ax1.set_title('Phase Portrait: Unstable Spiral', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-3, 3)
ax1.set_ylim(-3, 3)
ax1.set_aspect('equal')

# 右图：时间演化
ax2 = axes[1]
x0 = [1, 0.5]
sol = odeint(system, x0, t, args=(A,))

ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='x(t) - Rabbits')
ax2.plot(t, sol[:, 1], 'r-', linewidth=2, label='y(t) - Wolves')
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Time Evolution: Oscillating Growth', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('linear_system_spiral.png', dpi=150, bbox_inches='tight')
plt.show()

# 验证矩阵指数
t_test = 1.0
x0 = np.array([1, 0.5])

# 方法 1：用 scipy 计算
exp_At = expm(A * t_test)
x_expm = exp_At @ x0

# 方法 2：用 ODE 求解器
sol_ode = odeint(system, x0, [0, t_test], args=(A,))[-1]

print(f"\n 在 t={t_test}时的解:")
print(f"  矩阵指数方法: {x_expm}")
print(f"  ODE 求解器:   {sol_ode}")
print(f"  误差: {np.linalg.norm(x_expm - sol_ode):.2e}")

相平面分析： 2 × 2 系统的分类

特征值决定一切

对于 2 × 2 实矩阵，根据特征值 的性质，相图有以下几种类型：

情况 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
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
95
96
97
98
99
100
101
102
103
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def plot_phase_portrait(A, ax, title, xlim=(-3, 3), ylim=(-3, 3)):
    """绘制相图"""
    def system(x, t):
        return A @ x
    
    # 网格点
    x_range = np.linspace(xlim[0], xlim[1], 20)
    y_range = np.linspace(ylim[0], ylim[1], 20)
    X, Y = np.meshgrid(x_range, y_range)
    
    # 计算方向场
    U = A[0, 0] * X + A[0, 1] * Y
    V = A[1, 0] * X + A[1, 1] * Y
    
    # 归一化
    M = np.sqrt(U**2 + V**2)
    M[M == 0] = 1
    U_norm = U / M
    V_norm = V / M
    
    # 绘制方向场
    ax.quiver(X, Y, U_norm, V_norm, M, cmap='coolwarm', alpha=0.5)
    
    # 绘制轨迹
    t = np.linspace(0, 10, 500)
    t_back = np.linspace(0, -10, 500)
    
    # 选择初值
    angles = np.linspace(0, 2*np.pi, 8, endpoint=False)
    r = 2.5
    for angle in angles:
        x0 = [r * np.cos(angle), r * np.sin(angle)]
        
        # 正向积分
        sol = odeint(system, x0, t)
        ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=1, alpha=0.7)
        
        # 反向积分
        sol_back = odeint(system, x0, t_back)
        ax.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=1, alpha=0.7)
    
    # 画特征向量
    eigenvalues, eigenvectors = np.linalg.eig(A)
    for i, (ev, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
        if np.isreal(ev) and np.isreal(vec).all():
            vec = vec.real
            scale = 2.5
            ax.arrow(0, 0, scale*vec[0], scale*vec[1], 
                    head_width=0.15, head_length=0.1, fc='red', ec='red',
                    linewidth=2)
            ax.arrow(0, 0, -scale*vec[0], -scale*vec[1], 
                    head_width=0.15, head_length=0.1, fc='red', ec='red',
                    linewidth=2)
    
    ax.plot(0, 0, 'ko', markersize=8)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('y', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    # 显示特征值
    ev_str = ', '.join([f'{ev:.2f}' for ev in eigenvalues])
    ax.text(0.05, 0.95, f'λ = {ev_str}', transform=ax.transAxes, 
            fontsize=9, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# 创建 6 种典型情况
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 1. 稳定节点 (λ 1, λ 2 < 0)
A1 = np.array([[-2, 0], [0, -1]])
plot_phase_portrait(A1, axes[0, 0], 'Stable Node')

# 2. 不稳定节点 (λ 1, λ 2 > 0)
A2 = np.array([[2, 0], [0, 1]])
plot_phase_portrait(A2, axes[0, 1], 'Unstable Node')

# 3. 鞍点 (λ 1 < 0 < λ 2)
A3 = np.array([[-1, 0], [0, 2]])
plot_phase_portrait(A3, axes[0, 2], 'Saddle Point')

# 4. 稳定焦点 (α < 0)
A4 = np.array([[-0.5, 1], [-1, -0.5]])
plot_phase_portrait(A4, axes[1, 0], 'Stable Spiral')

# 5. 不稳定焦点 (α > 0)
A5 = np.array([[0.5, 1], [-1, 0.5]])
plot_phase_portrait(A5, axes[1, 1], 'Unstable Spiral')

# 6. 中心 (α = 0)
A6 = np.array([[0, 1], [-1, 0]])
plot_phase_portrait(A6, axes[1, 2], 'Center')

plt.tight_layout()
plt.savefig('phase_portraits.png', dpi=150, bbox_inches='tight')
plt.show()

基解矩阵

定义

基解矩阵 是一个 矩阵，其列是 个线性无关的解：

性质： 1. 2. 对所有 成立 3. 通解为，其中 是常向量

与矩阵指数的关系

定理：如果，则。

更一般地，。

刘维尔公式

定理（ Liouville/Abel）：

对于常系数情况：

物理意义：相空间中体积元的演化率等于。

如果，体积收缩（耗散系统）； 如果，体积守恒（哈密顿系统）； 如果，体积膨胀。

重特征值与广义特征向量

问题

当特征值有重根时，可能没有足够的特征向量。例如： $$

A =

$$

特征值（重根），但只有一个特征向量。

广义特征向量

定义：如果 但，则 是 的阶广义特征向量。

对上面的例子： 取，则。

对应的解

对于$ k阶广义特征向量，对应的解包含多项式因子：$

对于上面的例子，两个独立解是： 

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

# 重特征值的例子
A = np.array([[3, 1], [0, 3]])

print("矩阵 A:")
print(A)
print("\n 特征值:", np.linalg.eigvals(A))

# 相图
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

def system(x, t):
    return A @ x

# 轨迹
ax1 = axes[0]
t = np.linspace(0, 1, 200)
t_back = np.linspace(0, -1, 200)

initial_conditions = [
    [1, 0.5], [-1, 0.5], [1, -0.5], [-1, -0.5],
    [0.1, 1], [0.1, -1], [-0.1, 1], [-0.1, -1]
]

for x0 in initial_conditions:
    sol = odeint(system, x0, t)
    sol_back = odeint(system, x0, t_back)
    ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=1.5, alpha=0.7)
    ax1.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=1.5, alpha=0.7)
    ax1.scatter(x0[0], x0[1], s=50, zorder=3)

# 特征向量方向
ax1.arrow(0, 0, 1, 0, head_width=0.1, head_length=0.05, fc='red', ec='red', linewidth=2)
ax1.arrow(0, 0, -1, 0, head_width=0.1, head_length=0.05, fc='red', ec='red', linewidth=2)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Degenerate Node (Repeated Eigenvalue λ=3)', fontsize=14, fontweight='bold')
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# 时间演化
ax2 = axes[1]
t_long = np.linspace(-1, 1, 500)
x0 = [0.1, 1]
sol_long = odeint(system, x0, t_long)

ax2.plot(t_long, sol_long[:, 0], 'b-', linewidth=2, label='x(t)')
ax2.plot(t_long, sol_long[:, 1], 'r-', linewidth=2, label='y(t)')

# 解析解对比
t_pos = t_long[t_long >= 0]
# x(t) = (c1 + c2*t) * e^(3t), y(t) = c2 * e^(3t)
# 初值 x(0)=0.1, y(0)=1 => c1=0.1, c2=1
c1, c2 = 0.1, 1
x_analytical = (c1 + c2*t_long) * np.exp(3*t_long)
y_analytical = c2 * np.exp(3*t_long)

ax2.plot(t_long, x_analytical, 'b--', linewidth=1.5, alpha=0.7, label='x(t) analytical')
ax2.plot(t_long, y_analytical, 'r--', linewidth=1.5, alpha=0.7, label='y(t) analytical')

ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Value', fontsize=12)
ax2.set_title('Time Evolution with Polynomial Growth', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-5, 10)

plt.tight_layout()
plt.savefig('repeated_eigenvalue.png', dpi=150, bbox_inches='tight')
plt.show()

非齐次系统：常数变易法

问题

求解非齐次系统：

常数变易法

思想：设解的形式为，其中 是基解矩阵。

代入方程：

由于： $

$$

最终解：

对于常系数情况（）：

这就是著名的Duhamel 公式。

例子：受迫谐振子

考虑受周期力驱动的弹簧系统： $$

x'' + _0^2 x = F_0 (t) $$

写成一阶系统（令$ y = x'）：$ 

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

# 受迫谐振子
omega0 = 1.0  # 固有频率
F0 = 0.5  # 驱动力幅度
gamma = 0.1  # 阻尼系数

def forced_oscillator(X, t, omega_drive):
    x, v = X
    dxdt = v
    dvdt = -omega0**2 * x - 2*gamma*v + F0 * np.cos(omega_drive * t)
    return [dxdt, dvdt]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 100, 5000)
X0 = [1, 0]

# 不同驱动频率
driving_freqs = [0.5, 1.0, 1.5, 2.0]

# 稳态振幅 vs 驱动频率（理论）
omega_range = np.linspace(0.1, 2.5, 200)
amplitude_theory = F0 / np.sqrt((omega0**2 - omega_range**2)**2 + (2*gamma*omega_range)**2)

ax1 = axes[0, 0]
ax1.plot(omega_range, amplitude_theory, 'b-', linewidth=2)
ax1.axvline(x=omega0, color='r', linestyle='--', label=f'ω₀ = {omega0}')
ax1.set_xlabel('Driving Frequency ω', fontsize=12)
ax1.set_ylabel('Steady-State Amplitude', fontsize=12)
ax1.set_title('Resonance Curve', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 时间演化：共振情况
ax2 = axes[0, 1]
omega_resonance = omega0
sol_res = odeint(forced_oscillator, X0, t, args=(omega_resonance,))
ax2.plot(t, sol_res[:, 0], 'b-', linewidth=1)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Displacement x', fontsize=12)
ax2.set_title(f'At Resonance: ω = ω₀ = {omega0}', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 时间演化：远离共振
ax3 = axes[1, 0]
omega_off = 2.0
sol_off = odeint(forced_oscillator, X0, t, args=(omega_off,))
ax3.plot(t, sol_off[:, 0], 'g-', linewidth=1)
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Displacement x', fontsize=12)
ax3.set_title(f'Off Resonance: ω = {omega_off}', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)

# 相空间
ax4 = axes[1, 1]
# 稳态时的相空间轨迹
t_steady = t[t > 50]
sol_res_steady = sol_res[t > 50]
sol_off_steady = sol_off[t > 50]

ax4.plot(sol_res_steady[:, 0], sol_res_steady[:, 1], 'b-', linewidth=1, 
         label=f'ω = {omega_resonance} (resonance)', alpha=0.7)
ax4.plot(sol_off_steady[:, 0], sol_off_steady[:, 1], 'g-', linewidth=1, 
         label=f'ω = {omega_off} (off resonance)', alpha=0.7)
ax4.set_xlabel('x', fontsize=12)
ax4.set_ylabel('v', fontsize=12)
ax4.set_title('Phase Space (Steady State)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_aspect('equal')

plt.tight_layout()
plt.savefig('forced_oscillator.png', dpi=150, bbox_inches='tight')
plt.show()

应用：耦合振子

问题设定

两个质量通过弹簧连接： 

1
2
墙 ----/\/\/\---- [m1] ----/\/\/\---- [m2] ----/\/\/\---- 墙
       k1              k2              k3

设 为两个质量相对于平衡位置的位移。

运动方程： $$

m_1 x_1'' = -k_1 x_1 + k_2(x_2 - x_1) $

$

简化（令$ m_1 = m_2 = 1， k_1 = k_2 = k_3 = 1）：$

x_1'' = -2x_1 + x_2 $

$

矩阵形式

令：

正则模（ Normal Modes）

特征值分析揭示系统的正则模——整个系统以特定频率同步振动的模式。

对于对称耦合系统，正则模有简单的物理解释： - 对称模：两个质量同方向运动（$ x_1 = x_2）反对称模：两个质量反方向运动（ x_1 = -x_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
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
95
96
97
98
99
100
101
102
103
104
105
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.linalg import eig

# 耦合振子系统
# 令 m1 = m2 = 1, k1 = k2 = k3 = 1

def coupled_oscillators(X, t):
    x1, v1, x2, v2 = X
    dx1 = v1
    dv1 = -2*x1 + x2
    dx2 = v2
    dv2 = x1 - 2*x2
    return [dx1, dv1, dx2, dv2]

# 系统矩阵
A = np.array([
    [0, 1, 0, 0],
    [-2, 0, 1, 0],
    [0, 0, 0, 1],
    [1, 0, -2, 0]
])

# 特征值分析
eigenvalues, eigenvectors = eig(A)
print("特征值:", eigenvalues)
print("\n 特征频率（虚部的绝对值）:")
freqs = np.abs(eigenvalues.imag)
unique_freqs = np.unique(np.round(freqs, 6))
print(unique_freqs[unique_freqs > 0])

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 30, 1000)

# 正则模 1：对称模 (x1(0) = x2(0) = 1)
ax1 = axes[0, 0]
X0_sym = [1, 0, 1, 0]
sol_sym = odeint(coupled_oscillators, X0_sym, t)
ax1.plot(t, sol_sym[:, 0], 'b-', linewidth=2, label='$ x_1(t)$')
ax1.plot(t, sol_sym[:, 2], 'r--', linewidth=2, label='$ x_2(t)$')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Displacement', fontsize=12)
ax1.set_title('Symmetric Mode: $ x_1(0) = x_2(0) = 1$', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 正则模 2：反对称模 (x1(0) = 1, x2(0) = -1)
ax2 = axes[0, 1]
X0_antisym = [1, 0, -1, 0]
sol_antisym = odeint(coupled_oscillators, X0_antisym, t)
ax2.plot(t, sol_antisym[:, 0], 'b-', linewidth=2, label='$ x_1(t)$')
ax2.plot(t, sol_antisym[:, 2], 'r--', linewidth=2, label='$ x_2(t)$')
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Displacement', fontsize=12)
ax2.set_title('Antisymmetric Mode: $ x_1(0) = 1, x_2(0) = -1$', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# 一般初值：拍频现象
ax3 = axes[1, 0]
X0_general = [1, 0, 0, 0]  # 只激发第一个质量
sol_general = odeint(coupled_oscillators, X0_general, t)
ax3.plot(t, sol_general[:, 0], 'b-', linewidth=1.5, label='$ x_1(t)$')
ax3.plot(t, sol_general[:, 2], 'r-', linewidth=1.5, label='$ x_2(t)$')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Displacement', fontsize=12)
ax3.set_title('General Initial Condition: Beat Phenomenon', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# 能量交换
ax4 = axes[1, 1]
# 动能 + 势能
KE1 = 0.5 * sol_general[:, 1]**2
KE2 = 0.5 * sol_general[:, 3]**2
PE1 = 0.5 * sol_general[:, 0]**2  # 只考虑边界弹簧
PE2 = 0.5 * sol_general[:, 2]**2
PE_coupling = 0.5 * (sol_general[:, 2] - sol_general[:, 0])**2

E1 = KE1 + PE1 + 0.5*PE_coupling
E2 = KE2 + PE2 + 0.5*PE_coupling
E_total = KE1 + KE2 + PE1 + PE2 + PE_coupling

ax4.plot(t, E1, 'b-', linewidth=1.5, label='Energy in mass 1')
ax4.plot(t, E2, 'r-', linewidth=1.5, label='Energy in mass 2')
ax4.plot(t, E_total, 'k--', linewidth=1.5, label='Total energy')
ax4.set_xlabel('Time', fontsize=12)
ax4.set_ylabel('Energy', fontsize=12)
ax4.set_title('Energy Transfer Between Oscillators', fontsize=14, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('coupled_oscillators.png', dpi=150, bbox_inches='tight')
plt.show()

# 打印正则频率
omega1 = 1  # 对称模频率 sqrt(k1/m) = sqrt(1)
omega2 = np.sqrt(3)  # 反对称模频率 sqrt((k1+2k2)/m) = sqrt(3)
print(f"\n 正则模频率:")
print(f"  对称模: ω₁ = {omega1:.4f}")
print(f"  反对称模: ω₂ = {omega2:.4f}")
print(f"  拍频: Δω = {omega2 - omega1:.4f}")

应用：电路网络

RLC 电路

考虑一个简单的 RLC 串联电路：

$$

L + RI + = V(t) $$

其中。

写成一阶系统（状态变量：电荷 和电流）：

特征值与电路行为

设（共振频率），（阻尼系数）。

特征值： - 欠阻尼（）：复特征值，振荡衰减 - 临界阻尼（）：重实特征值，最快达到稳态 - 过阻尼（）：两个实特征值，单调衰减

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

# RLC 电路响应
L = 1.0  # 电感
C = 1.0  # 电容
omega0 = 1 / np.sqrt(L * C)

def rlc_circuit(X, t, R, V_func):
    Q, I = X
    dQdt = I
    dIdt = (V_func(t) - R*I - Q/C) / L
    return [dQdt, dIdt]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 20, 1000)
X0 = [1, 0]  # 初始电荷为 1，电流为 0
V_func = lambda t: 0  # 无外加电压，研究自由响应

# 不同阻尼情况
R_values = [0.5, 2.0, 4.0]  # 欠阻尼、临界阻尼、过阻尼
labels = ['Underdamped', 'Critically damped', 'Overdamped']
colors = ['b', 'g', 'r']

# 电荷响应
ax1 = axes[0, 0]
for R, label, color in zip(R_values, labels, colors):
    gamma = R / (2*L)
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax1.plot(t, sol[:, 0], color=color, linewidth=2, 
             label=f'{label} (R={R}, γ/ω₀={gamma/omega0:.2f})')

ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Charge Q', fontsize=12)
ax1.set_title('RLC Circuit: Charge Response', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 电流响应
ax2 = axes[0, 1]
for R, label, color in zip(R_values, labels, colors):
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax2.plot(t, sol[:, 1], color=color, linewidth=2, label=f'{label}')

ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Current I', fontsize=12)
ax2.set_title('RLC Circuit: Current Response', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 相空间
ax3 = axes[1, 0]
for R, label, color in zip(R_values, labels, colors):
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax3.plot(sol[:, 0], sol[:, 1], color=color, linewidth=2, label=f'{label}')
    ax3.scatter(X0[0], X0[1], s=100, color=color, zorder=3)

ax3.plot(0, 0, 'ko', markersize=8)
ax3.set_xlabel('Charge Q', fontsize=12)
ax3.set_ylabel('Current I', fontsize=12)
ax3.set_title('Phase Space', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_aspect('equal')

# 受迫响应：频率响应
ax4 = axes[1, 1]
omega_range = np.linspace(0.1, 3, 200)
R = 0.5  # 欠阻尼

# 稳态振幅（理论）
# |H(j ω)| = 1 / sqrt((1 - (ω/ω 0)²)² + (2 γω/ω 0 ²)²)
gamma = R / (2*L)
amplitude = 1 / np.sqrt((1 - (omega_range/omega0)**2)**2 + (2*gamma*omega_range/omega0**2)**2)

ax4.plot(omega_range, amplitude, 'b-', linewidth=2)
ax4.axvline(x=omega0, color='r', linestyle='--', label=f'ω₀ = {omega0}')
ax4.set_xlabel('Driving Frequency ω', fontsize=12)
ax4.set_ylabel('Amplitude Response', fontsize=12)
ax4.set_title('Frequency Response (R=0.5)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_circuit.png', dpi=150, bbox_inches='tight')
plt.show()

高维系统与 Jordan 标准形

Jordan 标准形

当矩阵不可对角化时，可以化为Jordan 标准形： $$

A = PJP^{-1} $$

其中 是 Jordan 矩阵，由 Jordan 块组成。

一个 的 Jordan 块： $$

J_k() =

$$

Jordan 块的指数

$$

e^{J_k()t} = e^{t}

$$

关键观察：即使所有特征值有负实部，多项式因子可能导致短期增长，但长期仍然衰减。

稳定性初步

零解的稳定性

考虑 的零解。

定理： 1. 如果 的所有特征值的实部，则零解渐近稳定 2. 如果 有某个特征值的实部，则零解不稳定 3. 如果所有特征值的实部，且实部 的特征值都是简单的（代数重数等于几何重数），则零解稳定但不渐近稳定

这为下一章的详细稳定性理论做铺垫。

数值方法

矩阵指数的数值计算

计算 有多种方法：

方法 1： Pad é逼近（ scipy.linalg.expm）

最稳定的通用方法。

方法 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
import numpy as np
from scipy.linalg import expm, eig
import time

def matrix_exp_series(A, n_terms=50):
    """用级数计算矩阵指数"""
    result = np.eye(len(A))
    A_power = np.eye(len(A))
    factorial = 1.0
    
    for k in range(1, n_terms):
        A_power = A_power @ A
        factorial *= k
        result = result + A_power / factorial
    
    return result

def matrix_exp_diag(A):
    """用对角化计算矩阵指数"""
    eigenvalues, P = eig(A)
    exp_D = np.diag(np.exp(eigenvalues))
    return P @ exp_D @ np.linalg.inv(P)

# 测试
A = np.array([[1, 2], [3, 4]], dtype=float)

print("矩阵 A:")
print(A)

# 三种方法
exp_A_scipy = expm(A)
exp_A_series = matrix_exp_series(A, 100)
exp_A_diag = matrix_exp_diag(A)

print("\nexp(A) - scipy.linalg.expm:")
print(exp_A_scipy)

print("\nexp(A) - 级数展开:")
print(exp_A_series)

print("\nexp(A) - 对角化:")
print(np.real(exp_A_diag))

print("\n 误差（级数 vs scipy）:", np.linalg.norm(exp_A_series - exp_A_scipy))
print("误差（对角化 vs scipy）:", np.linalg.norm(np.real(exp_A_diag) - exp_A_scipy))

# 性能比较
n = 100
A_large = np.random.randn(n, n)

start = time.time()
for _ in range(10):
    expm(A_large)
print(f"\n{n}x{n}矩阵， scipy.expm: {(time.time()-start)/10*1000:.2f} ms")

start = time.time()
for _ in range(10):
    matrix_exp_series(A_large, 30)
print(f"{n}x{n}矩阵，级数展开: {(time.time()-start)/10*1000:.2f} ms")

ODE 求解器

对于大型系统或长时间积分，直接用矩阵指数可能不是最好的选择。

scipy.integrate.odeint和solve_ivp提供了多种自适应步长方法：

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

# 比较不同的 ODE 求解器
A = np.array([[-1, 100], [0, -1]])  # 刚性系统

def system(t, x):
    return A @ x

x0 = [1, 1]
t_span = (0, 10)
t_eval = np.linspace(0, 10, 1000)

# 不同方法
methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()

for ax, method in zip(axes, methods):
    try:
        sol = solve_ivp(system, t_span, x0, method=method, t_eval=t_eval)
        ax.plot(sol.t, sol.y[0], 'b-', linewidth=1.5, label='x ₁(t)')
        ax.plot(sol.t, sol.y[1], 'r-', linewidth=1.5, label='x ₂(t)')
        ax.set_title(f'{method} (nfev={sol.nfev})', fontsize=12, fontweight='bold')
    except Exception as e:
        ax.text(0.5, 0.5, f'Error: {e}', transform=ax.transAxes, ha='center')
        ax.set_title(method, fontsize=12)
    
    ax.set_xlabel('Time', fontsize=10)
    ax.set_ylabel('Value', fontsize=10)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('ode_solvers.png', dpi=150, bbox_inches='tight')
plt.show()

总结

本章要点

  1. 线性系统：，解为$(t) = e{At}0矩阵指数：定义： e^A = {k=0}{} A^k/k!$ - 计算：对角化、 Jordan 标准形、 Pad é逼近

  2. 相空间分析：

    • 特征值决定平衡点类型
    • 节点、焦点、鞍点、中心

  3. 基解矩阵：列为独立解的矩阵

  4. 非齐次系统：常数变易法、 Duhamel 公式

  5. 应用：耦合振子、电路、生态模型

下一章预告

在《稳定性理论》中，我们将： - 深入研究 Lyapunov 稳定性 - 分析非线性系统的局部行为 - 学习 Lyapunov 函数方法 - 探索分岔与混沌的萌芽

练习题

基础题

  1. 计算矩阵 的指数。

  2. 求解初值问题： 3. 判断以下系统原点的类型（节点/焦点/鞍点/中心）：

    • - - $A = 证明(e^A) = e^{(A)}$。

  3. 对于（重特征值），求基解矩阵。

进阶题

  1. 耦合谐振子：

    • 求系统的正则模频率
    • 证明能量守恒
    • 分析拍频周期

  2. RLC 电路：

    • 证明临界阻尼条件 - 计算品质因数 的物理意义
    • 设计一个带宽为 10 Hz 、中心频率为 1000 Hz 的滤波器

  3. 证明如果 的所有特征值有负实部，则。

  4. 化学反应：考虑反应：

  • 写成矩阵形式
  • 证明总物质守恒
  • 求稳态浓度

编程题

  1. 实现函数matrix_exp(A, t)，用三种方法计算 并比较精度。

  2. 编写程序绘制 2 × 2 系统的完整相图，包括：

    • 方向场
    • 特征向量
    • 多条轨迹
    • 自动判断并标注平衡点类型

  3. 模拟三体问题的简化版：三个弹簧连接的质量在一条线上振动。

    • 找出三个正则模
    • 可视化能量在三个质量之间的流动

  4. 实现一个电路模拟器，可以分析任意 RLC 网络的响应。

思考题

  1. 为什么矩阵指数$ e^{A+B} e^A e^B（当AB BA$）？给出一个具体的反例。

  2. 相空间中的轨迹为什么不能相交（除了在平衡点）？这与解的唯一性有什么关系？

  3. 如何从物理直觉理解"迹为零意味着相空间体积守恒"？

参考资料

  1. Hirsch, M. W., Smale, S., & Devaney, R. L. (2012). Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press.
  2. Strang, G. (2019). Linear Algebra and Learning from Data. Wellesley-Cambridge Press.
  3. Perko, L. (2001). Differential Equations and Dynamical Systems. Springer.
  4. MIT 18.03SC: Differential Equations, Unit III: Linear Systems.

  5. Moler, C., & Van Loan, C. (2003). Nineteen dubious ways to compute the exponential of a matrix. SIAM Review.

本文是《常微分方程的世界》系列的第 6 章。

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