物理エンジンで作る縦の短い動画と、その中の数学

刻みの式の違いは、長く計算したときに出る。跳ねるボールの動画と同じ条件(重力 300 px/s²・円の半径 380 px・反発係数 1)で、1 個のボールだけを 120 秒計算し、力学的エネルギーが始めからどれだけずれたかを測った。反発係数 1 なので、正しい答えではエネルギーは変わらない(位置エネルギーの 0 は円の中心の高さに置き、始めのエネルギーに対する割合で表した)。

計算のしかた	10 秒	30 秒	60 秒	120 秒
陽的オイラー法	+7.2%	+18.6%	+37.0%	+75.8%
半陰的オイラー法	−0.8%	+0.6%	−2.4%	−2.3%

陽的オイラー法では、エネルギーが増え続け、120 秒で 1.76 倍になった。ボールは跳ねるたびに高く上がっていく。半陰的オイラー法のずれは、刻みごとに見ても −6.0% から +1.5% の間にとどまり、増え続けない。

なぜこう分かれるのかは、空中の 1 刻みを計算すると分かる。下向きを正にとり、上下の速さを $v_y$ とすると、1 刻みでのエネルギー $E$ の変化は、どの刻みでも同じ大きさになる。

$$\begin{array}{rl}\text{陽的オイラー法:}{E}_{n+1}-E_n& =+\frac{1}{2}{g}^2\Delta t^2\\ \text{半陰的オイラー法:}{E}_{n+1}-E_n& =-\frac{1}{2}{g}^2\Delta t^2\end{array}$$

どちらも空中だけなら 120 秒で 40% ぶん動く大きさである。半陰的オイラー法では、代わりに少しだけ形の違う量が空中でぴったり保たれる。

$$\tilde{E}=E+\frac{1}{2}g\Delta t{v}_y$$

$E$ と $\tilde{E}$ の差は $v_y$ に比例する小さな量なので、$E$ は $\tilde{E}$ のまわりで揺れるだけになる。このように、本当のエネルギーにごく近い別の量を長い時間保つ性質は、シンプレクティックと呼ばれる積分法(半陰的オイラー法・ベレ法など)に共通する[2]。

ただし壁での跳ね返りは、この性質の外にある。跳ね返りは刻みの終わりにまとめて行うので、本当にぶつかった瞬間とずれる。この動画では、壁にめり込んだ分を鏡のように内側へ戻す処理(次の章)が、そのずれをおおよそ打ち消している。戻す処理を外して同じ計算をすると、半陰的オイラー法でもずれは 120 秒で −40% まで広がった。陽的オイラー法の +75.8% のうち、空中の分は +40.2% で、残りは壁の処理との組み合わせで増えた分である。

4 次のルンゲ・クッタ法はシンプレクティックではないが、誤差がとても小さい。二重振り子を $\Delta t=1/480$ 秒で、放してから動画の終わりまで(約 29 秒)計算したとき、エネルギーのずれは最大でも $2.1\times 10^{-7}$(相対)だった。

ガルトンボードの名は、イギリスのフランシス・ゴルトン(Galton)から来ている。ゴルトンは 1870 年代に、釘を互い違いに並べた板を作り、1889 年の著書『Natural Inheritance』で、玉の山が釣り鐘の形になる様子を説明した[3]。

釘に当たるたびに玉が右か左へ半々で分かれるなら、14 段を通ったあとの箱の番号 $k$(左から 0〜14)は二項分布に従う。係数は二項定理とパスカルの三角形の、いちばん上を 0 段目として 14 段目の数である。

$$P(k)=\left(\genfrac{}{}{0ex}{}{14}{k}\right){\left(\frac{1}{2}\right)}^{14}$$

平均と標準偏差は次のとおり。段の数を増やすと、この形は正規分布に近づく(ド・モアブル=ラプラスの定理)。

$$\mu =14\cdot \frac{1}{2}=7,\sigma =\sqrt{14\cdot \frac{1}{2}\cdot \frac{1}{2}}\approx 1.87$$

