一开始的思路是——手动写欧拉法迭代。当时写的“玩具级”数值积分代码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才发现原来我根本不需要自己写数值积分。痛点场景还原在 Manim 中做物理模拟动画常见的“手工数值积分”方案有三大痛点痛点 具体表现数值误差累积 欧拉法能量不守恒弹簧振子越振越离谱步长调参折磨 dt 设太小渲染慢设太大精度炸反复试错无法利用解析解 明明简谐振子有精确的 cos(ωt) 解析表达式却硬要用数值去逼近更关键的是Manim 的动画是基于 ValueTracker 的时间驱动的。我们真正需要的是一个函数 x(t)能根据当前时刻 t 精确返回小球的位置。这和 SymPy 的解析解简直是天作之合。Manim 动画的本质需要一个时间→位置的映射函数def get_position(t):# 如果这里能直接用解析解该多好return ???SymPy 的 dsolve 可以帮我们求出这个解析的 x(t)然后通过 lambdify 把它转成 NumPy 函数直接注入 ValueTracker 的更新逻辑中。全程无误差累积因为每一帧的位置都是独立计算的。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 0ode sp.Eq(m * sp.diff(x, t, 2) k * x, 0)代入初始条件求解设初始位移 x(0) 1初始速度 x’(0) 0initial_conditions {x.subs(t, 0): 1, # x(0) 1sp.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 → ω2x_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.0print(x_func(np.pi/4)) # cos(π/2) ≈ 0.0print(x_func(np.pi/2)) # cos(π) -1.0lambdify 的三个关键优势速度快底层转成 NumPy 运算支持向量化比 SymPy 的 .subs() 快几个数量级无缝衔接返回的就是标准 Python 可调用对象直接用在任何需要函数的地方支持数组输入可以一次计算整段时间的所有位置便于做轨迹预览2.3 扩展到阻尼振动和单摆对于更复杂的微分方程dsolve lambdify 的组合依然好用欠阻尼振动$ \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) 2betasp.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})代入参数后 lambdifyx_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.0omega 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 π/6eq1 sp.Eq(theta_expr.subs(t, 0), theta_0)条件2θ’(0) 0 → ω·B 0 → B 0eq2 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 spimport numpy as npclass 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)效果展示说明运行上面的脚本你会看到这样的动画流程观察重点小球的运动完全符合余弦规律在两端停顿、中间最快完美再现简谐振动弹簧长度随小球位置实时变化视觉上非常自然整个过程没有任何数值积分代码——你写的只是方程和初始条件剩下的交给 SymPy换成阻尼振动有多容易只需要改两行代码阻尼振动方程ode_damped sp.Eq(sp.diff(x_sym, t_sym, 2) 0.6sp.diff(x_sym, t_sym) 4x_sym, 0)dsolve 会自动处理lambdify 之后的用法完全一样小球的振幅会逐渐减小直到停在平衡位置。而这一切你不需要修改任何动画逻辑——因为 x_func 这个黑盒函数的内部实现已经被 SymPy 替换了。本期小结核心公式dsolve → lambdify → ValueTracker步骤 工具 作用定义微分方程 sp.Eq(sp.diff(…)) 用符号表达物理规律求解析解 sp.dsolve(ode, ics…) 得出 x(t) 的封闭表达式编译为快函数 sp.lambdify(t, expr, ‘numpy’) 符号表达式 → 每秒可调用千次的 NumPy 函数驱动 Manim 动画 ValueTracker add_updater 每帧把虚拟时间 t 喂给函数更新物体位置