模拟的完整视频

Manim

Manim（Mathematical Animation Engine）是一个开源的动画制作引擎，主要用于制作数学、科学和教育类的视频动画。它是用Python编写的，可以通过编程控制动画的每一个细节，从而创作出非常精确和复杂的数学可视化效果。

Manim最初是由Grant Sanderson（数学YouTuber 3Blue1Brown的创始人）开发的，用于制作3Blue1Brown视频中的数学动画。Manim提供了强大的工具，可以通过代码来创建几何图形、函数图像、动态转换等。

Manim 教学网址 Manim Community Edition

Manim

双摆

首先构建双摆的数学机制，其动力学方程由两个摆的运动方程组成，可以通过拉格朗日方程推导出。假设摆1的质量为 $m_1$ ,长度为 $L_1$ ,角度为 $\theta_1$ ,摆2的质量为 $m_2$ ,长度为 $L_2$ ,角度为 $\theta_2$ ,重力加速度为 $g$ ,系统的方程可以表示为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
from manim import *
import numpy as np
from scipy.integrate import solve_ivp

# ===================================================================
# 1. 物理参数和方程求解
# ===================================================================

# --- 物理常量 ---
g = 9.8067  # 重力加速度

# --- 摆的属性 ---
L1 = 1.6  # 杆1长度
L2 = 1.6  # 杆2长度
m1 = 1.0  # 摆球1质量
m2 = 1.0  # 摆球2质量
t_max = 60  # 模拟总时长 (秒)
dt = 1/1000   # 时间步长

# --- 初始条件 ---
# 角度 (弧度) 和角速度 (弧度/秒)
# 角度与竖直向下方向逆时针方向夹角为正角度
theta1_0 = PI  + 0.01
theta2_0 = PI  
omega1_0 = 0.0
omega2_0 = 0.0

def rhs_using_linear_system(t, y):
    """
    右端函数 f(t,y)，实现与参考代码同一套数学逻辑，
    但使用状态顺序 [theta1, omega1, theta2, omega2]。
    返回 [θ1', ω1', θ2', ω2'] = [ω1, dω1, ω2, dω2]
    """
    th1, om1, th2, om2 = y
    A = np.zeros((2, 2), dtype=np.float64)
    b = np.zeros(2, dtype=np.float64)

    c12 = np.cos(th1 - th2)
    s12 = np.sin(th1 - th2)

    A[0, 0] = L1 * (m1 + m2)
    A[0, 1] = m2 * L2 * c12
    A[1, 0] = L2 * c12
    A[1, 1] = L1

    b[0] = (m1 + m2) * g * np.cos(th1) - m2 * L2 * s12 * (om2 ** 2)
    b[1] = L1 * s12 * (om1 ** 2) + g * np.cos(th2)

    dom1, dom2 = np.linalg.solve(A, b)

    # 组装导数，顺序对应 y = [th1, om1, th2, om2]
    return np.array([om1, dom1, om2, dom2], dtype=np.float64)

def rk4_integrate(f, t0, t1, h, y0):
    """
    固定步长四阶RK（与参考实现等价）
    返回与原动画兼容的 (t_eval, Y)：
      - t_eval: np.arange(t0, t1, h)
      - Y.shape = (len(t_eval), len(y0))
    """
    # 为了和你动画里的索引严格对齐，这里不包含 t1 本身
    t_eval = np.arange(t0, t1, h, dtype=np.float64)
    n = t_eval.size
    y = np.zeros((n, len(y0)), dtype=np.float64)
    y[0, :] = y0

    for i in range(n - 1):
        t = t_eval[i]
        yi = y[i, :]

        k1 = f(t, yi)
        k2 = f(t + h/2.0, yi + (h/2.0) * k1)
        k3 = f(t + h/2.0, yi + (h/2.0) * k2)
        k4 = f(t + h,       yi + h * k3)

        y[i+1, :] = yi + (h/6.0) * (k1 + 2.0*(k2 + k3) + k4)

    return t_eval, y

# ===================================================================
# 2. Manim 动画设置（将图元与更新逻辑收敛到 Scene 内，避免全局副作用）
# ===================================================================

