不同阻尼比下的相图

import numpy as np
import matplotlib.pyplot as plt

def second_order_system_response(omega_n, zeta, t):
    """
    计算二阶系统的响应。
    """
    if zeta == 1:
        return np.exp(-omega_n * t) * (1 + omega_n * t)
    elif zeta > 1:
        s1 = -omega_n * (zeta + np.sqrt(zeta**2 - 1))
        s2 = -omega_n * (zeta - np.sqrt(zeta**2 - 1))
        return np.exp(s1 * t) / (2 * np.sqrt(zeta**2 - 1)) * ((s2 - s1) * np.exp((2 * np.sqrt(zeta**2 - 1)) * t) + s2 + s1)
    else:
        wd = omega_n * np.sqrt(1 - zeta**2)
        return np.exp(-zeta * omega_n * t) * (np.cos(wd * t) + zeta / np.sqrt(1 - zeta**2) * np.sin(wd * t))

def plot_phase_trajectories(omega_n=1):
    """
    绘制不同阻尼比下的二阶系统相轨迹，分为三个子图。
    """
    t = np.linspace(0, 10, 1000)
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # 子图 1：低阻尼比
    zetas_low = [0.2, 0.5,0.7]
    for zeta in zetas_low:
        x = second_order_system_response(omega_n, zeta, t)
        v = np.diff(x) / np.diff(t)
        axes[0].plot(x[:-1], v, label=f'ζ = {zeta}')
    axes[0].set_title('Low Damping Ratio')
    axes[0].set_xlabel('Position (x)')
    axes[0].set_ylabel('Velocity (v)')
    axes[0].legend()
    axes[0].grid(True)

    # 子图 2：临界阻尼和过阻尼
    zetas_mid = [ 1, 1.2,1.5]
    for zeta in zetas_mid:
        x = second_order_system_response(omega_n, zeta, t)
        v = np.diff(x) / np.diff(t)
        axes[1].plot(x[:-1], v, label=f'ζ = {zeta}')
    axes[1].set_title('Medium to High Damping Ratio')
    axes[1].set_xlabel('Position (x)')
    axes[1].set_ylabel('Velocity (v)')
    axes[1].legend()
    axes[1].grid(True)

    # 子图 3：无阻尼
    zeta = 0
    x = second_order_system_response(omega_n, zeta, t)
    v = np.diff(x) / np.diff(t)
    axes[2].plot(x[:-1], v, label=f'ζ = {zeta}')
    axes[2].set_title('No Damping')
    axes[2].set_xlabel('Position (x)')
    axes[2].set_ylabel('Velocity (v)')
    axes[2].legend()
    axes[2].grid(True)

    plt.tight_layout()
    plt.show()

plot_phase_trajectories()