x_{n + 1} = x_{n} + v_{n} Δ t, v_{n + 1} = v_{n} + a (x_{n}) Δ t

計算の裏側

オイラー法とルンゲ=クッタ法

Euler and Runge-Kutta methods

式は正しいのに、惑星は回るたびに外へ逃げていく。

運動方程式を時間の刻み Δt で進める 4 つの方法を、同じ惑星の軌道で比べる。エネルギーが増える理由、誤差の次数、物理エンジンが使う方法、Δt の決め方まで。

種類: 計算の裏側
公式: 6 本
人物: レオンハルト・オイラー、カール・ルンゲ、マルティン・ヴィルヘルム・クッタ、ルー・ベルレ

運動方程式万有引力の法則ケプラーの第三法則テイラー展開振り子の周期フックの法則

場面

逃げていく惑星

運動方程式と万有引力の法則を、10 行ほどのプログラムに書いてみる。太陽の周りに惑星を 1 つ置き、1 コマごとに、位置に「速度 × Δt」を足し、速度に「加速度 × Δt」を足す。Δt は 1 コマ(フレーム)ぶんの時間の刻みである。

式に誤りは無い。それでも惑星は、回るたびに少しずつ太陽から離れていく。軌道は渦巻きになり、何周かすると画面の外へ飛び去る。現実の惑星は、ほぼ同じ楕円を何十億年も回り続けているのに、である。

ゲームを作ると、同じことが別の形で起こる。ばねでつないだ物体が、何もしていないのに揺れを大きくして暴れだす。跳ねるボールが、跳ねるたびに高く上がる。どれも物理の法則のせいではなく、時間を刻んで進める計算の方法から来ている。

この記事では、運動方程式を Δt ずつ進める代表的な 4 つの方法、陽的オイラー法・半陰的オイラー法・ベルレ法・4 次のルンゲ=クッタ法(RK4)を、同じ惑星の軌道で比べる。違いがどこから来るのか、物理エンジンが実際にどれを使っているのか、Δt をどう決めればよいのかを順に考える。

考え方

1 歩ずつ進める

運動方程式 $F = ma$ は、加速度が「速度の変わる速さ」、速度が「位置の変わる速さ」であることを使って、動きを決める。天体が 3 つ以上になったり衝突が絡んだりすると、動きを式で書き下せることはほとんど無い。そこで計算機は、時間を $Δ t$ ごとに区切り、飛び飛びの時刻 $t_{0}, t_{0} + Δ t, t_{0} + 2Δ t, \dots$ での位置と速度を順に求める。これが常微分方程式の数値解法(数値積分)である。

いちばん素直な考え方は、短い間なら速度も加速度もほとんど変わらないとみなすことである。 $Δ t$ の間は速度が一定だとすれば、位置は「速度 × $Δ t$ 」だけ進む。これがオイラー法である。

どれだけずれるかは、テイラー展開で分かる。本当の位置は

x (t + Δ t) = x (t) + v (t) Δ t + \frac{1}{2} a (t) Δ t^{2} + \dots

であり、オイラー法は 2 項目までしか使わない。捨てた部分が 1 歩ごとの誤差(局所誤差)で、 $Δ t^{2}$ に比例する。決まった時間 $T$ を進むには $T /Δ t$ 歩かかるので、積もる誤差(大域誤差)は $Δ t^{2} \times T /Δ t$ 、つまり $Δ t$ に比例する。 $Δ t$ を半分にすれば誤差も半分になる。こういう方法を 1 次の方法と呼ぶ。

次数の高い方法は、1 歩の中で加速度を何度か調べ、展開のより先の項まで本当の解に合わせる。ベルレ法は 2 次で、 $Δ t$ を半分にすると誤差は 1/4 になる。RK4 は 1 歩の中で傾きを 4 回求めて重みをつけて平均し、 $Δ t^{4}$ の項まで合わせる。 $Δ t$ を半分にすると誤差は 1/16 になる。

ただし、次数だけで方法の良し悪しは決まらない。長い時間回し続けたとき、エネルギーが増え続けるのか、減り続けるのか、揺れながら同じところに留まるのか。ゲームや天体の計算では、こちらの性質が結果を大きく左右する。