ただし、物理エンジンの玉が釘 1 本ごとにきれいに半々で分かれるとは限らない。反発係数と空気の抵抗(Matter.js の frictionAir の値)を変えて 300 個ずつ落とし(乱数の種 1)、ばらつき(標準偏差・単位は箱の数)を測った。よく跳ねるほど山は広く、勢いを止めるほど狭くなった。

反発係数	空気の抵抗	ばらつき
0.3	0.03	3.23
0.2	0.04	2.16
0.2	0.045	1.89
0.15	0.045	1.68
0.2	0.05	1.55

動画では、理論の 1.87 にいちばん近かった反発係数 0.2・空気の抵抗 0.045 を使った。つまり、理論に合うように物理の値を選んでいる。乱数の種 3 で落とした 300 個は、左の箱から 0, 1, 2, 6, 22, 37, 57, 60, 52, 30, 18, 9, 4, 1, 1 個になった。平均は 6.99、ばらつきは 2.04 で、理論の 7 と 1.87 に近い。

いちばん右の箱にも 1 個入った。理想の板なら 1 個の玉がそこに入る確率は 16384 分の 1 で、300 個のうち 1 個でも両端に入る確率は約 3.6% である。ただし、この種は端に入る玉があるものを選んでいる。同じ値で種 1 では、ばらつきは 1.89 で、両端の 2 箱ずつには 1 個も入らなかった。300 個を 1 回落としただけでは、この板が理想より端に転がりやすいかどうかは言えない。

跳ねるボールの中心の決まりは「壁に当たるたびに、その場所から 1 個増える」である。新しいボールは、跳ね返った向きを 0.06〜0.22 ラジアン(乱数で決める)回した向きに飛ぶ。

ボール 1 個が、増える当たりを平均して $\tau$ 秒に 1 回起こすとすると、$N$ 個のボールは $\tau$ 秒でおよそ $N$ 個増える。増え方は今の数に比例する。

$$\begin{array}{rl}\frac{dN}{dt}& \approx \frac{N}{\tau }\\ N(t)& =e^{t/\tau }=2^{t/T},T=\tau \ln 2\end{array}$$

$e$ はネイピア数、$T$ は数が倍になる時間である。動画の記録から、8 個・16 個・…・1024 個になった時刻を取り、${\log }_{2}N$ と時間に直線を当てはめると、$T\approx 1.67$ 秒になった。1 個から 1024 個($2^{10}$)までは 17.6 秒で、$10T$ より少し長い。最初の 1 個から 8 個までに 5.9 秒(倍になるのに 1 回あたり約 2.0 秒)かかったためである。数が少ないうちは、当たる間隔のばらつきがそのまま出る。動画の最後のグラフで、個数の対数がほぼ 1 本の直線に並ぶのは、この式のためである。

決まりには細かい条件を 2 つ付けた。壁へ向かう速さが 120 px/s に満たない当たりでは増やさない。同じボールは 0.25 秒あける。同じ種で計算し直すと、2 つとも外したときは 64 個から 1024 個まで 1.1 秒(条件ありは 7.0 秒)、1 個から 1024 個まで 10.9 秒で増えた。120 px/s の条件だけを外すと 15.5 秒、0.25 秒の条件だけを外すと 17.6 秒のままだった。弱い当たりが続けて起きたとき、2 つの条件がそろって増えすぎを止めている。

本物の生き物の数は、こうして倍々に増え続けることはなく、やがて頭打ちになる(ロジスティック方程式)。動画では 1024 個で増えるのを止め、壁を開いて終わりにした。

二重振り子は、振り子の先にもう 1 本の振り子をつないだものである。1 本の振り子は、揺れが小さければ同じ周期で規則正しく揺れる。2 本つなぐと、動きは初めの角度のわずかな差で大きく変わる。

2 つのおもりの質量と 2 本の腕の長さ $L$ を同じにし(腕の重さは考えない)、角度は真下から測る。角度の差を $\Delta ={\theta }_1-{\theta }_2$、角速度を ${\omega }_1,{\omega }_2$ とすると、ラグランジュの方程式から角加速度は次の形になる。

