探索迷人的双摆模拟

简介

双摆是物理学和数学中的一个经典问题。它涉及两个相互连接的摆，第二个摆的运动受第一个摆的运动影响。在本实验中，我们将使用 Python 和 Matplotlib 来模拟双摆的运动。

虚拟机使用提示

虚拟机启动完成后，点击左上角切换到笔记本标签页，以访问 Jupyter Notebook 进行练习。

有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作验证无法自动化。

如果你在学习过程中遇到问题，随时向 Labby 提问。课程结束后提供反馈，我们会及时为你解决问题。

导入所需库

我们将首先导入模拟所需的库。

from collections import deque
import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin
import matplotlib.animation as animation

设置参数

接下来，我们将为模拟定义参数。这些参数包括重力加速度、每个摆的长度和质量，以及模拟的时间间隔。

G = 9.8  ## 重力加速度，单位为 m/s^2
L1 = 1.0  ## 摆 1 的长度，单位为 m
L2 = 1.0  ## 摆 2 的长度，单位为 m
L = L1 + L2  ## 组合摆的最大长度
M1 = 1.0  ## 摆 1 的质量，单位为 kg
M2 = 1.0  ## 摆 2 的质量，单位为 kg
t_stop = 2.5  ## 模拟多少秒
history_len = 500  ## 显示多少个轨迹点

定义导数函数

我们将定义一个函数，用于计算系统在任何给定时间的导数。该函数接受系统的当前状态（每个摆的角度和角速度），并返回这些值的导数。

def derivs(t, state):
    dydx = np.zeros_like(state)

    dydx[0] = state[1]

    delta = state[2] - state[0]
    den1 = (M1+M2) * L1 - M2 * L1 * cos(delta) * cos(delta)
    dydx[1] = ((M2 * L1 * state[1] * state[1] * sin(delta) * cos(delta)
                + M2 * G * sin(state[2]) * cos(delta)
                + M2 * L2 * state[3] * state[3] * sin(delta)
                - (M1+M2) * G * sin(state[0]))
               / den1)

    dydx[2] = state[3]

    den2 = (L2/L1) * den1
    dydx[3] = ((- M2 * L2 * state[3] * state[3] * sin(delta) * cos(delta)
                + (M1+M2) * G * sin(state[0]) * cos(delta)
                - (M1+M2) * L1 * state[1] * state[1] * sin(delta)
                - (M1+M2) * G * sin(state[2]))
               / den2)

    return dydx

设置初始条件并求解常微分方程

现在我们将为模拟设置初始条件。我们将设定每个摆的初始角度和角速度，以及模拟的时间间隔。然后，我们将使用欧拉方法求解常微分方程（ODE），以得到每个摆在每个时间步的位置和速度。

## 创建一个从 0 到 t_stop 的时间数组，采样间隔为 0.02 秒
dt = 0.01
t = np.arange(0, t_stop, dt)

## th1 和 th2 是初始角度（度）
## w10 和 w20 是初始角速度（度每秒）
th1 = 120.0
w1 = 0.0
th2 = -10.0
w2 = 0.0

## 初始状态
state = np.radians([th1, w1, th2, w2])

## 使用欧拉方法求解常微分方程
y = np.empty((len(t), 4))
y[0] = state
for i in range(1, len(t)):
    y[i] = y[i - 1] + derivs(t[i - 1], y[i - 1]) * dt

计算每个摆的位置

现在我们将使用每个摆在每个时间步的位置和速度来计算每个摆的 x 和 y 坐标。

x1 = L1*sin(y[:, 0])
y1 = -L1*cos(y[:, 0])

x2 = L2*sin(y[:, 2]) + x1
y2 = -L2*cos(y[:, 2]) + y1

设置绘图

现在我们将为模拟设置绘图。我们将创建一个图形，其 x 轴和 y 轴的范围等于摆的最大长度，将纵横比设置为相等，并添加网格。

fig = plt.figure(figsize=(5, 4))
ax = fig.add_subplot(autoscale_on=False, xlim=(-L, L), ylim=(-L, 1.))
ax.set_aspect('equal')
ax.grid()

定义动画函数

我们将定义一个函数来为双摆的运动制作动画。这个函数将接受当前的时间步，并使用它来更新每个摆的位置。它还将更新第二个摆的轨迹。

line, = ax.plot([], [], 'o-', lw=2)
trace, = ax.plot([], [], '.-', lw=1, ms=2)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
history_x, history_y = deque(maxlen=history_len), deque(maxlen=history_len)

def animate(i):
    thisx = [0, x1[i], x2[i]]
    thisy = [0, y1[i], y2[i]]

    if i == 0:
        history_x.clear()
        history_y.clear()

    history_x.appendleft(thisx[2])
    history_y.appendleft(thisy[2])

    line.set_data(thisx, thisy)
    trace.set_data(history_x, history_y)
    time_text.set_text(time_template % (i*dt))
    return line, trace, time_text

创建动画

现在我们将使用 Matplotlib 中的 FuncAnimation 函数来创建动画。

ani = animation.FuncAnimation(
    fig, animate, len(y), interval=dt*1000, blit=True)
plt.show()

总结

在这个实验中，我们使用 Python 和 Matplotlib 来模拟双摆的运动。我们为模拟设置了参数，定义了一个函数来计算系统的导数，使用欧拉方法对常微分方程进行积分，计算了每个时间步每个摆的位置，并创建了双摆运动的动画。