4. 物理シミュレーションの基礎

4.1 常微分方程式の数値解法

いよいよ本章から、物理シミュレーションの中核となる基礎事項を学んでいきます。

4.1.1 解析解と数値解

物理シミュレーションの多くは、事象を支配する微分方程式を解くことで行われます。微分方程式の解は大きく、解析解（厳密解）と数値解に分けることができます。

解析解は、微分方程式の解を積分計算などにより代数的に求めたものです。たとえば、ボールの自由落下を表す二階常微分方程式 \(m\frac{d^2y}{dt^2} = mg\) の解析解は、\(y = y_0 + v_0 t + \frac{1}{2} gt^2\) です。

このように解析解が得られる場合、ある時刻 \(t\) におけるボールの位置は即座に求まります。しかしながら、系が少しでも複雑になると、たちまち微分方程式を解析的に解くことが困難・不可能になってしまいます。そのような場合には、微分方程式を差分近似することで、数値解を求めます。

数学モデル（微分方程式）を離散化して数値モデル（差分方程式）とし、これを解くことを数値シミュレーションと呼びます。

4.1.2 微分方程式から差分方程式へ

実際に微分方程式を離散化して差分方程式に変換する過程を追っていきましょう。以下の、単純な一階常微分方程式を数値的に解くことを考えます。

\[ \frac{dy}{dt} = f(t, y)\tag{1}\]

\(t\) は時間、\(y\) は座標を表すものとします（すなわち \(f(t, y)\) は速度を表します）。まず、微分項は以下のように書き換えることができます。

\[ \frac{dy}{dt} = \lim_{h\rightarrow 0} \frac{y(t + h) - y(t)}{h}\tag{2}\]

すなわち、十分小さい \(h\) を考えると、

\[ \frac{dy}{dt} \approx \frac{y(t + h) - y(t)}{h}\tag{3}\]

となります。これが微分方程式の差分近似です。式\((3)\)を式\((1)\)に代入して変形すると、次式が得られます。

\[ y(t + h) = y(t) + hf(t, y)\tag{4}\]

ここで、\(h\) は「時間の刻み幅」を表しています。ある時点の時刻を \(t_n\)、座標を \(y_n\)、そこから \(h\) だけ経過した時刻を \(t_{n+1}\)、座標を \(y_{n+1}\)とすると、式\((4)\)は以下のように表されます。

\[ y_{n+1} = y_n + hf(t_n, y_n) \tag{5}\]

この漸化式は、ある時刻 \(t_n\) における座標と速度がわかっていれば、次の時刻 \(t_{n+1}\) における座標を得ることができることを意味します。すなわち、初期条件 \(y_0\) を与えれば、以降の任意の時刻における座標 \(y_n\) は逐次的に算出可能です。このように差分方程式を用いて微分方程式を数値的に解くことをオイラー法 (Euler method)と呼びます。

4.1.3 テイラー展開によるオイラー法の導出

オイラー法の精度と誤差を説明する上で、テイラー展開を用います。ある時刻 \(t\) における座標を \(y(t)\)、そこから時間 \(h\) だけ経過したときの座標を \(y(t + h)\) とし、 \(y(t + h)\) を \(t\) の周りでテイラー展開すると、以下のようになります。

\[ y(t + h) = y(t) + y'(t)h + \frac{1}{2}y''(t)h^2 + \frac{1}{6}y'''(t)h^3 + \cdots \tag{6}\]

