Manim物理模拟：别自己写欧拉了！

合集 - manim动画(43)

1.【manim动画教程】-- 安装2023-03-282.【manim动画教程】-- 基本图形2023-03-293.【manim动画教程】-- 坐标系2023-04-104.【manim动画教程】-- 文本样式2023-04-075.【manim动画教程】-- 文字和公式2023-04-046.【manim动画教程】-- 图形样式2023-03-317.【manim动画教程】--相机2023-04-198.【manim动画教程】--高级动画效果2023-04-149.【manim动画教程】--常用动画效果2023-04-1210.【manim】之滚动字幕2022-12-0611.【manim】之圆规动画2023-01-3112.【manim】之目录动画2023-02-2813.Manim：动画制作背后的魔法2025-05-1814.manim 动画效果总结2025-05-2015.manim变换效果总结2025-05-2116.Manim动画渲染：从代码到屏幕的幕后故事2025-05-2217.轻松掌握Manim的.animate语法：让动画编程更简单2025-05-2518.Manim中三种函数图像类的比较2025-06-0519.掌握ChangeSpeed类：让数学动画速度随心而动！2025-06-1520.把数学对象画出来：Manim Mobject类库速查手册2025-08-1621.「亲手画出动态数学」：让数学可视化触手可及的Manim入门课2025-11-0722.让你的动画“活”过来：Manim 节奏控制指南 (Rate Functions)2025-11-2323.Manim进阶：用背景图片让你的数学视频脱颖而出2025-11-2624.Manim v0.19.1 发布啦！三大新特性让动画制作更丝滑2025-12-0225.🚀 Manim CE v0.20.0 发布：动画构建更丝滑，随机性终于“可控”了！02-2126.ManimCE v0.20.1 发布：LaTeX 渲染修复与动画稳定性提升03-0227.探索视觉的边界：用 Manim 重现有趣的知觉错觉03-0728.用Manim实现动态交点计算--从一个动点问题说起04-1629.Manim如何在数学公式中完美显示中文？04-1930.如何灵活设置公式中各个部分的颜色？04-2231.如何让多个动画“齐步走”？04-2632.告别手动计算，SymPy 初识与 Manim 联动05-1633.让数学公式自动推导05-1734.切线的魔法：用 SymPy 和 Manim 轻松搞定导数动画05-1935.填充与积累：积分与面积的可视化05-2136.轨迹的蓝图：方程求解与交点计算05-23

37.Manim物理模拟：别自己写欧拉了！05-25

38.用 SymPy 解决 Manim 曲线绘制速度不均的问题06-0139.一次函数图像工厂：用 SymPy 自动生成 y=kx+b 对比动画06-0340.用SymPy自动计算抛物线求根、判别式与顶点06-0541.用SymPy自动求解三角形构造与全等条件验证06-0842.用SymPy自动求解追及问题的方程06-1043.用SymPy自动因式分解：从面积拼图到代数恒等式06-13

收起

做物理模拟动画时，我遇到过一个坑。

当时想做一个弹簧振子的 Manim 动画：一个小球连接在弹簧上，在平衡位置附近往复振动。

我一开始的思路是——手动写欧拉法迭代。

# 当时写的“玩具级”数值积分代码
x = 1.0   # 初始位移
v = 0.0   # 初始速度
dt = 0.02 # 时间步长
k = 2.0   # 弹簧劲度系数
m = 1.0   # 质量

positions = []
for i in range(300):
    a = -k/m * x          # 加速度
    v += a * dt           # 更新速度
    x += v * dt           # 更新位置
    positions.append(x)

这段代码跑起来之后，小球确实动起来了。但看了几秒之后——小球越振幅度越大，能量明显不守恒。

欧拉法的数值误差在逐帧累积，像个隐形的外力在不断推着小球。

我当然可以换龙格-库塔法，但那意味着更复杂的代码、更长的调试时间。我只是想做一个弹簧振子动画，为什么要从数值分析开始写起？

直到我开始用 SymPy 的 dsolve，才发现原来我根本不需要自己写数值积分。

1. 痛点场景还原

在 Manim 中做物理模拟动画，常见的“手工数值积分”方案有三大痛点：

痛点	具体表现
数值误差累积	欧拉法能量不守恒，弹簧振子越振越离谱
步长调参折磨	dt 设太小渲染慢，设太大精度炸，反复试错
无法利用解析解	明明简谐振子有精确的 cos(ωt) 解析表达式，却硬要用数值去逼近

更关键的是，Manim 的动画是基于 ValueTracker 的时间驱动的。

我们真正需要的是一个函数 x(t)，能根据当前时刻 t 精确返回小球的位置。

这和 SymPy 的解析解简直是天作之合。

# Manim 动画的本质：需要一个时间→位置的映射函数
def get_position(t):
    # 如果这里能直接用解析解，该多好！
    return ???  

