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

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

当时想做一个弹簧振子的 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中做物理模拟动画，常见的“手工数值积分”方案有三大痛点：

更关键的是，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{d2x}{dt2} + kx = 0 ，即 \frac{d2x}{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(t)=cos(ωt)，这就是我们熟知的简谐振动公式。

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

有了符号表达式 x(t)=cos(ωt)，下一步是把它变成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{d2x}{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