常微分方程（七）稳定性理论

一个系统受到微小扰动后会怎样？是回到原来的状态，还是偏离越来越远？这个问题的答案决定了桥梁是否会在风中倒塌、生态系统能否维持平衡、经济周期是否会失控。稳定性理论给了我们回答这些问题的数学工具。本章我们将从直觉出发，逐步建立 Lyapunov 稳定性理论，并看到它在工程、生物和物理中的广泛应用。

从一个倒立摆开始

想象一根竖直的木棍，你用手指托着它的底端，试图让它保持平衡。

两种平衡状态： 1. 木棍朝下（重心在支点下方）：轻轻推一下，它会摆动几下然后回到垂直位置——这是稳定的 2. 木棍朝上（重心在支点上方）：稍有偏差就会倒下——这是不稳定的

数学上，两种情况的平衡点都满足相同的条件（导数为零），但它们的本质完全不同。稳定性理论的任务就是区分这些情况。

稳定性的精确定义

Lyapunov 稳定性

考虑自治系统：

设 是平衡点，即。

定义（ Lyapunov 稳定）：平衡点 是稳定的，如果对任意，存在，使得当 时，对所有 有。

直觉：从平衡点附近出发，轨迹永远不会跑太远。

定义（渐近稳定）：如果 是稳定的，且存在 使得当 时，。

直觉：不仅不跑远，还会最终回到平衡点。

定义（不稳定）：如果 不是稳定的，则称它是不稳定的。

吸引域

定义：平衡点 的吸引域（ basin of attraction）是所有满足 的初值 的集合。

吸引域可以是整个空间（全局渐近稳定），也可以只是平衡点周围的一个区域（局部渐近稳定）。

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

# 可视化不同类型的稳定性
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 系统 1：稳定但不渐近稳定（简谐振子）
def harmonic(x, t):
    return [x[1], -x[0]]

ax1 = axes[0]
t = np.linspace(0, 20, 1000)
for r in [0.5, 1, 1.5, 2]:
    for theta0 in np.linspace(0, 2*np.pi, 4, endpoint=False):
        x0 = [r*np.cos(theta0), r*np.sin(theta0)]
        sol = odeint(harmonic, x0, t)
        ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.5, alpha=0.7)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Stable (not asymptotically)\nCenter', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# 系统 2：渐近稳定（阻尼振子）
def damped(x, t):
    return [x[1], -x[0] - 0.3*x[1]]

ax2 = axes[1]
for r in [0.5, 1, 1.5, 2]:
    for theta0 in np.linspace(0, 2*np.pi, 4, endpoint=False):
        x0 = [r*np.cos(theta0), r*np.sin(theta0)]
        sol = odeint(damped, x0, t)
        ax2.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)

ax2.plot(0, 0, 'ko', markersize=8)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('y', fontsize=12)
ax2.set_title('Asymptotically Stable\nStable Spiral', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')
ax2.grid(True, alpha=0.3)

# 系统 3：不稳定（鞍点）
def saddle(x, t):
    return [x[0], -x[1]]

ax3 = axes[2]
t_short = np.linspace(0, 2, 500)
t_back = np.linspace(0, -2, 500)

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

for x0 in initial_conditions:
    sol = odeint(saddle, x0, t_short)
    sol_back = odeint(saddle, x0, t_back)
    ax3.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)
    ax3.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=0.8, alpha=0.7)

# 稳定和不稳定流形
ax3.arrow(0, 0, 0, 1.5, head_width=0.1, head_length=0.1, fc='green', ec='green', linewidth=2)
ax3.arrow(0, 0, 0, -1.5, head_width=0.1, head_length=0.1, fc='green', ec='green', linewidth=2)
ax3.arrow(0, 0, 1.5, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', linewidth=2)
ax3.arrow(0, 0, -1.5, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', linewidth=2)

ax3.plot(0, 0, 'ko', markersize=8)
ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('y', fontsize=12)
ax3.set_title('Unstable\nSaddle Point', fontsize=12, fontweight='bold')
ax3.text(1.2, 0.2, 'Unstable', fontsize=10, color='red')
ax3.text(0.1, 1.2, 'Stable', fontsize=10, color='green')
ax3.set_aspect('equal')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-2, 2)
ax3.set_ylim(-2, 2)

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

线性化方法

基本思想

对于非线性系统，在平衡点附近进行 Taylor 展开：

由于，令：

其中 是Jacobi 矩阵。

Jacobi 矩阵

$$

A = D(^{*}) =

_{ = ^{*}} $$

Hartman-Grobman 定理

定理：如果 的所有特征值的实部都非零（双曲平衡点），则非线性系统在附近的定性行为与线性化系统 相同。

"定性行为相同"的精确含义：存在一个同胚（连续双射且逆也连续），将非线性系统的轨迹映射到线性系统的轨迹，保持时间方向。

重要推论： - 如果 的所有特征值有负实部→原点渐近稳定 - 如果 有某个特征值有正实部→原点不稳定

警告：如果 有纯虚特征值，线性化方法失效！需要高阶分析。

例子：阻尼单摆

单摆方程（无小角近似）：

写成一阶系统（, $ y = '$

x' = y $

$

平衡点： 且$ y = 0(n, 0) n = 0, , , $ Jacobi 矩阵： $$

A =

$$

在 处： $$

A =

$$

特征值： 如果，两个特征值都有负实部→渐近稳定（悬挂位置）

在 处： $$

A =

$$

特征值： 一正一负→鞍点→不稳定（倒立位置）

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

# 阻尼单摆
gamma = 0.3  # 阻尼系数
omega0 = 1.0  # 固有频率

def pendulum(X, t):
    x, y = X
    dxdt = y
    dydt = -gamma * y - omega0**2 * np.sin(x)
    return [dxdt, dydt]

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

# 相空间
ax1 = axes[0]
t = np.linspace(0, 50, 2000)

# 从不同初值出发的轨迹
for x0 in np.linspace(-3*np.pi, 3*np.pi, 15):
    for y0 in np.linspace(-3, 3, 5):
        X0 = [x0, y0]
        sol = odeint(pendulum, X0, t)
        ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.3, alpha=0.5)