SymPy 的 dsolve 可以帮我们求出这个解析的 x(t)，然后通过 lambdify 把它转成 NumPy 函数，直接注入 ValueTracker 的更新逻辑中。全程无误差累积，因为每一帧的位置都是独立计算的。

2. SymPy 解决方案

2.1 dsolve：符号求解微分方程

SymPy 的 dsolve 可以求解常微分方程（ODE）的解析解。以弹簧振子（简谐振动）为例：

运动方程：$ m\frac{d^2x}{dt2} + kx = 0 $，即$ \frac{d^2x}{dt2} + \omega_0^2 x = 0 $，其中$ \omega_0 = \sqrt{\frac{k}{m}} $是系统的固有角频率。

import sympy as sp

# 定义符号
t = sp.Symbol('t', real=True)           # 时间
x = sp.Function('x')(t)                 # 位移函数 x(t)
m, k = sp.symbols('m k', positive=True)  # 质量和劲度系数
omega = sp.sqrt(k/m)                     # 角频率

# 定义微分方程：m * x'' + k * x = 0
ode = sp.Eq(m * sp.diff(x, t, 2) + k * x, 0)

# 代入初始条件求解
# 设初始位移 x(0) = 1，初始速度 x'(0) = 0
initial_conditions = {
    x.subs(t, 0): 1,                # x(0) = 1
    sp.diff(x, t).subs(t, 0): 0     # x'(0) = 0
}

# dsolve 求解！
solution = sp.dsolve(ode, ics=initial_conditions)
print(solution)  # 输出：x(t) = cos(√(k/m)·t)

SymPy 直接给出了解析解：$x\left(t\right)=\cos \left(\omega t\right)$，这就是我们熟知的简谐振动公式。

2.2 lambdify：符号表达式 → 高性能数值函数

有了符号表达式 $x\left(t\right)=\cos \left(\omega t\right)$，下一步是把它变成 Manim 可以每秒调用几十次的快速函数。lambdify 就是干这个的：

import numpy as np

# 提取解表达式右侧的 x(t)
x_expr = solution.rhs  # cos(sqrt(k/m)*t)

# 代入具体参数值：k=4, m=1 → ω=2
x_expr_sub = x_expr.subs({k: 4, m: 1})  # cos(2*t)

# lambdify：将 SymPy 表达式编译为 NumPy 函数
# 参数是 t，返回 NumPy 数组
x_func = sp.lambdify(t, x_expr_sub, modules='numpy')

# 现在 x_func 可以像普通 NumPy 函数一样使用
print(x_func(0))      # 1.0
print(x_func(np.pi/4)) # cos(π/2) ≈ 0.0
print(x_func(np.pi/2)) # cos(π) = -1.0

lambdify 的三个关键优势：

1. 速度快：底层转成 NumPy 运算，支持向量化，比 SymPy 的 .subs() 快几个数量级
2. 无缝衔接：返回的就是标准 Python 可调用对象，直接用在任何需要函数的地方
3. 支持数组输入：可以一次计算整段时间的所有位置，便于做轨迹预览

2.3 扩展到阻尼振动和单摆

对于更复杂的微分方程，dsolve + lambdify 的组合依然好用：

欠阻尼振动（$ \frac{d^2x}{dt2} + 2\beta\frac{dx}{dt} + \omega_0^2 x = 0 $）：

beta, omega0 = sp.symbols('beta omega0', positive=True)

# 欠阻尼方程
ode_damped = sp.Eq(sp.diff(x, t, 2) + 2*beta*sp.diff(x, t) + omega0**2*x, 0)
sol_damped = sp.dsolve(ode_damped, ics={
    x.subs(t, 0): 1,
    sp.diff(x, t).subs(t, 0): 0
})

# 代入参数后 lambdify
x_damped_expr = sol_damped.rhs.subs({beta: 0.3, omega0: 2})
x_damped_func = sp.lambdify(t, x_damped_expr, 'numpy')

小角度单摆（$ \theta'' + \frac{g}{L}\sin(\theta) = 0 $，小角度近似$ \sin(\theta) \approx \theta $）：

g_val, L_val = 9.8, 2.0
omega = sp.sqrt(g_val/L_val)  # 角频率 ω = √(g/L)

# ========== 直接写出通解形式 ==========
# 小角度近似后的方程：θ'' + ω²θ = 0
# 通解：θ(t) = A·cos(ωt) + B·sin(ωt)

A, B = sp.symbols('A B')  # 待定系数
theta_expr = A * sp.cos(omega * t) + B * sp.sin(omega * t)

# ========== 手动代入初始条件 ==========
theta_0 = sp.pi / 6          # 初始角度 30°
dtheta_expr = sp.diff(theta_expr, t)  # 角速度表达式