式

4 つの方法の式

運動方程式を、位置と速度の 2 つの式に分ける。 $a (x)$ は位置 $x$ での加速度(力 ÷ 質量)である。
$\frac{d x}{d t} = v, \frac{d v}{d t} = a (x)$
陽的オイラー法。今の速度で位置を、今の加速度で速度を進める。
$x_{n + 1} = x_{n} + v_{n} Δ t, v_{n + 1} = v_{n} + a (x_{n}) Δ t$
半陰的オイラー法(シンプレクティック・オイラー法)。先に速度を進め、その新しい速度で位置を進める。手順 2 の順番を入れ替えただけである。
$v_{n + 1} = v_{n} + a (x_{n}) Δ t, x_{n + 1} = x_{n} + v_{n + 1} Δ t$
ベルレ法(速度ベルレ法)。速度を半歩、位置を 1 歩進め、新しい位置での加速度で速度の残りの半歩を進める。最後に求めた加速度は次の歩でそのまま使えるので、力の計算は 1 歩に 1 回で済む。
$v_{n + \frac{1}{2}} = v_{n} + \frac{1}{2} a (x_{n}) Δ t, x_{n + 1} = x_{n} + v_{n + \frac{1}{2}} Δ t, v_{n + 1} = v_{n + \frac{1}{2}} + \frac{1}{2} a (x_{n + 1}) Δ t$
4 次のルンゲ=クッタ法(RK4)。位置と速度をまとめて $y = (x, v)$ 、その変わる速さを $f (y) = (v, a (x))$ と書き、1 歩の初め・半ば(2 回)・終わりで傾きを求める。
$k_{1} = f (y_{n}), k_{2} = f (y_{n} + \frac{Δ t}{2} k_{1}), k_{3} = f (y_{n} + \frac{Δ t}{2} k_{2}), k_{4} = f (y_{n} + Δ t k_{3})$
4 つの傾きを 1 : 2 : 2 : 1 の重みで平均して進める。力の計算は 1 歩に 4 回要る。
$y_{n + 1} = y_{n} + \frac{Δ t}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$
ばね( $a = - ω^{2} x$ )を陽的オイラー法で 1 歩進め、質量 1 あたりのエネルギー $E = \frac{1}{2} (v^{2} + ω^{2} x^{2})$ を計算する。 $xv$ の項が打ち消し合い、どんな状態からでもエネルギーは同じ割合で増える。
$E_{n + 1} = (1 + ω^{2} Δ t^{2}) E_{n}$
同じばねを半陰的オイラー法で進めると、エネルギーそのものは揺れるが、それを少しずらした次の量が正確に保たれる( $C$ は一定)。だから揺れが大きくなり続けることは無い。
$v_{n}^{2} + ω^{2} x_{n}^{2} - ω^{2} Δ t x_{n} v_{n} = C$

手順 7 と 8 は、1 歩の式を代入して展開すれば確かめられる。

図

1 歩の図

円の軌道の上の点 P から 1 歩進める。陽的オイラー法は接線の向きにまっすぐ進むので、中心からの距離は √(r² + (vΔt)²) になり、必ず円の外に出る。半陰的オイラー法は先に速度を中心へ曲げてから進むので、この歩では円の内側に入る。

陽的オイラー法の 1 歩は、円の接線に沿った直線である。直角をはさむ 2 辺が $r$ と $v Δ t$ の直角三角形の斜辺になるので、三平方の定理から中心までの距離は $r^{2} + (v Δ t)^{2}$ となり、必ず $r$ より長い。次の歩も、その次の歩も外へずれ、ずれは打ち消されずに積もる。これが渦巻きの正体である。半陰的オイラー法では、ある歩では内側へ、別の歩では外側へとずれが行き来する。

シミュレーション

4 つの方法を同時に動かす

