比一比三个积分器的能量

在 LabHub 中继续学习

目标

使用三种方法离散化并运行同一物理定律，让差异在图像中显现。完成本实验后，你将能够有依据地说明物理引擎为什么选择特定积分器。

为什么重要

物理引擎的外观是碰撞和摩擦，但底层每一步都在运行的是积分器。而积分方法不只是产生误差，它还会改变系统的性质。使用显式欧拉法模拟弹簧时，能量每一步都会精确增加到 (1 + dt²) 倍，在几秒内发散。不是力计算错误，而是顺序错误。

只交换两行顺序的半隐式欧拉法，在相同成本下却很稳定。这正是游戏物理引擎使用它的原因；亲自提取数字并绘制图表后，这一事实便不会忘记。因此，本实验将结果保存为 CSV，并叠加绘制到一张 PNG 中。

步骤

在 /root/integrator 中放置工具箱。
/root/integrator/out/euler.csv——显式欧拉法。
/root/integrator/out/semi.csv——半隐式欧拉法。
/root/integrator/out/verlet.csv——速度 Verlet。
/root/integrator/out/energy.png——将三条能量曲线绘制在一张图中。
/root/integrator/out/06-dt.txt——改变 dt 后的增长率。
/root/integrator/out/damped.csv 和 damped.png——以力的形式加入阻尼。

参考

使用 nohup python3 -m http.server 8080 -d /root/integrator/out & 提供图像，并在网页预览中打开 http://localhost:8080/energy.png。
CSV 包含一行表头和 401 行数据。不要遗漏第 0 步（初始状态）。
常见错误一：将显式欧拉法拆成两行编写，却在速度计算中使用了新 x。此时它就变成了半隐式欧拉法，两个 CSV 会完全相同。
常见错误二：在 Verlet 的循环前没有先计算初始加速度。这样第一步就会错误。

放置绘图工具箱

将示例中的 /root/integrator/gfxlib.py 原样保存，并使用 /root/integrator/check.py 绘制测试图案，生成 /root/integrator/out/00-check.png。图案是在 64x64 黑色背景上，从 (0,0) 到 (63,63) 绘制白色（255,255,255）对角线，然后在其上从 (0,32) 到 (63,32) 绘制红色（255,0,0）水平线。

从本实验开始，不再重新实现 PNG 编码器。我们将第一个实验中手工创建的相同代码作为工具提供——因为这里要学习的不是文件格式。

实验 Pod 没有卷，因此前一个实验创建的文件不会保留。所以每个实验都从重新放置工具箱开始。

创建 Canvas(w, h, bg)，使用 line(x0, y0, x1, y1, rgb) 绘制两条线，再用 write_png(path) 保存。必须后画水平线，这样交点 (32,32) 才会是红色。

本实验用它绘制能量曲线。只看数字很难看出发散。

显式欧拉法

使用 /root/integrator/euler.py 模拟弹簧（质量 1、刚度 1、加速度 a = -x）：从 x=1.0、v=0.0 开始，以 dt=0.05 运行 400 步，并在 /root/integrator/out/euler.csv 中写入 step,x,v,energy 表头和 401 行数据。能量为 0.5*v*v + 0.5*x*x。

显式欧拉法只使用单步开始时的值。可以理解为同时更新位置和速度。

x, v = x + dt * v, v - dt * x    # 오른쪽은 둘 다 옛 값

写成一行就不容易出错。拆成两行时，很容易在 v 的计算中使用新 x，那样就会变成半隐式欧拉法。

还必须将第 0 步（初始状态）写成一行，才能得到 401 行。最终能量约为初始值的 2.7 倍才是正常结果——这正是本步骤要观察的现象。

半隐式欧拉法

使用 /root/integrator/semi.py，在相同条件下按先更新速度，再用新速度更新位置的方式运行，并创建 /root/integrator/out/semi.csv。格式与前一步相同。

