二阶微分方程 欠阻尼

import numpy as np
import matplotlib.pyplot as plt

# 中文正常显示配置
plt.rcParams["font.family"] = ["SimHei", "Microsoft YaHei"]
plt.rcParams["axes.unicode_minus"] = False

# 欠阻尼参数
zeta = 0.3
omega_d = np.sqrt(1 - zeta ** 2)
t = np.linspace(0, 20, 2000)

# 1行3列子图布局
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6))

# ========== 子图1：固定相位，改变振幅 ==========
phi = 0
amp_list = [0.5, 1, 1.5, 2]
for A in amp_list:
    x = A * np.exp(-zeta * t) * np.cos(omega_d * t + phi)
    ax1.plot(t, x, lw=1.7, label=f'振幅 $A={A}$')
ax1.set_xlabel('$t$')
ax1.set_ylabel('$x(t)$')
ax1.set_title(r'$\ddot{x}+0.6\dot{x}+x=0$ 固定相位、不同振幅')
ax1.grid(alpha=0.3)
ax1.legend()

# ========== 子图2：固定振幅，改变相位 ==========
A = 1.5
phi_list = [0, np.pi/4, np.pi/2, np.pi]
for phi in phi_list:
    x = A * np.exp(-zeta * t) * np.cos(omega_d * t + phi)
    ax2.plot(t, x, lw=1.7, label=f'相位 $\phi={phi:.2f}$')
ax2.set_xlabel('$t$')
ax2.set_ylabel('$x(t)$')
ax2.set_title(r'$\ddot{x}+0.6\dot{x}+x=0$ 固定振幅、不同相位')
ax2.grid(alpha=0.3)
ax2.legend()

# ========== 子图3：带阻尼相平面轨迹（螺旋收敛到原点） ==========
theta = np.linspace(0, 12 * np.pi, 3000)
r_list = [0.5, 1, 1.5, 2]
for r0 in r_list:
    x = r0 * np.exp(-zeta * theta / omega_d) * np.cos(theta)
    dx = -r0 * np.exp(-zeta * theta / omega_d) * (zeta * np.cos(theta) + omega_d * np.sin(theta))
    ax3.plot(x, dx, lw=1.8, label=f'初始振幅 $r_0={r0}$')
ax3.set_xlabel('$x$')
ax3.set_ylabel(r'$\dot{x}$')
ax3.set_title(r'欠阻尼系统相平面轨迹（螺旋收敛至原点）')
ax3.grid(alpha=0.3)
ax3.axis('equal')
ax3.legend()

# 自动调整间距防止文字重叠
plt.tight_layout()
plt.show()