太陽の周りの楕円軌道(長半径 1 au、離心率 0.3、1 周 1 年)を、4 つの方法で同じ位置と速度から同時に進める。刻み $Δ t = 10$ 日では、陽的オイラー法は 1 周でエネルギーが 73 % 増えて外へ逃げ、半陰的オイラー法とベルレ法は軌道の向きが少しずつ回りながらも大きさを保ち、RK4 は点で描いた正しい軌道にほぼ重なる。「エネルギー(100 周)」では、半陰的オイラー法が -6 % から +9 %、ベルレ法が -1.5 % から 0 の間で揺れ続けるのに対し、RK4 は 1 周に約 0.024 % ずつ減り続ける(方法を消すと目盛りが細かくなる)。「誤差と Δt」では、 $Δ t$ を半分にしたときにずれが 1/2・1/4・1/16 に減る違いが、線の傾きとして見える。

考え方

違いが出る理由

確かめられていること。 陽的オイラー法がエネルギーを増やすのは、上の図のとおり、1 歩が接線に沿って外へずれるからである。ばねでは手順 7 のとおり、1 歩ごとにエネルギーがちょうど $1 + ω^{2} Δ t^{2}$ 倍になる。1 秒に 60 コマ( $Δ t = 1/60$ 秒)で周期 0.5 秒のばねを動かすと、1 歩で 1.044 倍、1 秒で約 13 倍になる。何もしていないのにばねが暴れだすのは、このためである。

半陰的オイラー法とベルレ法は、シンプレクティックと呼ばれる性質を持つ。位置と速度を座標にとった平面(位相空間)で、1 歩ぶんの写像が面積を変えない。ばねなら手順 8 の量が正確に保たれる。力がなめらかなら一般の場合でも、本当のエネルギーとのずれが $Δ t$ の次数の程度(半陰的オイラー法は $Δ t$ 、ベルレ法は $Δ t^{2}$ )に留まり、非常に長い時間にわたって増えていかないことが証明されている(Hairer・Lubich・Wanner の総説、2003 年)。中心に向かう力だけで動く惑星では、どちらの方法も角運動量を正確に保つ。軌道の大きさが変わらず、向きだけが少しずつ回るのはこのためである。

RK4 はシンプレクティックではない。同じ総説は、陽的なルンゲ=クッタ法では、エネルギーも角運動量も、保たれるはずの量がどれも一定の値から離れていくと書いている。1 周あたりのずれはとても小さいが、決まった向きに積もる。このシミュレーションでは、惑星が少しずつ内側へ落ちていく向きだった。

ここからは推測を含む。 ゲームの物理エンジンが半陰的オイラー法やベルレ法を選ぶのは、1 歩に力の計算が 1 回で済む軽さと、長く回してもエネルギーが増え続けない性質の両方を持つからだと考えられる。RK4 の 4 次という精度は、力がなめらかに変わることを前提にしている。衝突のように力がとぎれる場面が多いゲームでは、その前提が崩れ、4 回ぶんの計算に見合う利点を出しにくい。

歴史

研究の歴史

1768 年、レオンハルト・オイラーはペテルブルクで出版した『積分学教程』第 1 巻の第 2 部第 7 章「微分方程式の近似による積分について」で、この方法を示した。解を小さな刻みで順に進め、各刻みの間は変化の速さを刻みの初めの値のまま一定とみなす。

1895 年、ドイツの数学者カール・ルンゲは『マテマティッシェ・アナーレン』第 46 巻の論文で、オイラーの方法を、積分の中点公式や台形公式にならって改良した。1 歩の中で仮の値を先に求め、そこでの傾きも使う。1900 年にはカール・ホインが 4 次までの条件を調べて 3 次の方法を示した。1901 年、マルティン・ヴィルヘルム・クッタは、ミュンヘンでの学位論文をもとにした論文で分析を 5 次まで進め、4 次の方法をすべて分類して、いま古典的なルンゲ=クッタ法と呼ばれる方法を示した(ブッチャーの総説、1996 年)。

ベルレ法の名前は、フランスの物理学者ルー・ベルレが 1967 年に『フィジカル・レビュー』に発表した論文から来ている。ベルレは、レナード=ジョーンズのポテンシャルで引き合う 864 個の分子の運動を計算機で解き、アルゴンの性質と比べた。分子の動きを 1 つずつ追う分子動力学では、この方法が最も広く使われるようになった。

