常微分方程（八）非线性系统与相图

生态系统中狐狸和兔子的数量为何会周期性波动？为什么心脏能保持稳定的节律跳动？神经元如何在"静息"和"放电"之间切换？这些看似无关的问题，背后都隐藏着非线性动力学的秘密。当我们从线性系统跨入非线性世界，数学的面貌将发生根本性改变——叠加原理失效、解可能不唯一、小扰动可能引发巨大变化。但正是这些"麻烦"，让非线性系统能够描述自然界最丰富、最迷人的现象。

从线性到非线性：思维的跃迁

在前面的章节中，我们系统学习了线性微分方程。线性系统有一个美妙的性质：叠加原理——如果和都是解，那么也是解。这个性质让线性系统变得"可驯服"：我们可以把复杂问题分解成简单部分，各个击破后再组合。

但现实世界很少是线性的。考虑一个简单的例子：人口增长。如果增长率恒定，我们得到指数增长，解是。但这预言人口会无限增长——显然荒谬！真实世界有资源限制，增长率会随人口增加而下降。这就引出了著名的 Logistic 方程：

这里是环境承载力。当很小时，方程近似为（指数增长）；当接近时，增长率趋近于零。这个方程是非线性的——右边出现了项。

非线性系统的特点： 1. 叠加原理失效：两个解相加不再是解 2. 解的行为可能极其复杂：周期振荡、准周期、混沌... 3. 对初值敏感：微小的初始差异可能导致完全不同的结果（这是混沌的种子，下一章详述） 4. 通常没有解析解：必须依赖定性分析和数值方法

相空间：可视化系统状态

状态空间的概念

考虑一个二阶常微分方程，比如带阻尼的弹簧： $$

m + b + kx = 0 $$

我们可以引入速度 $ v = $，把它改写成一阶系统：$ $

这是一个关于状态向量 的一阶方程组。所有可能的取值构成的空间叫做相空间（ phase space）或状态空间。

关键洞察：系统在任意时刻的状态由相空间中的一个点完全确定。随着时间演化，这个点在相空间中移动，画出一条轨迹（ trajectory）或轨线（ orbit）。

相图：系统行为的全景图

相图（ phase portrait）是相空间中轨迹的集合。它展示了不同初始条件下系统如何演化，是理解非线性系统最强大的工具之一。

让我们用 Python 绘制阻尼弹簧的相图：

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

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

def damped_oscillator(state, t, m, b, k):
    """阻尼振子的状态方程"""
    x, v = state
    dxdt = v
    dvdt = -(k/m) * x - (b/m) * v
    return [dxdt, dvdt]

# 参数：质量、阻尼系数、弹簧常数
m, b, k = 1.0, 0.3, 1.0

# 创建相图
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 不同阻尼情况
damping_cases = [
    (0.0, 'Undamped (b=0)'),
    (0.3, 'Underdamped (b=0.3)'),
    (2.0, 'Overdamped (b=2.0)')
]

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

