Import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

1. 定义系统参数：二阶阻尼系统

zeta = 0.4 # 阻尼系数 (0 < zeta < 1，欠阻尼状态，呈现优美的螺旋收敛)
omega_n = 2.0 # 自然频率

模拟时间序列

t = np.linspace(0, 10, 400)

初始扰动点 (x0, v0)

x0, v0 = 3.0, 2.0

求解二阶微分方程的解析解（相空间轨迹：位移 x 与速度 v）

衰减系数

alpha = zeta * omega_n
omega_d = omega_n * np.sqrt(1 - zeta**2)

计算随时间变化的位移 x(t) 和速度 v(t)

x(t) = e^(-alpha*t) * (C1*cos(omega_d*t) + C2*sin(omega_d*t))

C1 = x0
C2 = (v0 + alpha * x0) / omega_d
x = np.exp(-alpha * t) * (C1 * np.cos(omega_d * t) + C2 * np.sin(omega_d * t))
v = np.exp(-alpha * t) * (-alpha * C1 * np.cos(omega_d * t) - alpha * C2 * np.sin(omega_d * t)
- C1 * omega_d * np.sin(omega_d * t) + C2 * omega_d * np.cos(omega_d * t))

2. 绘制相空间图 (Phase Portrait) 与李雅普诺夫稳定域

fig, ax = plt.subplots(figsize=(8, 8), facecolor='#0d1117')
ax.set_facecolor('#0d1117')

绘制李雅普诺夫稳定域的同心椭圆（代表能量等高线 V(x,v) = C）

theta = np.linspace(0, 2*np.pi, 200)
for r in np.linspace(0.5, 4.0, 8):
# 简化的能量椭圆边界
ax.plot(r * np.cos(theta), r * np.sin(theta) / omega_n,
color='#1f6feb', alpha=0.3, linestyle='--', linewidth=1)

动图元素初始化

line, = ax.plot([], [], color='#58a6ff', linewidth=2.5, label='系统演化轨迹 (Trajectory)')
point, = ax.plot([], [], marker='o', color='#f85149', markersize=10, label='当前状态 (Perturbed State)')
origin, = ax.plot([0], [0], marker='*', color='#3fb950', markersize=15, label='平衡态 (Equilibrium)')

图表美化

ax.set_xlim(-4.5, 4.5)
ax.set_ylim(-4.5, 4.5)
ax.axhline(0, color='#30363d', linewidth=1)
ax.axvline(0, color='#30363d', linewidth=1)
ax.set_title("Lyapunov Stability Region & Damping Convergence (ζ = {})".format(zeta),
color='#c9d1d9', fontsize=14, pad=15)
ax.set_xlabel("Displacement x (位移)", color='#8b949e', fontsize=12)
ax.set_ylabel("Velocity v (速度)", color='#8b949e', fontsize=12)
ax.tick_params(colors='#8b949e')
for spine in ax.spines.values():
spine.set_edgecolor('#30363d')
ax.legend(facecolor='#161b22', edgecolor='#30363d', labelcolor='#c9d1d9', loc='upper right')

3. 动画更新函数

def init():
line.set_data([], [])
point.set_data([], [])
return line, point

def update(frame):
# 轨迹动态延伸
line.set_data(x[:frame], v[:frame])
# 当前扰动点位置
point.set_data([x[frame]], [v[frame]])
return line, point

生成动画

ani = animation.FuncAnimation(fig, update, frames=len(t), init_func=init, interval=30, blit=True)

保存为彩色的动图文件

ani.save('lyapunov_damping.gif', writer='pillow', fps=30)
plt.show()