class DoublePendulumAnimation(Scene):
    def construct(self):
        t_span = [0.0, t_max]
        t_eval = np.arange(0.0, t_max, dt, dtype=np.float64)
        y0 = np.array([-theta1_0+PI/2, omega1_0, -theta2_0+PI/2, omega2_0], dtype=np.float64)

        _t, _Y = rk4_integrate(rhs_using_linear_system, t_span[0], t_span[1], dt, y0)
        theta1_sol = _Y[:, 0]
        theta2_sol = _Y[:, 2]

        origin = Dot(point=DOWN * 0.5)
        title = Tex("Double Pendulum Simulation").to_edge(UP)

        initial_theta1 = theta1_sol[0]
        initial_theta2 = theta2_sol[0]

        p1_initial_pos = origin.get_center() + np.array([
            L1 * np.cos(initial_theta1),
            -L1 * np.sin(initial_theta1),
            0
        ])

        p2_initial_pos = p1_initial_pos + np.array([
            L2 * np.cos(initial_theta2),
            -L2 * np.sin(initial_theta2),
            0
        ])

        # --- 创建摆的组件 (Line 和 Dot) ---
        bob1 = Dot(point=p1_initial_pos, color=BLUE)
        bob2 = Dot(point=p2_initial_pos, color=GREEN)
        rod1 = Line(origin.get_center(), bob1.get_center(), stroke_width=3, color=GRAY)
        rod2 = Line(bob1.get_center(), bob2.get_center(), stroke_width=3, color=GRAY)

        pendulum_group = VGroup(rod1, rod2, bob1, bob2)

        # --- 追踪路径 ---
        trace = TracedPath(bob2.get_center, stroke_width=3, stroke_color=GREEN, stroke_opacity=0.8)

        # --- 时间驱动器与更新函数 ---
        time = ValueTracker(0)

        def update_pendulum(_):
            current_t = time.get_value()
            max_index = len(t_eval) - 1
            index = min(int(current_t / dt), max_index)

            theta1 = theta1_sol[index]
            theta2 = theta2_sol[index]

            p1_x = origin.get_x() + L1 * np.cos(theta1)
            p1_y = origin.get_y() - L1 * np.sin(theta1)
            p2_x = p1_x + L2 * np.cos(theta2)
            p2_y = p1_y - L2 * np.sin(theta2)

            bob1.move_to([p1_x, p1_y, 0])
            bob2.move_to([p2_x, p2_y, 0])
            rod1.put_start_and_end_on(origin.get_center(), bob1.get_center())
            rod2.put_start_and_end_on(bob1.get_center(), bob2.get_center())

        pendulum_group.add_updater(update_pendulum)

        # --- 播放动画 ---
        self.add(title, origin, trace, pendulum_group)
        self.play(time.animate.set_value(t_max), run_time=t_max, rate_func=linear)
        self.wait()

三摆

为处理更多元素的模拟，采用新的思路，用拉格朗日力学的矩阵形式建模三摆（杆无质量、端点为点质量），把动力学写成

在每个时间步先解线性方程得到 $\ddot{\boldsymbol{\theta}}=M^{-1}[-\mathbf{C}-\mathbf{G}]$ 再用 固定步长四阶 Runge–Kutta（RK4） 做时间积分。

各项是怎样构造的

质量矩阵 $M(\theta)$（点质量、无杆质量模型）
对角元 非对角元 （见 mass_matrix(theta)）

重力项 $\mathbf G(\theta)$
取势能 因而在代码约定的坐标与号约定下得到 （见 gravity_vector(theta)）

科氏/离心项 $\mathbf C(\theta,\dot\theta)$
采用数值 Christoffel 形式：先对 $M(\theta)$ 做有限差分近似 $\partial M/\partial \theta_k$（步长 $\varepsilon=10^{-7}$ ），再按

组装（见 coriolis_vector(...)

状态方程与求解
以状态 $y=[\theta_1,\omega_1,\theta_2,\omega_2,\theta_3,\omega_3]$，在 triple_rhs 中先算 $M,C,G$，解 $\;M\,\ddot{\boldsymbol\theta}=-(\mathbf C+\mathbf G)$，输出 $[\,\omega_1,\ddot\theta_1,\omega_2,\ddot\theta_2,\omega_3,\ddot\theta_3\,]$。

时间推进
用的固定步长 RK4 积分器 rk4_integrate(...)，步长 dt=1e-3，全程无自适应步长与误差控制。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
from manim import *
import numpy as np

# ===================================================================
# 1. 三摆动力学与数值积分（仿 double_pendulum_2.py 结构）
# ===================================================================

# --- 物理常量与参数 ---
g = 9.8067

L1 = 1.0
L2 = 1.0
L3 = 1.0