方法そのものは、もっと古い。Hairer らの総説によれば、ノルウェーのカール・ストルマー(Carl Størmer)は 1907 年、地球の磁場の中を進むオーロラの荷電粒子の計算に、これを高い次数にしたものを使った。後に科学史に関心を持ったベルレ自身が、1792 年のドランブルの対数表と天文表の計算に同じ方法を見つけている。さらに、ニュートンが 1687 年の『プリンキピア』でケプラーの第二法則を示した議論は、「まっすぐ進み、途中で速度を中心へ一度に曲げ、またまっすぐ進む」という、ベルレ法と同じ形をしている。物理学者ファインマンが 1964 年のコーネル大学での講演でこの議論を紹介している。

使いどころ

実務での使いどころ

物理エンジンの中身。 2 次元の物理エンジン Box2D のソースコード(v2.4)では、まず重力と力から求めた加速度で速度を更新し、次にその新しい速度で位置を進めている。半陰的オイラー法の順番である。v3 では 1 コマを細かいサブステップに分け、それぞれで同じ順に速度と位置を進める。文書はサブステップの数を 4 にすることを勧めている。JavaScript の物理エンジン Matter.js は、前のコマとの位置の差から速度を求めるベルレ法の形で、ソースコードにも「ベルレ積分で速度を更新する」と注記がある。ゲーム開発者のグレン・フィードラーは、2004 年の記事「Integration Basics」の初めの版では RK4 を勧めていた。のちに書き直した今の版では、市販のゲーム用物理エンジンの多くは半陰的オイラー法を使っていると書いている。なお、このサイトの万有引力の法則のシミュレーションは速度ベルレ法で、振り子の周期のシミュレーションは RK4 で 1 コマを 8 つに刻んで動かしている。

Δt の決め方。

固定する。画面のコマの時間は機械や負荷で揺れるので、物理の $Δ t$ は一定にし、たまった時間を固定の $Δ t$ で何回かに分けて進める(フィードラー「Fix Your Timestep!」)。Box2D の文書は、刻みを 1/30 秒以下にし、1/60 秒ならたいてい質の高い計算になると書いている。
いちばん速い振動に合わせる。角振動数 $ω$ のばねを半陰的オイラー法かベルレ法で動かすとき、発散しない条件は $ω Δ t < 2$ 、RK4 では $ω Δ t < 22 \approx 2.83$ である。陽的オイラー法には、ばねが発散しない $Δ t$ が無い。1 秒に 60 コマなら、半陰的オイラー法で扱えるのは $ω < 120$ rad/s、周期が 0.052 秒より長いばねまでになる。発散しなくても正確とは限らない。ベルレ法で 1 周期を 10 歩で刻むと振動数が約 1.7 % ずれ、20 歩なら約 0.4 %、40 歩なら約 0.1 % まで減る。
衝突は別に考える。1 コマで進む距離が、物体の半径と壁の厚さの半分の和より大きいと、壁をすり抜けることがある。これは衝突判定とすり抜けで扱う。

方法の選び方の目安。 毎コマ長く回し続けるゲームの物理には、半陰的オイラー法かベルレ法。空気の抵抗を受ける弾道を数秒先まで求めるような、短い時間を正確に計算したい場面には RK4。何百周も回す軌道には、シンプレクティックな方法が向く。同じ計算量で比べると違いがはっきりする。シミュレーションの軌道を、1 周 37 歩(約 10 日)の RK4 と、力の計算の回数が同じになる 1 周 148 歩(約 2.5 日)のベルレ法で進めた。1 周目のエネルギーのずれは RK4 が 0.022 %、ベルレ法が最大 0.097 % で、RK4 のほうが小さい。ところが RK4 のずれは 5 周で 0.11 % と並び、100 周では 2.3 % に広がる。ベルレ法は 100 周たっても 0.097 % のままである。なお、ゲームの中で弾の行き先を予測線で描くなら、本体の物理と同じ方法・同じ $Δ t$ で計算する。方法が違うと、線と実際の動きがずれる。

