三つの積分器のエネルギーを比べる

目標

同じ物理法則を、3つの方法で離散化して回してみて、その違いが、絵として表れるようにします。このラボが終われば、物理エンジンがなぜ特定の積分器を選ぶのか、根拠を持って語れるようになります。

なぜ重要なのか

物理エンジンの見た目は、衝突と摩擦ですが、その下で毎ステップ回っているのは、積分器です。そして、積分の方法は、誤差を出すだけでなく、系の性質を変えます。陽的オイラー法でばねを回すと、エネルギーが毎ステップ、ちょうど(1 + dt²)倍ずつ増えて、数秒で発散します。力の計算が間違っているのではなく、順序が間違っているのです。

2行の順序だけを変えた半陰的オイラー法は、同じコストで安定します。ゲームの物理エンジンがこれを使う理由で、この事実は、数字を自分で取り出して、グラフで見れば、忘れません。そのため、このラボは、結果をCSVに残し、1枚のPNGに重ねて描きます。

ステップ

/root/integratorに、ツールボックスを置きます。
/root/integrator/out/euler.csvに、陽的オイラー法の結果を書きます。
/root/integrator/out/semi.csvに、半陰的オイラー法の結果を書きます。
/root/integrator/out/verlet.csvに、速度ヴェルレ法の結果を書きます。
/root/integrator/out/energy.pngに、3つのエネルギー曲線を、1枚に描きます。
/root/integrator/out/06-dt.txtに、dtを変えた増加率を書きます。
/root/integrator/out/damped.csvとdamped.pngに、減衰を力として入れた結果を書きます。

参考

絵は、nohup python3 -m http.server 8080 -d /root/integrator/out &で立てて、Webプレビューでhttp://localhost:8080/energy.pngを開きます。
CSVは、ヘッダー1行 + 401行です。0ステップ目(初期状態)を抜かさないでください。
よくあるミスの1つ目は、陽的オイラー法を2行に分けて書くうちに、新しいxを速度の計算に使ってしまうことです。その瞬間、それは半陰的オイラー法になり、2つのCSVがまったく同じになります。
よくあるミスの2つ目は、ヴェルレ法で、ループの前に初期の加速度を求めておかないことです。最初のステップが間違います。

描画のツールボックスを置く

/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)で2本の線を引いたあと、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です。

陽的オイラー法は、1ステップの開始時点の値だけを使います。位置と速度を同時に更新すると考えればよいです。

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

一行で書けば、間違えません。2行に分けて書くと、新しいxをvの計算に使ってしまいやすく、そうすると、それは半陰的オイラー法になります。

0ステップ目(初期状態)も1行として書いて、はじめて401行になります。最後のエネルギーが、最初の2.7倍ほどになれば正常です。それが、このステップで見たいことです。

半陰的オイラー法

/root/integrator/semi.pyで、同じ条件を、速度を先に更新して、その新しい速度で位置を更新する方式で回して、/root/integrator/out/semi.csvを作ってください。形式は、前のステップと同じです。

順序だけが変わります。

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

最初のステップの値で、2つの方法を区別できます。陽的オイラー法は、xがそのまま1.0ですが、半陰的は、1 - dt*dt = 0.9975になります。この一桁が、2つの方法の違いのすべてです。

エネルギーが、小さな幅で上下するだけで、増えなければ正常です。計算コストは、陽的オイラー法とまったく同じなのに、結果は、まったく違います。

速度ヴェルレ法

/root/integrator/verlet.pyで、同じ条件を速度ヴェルレ法で回して、/root/integrator/out/verlet.csvを作ってください。1ステップは、x += v*dt + 0.5*a*dt*dt、a_new = -x、v += 0.5*(a + a_new)*dtで、aをa_newに入れ替えて、次のステップへ進みます。

加速度を、ステップごとに2回使います。ステップの開始時点のaと、位置を動かしたあとのa_newを平均して、速度を更新することが核心です。

ループに入る前に、a = -xで初期の加速度を求めておく必要があります。これを抜かすと、最初のステップが間違います。

最初のステップのxは、1 - 0.5*0.0025 = 0.99875になります。3つの方法が、最初のステップから互いに異なる値を出すので、ここで確認できます。

エネルギーが、ほとんどそのまま保たれれば正常です。400ステップ後も、最初の1000分の1以内に収まっています。

3本の曲線を1枚に

/root/integrator/plot.pyで、3つの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)、ヴェルレ法が(80,120,255)です。描く順序も、この順序です。

y軸は、エネルギー0を画面の一番下(255)に、エネルギー2.0を一番上(0)に置くものです。最初のエネルギー0.5は、255 - int(0.5/2*255) = 192のあたりに打たれます。

線をつなぐ必要はなく、点だけを打っても構いません。401個の点が256個の列に入るので、自然につながって見えます。

CSVは、open(path).read().splitlines()で読んで、最初の行(ヘッダー)を飛ばしたあと、カンマで切れば済みます。

赤い曲線だけが上に跳ね上がり、緑と青は、下のほうに張り付いている必要があります。その絵が、このラボの結論です。

dtを変えるとどれだけ悪くなるか

陽的オイラー法で、同じシミュレーション時間T=20を、dt=0.01(2000ステップ)、dt=0.05(400ステップ)、dt=0.2(100ステップ)で回して、最後のエネルギーが最初の何倍かを、/root/integrator/out/06-dt.txtに、dt001=、dt005=、dt02=の3行で書いてください(小数点以下6桁)。

ステップ数は、T / dtです。3つの場合とも、同じ20秒を真似ていますが、結果はまったく違います。

この系には、閉じた式があります。1ステップごとに、エネルギーがちょうど(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))を作ってください。

1ステップは、次のとおりです。

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

エネルギーが、一度も増えずに減る必要があります。400ステップ後には、最初の1パーセント未満まで下がります。

速度に毎ステップ0.99を掛けるという、よくある方式と、比べてみてください。その方式は、簡単ですが、減衰量がdtによって変わるので、フレームレートが異なる機械で、物体が異なる動きをします。力として入れれば、dtが変わっても、同じ物理を真似ます。