m1 = 1.0
m2 = 1.0
m3 = 1.0

t_max = 60.0
dt = 1.0 / 1000.0

# --- 初始条件 ---
# 角度 (弧度) 和角速度 (弧度/秒)
# 角度与竖直向下方向逆时针方向夹角为正角度
theta1_0 = PI 
theta2_0 = PI
theta3_0 = PI + 0.01
omega1_0 = 0.0
omega2_0 = 0.0
omega3_0 = 0.0


# -------------------------------------------------------------------
# 质量矩阵 M(θ)，点质量、无杆质量模型
# M_ij = sum_{k=max(i,j)}^3 m_k * l_i * l_j * cos(θ_i - θ_j), i ≠ j
# M_ii = sum_{k=i}^3 m_k * l_i^2
# -------------------------------------------------------------------
def mass_matrix(theta):
    th1, th2, th3 = float(theta[0]), float(theta[1]), float(theta[2])
    l = [L1, L2, L3]
    m = [m1, m2, m3]

    def tail_mass(i):
        return sum(m[i:])

    M = np.zeros((3, 3), dtype=np.float64)
    # Diagonal terms
    M[0, 0] = tail_mass(0) * (l[0] ** 2)
    M[1, 1] = tail_mass(1) * (l[1] ** 2)
    M[2, 2] = tail_mass(2) * (l[2] ** 2)

    # Off-diagonal terms
    th = [th1, th2, th3]
    for i in range(3):
        for j in range(i + 1, 3):
            cij = sum(m[j:]) * l[i] * l[j] * np.cos(th[i] - th[j])
            M[i, j] = cij
            M[j, i] = cij

    return M


# -------------------------------------------------------------------
# 重力项 G(θ)
# G_i = (sum_{k=i}^3 m_k) * g * l_i * sin(θ_i)
# 坐标系与动画一致：x=Σ l_i cosθ_i, y=-Σ l_i sinθ_i
# -------------------------------------------------------------------
def gravity_vector(theta):
    th1, th2, th3 = float(theta[0]), float(theta[1]), float(theta[2])
    l = [L1, L2, L3]
    m = [m1, m2, m3]

    def tail_mass(i):
        return sum(m[i:])

    G = np.zeros(3, dtype=np.float64)
    # 向下的重力：U = \sum m g y = -\sum m g l sin(θ)，故 G_i = ∂U/∂θ_i = -\sum m g l cos(θ_i)
    G[0] = -tail_mass(0) * g * l[0] * np.cos(th1)
    G[1] = -tail_mass(1) * g * l[1] * np.cos(th2)
    G[2] = -tail_mass(2) * g * l[2] * np.cos(th3)
    return G


# -------------------------------------------------------------------
# 科氏/离心项 C(θ, θdot)θdot = Γ(θ, θdot)
# 使用数值 Christoffel 形式：
#   C_i = Σ_j Σ_k 0.5*(∂M_ij/∂θ_k + ∂M_ik/∂θ_j - ∂M_jk/∂θ_i) θ̇_j θ̇_k
# 通过数值微分 dM/dθ 计算
# -------------------------------------------------------------------
def coriolis_vector(theta, theta_dot, eps=1e-7):
    th = np.asarray(theta, dtype=np.float64)
    dth = np.asarray(theta_dot, dtype=np.float64)

    M0 = mass_matrix(th)
    dM = [None, None, None]
    for k in range(3):
        th_perturb = th.copy()
        th_perturb[k] += eps
        Mk = mass_matrix(th_perturb)
        dM[k] = (Mk - M0) / eps

    C = np.zeros(3, dtype=np.float64)
    for i in range(3):
        s = 0.0
        for j in range(3):
            for k in range(3):
                term = 0.5 * (dM[k][i, j] + dM[j][i, k] - dM[i][j, k])
                s += term * dth[j] * dth[k]
        C[i] = s
    return C


# 右端函数：y = [θ1, ω1, θ2, ω2, θ3, ω3] -> y' = [ω1, α1, ω2, α2, ω3, α3]
def triple_rhs(t, y):
    th1, om1, th2, om2, th3, om3 = y
    theta = np.array([th1, th2, th3], dtype=np.float64)
    theta_dot = np.array([om1, om2, om3], dtype=np.float64)

    M = mass_matrix(theta)
    C = coriolis_vector(theta, theta_dot)
    Gv = gravity_vector(theta)

    # M * theta_ddot + C + G = 0
    rhs = -(C + Gv)
    theta_ddot = np.linalg.solve(M, rhs)

    return np.array([om1, theta_ddot[0], om2, theta_ddot[1], om3, theta_ddot[2]], dtype=np.float64)