跳ねるボール。 床で跳ねるボールを、反発係数 1(エネルギーを失わない)・1 秒に 60 コマ・重力加速度 9.8 m/s² で、高さ 1 m から落として計算した。陽的オイラー法では、空中にいる間も 1 歩ごとにエネルギーが増え、最高点は 1 秒に $g Δ t /2 \approx 8$ cm ずつ上がる。床に食い込んだら床の上に戻して速度を逆向きにするだけの処理と組み合わせると、20 秒後の最高点は、陽的オイラー法では約 3.3 m まで上がり、半陰的オイラー法では逆に 0.2 m 以下まで下がった。床に戻す処理が、食い込んだ深さに応じてそのたびにエネルギーを変えるからで、その深さは床に触れるときの刻みの位置で決まる。重力加速度を 9.81 m/s² にしただけで、半陰的オイラー法の最高点は 0.89 m で止まった。コマの途中で床に触れた時刻まで考えて折り返すと、半陰的オイラー法は 1 m 前後を保った。積分の方法と衝突の処理は、組にして確かめる必要がある。

縦の動画で見る。 ラボの縦の短い動画のうち、跳ねるたびにボールが増える動画は半陰的オイラー法で、二重振り子の動画は 4 次のルンゲ=クッタ法で計算した。二重振り子は、初めの角度を 1000 万分の 1 ラジアンだけ変えた 2 つを同じ方法・同じ $Δ t$ で進めている。

跳ねるたびに 1 個増えるボール(半陰的オイラー法・乱数の種 1)声 VOICEVOX:玄野武宏・効果音と音楽は合成音

二重振り子 1000 万分の 1 の差はいつばらける？(4 次のルンゲ=クッタ法)声 VOICEVOX:玄野武宏・効果音と音楽は合成音

確かめ方。 $Δ t$ を半分にして、同じ計算をもう一度行う。正確な答え(か、 $Δ t$ をずっと小さくした計算)とのずれが 1/2 になれば 1 次、1/4 なら 2 次、1/16 なら 4 次の方法が期待どおりに動いている。減り方が合わなければ、衝突や摩擦のように力がとぎれる部分か、プログラムの誤りを疑う。二重振り子のように初めのわずかな違いが急に広がる系(カオス)では、方法や $Δ t$ を変えるとしばらく後には別の動きになるので、軌道ではなく、エネルギーのような保存量のずれで計算の良し悪しを見る。

次にやること。 自分のゲームやツールの物理の更新が「速度を先に、位置を後に」の順になっているかを確かめる。位置を先に進めているなら陽的オイラー法で、2 行の順番を入れ替えるだけで半陰的オイラー法になる。20 行ほどのプログラムでばねを 4 つの方法で動かし、エネルギーを並べて表示してみると、この記事の数字を自分の手で確かめられる。

関連する読みもの。 速度が時刻だけで決まるとき、オイラー法で足していく和は、Σ と ∫ の読み方の左端の長方形の和と同じになる。たくさんの振動子を同時に刻んで進める例は、ホタルとメトロノームはなぜそろうのかにある。

参考文献

L. Euler, Institutionum calculi integralis, Vol. 1, St. Petersburg, 1768, Sect. II, Ch. 7(Eneström 番号 E342、Opera Omnia I-11)
C. Runge, "Ueber die numerische Auflösung von Differentialgleichungen", Mathematische Annalen 46 (1895), 167-178.
W. Kutta, "Beitrag zur näherungsweisen Integration totaler Differentialgleichungen", Zeitschrift für Mathematik und Physik 46 (1901), 435-453.
L. Verlet, "Computer 'Experiments' on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules", Physical Review 159 (1967), 98-103.
J. C. Butcher, "A history of Runge-Kutta methods", Applied Numerical Mathematics 20 (1996), 247-260.
E. Hairer, C. Lubich, G. Wanner, "Geometric numerical integration illustrated by the Störmer-Verlet method", Acta Numerica 12 (2003), 399-450.
E. Catto, Box2D documentation "Simulation"(box2d.org)と、ソースコード v2.4.1 の src/dynamics/b2_island.cpp・v3 の src/solver.c
Matter.js のソースコード src/body/Body.js(Body.update)
G. Fiedler, "Integration Basics"(2004・のちに改訂。初めの版は web.archive.org に保存がある)・"Fix Your Timestep!"(2004), Gaffer On Games.