for ax, (b_val, title) in zip(axes, damping_cases):
    # 多个初始条件
    initial_conditions = [
        (2, 0), (0, 2), (-2, 0), (0, -2),
        (1, 1), (-1, -1), (1, -1), (-1, 1),
        (3, 0), (0, 3)
    ]
    
    for x0, v0 in initial_conditions:
        sol = odeint(damped_oscillator, [x0, v0], t, args=(m, b_val, k))
        ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.8, alpha=0.7)
        ax.plot(x0, v0, 'go', markersize=4)  # 起点
    
    # 绘制向量场
    x_range = np.linspace(-3.5, 3.5, 15)
    v_range = np.linspace(-3.5, 3.5, 15)
    X, V = np.meshgrid(x_range, v_range)
    dX = V
    dV = -(k/m) * X - (b_val/m) * V
    
    # 归一化箭头
    M = np.sqrt(dX**2 + dV**2)
    M[M == 0] = 1
    ax.quiver(X, V, dX/M, dV/M, M, cmap='Reds', alpha=0.4)
    
    ax.set_xlabel('Position x', fontsize=11)
    ax.set_ylabel('Velocity v', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.set_xlim(-4, 4)
    ax.set_ylim(-4, 4)
    ax.grid(True, alpha=0.3)
    ax.set_aspect('equal')
    ax.plot(0, 0, 'ro', markersize=8)  # 原点/平衡点

plt.tight_layout()
plt.savefig('phase_portraits_damping.png', dpi=150, bbox_inches='tight')
plt.close()

观察三种情况： - 无阻尼（$ b=0 $）：轨迹是闭合椭圆，系统永远振荡欠阻尼（$ b$ 较小）：轨迹螺旋收敛到原点 - 过阻尼（较大）：轨迹直接收敛到原点，不振荡

平衡点：系统的静止状态

平衡点的定义

对于自治系统（不显含时间的系统）：

如果，则是平衡点（ equilibrium point）或不动点（ fixed point）。

物理意义：如果系统处于平衡点，它将永远停留在那里——因为所有导数都是零，没有变化的驱动力。

平衡点的分类：线性化方法

在平衡点附近，我们可以用泰勒展开近似系统行为。设，其中是小偏离：

这里是雅可比矩阵（ Jacobian matrix）：

平衡点的性质由雅可比矩阵的特征值决定。对于二维系统，设特征值为：

特征值类型	平衡点类型	稳定性
（实数）	稳定结点（ stable node）	渐近稳定
（实数）	不稳定结点（ unstable node）	不稳定
	鞍点（ saddle）	不稳定
（复数）	稳定焦点（ stable spiral）	渐近稳定
（复数）	不稳定焦点（ unstable spiral）	不稳定
（纯虚数）	中心（ center）	临界稳定

Python 实现：自动分类平衡点

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

def classify_equilibrium(J):
    """根据雅可比矩阵的特征值分类平衡点"""
    eigenvalues = np.linalg.eigvals(J)
    
    real_parts = np.real(eigenvalues)
    imag_parts = np.imag(eigenvalues)
    
    # 判断是否有虚部
    has_imaginary = np.any(np.abs(imag_parts) > 1e-10)
    
    if has_imaginary:
        # 复特征值情况
        if real_parts[0] < -1e-10:
            return "Stable Spiral (稳定焦点)", eigenvalues
        elif real_parts[0] > 1e-10:
            return "Unstable Spiral (不稳定焦点)", eigenvalues
        else:
            return "Center (中心)", eigenvalues
    else:
        # 实特征值情况
        if np.all(real_parts < -1e-10):
            return "Stable Node (稳定结点)", eigenvalues
        elif np.all(real_parts > 1e-10):
            return "Unstable Node (不稳定结点)", eigenvalues
        elif real_parts[0] * real_parts[1] < 0:
            return "Saddle Point (鞍点)", eigenvalues
        else:
            return "Degenerate Case (退化情况)", eigenvalues

# 示例：分析不同系统
def example_system1(state):
    """线性系统: dx/dt = -x + y, dy/dt = -x - y"""
    x, y = state
    return [-x + y, -x - y]

def jacobian_system1(state):
    """雅可比矩阵（对于线性系统，与状态无关）"""
    return np.array([[-1, 1], [-1, -1]])

J = jacobian_system1([0, 0])
eq_type, eigs = classify_equilibrium(J)
print(f"System 1: {eq_type}")
print(f"Eigenvalues: {eigs}")

经典案例： Lotka-Volterra 捕食者-猎物模型

模型背景

1925 年，意大利数学家 Vito Volterra 研究亚得里亚海的鱼类种群时，提出了著名的捕食者-猎物模型。同时期，美国数学家 Alfred Lotka 独立得出了相同的方程。

想象一个简化的生态系统，只有兔子（猎物）和狐狸（捕食者）：

兔子：没有狐狸时，兔子会指数增长（食物充足）
狐狸：没有兔子时，狐狸会饿死（指数衰减）
相遇：狐狸吃兔子 → 兔子减少，狐狸增多

数学模型

设为猎物（兔子）数量，为捕食者（狐狸）数量：

参数解释： - ：猎物的自然增长率 - ：捕食率（每次相遇时猎物被吃的概率） - ：捕食者的自然死亡率 - ：捕食效率（吃掉猎物后转化为捕食者繁殖的效率）

平衡点分析

令和：

平衡点 1：— 物种灭绝

平衡点 2：— 共存平衡

在共存平衡点处计算雅可比矩阵：

特征值：纯虚数特征值！这意味着共存平衡点是中心，周围的轨迹是闭合曲线——周期解！

Python 模拟与相图

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

def lotka_volterra(state, t, alpha, beta, delta, gamma):
    """Lotka-Volterra 捕食者-猎物模型"""
    x, y = state  # x: 猎物, y: 捕食者
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# 参数设置
alpha = 1.0   # 兔子自然增长率
beta = 0.1    # 捕食率
delta = 0.075 # 捕食效率
gamma = 0.5   # 狐狸自然死亡率

# 计算平衡点
x_eq = gamma / delta
y_eq = alpha / beta
print(f"共存平衡点: ({x_eq:.2f}, {y_eq:.2f})")

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

# 左图：时间演化
t = np.linspace(0, 100, 2000)
x0, y0 = 10, 5  # 初始条件

sol = odeint(lotka_volterra, [x0, y0], t, args=(alpha, beta, delta, gamma))
x_sol, y_sol = sol[:, 0], sol[:, 1]

axes[0].plot(t, x_sol, 'b-', linewidth=2, label='Prey (Rabbits)')
axes[0].plot(t, y_sol, 'r-', linewidth=2, label='Predator (Foxes)')
axes[0].axhline(y=x_eq, color='b', linestyle='--', alpha=0.5)
axes[0].axhline(y=y_eq, color='r', linestyle='--', alpha=0.5)
axes[0].set_xlabel('Time', fontsize=12)
axes[0].set_ylabel('Population', fontsize=12)
axes[0].set_title('Population Dynamics Over Time', fontsize=13, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# 右图：相图
ax = axes[1]

# 多条轨迹
initial_conditions = [
    (5, 5), (10, 5), (15, 5), (20, 5),
    (10, 3), (10, 8), (10, 12)
]

for x0, y0 in initial_conditions:
    sol = odeint(lotka_volterra, [x0, y0], t, args=(alpha, beta, delta, gamma))
    ax.plot(sol[:, 0], sol[:, 1], linewidth=1.5, alpha=0.7)
    ax.plot(x0, y0, 'go', markersize=5)

# 向量场
x_range = np.linspace(0.5, 25, 15)
y_range = np.linspace(0.5, 18, 15)
X, Y = np.meshgrid(x_range, y_range)
dX = alpha * X - beta * X * Y
dY = delta * X * Y - gamma * Y
M = np.sqrt(dX**2 + dY**2)
M[M == 0] = 1
ax.quiver(X, Y, dX/M, dY/M, M, cmap='Greys', alpha=0.4)

# 标记平衡点
ax.plot(x_eq, y_eq, 'r*', markersize=15, label=f'Equilibrium ({x_eq:.1f}, {y_eq:.1f})')
ax.plot(0, 0, 'ko', markersize=8)

ax.set_xlabel('Prey Population x', fontsize=12)
ax.set_ylabel('Predator Population y', fontsize=12)
ax.set_title('Phase Portrait: Lotka-Volterra Model', fontsize=13, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 28)
ax.set_ylim(0, 20)

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

模型的生态学意义

从相图可以清晰看到种群振荡： 1. 兔子多 → 狐狸有充足食物 → 狐狸增加 2. 狐狸多 → 兔子被大量捕食 → 兔子减少 3. 兔子少 → 狐狸挨饿 → 狐狸减少 4. 狐狸少 → 兔子少受捕食 → 兔子增加 5. 回到 1，循环往复

这解释了自然界观察到的种群周期性波动！

模型的局限性

结构不稳定：任何参数的微小改变都会破坏闭合轨道
忽略了环境承载力：兔子可以无限增长（无狐狸时）
忽略了空间因素：所有个体均匀混合
忽略了时间延迟：繁殖需要时间

改进模型见下一节。

竞争模型：两个物种争夺资源

竞争排斥原理

两个物种竞争同一资源会发生什么？苏联生态学家 G.F. Gause 在 1930 年代用草履虫做实验发现：两个生态位完全重叠的物种不能长期共存。

数学模型

设两个物种的种群数量为和，它们各自服从 Logistic 增长，同时相互竞争：

参数解释： - ：各自的内禀增长率 - ：各自的环境承载力 - ：物种 2 对物种 1 的竞争系数（ 1 个 y 对 x 的影响相当于个 x） - ：物种 1 对物种 2 的竞争系数

平衡点与零等倾线

零等倾线（ nullcline）是相空间中导数为零的曲线： - $ x $零等倾线：$ = 0 x = 0 $或$ x + {12}y = K_1 y $零等倾线：$ = 0 y = 0 $或$ y + {21}x = K_2$ 两条非零等倾线的交点就是共存平衡点。

Python 分析与可视化

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

def competition_model(state, t, r1, r2, K1, K2, a12, a21):
    """两物种竞争模型"""
    x, y = state
    dxdt = r1 * x * (1 - (x + a12 * y) / K1)
    dydt = r2 * y * (1 - (y + a21 * x) / K2)
    return [dxdt, dydt]

# 四种竞争结果的参数设置
scenarios = {
    'Species 1 wins': {'r1': 1, 'r2': 1, 'K1': 100, 'K2': 80, 'a12': 0.5, 'a21': 1.5},
    'Species 2 wins': {'r1': 1, 'r2': 1, 'K1': 80, 'K2': 100, 'a12': 1.5, 'a21': 0.5},
    'Coexistence': {'r1': 1, 'r2': 1, 'K1': 100, 'K2': 100, 'a12': 0.5, 'a21': 0.5},
    'Bistability': {'r1': 1, 'r2': 1, 'K1': 100, 'K2': 100, 'a12': 1.5, 'a21': 1.5}
}

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

for ax, (scenario_name, params) in zip(axes, scenarios.items()):
    r1, r2 = params['r1'], params['r2']
    K1, K2 = params['K1'], params['K2']
    a12, a21 = params['a12'], params['a21']
    
    # 绘制零等倾线
    x_line = np.linspace(0, 150, 100)
    
    # x-nullcline: x + a12*y = K1 => y = (K1 - x) / a12
    y_nullcline_x = (K1 - x_line) / a12
    y_nullcline_x = np.clip(y_nullcline_x, 0, 150)
    
    # y-nullcline: y + a21*x = K2 => y = K2 - a21*x
    y_nullcline_y = K2 - a21 * x_line
    y_nullcline_y = np.clip(y_nullcline_y, 0, 150)
    
    ax.plot(x_line, y_nullcline_x, 'b-', linewidth=2, label='x-nullcline')
    ax.plot(x_line, y_nullcline_y, 'r-', linewidth=2, label='y-nullcline')
    
    # 绘制轨迹
    t = np.linspace(0, 100, 1000)
    initial_conditions = [
        (10, 10), (80, 10), (10, 80), (50, 50),
        (20, 60), (60, 20), (90, 90)
    ]
    
    for x0, y0 in initial_conditions:
        sol = odeint(competition_model, [x0, y0], t, 
                     args=(r1, r2, K1, K2, a12, a21))
        ax.plot(sol[:, 0], sol[:, 1], 'g-', linewidth=0.8, alpha=0.6)
        ax.plot(x0, y0, 'ko', markersize=4)
        ax.plot(sol[-1, 0], sol[-1, 1], 'r*', markersize=8)
    
    ax.set_xlabel('Species 1 (x)', fontsize=11)
    ax.set_ylabel('Species 2 (y)', fontsize=11)
    ax.set_title(scenario_name, fontsize=12, fontweight='bold')
    ax.legend(fontsize=9, loc='upper right')
    ax.set_xlim(0, 120)
    ax.set_ylim(0, 120)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('competition_model.png', dpi=150, bbox_inches='tight')
plt.close()

四种竞争结局

根据参数不同，竞争可能有四种结果：

物种 1 获胜：物种 2 灭绝
物种 2 获胜：物种 1 灭绝
共存：两物种达到稳定平衡
双稳态：谁先到达足够数量谁获胜（历史依赖）

判别条件： - 若且：物种 1 获胜 - 若且：物种 2 获胜 - 若且：稳定共存 - 若且：双稳态

Van der Pol 振子：自激振荡

从收音机到心脏

1920 年代，荷兰工程师 Balthasar van der Pol 研究收音机电路中的电子管时，发现了一种奇特的振荡器。它的独特之处在于：无论从什么初始状态出发，系统最终都会收敛到同一个周期振荡。

这种行为叫做极限环（ limit cycle）——一个孤立的闭合轨道，周围的轨迹要么螺旋进入，要么螺旋离开。

Van der Pol 方程：

改写成一阶系统：

参数控制非线性强度： - 当时，：负阻尼，系统获得能量 - 当时，：正阻尼，系统损失能量

这种"智能阻尼"使系统自动维持稳定振荡——正是心脏起搏器的工作原理！

Python 绘制极限环

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

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

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
mu_values = [0.5, 1.0, 3.0]
t = np.linspace(0, 50, 5000)

for ax, mu in zip(axes, mu_values):
    # 不同初始条件的轨迹
    initial_conditions = [
        (0.1, 0), (0.5, 0), (3, 0), (0, 3),
        (-2, -2), (2, 2), (4, 0)
    ]
    
    for x0, y0 in initial_conditions:
        sol = odeint(van_der_pol, [x0, y0], t, args=(mu,))
        ax.plot(sol[:, 0], sol[:, 1], linewidth=0.8, alpha=0.7)
        ax.plot(x0, y0, 'go', markersize=5)
    
    # 标记原点（不稳定平衡点）
    ax.plot(0, 0, 'ro', markersize=8)
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y = dx/dt', fontsize=12)
    ax.set_title(f'Van der Pol Oscillator (μ = {mu})', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-5, 5)
    ax.set_ylim(-8, 8)
    ax.set_aspect('equal')

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

弛豫振荡

当很大时，振荡变成弛豫振荡（ relaxation oscillation）： - 大部分时间系统缓慢演化 - 然后突然快速跳跃到另一状态 - 再缓慢演化...

这描述了许多"充电-放电"型现象： - 神经元的动作电位 - 心脏节律 - 间歇泉喷发

李雅普诺夫稳定性理论

稳定性的精确定义

直觉上，稳定的平衡点是"抗干扰"的——轻轻推一下，系统会回到平衡点或停留在附近。

李雅普诺夫稳定：对于任意，存在，使得只要初值满足，则对所有都有。

渐近稳定：李雅普诺夫稳定 +

李雅普诺夫函数方法

不用求解微分方程就能判断稳定性！核心思想：构造一个能量函数 ，如果它沿轨迹单调递减，系统就是稳定的。

李雅普诺夫第二方法：设是正定函数（，对）。计算沿轨迹的导数： - 若（负半定）：李雅普诺夫稳定 - 若（负定）：渐近稳定

例子：阻尼摆

考虑带阻尼的单摆：

改写为系统：，平衡点：构造李雅普诺夫函数（总能量）： $$

V(, ) = ^2 + (1 - ) $$

这是正定的（动能 + 势能，以最低点为零点）。

沿轨迹的导数：

因此原点是稳定的！当时（有阻尼），可以进一步证明渐近稳定。

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

def damped_pendulum(state, t, b):
    """阻尼单摆"""
    theta, omega = state
    dthetadt = omega
    domegadt = -b * omega - np.sin(theta)
    return [dthetadt, domegadt]

# 绘制李雅普诺夫函数的等高线和轨迹
fig, ax = plt.subplots(figsize=(10, 8))

b = 0.3  # 阻尼系数

# 李雅普诺夫函数 V = omega^2/2 + (1 - cos(theta))
theta_range = np.linspace(-np.pi, np.pi, 100)
omega_range = np.linspace(-3, 3, 100)
Theta, Omega = np.meshgrid(theta_range, omega_range)
V = 0.5 * Omega**2 + (1 - np.cos(Theta))

# 等高线
contour = ax.contour(Theta, Omega, V, levels=15, cmap='coolwarm', alpha=0.6)
ax.clabel(contour, inline=True, fontsize=8)

# 轨迹
t = np.linspace(0, 30, 1000)
initial_conditions = [
    (2.5, 0), (2, 1), (1, 2), (-2, 1), (-1, -2), (0.5, 0.5)
]

for theta0, omega0 in initial_conditions:
    sol = odeint(damped_pendulum, [theta0, omega0], t, args=(b,))
    ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=1.5, alpha=0.7)
    ax.plot(theta0, omega0, 'go', markersize=6)

ax.plot(0, 0, 'r*', markersize=15, label='Stable equilibrium')
ax.set_xlabel(r'$\theta$(angle)', fontsize=12)
ax.set_ylabel(r'$\omega$(angular velocity)', fontsize=12)
ax.set_title(f'Damped Pendulum Phase Portrait (b={b})\nContours: Lyapunov function V', 
             fontsize=12, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_xlim(-np.pi, np.pi)
ax.set_ylim(-3, 3)

plt.tight_layout()
plt.savefig('lyapunov_pendulum.png', dpi=150, bbox_inches='tight')
plt.close()

梯度系统与哈密顿系统

梯度系统：永远下坡

如果系统可以写成：

其中是某个势函数，这就是梯度系统。

性质： - 轨迹总是沿着下降最快的方向 - 不存在周期解（不可能绕圈） - 所有轨迹最终收敛到的极小点

机器学习中的梯度下降就是梯度系统的离散版本！

哈密顿系统：能量守恒

如果系统可以写成：

其中是哈密顿量（通常是总能量），这就是哈密顿系统。

性质： - 沿轨迹守恒： - 轨迹是的等高线 - 相空间体积守恒（刘维尔定理）

经典力学中的无阻尼系统都是哈密顿系统。

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

def simple_pendulum(state, t):
    """无阻尼单摆（哈密顿系统）"""
    theta, p = state  # p = d(theta)/dt
    # H = p^2/2 + (1 - cos(theta))
    dthetadt = p  # = dH/dp
    dpdt = -np.sin(theta)  # = -dH/dtheta
    return [dthetadt, dpdt]

fig, ax = plt.subplots(figsize=(10, 8))

# 哈密顿量的等高线
theta_range = np.linspace(-2*np.pi, 2*np.pi, 200)
p_range = np.linspace(-3, 3, 200)
Theta, P = np.meshgrid(theta_range, p_range)
H = 0.5 * P**2 + (1 - np.cos(Theta))

# 绘制等高线（也是轨迹！）
contour = ax.contour(Theta, P, H, levels=20, cmap='viridis')
ax.clabel(contour, inline=True, fontsize=8)

# 标记分界线（ separatrix） H = 2
separatrix = ax.contour(Theta, P, H, levels=[2], colors='red', linewidths=2)

# 标记平衡点
ax.plot(0, 0, 'go', markersize=10, label='Stable (oscillation center)')
ax.plot(np.pi, 0, 'ro', markersize=10, label='Unstable (inverted)')
ax.plot(-np.pi, 0, 'ro', markersize=10)

ax.set_xlabel(r'$\theta$', fontsize=12)
ax.set_ylabel(r'$ p = d\theta/dt$', fontsize=12)
ax.set_title('Phase Portrait of Simple Pendulum (Hamiltonian System)\nContours = Energy levels = Trajectories', 
             fontsize=12, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_xlim(-2*np.pi, 2*np.pi)
ax.set_ylim(-3.5, 3.5)

plt.tight_layout()
plt.savefig('hamiltonian_pendulum.png', dpi=150, bbox_inches='tight')
plt.close()

全局分析： Poincar é-Bendixson 定理

极限集

当时，轨迹会趋向什么？可能的极限集（ limit set）包括： - 平衡点 - 极限环 - 混沌吸引子（三维以上才可能，见下一章）

Poincar é-Bendixson 定理

对于二维连续系统，有一个惊人的结论：

定理：设是平面上的有界闭区域，轨迹在内且不趋向平衡点。则的极限集是一个闭轨道（周期解）。

推论：二维连续系统不可能出现混沌！

这就是为什么研究混沌必须至少三维（如洛伦兹系统）。

判断极限环存在

Bendixson 判据：如果在某区域内不变号（恒正或恒负），则该区域内不存在闭轨道。

Dulac 判据： Bendixson 判据的推广，允许乘以一个正因子。

生活中的非线性系统

案例 1：心脏节律

心脏的窦房结是天然的起搏器，它的电活动可以用 FitzHugh-Nagumo 模型描述——一个简化的 Van der Pol 型方程。

正常情况下，系统有稳定的极限环（正常心跳）。某些病理条件下，极限环可能消失或变形，导致心律失常。

案例 2：传染病传播

SIR 模型（见后续章节）是非线性系统。（基本传染数）决定了疾病能否爆发——这是一个分岔现象，下下章详述。

案例 3：气候系统

全球气候是典型的高维非线性系统。小冰期、温室效应、气候临界点都可以用非线性动力学理解。

案例 4：神经网络

人工神经网络的训练过程可以看作梯度系统在高维空间的演化。局部极小点、鞍点、平坦区域都影响训练效率。

数值方法简述

常用方法

欧拉法：最简单，但精度低 2. 改进欧拉法（ Heun）：二阶精度
四阶 Runge-Kutta（ RK4）：最常用，精度和效率平衡好

4. 自适应步长方法：如 RK45 、 DOPRI5

Python 实现 RK4

import numpy as np

def rk4_step(f, x, t, h):
    """单步 RK4"""
    k1 = h * np.array(f(x, t))
    k2 = h * np.array(f(x + k1/2, t + h/2))
    k3 = h * np.array(f(x + k2/2, t + h/2))
    k4 = h * np.array(f(x + k3, t + h))
    return x + (k1 + 2*k2 + 2*k3 + k4) / 6

def rk4_solve(f, x0, t_span, h):
    """用 RK4 求解 ODE 系统"""
    t_start, t_end = t_span
    t_values = [t_start]
    x_values = [np.array(x0)]
    
    t = t_start
    x = np.array(x0)
    
    while t < t_end:
        if t + h > t_end:
            h = t_end - t
        x = rk4_step(f, x, t, h)
        t = t + h
        t_values.append(t)
        x_values.append(x.copy())
    
    return np.array(t_values), np.array(x_values)

# 测试： Lotka-Volterra
def lotka_volterra_f(state, t):
    x, y = state
    alpha, beta, delta, gamma = 1.0, 0.1, 0.075, 0.5
    return [alpha*x - beta*x*y, delta*x*y - gamma*y]

t_vals, sol = rk4_solve(lotka_volterra_f, [10, 5], (0, 100), 0.01)
print(f"Final state: x = {sol[-1, 0]:.2f}, y = {sol[-1, 1]:.2f}")

总结

本章我们深入探索了非线性系统的核心概念：

相空间与相图：用几何直觉理解系统行为
平衡点分类：通过雅可比矩阵的特征值判断稳定性
经典模型：
- Lotka-Volterra（周期振荡）
- 竞争模型（四种结局）
- Van der Pol（极限环）
李雅普诺夫方法：不解方程判断稳定性
特殊结构：梯度系统、哈密顿系统
全局理论： Poincar é-Bendixson 定理

非线性世界比线性世界丰富得多：周期、准周期、极限环...而这还只是冰山一角。下一章，我们将见证非线性系统最惊人的行为——混沌。

练习题

概念题

解释为什么线性系统的叠加原理不适用于非线性系统。给出一个具体的反例。
什么是相空间？为什么它对分析动力系统很重要？
区分李雅普诺夫稳定和渐近稳定。给出各自的例子。
为什么二维连续系统不能出现混沌？

计算题

对于系统

找出所有平衡点
计算每个平衡点的雅可比矩阵
分类每个平衡点的类型

证明是系统，在原点的李雅普诺夫函数，并判断原点的稳定性。
对于 Lotka-Volterra 系统，证明函数沿轨迹是常数（首次积分）。
对于 Van der Pol 方程（），计算极限环的近似周期（可用数值方法）。

分析题

考虑带密度依赖死亡的捕食者-猎物模型：

分析额外的项如何影响系统的定性行为。

在竞争模型中，如果两个物种的参数完全对称（$ r_1 = r_2 $，$ K_1 = K_2 $，$ {12} = {21}$），系统会如何演化？这在生态学上意味着什么？

编程题

编写 Python 程序绘制下列系统的相图，并标出所有平衡点和它们的类型： 12. 实现一个交互式程序，让用户调整 Lotka-Volterra 模型的参数，实时观察相图和时间序列的变化。
用数值方法验证 Poincar é-Bendixson 定理：对于一个有界轨迹不趋向平衡点的二维系统，证明它确实收敛到周期轨道。
比较欧拉法、改进欧拉法和 RK4 在求解 Van der Pol 方程时的精度和效率。用不同步长测试，绘制误差-步长关系图。

参考文献

Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Westview Press.
Perko, L. (2001). Differential Equations and Dynamical Systems. Springer.
Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.
Murray, J. D. (2002). Mathematical Biology I: An Introduction. Springer.
Hirsch, M. W., Smale, S., & Devaney, R. L. (2012). Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press.

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

从线性到非线性：思维的跃迁

相空间：可视化系统状态

状态空间的概念

相图：系统行为的全景图

平衡点：系统的静止状态

平衡点的定义

平衡点的分类：线性化方法

Python 实现：自动分类平衡点

经典案例： Lotka-Volterra 捕食者-猎物模型

模型背景

数学模型

平衡点分析

Python 模拟与相图

模型的生态学意义

模型的局限性

竞争模型：两个物种争夺资源

竞争排斥原理

数学模型

平衡点与零等倾线

Python 分析与可视化

四种竞争结局

Van der Pol 振子：自激振荡

从收音机到心脏

Python 绘制极限环

弛豫振荡

李雅普诺夫稳定性理论

稳定性的精确定义

李雅普诺夫函数方法

例子：阻尼摆

梯度系统与哈密顿系统

梯度系统：永远下坡

哈密顿系统：能量守恒

全局分析： Poincar é-Bendixson 定理

极限集

Poincar é-Bendixson 定理

判断极限环存在

生活中的非线性系统

案例 1：心脏节律

案例 2：传染病传播

案例 3：气候系统

案例 4：神经网络

数值方法简述

常用方法

Python 实现 RK4

总结

练习题

概念题

计算题

分析题

编程题

参考文献

Hirsch, M. W., Smale, S., & Devaney, R. L. (2012). Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press.