# 固定步长 RK4
def rk4_integrate(f, t0, t1, h, y0):
    t_eval = np.arange(t0, t1, h, dtype=np.float64)
    n = t_eval.size
    y = np.zeros((n, len(y0)), dtype=np.float64)
    y[0, :] = y0

    for i in range(n - 1):
        t = t_eval[i]
        yi = y[i, :]
        k1 = f(t, yi)
        k2 = f(t + h / 2.0, yi + (h / 2.0) * k1)
        k3 = f(t + h / 2.0, yi + (h / 2.0) * k2)
        k4 = f(t + h, yi + h * k3)
        y[i + 1, :] = yi + (h / 6.0) * (k1 + 2.0 * (k2 + k3) + k4)

    return t_eval, y


# ===================================================================
# 2. Manim 动画
# ===================================================================

class TriplePendulumAnimation(Scene):
    def construct(self):
        t_span = [0.0, t_max]
        t_eval = np.arange(0.0, t_max, dt, dtype=np.float64)
        y0 = np.array([-theta1_0+PI/2, omega1_0, -theta2_0+PI/2, omega2_0, -theta3_0+PI/2, omega3_0], dtype=np.float64)

        _t, _Y = rk4_integrate(triple_rhs, t_span[0], t_span[1], dt, y0)
        theta1_sol = _Y[:, 0]
        theta2_sol = _Y[:, 2]
        theta3_sol = _Y[:, 4]

        # --- 基础对象 ---
        origin = Dot(point=DOWN * 0.5)
        title = Tex("Triple Pendulum Simulation").to_edge(UP)

        th1_0 = theta1_sol[0]
        th2_0 = theta2_sol[0]
        th3_0 = theta3_sol[0]

        p1_initial = origin.get_center() + np.array([
            L1 * np.cos(th1_0),
            -L1 * np.sin(th1_0),
            0.0,
        ])
        p2_initial = p1_initial + np.array([
            L2 * np.cos(th2_0),
            -L2 * np.sin(th2_0),
            0.0,
        ])
        p3_initial = p2_initial + np.array([
            L3 * np.cos(th3_0),
            -L3 * np.sin(th3_0),
            0.0,
        ])

        bob1 = Dot(point=p1_initial, color=BLUE)
        bob2 = Dot(point=p2_initial, color=GREEN)
        bob3 = Dot(point=p3_initial, color=RED)

        rod1 = Line(origin.get_center(), bob1.get_center(), stroke_width=3, color=GRAY)
        rod2 = Line(bob1.get_center(), bob2.get_center(), stroke_width=3, color=GRAY)
        rod3 = Line(bob2.get_center(), bob3.get_center(), stroke_width=3, color=GRAY)

        group = VGroup(rod1, rod2, rod3, bob1, bob2, bob3)

        trace = TracedPath(bob3.get_center, stroke_width=3, stroke_color=RED, stroke_opacity=0.8)

        time = ValueTracker(0.0)

        def update_system(_):
            current_t = time.get_value()
            max_index = len(t_eval) - 1
            index = min(int(current_t / dt), max_index)

            th1 = theta1_sol[index]
            th2 = theta2_sol[index]
            th3 = theta3_sol[index]

            p1_x = origin.get_x() + L1 * np.cos(th1)
            p1_y = origin.get_y() - L1 * np.sin(th1)
            p2_x = p1_x + L2 * np.cos(th2)
            p2_y = p1_y - L2 * np.sin(th2)
            p3_x = p2_x + L3 * np.cos(th3)
            p3_y = p2_y - L3 * np.sin(th3)

            bob1.move_to([p1_x, p1_y, 0])
            bob2.move_to([p2_x, p2_y, 0])
            bob3.move_to([p3_x, p3_y, 0])

            rod1.put_start_and_end_on(origin.get_center(), bob1.get_center())
            rod2.put_start_and_end_on(bob1.get_center(), bob2.get_center())
            rod3.put_start_and_end_on(bob2.get_center(), bob3.get_center())

        group.add_updater(update_system)

        self.add(title, origin, trace, group)
        self.play(time.animate.set_value(t_max), run_time=t_max, rate_func=linear)
        self.wait()

后记

若要推广到更多维度的摆，采用三摆的方案即可