做物理模拟动画时我遇到过一个坑。当时想做一个弹簧振子的 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{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, realTrue) # 时间 x sp.Function(x)(t) # 位移函数 x(t) m, k sp.symbols(m k, positiveTrue) # 质量和劲度系数 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, icsinitial_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) # 代入具体参数值k4, m1 → ω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, modulesnumpy) # 现在 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.0lambdify 的三个关键优势速度快底层转成 NumPy 运算支持向量化比 SymPy 的.subs()快几个数量级无缝衔接返回的就是标准 Python 可调用对象直接用在任何需要函数的地方支持数组输入可以一次计算整段时间的所有位置便于做轨迹预览2.3 扩展到阻尼振动和单摆对于更复杂的微分方程dsolvelambdify的组合依然好用欠阻尼振动$ \frac{d2x}{dt2} 2\beta\frac{dx}{dt} \omega_0^2 x 0 $beta, omega0 sp.symbols(beta omega0, positiveTrue) # 欠阻尼方程 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(fA {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, realTrue) m_sym, k_sym sp.symbols(m k, positiveTrue) # 弹簧振子方程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]) # 代入参数m1, k4 → ω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], length8).shift(DOWN * 1.5) self.play(Create(number_line)) fixed_point number_line.number_to_point(-3) UP * 0.8 ball Circle(radius0.2, colorRED, fill_opacity0.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_time5, rate_funclinear, ) 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 动画ValueTrackeradd_updater每帧把虚拟时间t喂给函数更新物体位