$$\begin{array}{rl}{\ddot{\theta }}_1& =\frac{N_1}{L(3-\cos 2\Delta )}\\ {\ddot{\theta }}_2& =\frac{N_2}{L(3-\cos 2\Delta )}\end{array}$$ $$\begin{array}{rl}{N}_1& =-3g\sin {\theta }_1-g\sin ({\theta }_1-2{\theta }_2)\\ & -2\sin \Delta ({\omega }_2^{2}L+{\omega }_1^{2}L\cos \Delta )\\ N_2& =2\sin \Delta (2{\omega }_1^{2}L+2g\cos {\theta }_1\\ & +{\omega }_2^{2}L\cos \Delta )\end{array}$$

動画の 7 本は、上の腕 1.9 ラジアン・下の腕 2.9 ラジアンから放した。違うのは下の腕の角度だけで、隣どうし $10^{-7}$ ラジアン(1000 万分の 1)ずつずらした。この記事の「差」は、7 本の下の腕の角度のうち、いちばん離れた 2 本の差である(1 回転の違いは数えない)。始めは $6\times 10^{-7}$ ラジアンである。

カオスでは、差がおよそ一定の割合で増えていく。

$$\delta (t)\approx {\delta }_0{e}^{\lambda t}$$

$\lambda$ はリアプノフ指数と呼ばれ、正なら近い 2 つの動きが指数関数的に離れる。差の対数と時間に、放して 0.5 秒後から差が 0.1 ラジアンになるまでの区間で直線を当てはめると、$\lambda \approx 0.81$ /秒だった。差は約 2.83 秒ごとに 10 倍になる。差が 0.01 ラジアンになったのは放してから 10.5 秒後、1.5 ラジアンで 13.4 秒後である。この $\lambda$ は、この初めの角度と短い区間で測った目安の値である。カオスの振り子の実験と計算は、Shinbrot らの論文がくわしい[4]。近い 2 つの初めの値の差が指数関数的に広がる点は、ロジスティック写像のカオスと同じである。

同じことは計算の誤差にも起きる。刻みを半分の 1/960 秒にして計算し直すと、1 本ごとの下の腕の角度は、放して 10 秒で $7\times 10^{-6}$ ラジアン、20 秒で 0.7 ラジアン変わった。誤差も $e^{\lambda t}$ で増えるため、放して 20 秒を過ぎた動きは、細かいところまでは正しくない。それでも、ばらけるまでの時刻(10.5 秒・13.4 秒)と $\lambda$ は刻みを半分にしても同じだった。刻みの誤差は 7 本にほぼ同じように乗るので、差の広がり方は追える。

動画の音は、7 本の下の腕の角速度の大きさに比例して高くした音である。そろっている間は 1 つの音に聞こえ、ばらけると 7 つの音に分かれる。

コマは Node.js の canvas(@napi-rs/canvas)で描き、色は黒 #0a0a0a と白 #f6f6f4 だけにした。数式は KaTeX で組み、Chrome で透明な画像にして重ねた。描いたコマを 1 枚ずつ ffmpeg に渡し、H.264(縦 1080×1920・30 コマ/秒)と AAC の mp4 にした。音の大きさは -14 LUFS にそろえた(実測 -15.4〜-14.8 LUFS)。

縦の動画は、上と下に再生の操作やアカウント名が重なる。重なる範囲は資料によって数字が違う(上 108〜269 px・下 320〜672 px)。この動画では、文字を上 250 px・下 420 px・左 80 px・右 150 px より内側に置いた。

効果音と BGM は、すべてプログラムで合成した音である。ガラス玉の当たる音は 2 つの高い音と短い雑音、釘の音ははじく音を 5 音の音階で鳴らした。素材サイトの音を使うときは、商用利用とクレジットの決まりがサイトごとに違う。比べた表は、サービス・システム大全の効果音の素材サイト 12 選にある。