Python中用梯形法则绘制摆运动：θ、ω随时间变化图绘制遇阻

嘿，我看到你尝试用梯形法则求解无摩擦、无驱动力的线性摆微分方程，还想绘制角度θ和角速度ω随时间t的变化曲线，但卡在了图像生成这一步。先帮你梳理下问题所在，再给出完整的可运行代码吧！

问题核心梳理 #

你给出的代码片段只完成了初始参数定义，缺少梯形法则的迭代计算逻辑、数据存储数组以及最终的绘图代码，这应该是图像无法生成的主要原因。另外，我们需要先把二阶微分方程拆成一阶方程组，才能用梯形法则求解。

线性摆的一阶微分方程组（小角度近似）

无摩擦无驱动力的线性摆，运动方程为：
$$\ddot{\theta} = -\frac{g}{L}\theta$$
拆成两个一阶方程：

• $\frac{d\theta}{dt} = \omega$
• $\frac{d\omega}{dt} = -\omega_0^2 \theta$，其中$\omega_0 = \sqrt{g/L}$，这里我们取$g=9.8\mathrm{m/s^{2}$，摆长$L=1\mathrm{m}$，所以$\omega_0}2=9.8$

梯形法则的迭代逻辑

对于一阶方程组$\mathbf{y}' = \mathbf{f}(t, \mathbf{y})$，梯形法则的更新公式需要解线性方程组：
$$\mathbf{y}{n+1} = \mathbf{y}n + \frac{h}{2}\left[\mathbf{f}(t_n, \mathbf{y}n) + \mathbf{f}(t{n+1}, \mathbf{y}{n+1})\right]$$
展开后针对我们的方程组，可以推导出$\theta{n+1}$和$\omega_{n+1}$的简化求解公式，避免直接解矩阵。

完整可运行代码 #

import matplotlib.pylab as plt
import math

# 初始化参数
theta = 0.2  # 初始角度（弧度）
omega = 0.0  # 初始角速度
t = 0.0      # 初始时间
h = 0.01     # 时间步长
t_total = 10 # 总模拟时长
omega0_sq = 9.8  # (g/L)，g=9.8, L=1

# 存储数据的数组
t_list = [t]
theta_list = [theta]
omega_list = [omega]

# 梯形法则迭代计算
while t < t_total:
    # 计算梯形法则的中间项
    theta_temp = theta + h * omega / 2
    omega_temp = omega - h * omega0_sq * theta / 2
    
    # 求解下一个时间步的theta和omega
    denominator = 1 + (h**2 * omega0_sq) / 4
    theta_next = (theta_temp - (h**2 * omega0_sq * theta_temp) / 4 + h * omega_temp / 2) / denominator
    omega_next = (omega_temp - h * omega0_sq * theta_temp) / denominator
    
    # 更新参数
    theta = theta_next
    omega = omega_next
    t += h
    
    # 存储数据
    t_list.append(t)
    theta_list.append(theta)
    omega_list.append(omega)

# 绘制图像
plt.figure(figsize=(12, 6))

# 角度θ随时间t的变化
plt.subplot(1, 2, 1)
plt.plot(t_list, theta_list, label=r'$\theta(t)$')
plt.xlabel('时间 t (s)')
plt.ylabel('角度 θ (rad)')
plt.title('线性摆角度随时间变化')
plt.grid(True)
plt.legend()

# 角速度ω随时间t的变化
plt.subplot(1, 2, 2)
plt.plot(t_list, omega_list, label=r'$\omega(t)$', color='orange')
plt.xlabel('时间 t (s)')
plt.ylabel('角速度 ω (rad/s)')
plt.title('线性摆角速度随时间变化')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

非线性摆的修改方法 #

如果要模拟非线性摆（不做小角度近似），只需要把角速度的更新公式里的$\theta$换成$\sin\theta$即可，对应的迭代部分修改为：

# 非线性摆的梯形法则迭代（替换原迭代块）
while t < t_total:
    theta_temp = theta + h * omega / 2
    omega_temp = omega - h * omega0_sq * math.sin(theta) / 2
    
    # 非线性情况下用一次迭代近似求解（也可以用更精确的牛顿迭代）
    theta_next = theta + h * (omega + omega_temp) / 2
    omega_next = omega - h * omega0_sq * (math.sin(theta) + math.sin(theta_next)) / 2
    
    theta = theta_next
    omega = omega_next
    t += h
    
    t_list.append(t)
    theta_list.append(theta)
    omega_list.append(omega)

运行上面的代码后，你就能看到线性摆的θ-t和ω-t的正弦/余弦曲线了，非线性摆的曲线会因为大角度出现明显的非线性特征（不再是纯正弦）。

内容的提问来源于stack exchange，提问作者Aran G