3.2 3.2 相平面与稳定性判据

3.2 相平面与稳定性判据

本节摘要：动力系统分析的核心问题不是算出每条轨迹，而是定性回答"系统最终去哪、在那个地方站得稳不稳"。本节以捕食者-猎物系统为例讲相平面的读法，用线性化与特征值给平衡点分型，并介绍李雅普诺夫直接法——不解方程也能证明稳定性的"能量函数"路线。

目标能力

阅读完本节，你应当能够：

1. 在相平面上识别零斜线、轨线方向与平衡点；
2. 对平衡点做线性化，用雅可比矩阵特征值分类（结点、鞍点、焦点）；
3. 理解李雅普诺夫函数的构造思想并验证一个实例。

相平面：把时间轴藏起来的读图法

生态学经典的捕食者-猎物模型：猎物数量增长率为"自然增长减捕食"，捕食者增长率为"捕食获益减自然死亡"。两个状态量（猎物 x、捕食者 y）构成一张平面，每个点上都有一支箭头指明系统的演化方向，箭头场就是相平面。

读图三步：先找零斜线（x 的变化率为零的曲线、y 的变化率为零的曲线），两族零斜线的交点就是平衡点；再看轨线在箭头场中的走向；最后盯住平衡点附近轨线是汇聚、发散还是打转。捕食者-猎物系统的相平面是一圈圈闭合轨道——猎物多了捕食者涨、捕食者多了猎物跌、跌完捕食者饿死、猎物再起。闭合轨道意味着周期振荡，不用解方程就能从图上读出来。

import numpy as np
from scipy.integrate import solve_ivp

def lotka(t, z, a=1.0, b=0.1, c=1.5, d=0.75):
    """猎物 z0 与捕食者 z1 的耦合方程"""
    x, y = z
    return [a*x - b*x*y, c*b*x*y - d*y]

# 从三个不同初值出发，看轨线如何绕行
for z0 in [(4, 1), (6, 2), (2.5, 1.5)]:
    sol = solve_ivp(lotka, [0, 12], z0, rtol=1e-8, dense_output=True)
    tt = np.linspace(0, 12, 400)
    traj = sol.sol(tt)
    period_guess = tt[np.argmax(traj[0][:100])] * 2
    print(f"初值 {z0}: 猎物峰值 {traj[0].max():.2f}, "
          f"捕食者峰值 {traj[1].max():.2f}")

# 零斜线交点即平衡点
x_star = 0.75 / (1.5 * 0.1)
y_star = 1.0 / 0.1
print(f"正平衡点: 猎物 {x_star:.2f}, 捕食者 {y_star:.2f}")
​

三组初值各自绕出不同的闭合圈：离平衡点越远圈越大、振荡越猛。这个结论有管理含义——若渔业部门想把系统"稳定"在某个产量水平，任何把它推离平衡点的收割政策都会换来更剧烈的振荡。

线性化：平衡点分型

平衡点附近把方程做一阶泰勒展开，得到雅可比矩阵，平衡点的"性格"由它的特征值决定：两个都为负实数是稳定结点（轨线直奔而去）；一正一负是鞍点（沿一个方向吸引、另一个方向排斥）；实部为负的共轭复数是稳定焦点（旋进）；实部为零的纯虚数是中心（闭合轨道，临界情形）。捕食者-猎物的正平衡点线性化后特征值是纯虚数——中心，与闭合轨道的图景一致。

import numpy as np

def jacobian_lotka(x, y, a=1.0, b=0.1, cb=0.15, d=0.75):
    """在点 处的雅可比矩阵"""
    return np.array([
        [a - b*y,        -b*x],
        [cb*y,  cb*x - d],
    ])

x_star, y_star = 0.75 / 0.15, 1.0 / 0.1
J = jacobian_lotka(x_star, y_star)
eig = np.linalg.eigvals(J)
print("雅可比特征值:", np.round(eig, 4))
print("实部最大者:", np.round(eig.real.max(), 6))

# 加一点猎物自限项后重新线性化：结果会从中心变成稳定焦点
def jacobian_selflimit(x, y, a=1.0, b=0.1, cb=0.15, d=0.75, K=20.0):
    return np.array([
        [a - 2*a*x/K - b*y,          -b*x],
        [cb*y,               cb*x - d],
    ])

# 自限情形的平衡点需数值求，这里取其近似位置演示判据用法
J2 = jacobian_selflimit(4.5, 2.0)
print("自限情形特征值:", np.round(np.linalg.eigvals(J2), 4))
​

线性化判据的适用边界要说清：它只保证在平衡点足够小的邻域内有效。中心是临界情形——非线性项可能把闭合轨道变成缓慢旋进或旋出，特征值判据在这里"不作证"。所以工程结论不能只报特征值，要在相平面上数值积分几条轨线做交叉验证——上一段代码正是这个用途。

李雅普诺夫直接法提供了第三条路：找一个"能量样"的正定函数 V，如果它沿系统轨线的导数处处非正，能量只降不升，平衡点就稳定。它不解方程、不线性化，适用范围比前两条路都宽，代价是 V 的构造靠功夫。对阻尼摆，取机械能（动能加势能）即得：阻尼使能量导数为负，摆最终停在最低点。

import numpy as np

def pendulum(t, s, g=9.81, L=1.0, c=0.3):
    """阻尼单摆：角度 s0 与角速度 s1"""
    theta, omega = s
    return [omega, -g/L*np.sin(theta) - c*omega]

def energy(s, g=9.81, L=1.0, m=1.0):
    """李雅普诺夫函数：机械能（动能+势能）"""
    theta, omega = s
    return 0.5*m*(L*omega)**2 + m*g*L*(1 - np.cos(theta))

# 沿轨线采样能量：应单调下降
from scipy.integrate import solve_ivp
sol = solve_ivp(pendulum, [0, 15], [2.5, 0.0], rtol=1e-9, dense_output=True)
tt = np.linspace(0, 15, 60)
E = [energy(sol.sol(t)) for t in tt]
drops = np.diff(E)
print("能量沿轨线单调不增:", bool((drops <= 1e-9).all()))
print(f"初始能量 {E[0]:.3f} → 终末能量 {E[-1]:.3f}")
​

能量从初值一路降到零附近且全程单调——这就是李雅普诺夫思想的数值现身：稳定性证明被转化为一个标量函数的单调性检查。在现代控制与优化（第 6、8 章）里，这一思想反复出现：神经网络的损失函数、优化算法的下降函数，本质都是李雅普诺夫函数的亲戚。

一个常见误区是只看线性化就宣布全局稳定。线性化只管邻域，全局行为必须配合相平面数值实验或李雅普诺夫论证——"鞍点附近轨线先被吸引再被甩出"的误判，多半源于把局部结论当全局用。

💡 关键直觉：雅可比的特征值是"平衡点的体温计"。它不能告诉你病人整体健康状况（全局动力学），但能立刻告诉你局部有没有急性炎症（不稳定方向）。查完体温还要做全套检查，两者顺序不能反。

平衡点分型速查

本节要点回顾

• 相平面读图三步：零斜线找平衡点、看轨线走向、判断汇散旋；
• 闭合轨道意味着周期振荡，图形结论先于数值细节；
• 雅可比特征值给平衡点分型：结点、鞍点、焦点、中心各有判据；
• 临界情形不作证：中心附近要靠非线性分析或数值实验收尾；
• 李雅普诺夫函数把稳定性证明化成标量单调性，是能量思想的推广。