# 条件1：θ(0) = π/6  →  A = π/6
eq1 = sp.Eq(theta_expr.subs(t, 0), theta_0)

# 条件2：θ'(0) = 0  →  ω·B = 0  →  B = 0
eq2 = sp.Eq(dtheta_expr.subs(t, 0), 0)

# 求解待定系数
solution = sp.solve([eq1, eq2], [A, B])
print(f"A = {solution[A]}, B = {solution[B]}")

# 代入得到精确解
theta_final = theta_expr.subs(solution)
print(f"θ(t) = {theta_final}")
# 输出：θ(t) = pi*cos(√(g/L)*t)/6

# ========== lambdify 转为 NumPy 函数 ==========
theta_func = sp.lambdify(t, theta_final, 'numpy')

3. Manim 联动实战

下面是一个弹簧振子动画的核心代码示例。

核心思路是：用 SymPy 解出 x(t)，用 lambdify 转成快速函数，再通过 ValueTracker 的追踪器模式驱动小球运动。

from manim import *
import sympy as sp
import numpy as np

class SpringOscillator(Scene):
    def construct(self):
        # ========== SymPy 符号求解微分方程 ==========
        t_sym = sp.Symbol("t", real=True)
        m_sym, k_sym = sp.symbols("m k", positive=True)

        # 弹簧振子方程：m·x'' + k·x = 0  →  通解 x(t) = A·cos(ωt) + B·sin(ωt)
        A, B = sp.symbols("A B")
        omega = sp.sqrt(k_sym / m_sym)
        x_expr = A * sp.cos(omega * t_sym) + B * sp.sin(omega * t_sym)

        # 手动代入初始条件：x(0)=2, x'(0)=0
        eq1 = sp.Eq(x_expr.subs(t_sym, 0), 2)              # A = 2
        eq2 = sp.Eq(sp.diff(x_expr, t_sym).subs(t_sym, 0), 0)  # ω·B = 0
        sol_const = sp.solve([eq1, eq2], [A, B])

        # 代入参数：m=1, k=4 → ω=2 → x(t) = 2cos(2t)
        x_final = x_expr.subs(sol_const).subs({m_sym: 1, k_sym: 4})
        x_func = sp.lambdify(t_sym, x_final, "numpy")  # 编译为 NumPy 函数

        # ========== Manim 场景搭建 ==========
        number_line = NumberLine(x_range=[-3, 3, 1], length=8).shift(DOWN * 1.5)
        self.play(Create(number_line))

        fixed_point = number_line.number_to_point(-3) + UP * 0.8
        ball = Circle(radius=0.2, color=RED, fill_opacity=0.8)
        ball.move_to(number_line.number_to_point(x_func(0)) + UP * 0.8)

        self.play(DrawBorderThenFill(ball))

        # ========== ValueTracker 驱动动画（核心机制）==========
        time_tracker = ValueTracker(0)  # 虚拟时钟

        # 更新器：每帧用 SymPy 解析解计算精确位置
        ball.add_updater(lambda m: m.set_x(
            number_line.number_to_point(
                x_func(time_tracker.get_value())  # x(t) 精确值！
            )[0]
        ))

        # 播放：虚拟时间从 0 流逝到 5 秒
        self.play(
            time_tracker.animate.set_value(5),
            run_time=5,
            rate_func=linear,
        )
        self.wait(1)

4. 效果展示说明

运行上面的脚本，你会看到这样的动画流程：

观察重点：

• 小球的运动完全符合余弦规律：在两端停顿、中间最快，完美再现简谐振动
• 弹簧长度随小球位置实时变化，视觉上非常自然
• 整个过程没有任何数值积分代码——你写的只是方程和初始条件，剩下的交给 SymPy

换成阻尼振动有多容易？

只需要改两行代码：

# 阻尼振动方程
ode_damped = sp.Eq(
    sp.diff(x_sym, t_sym, 2) + 0.6*sp.diff(x_sym, t_sym) + 4*x_sym, 0
)
# dsolve 会自动处理，lambdify 之后的用法完全一样

小球的振幅会逐渐减小，直到停在平衡位置。而这一切，你不需要修改任何动画逻辑——因为 x_func 这个黑盒函数的内部实现已经被 SymPy 替换了。

5. 本期小结

核心公式：dsolve → lambdify → ValueTracker

步骤	工具	作用
1. 定义微分方程	sp.Eq(sp.diff(...))	用符号表达物理规律
2. 求解析解	sp.dsolve(ode, ics=...)	得出 x(t) 的封闭表达式
3. 编译为快函数	sp.lambdify(t, expr, 'numpy')	符号表达式 → 每秒可调用千次的 NumPy 函数
4. 驱动 Manim 动画	ValueTracker + add_updater	每帧把虚拟时间 t 喂给函数，更新物体位置