只改变顺序。

v = v - dt * x        # 속도 먼저
x = x + dt * v        # 방금 구한 새 속도를 쓴다

可以通过第一步的值区分两种方法。显式欧拉法中的 x 仍为 1.0，而半隐式方法会得到 1 - dt*dt = 0.9975。这一处就是两种方法的全部差异。

能量应该只在一个小范围内上下波动，而不会持续增加。计算成本与显式欧拉法完全相同，结果却截然不同。

速度 Verlet

使用 /root/integrator/verlet.py，在相同条件下采用速度 Verlet 运行，并创建 /root/integrator/out/verlet.csv。一步为 x += v*dt + 0.5*a*dt*dt、a_new = -x、v += 0.5*(a + a_new)*dt，然后将 a 替换为 a_new，进入下一步。

每一步会使用加速度两次。关键是对步骤开始时的 a 和移动位置后的 a_new 求平均，再更新速度。

进入循环前，必须先用 a = -x 计算初始加速度。遗漏它会导致第一步错误。

第一步的 x 应为 1 - 0.5*0.0025 = 0.99875。三种方法从第一步开始就产生不同值，因此可以在此确认。

能量应几乎保持不变——即使 400 步后，也应位于初始值的千分之一以内。

将三条曲线绘制在一张图中

使用 /root/integrator/plot.py，将三个 CSV 的能量绘制到 /root/integrator/out/energy.png（256x256，黑色背景）。第 i 步的点为 px = int(i*255/400)、py = 255 - int(E/2.0*255)（超出 0~255 时裁剪），颜色分别为：欧拉法 (255,80,80)、半隐式 (80,255,80)、Verlet (80,120,255)。绘制顺序也按此顺序。

y 轴将能量 0 放在屏幕最下方（255），将能量 2.0 放在最上方（0）。初始能量 0.5 会绘制在 255 - int(0.5/2*255) = 192 附近。

无需连线，只绘制点也可以。401 个点会落在 256 列中，自然看起来连续。

可以使用 open(path).read().splitlines() 读取 CSV，跳过第一行（表头）后按逗号拆分。

只有红色曲线应向上急升，绿色和蓝色应贴近下方。这张图就是本实验的结论。

改变 dt 后会恶化多少

使用显式欧拉法，在相同模拟时间 T=20 下，分别以 dt=0.01（2000 步）、dt=0.05（400 步）、dt=0.2（100 步）运行，并在 /root/integrator/out/06-dt.txt 中以 dt001=、dt005=、dt02= 三行写入最终能量是初始值的多少倍（保留小数点后六位）。

步数为 T / dt。三种情况都模拟相同的 20 秒，但结果完全不同。

对于此系统存在闭式公式。每一步能量都精确变为 (1 + dt²) 倍，因此 n 步后为 (1 + dt²)^n 倍。请确认直接运行的值是否与此公式一致——一致就是实现正确的有力证据。

dt 只增加了 4 倍，增长率却增加数十倍。因此，物理引擎会选择多运行几步，而不是增大 dt。

以力的形式加入阻尼

在半隐式欧拉法中以力的形式加入阻尼（a = -x - 0.3*v，其中 v 为更新前的值），以 dt=0.05 运行 400 步，并创建 /root/integrator/out/damped.csv（格式与前面相同）和 /root/integrator/out/damped.png（能量曲线，坐标转换与第 5 步相同，颜色为 (255,200,80)）。

一步如下。

a = -x - 0.3 * v      # v 는 아직 갱신하기 전 값
v = v + dt * a
x = x + dt * v

能量必须持续下降，一次也不能增加。400 步后应降到初始值的 1% 以下。

请与每一步将速度乘以 0.99 的常见方式比较。该方式简单，但阻尼量会随 dt 改变，导致帧率不同的机器上物体运动不同。以力的形式加入时，即使 dt 不同，也能模拟相同的物理现象。