# 标记平衡点
for n in range(-2, 3):
    if n % 2 == 0:  # 稳定平衡点
        ax1.plot(n*np.pi, 0, 'go', markersize=10, label='Stable' if n == 0 else '')
    else:  # 不稳定平衡点（鞍点）
        ax1.plot(n*np.pi, 0, 'rx', markersize=10, mew=2, label='Unstable' if n == 1 else '')

ax1.set_xlabel('θ (angle)', fontsize=12)
ax1.set_ylabel('ω (angular velocity)', fontsize=12)
ax1.set_title('Phase Portrait of Damped Pendulum', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-3*np.pi, 3*np.pi)
ax1.set_ylim(-4, 4)

# 时间演化
ax2 = axes[1]
initial_conditions = [
    [0.5, 0, 'Near stable eq'],
    [np.pi - 0.3, 0, 'Near unstable eq'],
    [0, 2, 'Large velocity']
]

for x0, y0, label in initial_conditions:
    sol = odeint(pendulum, [x0, y0], t)
    ax2.plot(t, sol[:, 0], linewidth=2, label=label)

ax2.axhline(y=0, color='g', linestyle='--', alpha=0.5)
ax2.axhline(y=np.pi, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('θ (angle)', fontsize=12)
ax2.set_title('Time Evolution', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

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

# 计算特征值
print("平衡点分析:")
for eq_name, x_eq in [("悬挂位置 (0, 0)", 0), ("倒立位置 (π, 0)", np.pi)]:
    A = np.array([[0, 1], [-omega0**2 * np.cos(x_eq), -gamma]])
    eigenvalues = np.linalg.eigvals(A)
    print(f"\n{eq_name}:")
    print(f"  Jacobi 矩阵:\n{A}")
    print(f"  特征值: {eigenvalues}")
    if all(ev.real < 0 for ev in eigenvalues):
        print(f"  结论: 渐近稳定")
    elif any(ev.real > 0 for ev in eigenvalues):
        print(f"  结论: 不稳定")

Lyapunov 直接法

能量的启发

物理系统的稳定性通常与能量有关： - 能量最小的状态是稳定的 - 如果能量持续减少，系统会趋向稳定状态

Lyapunov 直接法将这个思想数学化。

Lyapunov 函数

定义：对于系统，如果存在函数 满足： 1. $V(^) = 0$2. 对所有$ ^$3. 则 称为Lyapunov 函数。

稳定性定理

定理（ Lyapunov 稳定性定理）： 1. 如果存在 Lyapunov 函数且，则稳定 2. 如果 对所有，则渐近稳定 3. 如果存在函数 使得，则不稳定

直觉： 就像"广义能量"，如果它总是减少，系统最终会到达能量最低点（平衡点）。

如何构造 Lyapunov 函数？

这是一门艺术，没有万能公式。但有一些常用技巧：

技巧 1：物理能量

对于力学系统，总能量（动能+势能）常常是好的候选。

技巧 2：二次型

对于线性系统，尝试，其中 是正定矩阵。

如果能找到 使得（ 正定），则→渐近稳定。

这称为Lyapunov 方程。

技巧 3：试探法

从简单形式开始（如），检验 的符号。

例子： Van der Pol 振子的稳定性

Van der Pol 方程： $$

x'' - (1 - x^2)x' + x = 0 $$

这是描述电子管振荡器的方程，有非线性阻尼。

写成一阶系统： $$

x' = y $

$

平衡点： Jacobi 矩阵： $$

A =

$$

特征值： 当 时，→不稳定！

但这里有个有趣的现象：尽管原点不稳定，系统存在一个稳定极限环。

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

# Van der Pol 振子
def van_der_pol(X, t, mu):
    x, y = X
    dxdt = y
    dydt = mu * (1 - x**2) * y - x
    return [dxdt, dydt]

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

t = np.linspace(0, 50, 2000)

# 不同μ值的相图
mu_values = [0.5, 1.0, 3.0]

for ax, mu in zip(axes, mu_values):
    # 从内部出发
    for r in [0.1, 0.3, 0.5]:
        X0 = [r, 0]
        sol = odeint(van_der_pol, X0, t, args=(mu,))
        ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)
    
    # 从外部出发
    for r in [3, 4, 5]:
        X0 = [r, 0]
        sol = odeint(van_der_pol, X0, t, args=(mu,))
        ax.plot(sol[:, 0], sol[:, 1], 'r-', linewidth=0.8, alpha=0.7)
    
    ax.plot(0, 0, 'ko', markersize=8)
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(f'Van der Pol Oscillator (μ = {mu})', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')

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

# 时间演化展示极限环
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

mu = 1.0
t_long = np.linspace(0, 100, 5000)

# 从内部出发
ax1 = axes[0]
X0_in = [0.1, 0]
sol_in = odeint(van_der_pol, X0_in, t_long, args=(mu,))
ax1.plot(t_long, sol_in[:, 0], 'b-', linewidth=1)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('x(t)', fontsize=12)
ax1.set_title('Starting Inside Limit Cycle', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 从外部出发
ax2 = axes[1]
X0_out = [4, 0]
sol_out = odeint(van_der_pol, X0_out, t_long, args=(mu,))
ax2.plot(t_long, sol_out[:, 0], 'r-', linewidth=1)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('x(t)', fontsize=12)
ax2.set_title('Starting Outside Limit Cycle', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

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

LaSalle 不变原理

问题

当 但不严格小于零时，标准 Lyapunov 定理只能保证稳定，不能保证渐近稳定。

LaSalle 不变原理

定理：设 是 Lyapunov 函数，。定义集合。

设 是包含在 中的最大不变集（即：如果轨迹从 出发，它永远留在 中）。

则所有有界轨迹趋向。

应用：如果 只包含平衡点，则 渐近稳定。

例子

考虑带摩擦的单摆：

取能量函数： $$

V = ^2 + (1 - ) $$

计算： 。

在 上， 是常数，所以。代入方程：→。

对于，这确实是 中的不变集。所以 渐近稳定。

全局稳定性

全局渐近稳定

定义：如果平衡点 的吸引域是整个空间，则称它是全局渐近稳定的。

定理：如果： 1. 对所有2. 当（径向无界） 3. 对所有 则 全局渐近稳定。

例子： Logistic 方程

平衡点： 和。

对于，取。

当 且 时，。所以在 上， 是渐近稳定的。

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

# Logistic 方程
r = 1.0  # 增长率
K = 100  # 环境容纳量

def logistic(N, t):
    return r * N * (1 - N/K)

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

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

# 不同初值的解
ax1 = axes[0]
N0_values = [10, 30, 50, 80, 120, 150, 200]

for N0 in N0_values:
    sol = odeint(logistic, N0, t)
    ax1.plot(t, sol, linewidth=2, label=f'N ₀ = {N0}')

ax1.axhline(y=K, color='r', linestyle='--', linewidth=2, label=f'K = {K}')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Population N', fontsize=12)
ax1.set_title('Logistic Growth: Global Stability of K', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# Lyapunov 函数等高线
ax2 = axes[1]
N_range = np.linspace(1, 200, 200)
V = 0.5 * (N_range - K)**2

ax2.plot(N_range, V, 'b-', linewidth=2)
ax2.axvline(x=K, color='r', linestyle='--', linewidth=2, label=f'N = K = {K}')
ax2.set_xlabel('Population N', fontsize=12)
ax2.set_ylabel('V(N) = ½(N - K)²', fontsize=12)
ax2.set_title('Lyapunov Function', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# 在同一图上画 dV/dt
ax2_twin = ax2.twinx()
dV = -(r/K) * N_range * (N_range - K)**2
ax2_twin.plot(N_range, dV, 'g-', linewidth=2, alpha=0.7)
ax2_twin.set_ylabel('$\\dot{V}$', fontsize=12, color='g')
ax2_twin.tick_params(axis='y', labelcolor='g')

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

线性系统的 Lyapunov 分析

Lyapunov 方程

对于线性系统，设。

定理：如果 所有特征值有负实部，则对任意正定矩阵，存在唯一正定矩阵 满足： $$

A^T P + PA = -Q $$

这称为连续 Lyapunov 方程。

推论：如果能找到 使得 负定，则系统渐近稳定。

数值求解

import numpy as np
from scipy.linalg import solve_continuous_lyapunov
import matplotlib.pyplot as plt

# 线性系统的 Lyapunov 分析
A = np.array([[-2, 1], 
              [-1, -3]])

# 检查稳定性（特征值）
eigenvalues = np.linalg.eigvals(A)
print("特征值:", eigenvalues)
print("所有特征值实部 < 0:", all(ev.real < 0 for ev in eigenvalues))

# 求解 Lyapunov 方程 A'P + PA = -Q
Q = np.eye(2)  # 取单位矩阵
P = solve_continuous_lyapunov(A.T, -Q)

print("\nQ =")
print(Q)
print("\nP =")
print(P)

# 验证 P 是正定的
eigenvalues_P = np.linalg.eigvals(P)
print("\nP 的特征值:", eigenvalues_P)
print("P 正定:", all(ev > 0 for ev in eigenvalues_P))

# 可视化 Lyapunov 函数
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 等高线图
ax1 = axes[0]
x1 = np.linspace(-2, 2, 100)
x2 = np.linspace(-2, 2, 100)
X1, X2 = np.meshgrid(x1, x2)

# V(x) = x'Px
V = P[0,0]*X1**2 + (P[0,1]+P[1,0])*X1*X2 + P[1,1]*X2**2

contour = ax1.contour(X1, X2, V, levels=20, cmap='viridis')
ax1.clabel(contour, inline=True, fontsize=8)

# 画轨迹
from scipy.integrate import odeint

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

t = np.linspace(0, 5, 500)
for angle in np.linspace(0, 2*np.pi, 8, endpoint=False):
    x0 = [1.5*np.cos(angle), 1.5*np.sin(angle)]
    sol = odeint(linear_system, x0, t)
    ax1.plot(sol[:, 0], sol[:, 1], 'r-', linewidth=1, alpha=0.7)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('$ x_1$', fontsize=12)
ax1.set_ylabel('$ x_2$', fontsize=12)
ax1.set_title('Lyapunov Function V(x) = x\'Px\nwith Trajectories', fontsize=14, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# V 沿轨迹的变化
ax2 = axes[1]
x0 = [1.5, 1.0]
sol = odeint(linear_system, x0, t)

V_traj = []
for state in sol:
    V_traj.append(state @ P @ state)

ax2.plot(t, V_traj, 'b-', linewidth=2)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('V(x(t))', fontsize=12)
ax2.set_title('Lyapunov Function Decreases Along Trajectory', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

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

二维系统的分类

完整分类

对于二维线性系统，根据特征值 分类：

类型	特征值条件	稳定性	相图
稳定节点	（实）	渐近稳定	轨迹趋向原点
不稳定节点	（实）	不稳定	轨迹远离原点
鞍点	（实）	不稳定	鞍型轨迹
稳定焦点	，（复）	渐近稳定	螺旋趋向原点
不稳定焦点	，（复）	不稳定	螺旋远离原点
中心	，（纯虚）	稳定（非渐近）	闭合轨道
退化节点	，缺特征向量	取决于符号	星形/退化

用迹和行列式分类

设，。

判别式： - ：鞍点 - ，，：稳定节点 - ，，：不稳定节点 - ，，：稳定焦点 - ，，：不稳定焦点 - ，：中心 - ：非孤立平衡点（一条直线）

import numpy as np
import matplotlib.pyplot as plt

# 迹-行列式平面的分类图
fig, ax = plt.subplots(1, 1, figsize=(10, 8))

tau = np.linspace(-4, 4, 400)
delta = np.linspace(-4, 4, 400)
Tau, Delta = np.meshgrid(tau, delta)

# 分界线
# 1. Delta = 0 (x 轴)
ax.axhline(y=0, color='black', linewidth=2)

# 2. tau = 0 (y 轴，在 Delta > 0 区域)
ax.axvline(x=0, color='black', linewidth=2)

# 3. 判别式 D = tau^2 - 4*Delta = 0，即 Delta = tau^2/4
tau_curve = np.linspace(-4, 4, 200)
delta_curve = tau_curve**2 / 4
ax.plot(tau_curve, delta_curve, 'k-', linewidth=2)

# 填充区域
# 稳定节点： Delta > 0, tau < 0, tau^2 > 4*Delta
ax.fill_between(tau[tau < 0], 0, tau[tau < 0]**2/4, 
                alpha=0.3, color='blue', label='Stable Node')

# 稳定焦点： Delta > tau^2/4, tau < 0
tau_neg = tau[tau < 0]
ax.fill_between(tau_neg, tau_neg**2/4, 4, 
                alpha=0.3, color='cyan', label='Stable Spiral')

# 不稳定节点： Delta > 0, tau > 0, tau^2 > 4*Delta
ax.fill_between(tau[tau > 0], 0, tau[tau > 0]**2/4, 
                alpha=0.3, color='red', label='Unstable Node')

# 不稳定焦点： Delta > tau^2/4, tau > 0
tau_pos = tau[tau > 0]
ax.fill_between(tau_pos, tau_pos**2/4, 4, 
                alpha=0.3, color='orange', label='Unstable Spiral')

# 鞍点： Delta < 0
ax.fill_between(tau, -4, 0, alpha=0.3, color='yellow', label='Saddle')

# 中心（ tau = 0, Delta > 0）用线标记
ax.plot([0, 0], [0, 4], 'g-', linewidth=3, label='Center (τ=0, Δ>0)')

# 标注
ax.text(-2.5, 2, 'Stable\nSpiral', fontsize=12, ha='center', fontweight='bold')
ax.text(-2.5, 0.3, 'Stable\nNode', fontsize=12, ha='center', fontweight='bold')
ax.text(2.5, 2, 'Unstable\nSpiral', fontsize=12, ha='center', fontweight='bold')
ax.text(2.5, 0.3, 'Unstable\nNode', fontsize=12, ha='center', fontweight='bold')
ax.text(0, -2, 'Saddle', fontsize=12, ha='center', fontweight='bold')

ax.set_xlabel('τ = tr(A)', fontsize=14)
ax.set_ylabel('Δ = det(A)', fontsize=14)
ax.set_title('Classification of 2D Linear Systems', fontsize=16, fontweight='bold')
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.legend(loc='upper left', fontsize=10)
ax.grid(True, alpha=0.3)

# 添加数学关系
ax.text(2.8, 3.5, '$D = \\tau^2 - 4\\Delta = 0$', fontsize=10)
ax.annotate('', xy=(2.5, 1.56), xytext=(2.8, 3.3),
            arrowprops=dict(arrowstyle='->', color='black'))

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

分岔理论初步

什么是分岔？

当系统参数变化时，平衡点的数量或稳定性可能发生突变。这种突变称为分岔（ bifurcation）。

鞍点-节点分岔

考虑一维系统： - 当$ r < 0 x = r = 0 x = 0 r > 0$：没有平衡点

在 处，两个平衡点"碰撞并消失"。

跨临界分岔

平衡点： 和 - 当$ r < 0 x = 0 x = r r > 0 x = 0 x = r$ 稳定

在 处，两个平衡点"交换稳定性"。

叉形分岔（ Pitchfork）

超临界叉形分岔： - 当$ r < 0 x = 0 r > 0 x = 0 x = $ 亚临界叉形分岔：

分岔方向相反，可能导致"跳跃"行为。

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

# 三种分岔类型
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 鞍点-节点分岔
ax1 = axes[0]
r_range = np.linspace(-2, 2, 200)

# 画平衡点
r_neg = r_range[r_range < 0]
x_eq_pos = np.sqrt(-r_neg)
x_eq_neg = -np.sqrt(-r_neg)

ax1.plot(r_neg, x_eq_pos, 'b-', linewidth=2, label='Stable')
ax1.plot(r_neg, x_eq_neg, 'r--', linewidth=2, label='Unstable')
ax1.plot(0, 0, 'ko', markersize=8)

ax1.set_xlabel('r', fontsize=12)
ax1.set_ylabel('x*', fontsize=12)
ax1.set_title('Saddle-Node Bifurcation\n$\\dot{x} = r + x^2$', fontsize=12, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.axhline(y=0, color='k', linewidth=0.5)
ax1.axvline(x=0, color='k', linewidth=0.5)

# 跨临界分岔
ax2 = axes[1]
ax2.plot(r_range, np.zeros_like(r_range), 'b-', linewidth=2, label='x* = 0')
ax2.plot(r_range, r_range, 'r-', linewidth=2, label='x* = r')

# 稳定性标记
ax2.plot(r_range[r_range < 0], np.zeros_like(r_range[r_range < 0]), 'b-', linewidth=3)
ax2.plot(r_range[r_range > 0], np.zeros_like(r_range[r_range > 0]), 'b--', linewidth=3)
ax2.plot(r_range[r_range < 0], r_range[r_range < 0], 'r--', linewidth=3)
ax2.plot(r_range[r_range > 0], r_range[r_range > 0], 'r-', linewidth=3)

ax2.plot(0, 0, 'ko', markersize=8)

ax2.set_xlabel('r', fontsize=12)
ax2.set_ylabel('x*', fontsize=12)
ax2.set_title('Transcritical Bifurcation\n$\\dot{x} = rx - x^2$', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(-2, 2)
ax2.set_ylim(-2, 2)
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=0, color='k', linewidth=0.5)
ax2.text(-1.5, -1.5, 'Unstable', fontsize=10, color='red')
ax2.text(1.5, 1.5, 'Stable', fontsize=10, color='red')
ax2.text(-1.5, 0.2, 'Stable', fontsize=10, color='blue')
ax2.text(1.5, 0.2, 'Unstable', fontsize=10, color='blue')

# 叉形分岔（超临界）
ax3 = axes[2]
r_pos = r_range[r_range >= 0]

ax3.plot(r_range, np.zeros_like(r_range), 'k-', linewidth=2)
ax3.plot(r_range[r_range <= 0], np.zeros_like(r_range[r_range <= 0]), 'b-', linewidth=3)
ax3.plot(r_range[r_range >= 0], np.zeros_like(r_range[r_range >= 0]), 'r--', linewidth=3)

# 新的稳定分支
x_branch = np.sqrt(r_pos)
ax3.plot(r_pos, x_branch, 'b-', linewidth=3)
ax3.plot(r_pos, -x_branch, 'b-', linewidth=3)

ax3.plot(0, 0, 'ko', markersize=8)

ax3.set_xlabel('r', fontsize=12)
ax3.set_ylabel('x*', fontsize=12)
ax3.set_title('Supercritical Pitchfork\n$\\dot{x} = rx - x^3$', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim(-2, 2)
ax3.set_ylim(-2, 2)
ax3.axhline(y=0, color='k', linewidth=0.5)
ax3.axvline(x=0, color='k', linewidth=0.5)
ax3.text(-1, 0.2, 'Stable', fontsize=10, color='blue')
ax3.text(1, 0.2, 'Unstable', fontsize=10, color='red')
ax3.text(1.5, 1.2, 'Stable', fontsize=10, color='blue')

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

Hopf 分岔

极限环的诞生

当复特征值的实部从负变正时，会发生Hopf 分岔：平衡点从稳定焦点变为不稳定焦点，同时产生一个极限环。

超临界 Hopf 分岔：极限环是稳定的 亚临界 Hopf 分岔：极限环是不稳定的

例子

$

在极坐标 下： $

- 当：原点是稳定焦点 - 当：原点是不稳定焦点，存在稳定极限环

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

# Hopf 分岔示例
def hopf_system(X, t, mu):
    x, y = X
    r2 = x**2 + y**2
    dxdt = mu*x - y - x*r2
    dydt = x + mu*y - y*r2
    return [dxdt, dydt]

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

mu_values = [-0.5, 0, 0.5]
t = np.linspace(0, 50, 2000)

for i, mu in enumerate(mu_values):
    ax1 = axes[0, i]
    ax2 = axes[1, i]
    
    # 相图
    # 从不同初值出发
    for r0 in [0.3, 0.6, 1.0, 1.5]:
        for theta0 in [0, np.pi/2, np.pi, 3*np.pi/2]:
            x0 = [r0*np.cos(theta0), r0*np.sin(theta0)]
            sol = odeint(hopf_system, x0, t, args=(mu,))
            ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.5, alpha=0.7)
    
    # 画极限环（如果存在）
    if mu > 0:
        theta = np.linspace(0, 2*np.pi, 100)
        r_lim = np.sqrt(mu)
        ax1.plot(r_lim*np.cos(theta), r_lim*np.sin(theta), 'r-', linewidth=2, label='Limit cycle')
    
    ax1.plot(0, 0, 'ko', markersize=8)
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title(f'μ = {mu}', fontsize=12, fontweight='bold')
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-2, 2)
    ax1.set_ylim(-2, 2)
    if mu > 0:
        ax1.legend(fontsize=10)
    
    # 时间演化
    x0 = [1.5, 0]
    sol = odeint(hopf_system, x0, t, args=(mu,))
    r_traj = np.sqrt(sol[:, 0]**2 + sol[:, 1]**2)
    
    ax2.plot(t, r_traj, 'b-', linewidth=1.5)
    if mu > 0:
        ax2.axhline(y=np.sqrt(mu), color='r', linestyle='--', linewidth=2, label=f'$ r = \\sqrt{{ \\mu }} = {np.sqrt(mu):.2f}$')
        ax2.legend(fontsize=10)
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('r(t)', fontsize=12)
    ax2.set_title(f'Radius vs Time (μ = {mu})', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)

plt.suptitle('Supercritical Hopf Bifurcation', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('hopf_bifurcation.png', dpi=150, bbox_inches='tight')
plt.show()

# 分岔图
fig, ax = plt.subplots(figsize=(10, 6))

mu_range = np.linspace(-1, 1, 200)

# 原点的稳定性
ax.plot(mu_range[mu_range <= 0], np.zeros(np.sum(mu_range <= 0)), 'b-', linewidth=3, label='Stable eq.')
ax.plot(mu_range[mu_range >= 0], np.zeros(np.sum(mu_range >= 0)), 'r--', linewidth=3, label='Unstable eq.')

# 极限环振幅
mu_pos = mu_range[mu_range > 0]
ax.plot(mu_pos, np.sqrt(mu_pos), 'g-', linewidth=3, label='Stable limit cycle')
ax.plot(mu_pos, -np.sqrt(mu_pos), 'g-', linewidth=3)

ax.axvline(x=0, color='k', linestyle=':', linewidth=1)
ax.set_xlabel('μ (bifurcation parameter)', fontsize=14)
ax.set_ylabel('Amplitude', fontsize=14)
ax.set_title('Hopf Bifurcation Diagram', fontsize=16, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)

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

应用：生态系统的稳定性

Lotka-Volterra 捕食者-猎物模型

$

平衡点： 1. ：两种都灭绝 2. ：共存

在共存平衡点处的 Jacobi 矩阵： $$

A =

$$

特征值：（纯虚）

结论：共存平衡点是中心（非双曲），线性化分析失效。实际上，系统有无穷多个闭合轨道，种群周期性振荡。

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

# Lotka-Volterra 模型
a = 1.0  # 猎物增长率
b = 0.5  # 捕食率
c = 0.5  # 捕食者死亡率
d = 0.2  # 转化效率

def lotka_volterra(X, t):
    x, y = X
    dxdt = a*x - b*x*y
    dydt = -c*y + d*x*y
    return [dxdt, dydt]

# 平衡点
x_eq = c/d
y_eq = a/b

print(f"共存平衡点: ({x_eq:.2f}, {y_eq:.2f})")

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

# 相空间
ax1 = axes[0]
t = np.linspace(0, 50, 2000)

# 不同初值的轨迹
initial_conditions = [
    [1, 1], [2, 1], [3, 1], [4, 1], [5, 1],
    [1, 2], [1, 3], [2, 3], [3, 2]
]

for x0 in initial_conditions:
    sol = odeint(lotka_volterra, x0, t)
    ax1.plot(sol[:, 0], sol[:, 1], linewidth=1, alpha=0.7)

ax1.plot(x_eq, y_eq, 'ro', markersize=10, label=f'Equilibrium ({x_eq:.1f}, {y_eq:.1f})')
ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('Prey x', fontsize=12)
ax1.set_ylabel('Predator y', fontsize=12)
ax1.set_title('Lotka-Volterra: Closed Orbits (Center)', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 时间演化
ax2 = axes[1]
x0 = [4, 1]
sol = odeint(lotka_volterra, x0, t)

ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='Prey x(t)')
ax2.plot(t, sol[:, 1], 'r-', linewidth=2, label='Predator y(t)')
ax2.axhline(y=x_eq, color='b', linestyle='--', alpha=0.5)
ax2.axhline(y=y_eq, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Population Dynamics: Periodic Oscillations', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

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

# 能量函数（积分常数）
# V = d*x - c*ln(x) + b*y - a*ln(y) = const
fig, ax = plt.subplots(figsize=(10, 8))

x_range = np.linspace(0.1, 8, 200)
y_range = np.linspace(0.1, 6, 200)
X, Y = np.meshgrid(x_range, y_range)

V = d*X - c*np.log(X) + b*Y - a*np.log(Y)

contour = ax.contour(X, Y, V, levels=30, cmap='viridis')
ax.clabel(contour, inline=True, fontsize=8)

ax.plot(x_eq, y_eq, 'ro', markersize=10)
ax.set_xlabel('Prey x', fontsize=12)
ax.set_ylabel('Predator y', fontsize=12)
ax.set_title('Lotka-Volterra: Conserved Quantity (Level Sets)', fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)

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

应用：控制系统设计

反馈稳定化

考虑不稳定系统，如何通过反馈使其稳定？

状态反馈： 闭环系统： 目标：选择 使得 的所有特征值有负实部。

极点配置定理：如果 可控，则可以任意指定 的特征值。

例子：倒立摆控制

倒立摆线性化方程：

其中 是施加在小车上的力。

import numpy as np
from scipy.linalg import solve_continuous_are
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 倒立摆参数
M = 1.0  # 小车质量
m = 0.1  # 摆杆质量
l = 0.5  # 摆杆长度
g = 9.81  # 重力加速度

# 状态空间矩阵（线性化在θ=0 附近）
A = np.array([
    [0, 1, 0, 0],
    [0, 0, -m*g/M, 0],
    [0, 0, 0, 1],
    [0, 0, (M+m)*g/(M*l), 0]
])

B = np.array([[0], [1/M], [0], [-1/(M*l)]])

print("开环特征值:", np.linalg.eigvals(A))
print("系统不稳定（有正特征值）")

# LQR 设计
Q = np.diag([10, 1, 100, 1])  # 状态权重
R = np.array([[0.1]])  # 控制权重

# 求解 Riccati 方程
P = solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print("\n 反馈增益 K:", K.flatten())
print("闭环特征值:", np.linalg.eigvals(A - B @ K))
print("系统稳定（所有特征值实部<0）")

# 仿真
def inverted_pendulum(X, t, K, A, B):
    u = -K @ X
    return (A @ X + B @ u).flatten()

def inverted_pendulum_open(X, t, A):
    return (A @ X).flatten()

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

t = np.linspace(0, 5, 500)
X0 = [0, 0, 0.2, 0]  # 初始角度偏离 0.2 rad

# 开环响应
sol_open = odeint(inverted_pendulum_open, X0, t, args=(A,))

# 闭环响应
sol_closed = odeint(inverted_pendulum, X0, t, args=(K, A, B))

# 角度
ax1 = axes[0, 0]
ax1.plot(t, sol_open[:, 2] * 180/np.pi, 'r--', linewidth=2, label='Open-loop')
ax1.plot(t, sol_closed[:, 2] * 180/np.pi, 'b-', linewidth=2, label='Closed-loop')
ax1.set_xlabel('Time (s)', fontsize=12)
ax1.set_ylabel('Angle θ (degrees)', fontsize=12)
ax1.set_title('Pendulum Angle', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 小车位置
ax2 = axes[0, 1]
ax2.plot(t, sol_open[:, 0], 'r--', linewidth=2, label='Open-loop')
ax2.plot(t, sol_closed[:, 0], 'b-', linewidth=2, label='Closed-loop')
ax2.set_xlabel('Time (s)', fontsize=12)
ax2.set_ylabel('Position x (m)', fontsize=12)
ax2.set_title('Cart Position', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# 角速度
ax3 = axes[1, 0]
ax3.plot(t, sol_open[:, 3] * 180/np.pi, 'r--', linewidth=2, label='Open-loop')
ax3.plot(t, sol_closed[:, 3] * 180/np.pi, 'b-', linewidth=2, label='Closed-loop')
ax3.set_xlabel('Time (s)', fontsize=12)
ax3.set_ylabel('Angular velocity (deg/s)', fontsize=12)
ax3.set_title('Angular Velocity', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# 控制输入
ax4 = axes[1, 1]
u_closed = np.array([-K @ sol_closed[i] for i in range(len(t))]).flatten()
ax4.plot(t, u_closed, 'b-', linewidth=2)
ax4.set_xlabel('Time (s)', fontsize=12)
ax4.set_ylabel('Control Force u (N)', fontsize=12)
ax4.set_title('Control Input', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.suptitle('Inverted Pendulum Stabilization via LQR', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('inverted_pendulum.png', dpi=150, bbox_inches='tight')
plt.show()

总结

本章要点

1. Lyapunov 稳定性：
   • 稳定：不跑远
   • 渐近稳定：最终回到平衡点
   • 不稳定：会跑远
2. 线性化方法：
   • Jacobi 矩阵
   • 特征值决定稳定性
   • Hartman-Grobman 定理
3. Lyapunov 直接法：
   • Lyapunov 函数（广义能量）
   • LaSalle 不变原理
4. 分岔理论：
   • 鞍点-节点、跨临界、叉形分岔
   • Hopf 分岔与极限环
5. 应用：
   • 生态系统、电路、控制

下一章预告

在《数值方法》中，我们将： - 学习各种 ODE 数值求解算法 - 理解步长选择和误差控制 - 处理刚性方程 - 实现自适应方法

练习题

基础题

1. 判断以下系统原点的稳定性：

   • , - , （在原点附近）

2. 对于，用 Lyapunov 函数 分析原点的稳定性。

3. 找出方程 的所有分岔点，并画出分岔图。

4. 对于阻尼单摆，证明总能量 是 Lyapunov 函数。

5. 计算以下矩阵的迹和行列式，并判断平衡点类型： $$

A =

$$

进阶题

6. 双稳态系统：分析 的所有平衡点及其稳定性。画出- 图和轨迹。

7. 滞后现象：对于，

   • 画出不同 值下的- 图
   • 找出分岔点
   • 解释为什么会出现滞后

8. 证明如果 是 Lyapunov 函数且，则所有有界轨迹趋向平衡点集。

9. 竞争排斥： Lotka-Volterra 竞争模型：

分析当 和 时的稳定性。

编程题

10. 编写程序自动判断 2D 线性系统的平衡点类型（给定矩阵）。

11. 实现 Lyapunov 方程求解器，并用它验证线性系统的稳定性。

12. 模拟 Van der Pol 振子，制作动画展示极限环的形成。

13. 设计一个倒立摆的 PID 控制器，并与 LQR 控制器比较性能。

思考题

14. 为什么 Hartman-Grobman 定理要求特征值实部非零？举一个反例说明中心情况下线性化可能失效。

15. 物理系统的能量总是好的 Lyapunov 函数候选吗？什么情况下不是？

16. 分岔点附近的"临界慢化"现象（ critical slowing down）是什么？它有什么实际意义？

参考资料

1. Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos. CRC Press.
2. Khalil, H. K. (2002). Nonlinear Systems. Prentice Hall.
3. Perko, L. (2001). Differential Equations and Dynamical Systems. Springer.
4. Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.

5. Slotine, J.-J. E., & Li, W. (1991). Applied Nonlinear Control. Prentice Hall.

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