式\((6)\)において \(h\) が 1 次の項までを見ると、\(y' = f(t, y)\) であることから、オイラー法の式と一致します。このため、オイラー法の精度は 1 次であり、誤差は \(\mathcal{O}(h^2)\) であることがわかります。

十分小さい \(h\) を選べば問題なくシミュレーションすることも可能ですが、オイラー法は一般に精度が悪く、実際の現場では 2 次以上の精度を持つ修正オイラー法（ホイン法）、ベルレ法、ルンゲ・クッタ法などが用いられます。

高次の精度を持つ手法について学ぶ前に、まずはオイラー法の実装からはじめていきましょう。

4.2 ボールの運動

ボールの運動を題材に、オイラー法による物理シミュレーションを実装していきましょう。

4.2.1 ボールの運動方程式

実は、3 章のアニメーションで、オイラー法と同等のシミュレーションは実装済みです。摩擦のない床を水平方向に速度 \(v_x\) で等速運動するボールの運動方程式は、以下のようになります。

\[ \frac{dx}{dt} = v_x \tag{7}\]

オイラー法を適用するため、前節の手順で式(7)を差分近似すると、以下の式が得られます。

\[ x_{n + 1} = x_n + h \cdot v_x\tag{8}\]

式\((8)\) において \(h=1\) とすると \(x_{n + 1} = x_n + v_x\) となり、3 章のアニメーションで実装したプログラムとなります。 \(y\) 方向も同様に考え、時間刻み幅 \(h\) を定数（定数名: STEP）として導入したプログラムは次のようになります。

ここで、前回からプログラムに幾つか変更点を加えていることに注意して下さい。具体的な変更点を以下に記載します。

4.2.2 時間刻み幅の導入

前述のとおり、時間刻み幅の定数 STEP を導入し、初期値に 0.1 を与えています（変更してみましょう）。ボールの速度 vx, vy に時間刻み幅を乗じて座標 x, y に加算することで、次の時間における座標を計算します。

4.2.3 init() 関数の導入

シミュレーション空間の初期化を行うための関数 init() を導入しました。 ball オブジェクトの情報と経過時間を初期化する処理を記述しています。また、Canvas 要素やコンテキストの取得は初期化の際に一度だけ呼び出せばよいので、ここに移動させています。 init() は JavaScript 読込時に加えて、ユーザ操作により初期化の命令が与えられたときに呼び出します。

4.2.4 「1 ステップ進む」「リセット」ボタンの追加

シミュレーション実装の演習効率を向上させるために、二つのボタンを追加しました。

「1 ステップ進む」ボタンは、アニメーションを 1 フレーム分だけ進める機能を有します（stepOnce() 関数）。 boolean 型の変数 stepAnime をフラグ変数とし、一度だけ描画・座標の更新を行うようにしています。

「リセット」ボタンは、シミュレーション空間を初期化する機能を有します。 HTML のコードを見てもらうとわかりますが、先ほど導入した init() 関数を再利用しています。

4.2.5 時間情報の表示

Canvas 要素の上に div 要素（id="info"）を置き、経過時間の情報を表示するようにしました。 JavaScript コードの 1 行目で div 要素の取得、45 行目で innerHTML を参照して div 要素の中身を上書きしています。このとき、time.toFixed(2) とすることで、time の値を少数第 2 位まで表示させています。

経過時間は、71 行目で変数 time に STEP を加算代入することで計算しています。

4.3 重力の導入

等速直線運動であれば、上記のコードでオイラー法の実装が完了となりますが、重力の影響がある場合、物体は等加速度運動をします。

物体に加わる力（重力）は、ニュートンの運動方程式より以下のように表されます。

\[ F = m\frac{d^2y}{dt^2} = mg \tag{9}\]

ここで、\(m\) は質量、\(g\) は重力加速度です。式\((9)\)の右側の等式に注目すると、以下の等式が成り立ちます。

\[ \frac{d^2y}{dt^2} = g \tag{10}\]

式\((10)\)は、速度 \(v_y\) を導入することで、次のように分解することができます。

\[ \frac{dv_y}{dt} = g,\ \ \ \frac{dy}{dt} = v_y \tag{11}\]

さらに、各微分方程式を差分近似すると、次のようになります。

\[ y_{n + 1} = y_n + h\cdot v_{y_n} \tag{12}\]

\[ v_{n + 1} = v_n + h\cdot g \tag{13}\]

重力を考慮した数値シミュレーションをオイラー法を用いて実行するには、式\((12)(13)\) の差分方程式を利用すればよいわけです。実際に JavaScript で記述すると、以下のようになります。

// 位置の更新
y = y + STEP * vy;

// 速度の更新 (G: 重力加速度)
vy = vy + STEP \* G;

加算代入演算子を用いて以下のように書き直すこともできます。

// 位置の更新
y += STEP * vy;

// 速度の更新 (G: 重力加速度)
vy += STEP \* G;

これら更新式に基づき、重力の影響を考慮したシミュレーションのコード実装例を以下に示します。ここでは、重力加速度 \(G = 9.8\) と定義しています。

演習 8-2

オイラー法を用いたボールの運動シミュレーションを実装しなさい。