5. さまざまな物理シミュレーション
「4. 物理シミュレーションの基礎」では、単純なボールの運動シミュレーションを紹介しました。 本章では、オイラー法で解くことができる、さまざまな物理シミュレーションを紹介します。
5.1 ばねの振動・減衰
オイラー法を用いてばね振り子の運動を計算する方法を見ていきましょう。
5.1.1 ばねの振動
単純なばね振り子の運動方程式を考えていきます。 おもりの質量を \(m\)、ばね定数を \(k\)、自然長からの変位を \(y\) とすると、 おもりに加わるばねの弾性力は \(-ky\) であるため、以下の運動方程式が成り立ちます。
\[ F = m\frac{d^2y}{dt^2} = -ky \tag{1}\]
式\((1)\)の右側の等式に注目すると、以下の等式が成り立ちます。
\[ \frac{d^2y}{dt^2} = -\frac{k}{m} y \tag{2}\]
ボールの自由落下のときと同様に、速度 \(v\) を導入することで、式\((2)\)は次のように分解することができます。
\[\frac{dy}{dt} = v \tag{3}\]
\[\frac{dv}{dt} = -\frac{k}{m} y \tag{4}\]
さらに、各微分方程式を差分近似すると、次のようになります。
\[ y_{n + 1} = y_n + h\cdot v_n \tag{5}\]
\[ v_{n + 1} = v_n + h\cdot \left(-\frac{k}{m} y_n\right) \tag{6}\]
繰り返しになりますが、この差分方程式を逐次的に解くことが、オイラー法です。
5.1.2 ばね振り子の数値シミュレーション
前節で算出した差分方程式を用いて、ばね振り子の運動を計算してみましょう。 ボールの時と異なり、位置と速度が互いに影響しあっている点に、注意が必要です。
例えば、ボールと同様に以下のようにすると、更新後の座標が速度の計算に使われてしまいます。
これは、更新の順序を入れ替えても同じです。これではオイラー法にはなりません。
オイラー法を実装するには、更新前の座標と速度のデータも保持しておく必要があります。
更新前の座標と速度をそれぞれ yPrev、vyPrev とすると、オイラー法により座標と速度を計算するソースコードは次のようになります。
新しい座標と速度を算出したら、次のループではそれらを用いてまた新しい座標と速度を算出します。
そのため、以下のように y を yPrev に、vy を vyPrev に入れ直す「値の更新」を行う必要があります。
おもりの物体情報は、オブジェクトで管理することとしましょう。 今回利用する定数やオブジェクトをまとめると以下のようになります。
const K = 0.5; // バネ定数
const M = 2; // 質量
const Y_0 = 100; // 初期座標（自然長からの変位）
const VY_0 = 0; // 初速度
const SIDE = 25; // 物体の一辺の長さ
// おもりオブジェクトの初期化
let box = { y: Y_0, vy: VY_0, yPrev: Y_0, vyPrev: VY_0 };
ここでは正方形の形をしたおもりを想定しているため、一辺の長さを表す定数 SIDE を定義していますが、どのような形でも構いません。
ある座標点 (x, y) を中心に一辺が SIDE の正方形を描くコードは以下のようになります。
思い描いたとおりの位置に描画できるよう慣れておきましょう。
上記 box オブジェクトに対するばね振り子による速度・座標の計算および値の更新をまとめると、次のようになります。
// 速度・座標の計算
box.y = box.yPrev + STEP * box.vyPrev;
box.vy = box.vyPrev - STEP * (K / M) * box.yPrev;
// 値の更新
box.yPrev = box.y;
box.vyPrev = box.vy;
上述の説明をすべてまとめると、以下のようなコードになります。 Canvas 右側には簡単なグラフ描画をしたいので、左側に寄せて表示しています。 ばね部分はおもりから画面端までを一定回数折り返しながら結ぶ描写をしていますが、 単なる見栄えの問題なので今は深く理解しなくても問題ありません（描写しなくても構いません）。
ばねが振動するようすが確認できると思います。 しかしながら、だんだんと振幅が大きくなる、おかしな挙動を示しています。これは、オイラー法の誤差によるものです。
また、このコード例では、おもりの座標を yData 配列に記憶させ、順番に
演習 9-1
演習 9-1
オイラー法による垂直ばね振り子のシミュレーションを完成させなさい。
5.1.3 ばねの減衰
同様にして、減衰器付き水平ばね振り子の運動方程式を考えてみましょう。 おもりの質量を \(m\)、ばね定数を \(k\)、減衰器の減衰係数を \(c\)、自然長からの変位を \(y\) とすると、 おもりに加わるばねの弾性力は \(-ky\) であり、減衰器は速度に比例した抵抗力 \(-c\frac{dy}{dt}\) が発生するため、 以下の運動方程式が成り立ちます。
\[ m\frac{d^2y}{dt^2} = -ky - c\frac{dy}{dt} \tag{7}\]
両辺を \(m\) で割り、以下の形にすることで、数値的に解を求めることができます。
\[ \frac{d^2y}{dt^2} = -\frac{k}{m} y - \frac{c}{m}\frac{dy}{dt} \tag{8}\]
速度 \(v\) を導入することで、式\((8)\)は次のように分解することができます。
\[\frac{dy}{dt} = v \tag{9} \]
\[\frac{dv}{dt} = -\frac{k}{m} y - \frac{c}{m} v \tag{10}\]
さらに、各微分方程式を差分近似すると、次のようになります。
\[ y_{n + 1} = y_n + h\cdot v_n \tag{11}\]
\[ v_{n + 1} = v_n + h\cdot \left(-\frac{k}{m} y_n - \frac{c}{m} v_n \right) \tag{12}\]
コードを記述して、動きを確認してみましょう。
演習 9-2
演習 9-2
垂直ばね振り子にダンパを実装しなさい。また、Cnavas に質点の変位を表すグラフを描画しなさい。
5.1.4 解析解との比較
オイラー法で計算した座標と解析解（厳密解）の結果を比較してみましょう。 単振動の解析解は、初期位置を \(y_0\)、初速度を \(v_0\) とすると、以下のように求まります。
\[ y = y_0 \cos\left( \sqrt{\frac{k}{m}} t \right) + v_0 \sqrt{\frac{m}{k}} \sin\left( \sqrt{\frac{k}{m}} t \right)\tag{13}\]
これを JavaScript で書くと、少々長いですが以下のようになります。
Y_0 * Math.cos(Math.sqrt(K / M) * time) +
VY_0 * Math.sqrt(M / K) * Math.sin(Math.sqrt(K / M) * time);
オイラー法による解を青線、解析解による解を赤線で描画するコードは以下のようになります。
実行すると、解析解（赤い線）はきれいな単振動となっていることがわかります。
5.2 水中での物体の運動（浮力）
水中の物体には、重力とは逆の方向に浮力がはたらきます。 物体の質量を\(m\)、重力加速度を\(g\)、水の密度を\(\rho_0\)、水中に沈んでいる物体の体積を\(V\)とすると、以下の運動方程式が成立します。
\[F = m \frac{d^2 y}{dt^2} = mg - \rho_0 V \tag{1}g\]
両辺を \(m\) で割り、以下の形にすることで、数値的に解を求めることができます。
\[\frac{d^2 y}{dt^2} = \left(1 - \frac{\rho_0 V}{m}\right) g \tag{2}\]
速度 \(v\) を導入することで、式\((2)\)は次のように分解することができます。
\[\frac{dy}{dt} = v \tag{3} \]
\[\frac{dv}{dt} = \left(1 - \frac{\rho_0 V}{m}\right) g \tag{4}\]
さらに、各微分方程式を差分近似すると、次のようになります。
\[ y_{n + 1} = y_n + h\cdot v_n \tag{5}\]
\[ v_{n + 1} = v_n + h\cdot \left\{\left(1 - \frac{\rho_0 V}{m}\right) g \right\} \tag{6}\]
以下は、実装例です。
55 ～ 57 行目で、水中に沈んでいる物体の体積を計算するための、沈み込み量を計算しています。 オイラー法による計算であるため、誤差が徐々に大きくなっていくようすも確認できます。
5.3 振り子の運動
単振り子は、おもりの質量を \(m\)、棒の長さを \(l\)、棒と鉛直軸のなす角を \(\theta\)、重力加速度を \(g\) とすると、 \(\theta\) に関する以下の運動方程式が成り立ちます（棒の質量は考えないものとしています）。
\[F = ml \frac{d^2\theta}{dt^2} = -mg \sin \theta \tag{1}\]
両辺を \(ml\) で割り、以下の形にすることで、数値的に解を求めることができます。
\[\frac{d^2\theta}{dt^2} = -\frac{g}{l} \sin \theta \tag{2}\]
角速度 \(\omega\) を導入することで、式\((2)\)は次のように分解することができます。
\[\frac{d\theta}{dt} = \omega \tag{3} \]
\[\frac{d\omega}{dt} = -\frac{g}{l} \sin \theta \tag{4}\]
さらに、各微分方程式を差分近似すると、次のようになります。
\[ \theta_{n + 1} = \theta_n + h\cdot \omega_n \tag{5}\]
\[ \omega_{n + 1} = \omega_n + h\cdot \left( -\frac{g}{l} \sin \theta_n \right) \tag{6}\]
おもりの座標は、\((x, y) = (l\cos\theta, l\sin\theta)\) と求めることができます。
（座標系のとり方によっては、変換式を間に挟む必要があります。）
以下は、実装例です。
5.3.1 二重振り子の運動方程式
カオス的挙動を示すことで知られる二重振り子の運動方程式は、固定点からの振れ角を \(\theta_1, \theta_2\) とすると、ラグランジュ関数を用いて以下のように導かれます。
\[\ddot{\theta}_{1} = \frac{g\left(m_2\sin\theta_2\cos(\theta_1-\theta_2) - (m_1 + m_2)\sin\theta_1\right) - m_2\left(l_1 \dot{\theta_1}^2\cos(\theta_1-\theta_2) + l_2 \dot{\theta_2}^2\right)\sin(\theta_1-\theta_2)}{l_1\left(m_1 + m_2\sin^2(\theta_1-\theta_2)\right)}\]
\[\ddot{\theta}_2 = \frac{g(m_1 + m_2)\left(\sin\theta_1\cos(\theta_1-\theta_2) - \sin\theta_2\right) + \left( (m_1 + m_2)l_1 \dot{\theta_1}^2 + m_2 l_2\dot{\theta_2}^2\cos(\theta_1-\theta_2) \right)\sin(\theta_1-\theta_2)}{l_2 \left( m_1 + m_2\sin^2(\theta_1-\theta_2) \right)}\]
この微分方程式を解析的に解くことは実質的に不可能であり、数値シミュレーションの出番というわけです。 煩雑な式に見えますが、数値計算法を用いればこれまでと同様の手順で解くことが可能です。
5.4 RLC 回路における電流の減衰
RLC 回路において、抵抗値を \(R\)、コイルのリアクタンスを \(L\)、コンデンサの静電容量を \(C\) とすると、 電流値 \(i\) に関して、キルヒホッフの法則から以下の微分方程式を導くことができます。
\[L \frac{d^2i}{dt^2} + R\frac{di}{dt} + \frac{i}{C} = 0 \tag{1}\]
これを変形して以下の形にすることで、数値的に解を求めることができます。
\[\frac{d^2i}{dt^2} = -\frac{R}{L}\frac{di}{dt} - \frac{i}{LC} \tag{2}\]
電流変化量 \(j\) を導入することで、式\((2)\)は次のように分解することができます。
\[\frac{di}{dt} = j \tag{3} \]
\[\frac{dj}{dt} = -\frac{R}{L}\frac{di}{dt} - \frac{i}{LC} \tag{4}\]
さらに、各微分方程式を差分近似すると、次のようになります。
\[ i_{n + 1} = i_n + h\cdot j_n \tag{5}\]
\[ j_{n + 1} = j_n + h\cdot \left( -\frac{R}{L}j_n - \frac{i_n}{LC} \right) \tag{6}\]
以下は、実装例です。
scale() 関数を用いることで、描画する三角形の大きさを電流の大きさに依存させています。
5.5 磁場・電場中の荷電粒子の運動
電荷 \(q\) を帯びた荷電粒子は電場や磁場からクーロン力やローレンツ力を受けるので、 以下の運動方程式が成り立ちます。
\[\boldsymbol{F} = m \frac{d^2\boldsymbol{p}}{dt^2} = q(\boldsymbol{E} + \boldsymbol{v}\times\boldsymbol{B})\]
ここで、\(\boldsymbol{p}\) は位置ベクトル、\(\boldsymbol{v}\) は速度ベクトル、\(\boldsymbol{E}\) は電場ベクトル、\(\boldsymbol{B}\) は磁場ベクトルを表します。 以下のコードは、\(\boldsymbol{E} = \boldsymbol{0}, \boldsymbol{B} = (0, 0, B)\) としたときの実装例です。
5.6 太陽系の惑星の運動
惑星の運動に関しては、万有引力の法則とケプラーの法則から、以下の運動方程式を導出することができます。
\[\ddot{r} = r\dot{\theta}^2 - \frac{GM}{r^2}, \ \ \ \ddot{\theta} = -\frac{2\dot{r}\dot{\theta}}{r}\]
ここで \(r\) と \(\theta\) は極座標系における動径と偏角、\(G\) は万有引力定数、\(M\) は太陽の質量を表します。 極座標空間で上記微分方程式を解き、直交座標系に変換して可視化を行うコード例を以下に示します。
5.7 修正オイラー法
4.1 節で導入したオイラー法よりも、精度が高い方法を見ていきましょう。
ある微分方程式 $ \frac{dy}{dt} = f(t, y)$ をオイラー法で解くためには、以下のように離散化することを 4.1 節で紹介しました。
\[ y_{n+1} = y_n + hf(t_n, y_n) \tag{1}\]
これは、積分値 \(\int_{t}^{t+h} f(t, y) dt\) を、\(hf(t_n, y_n)\) という矩形の面積に近似していることを意味します。
積分区間の刻み幅 \(h\) が十分小さければある程度の精度が得られることがわかりますが、 4.1 節のテイラー展開からもわかるように、オイラー法では精度向上に限界があります。
離散化誤差を小さくする方法として、台形公式に基づく面積計算 が考えられます。 オイラー法では \(t\) における勾配 \(f(t,y)\) のみを用いて矩形の面積を計算していましたが、 さらに次の時間における \(f(t_{n+1}, y_{n+1})\) を導入し、台形の面積を計算するというわけです。 この考えをオイラー法に反映させると、次のようになります。
\[ y_{n+1} = y_n + \frac{h}{2}\left\{f(t_n, y_n) + f(t_{n+1}, y^*_{n+1})\right\}\tag{2}\]
しかしながら、ここで問題が発生します。 右辺の \(t_{n+1}\) は \(t_n + h\) と簡単に求めることができますが、 \(y^*_{n+1}\) は現時点では値が不明です。 そこで、\(y^*_{n+1}\) はオイラー法を用いることで求めます。 式をまとめると次のようになります。
\[ y_{n+1} = y_n + \frac{h}{2}\left\{f(t_n, y_n) + f(t_{n+1}, y_n + hf(t_n, y_n))\right\}\tag{3}\]
この差分方程式を解くことを、修正オイラー法 (Modified Euler's method) もしくはホイン法 (Heun's method) と呼びます。 修正オイラー法の精度は二次であることが知られています。
5.7.1 修正オイラー法の実装
垂直ばね振り子を題材に、修正オイラー法の実装例を見てみましょう。 式\((3)\) をそのままコード化することもできますが、可読性を考慮して式\((1)\) と式\((2)\) の部分に分けることにします。
再掲になりますが、単純なばね振り子の運動方程式は、 おもりの質量を \(m\)、ばね定数を \(k\)、自然長からの変位を \(y\) とすると、以下のようになります。
\[ m\frac{d^2y}{dt^2} = -ky \tag{4}\]
これをオイラー法では以下のように離散近似します。
\[ y_{n + 1} = y_n + h\cdot v_n \tag{5}\]
\[ v_{n + 1} = v_n + h\cdot \left(-\frac{k}{m} y_n\right) \tag{6}\]
式\((5)\)、式\((6)\) で求まる \(y_{n+1}, v_{n+1}\) を、修正オイラー法で用いる仮の数値 \(y^*_{n+1}, v^*_{n+1}\) であるとしましょう。 このとき、修正オイラー法による差分方程式は、式\((2)\) より以下のようになります。
\[ y_{n + 1} = y_n + \frac{h}{2}(v_n + v^*_{n+1})\tag{7}\]
\[ v_{n + 1} = v_n - \frac{h}{2} \left(\frac{k}{m} y_n + \frac{k}{m} y^*_{n+1} \right) \tag{8}\]
式\((5)\)～\((8)\) をプログラムにすると、次のようになります。
// 速度・座標の計算（修正オイラー法）
yEuler = yPrev + STEP * vyPrev;
vyEuler = vyPrev - STEP * (K / M) * yPrev;
y = yPrev + 0.5 * STEP * (vyPrev + vyEuler);
vy = vyPrev - 0.5 * STEP * ((K / M) * yPrev + (K / M) * yEuler);
ここでは、yEuler、vyEuler をオイラー法の計算結果を入れておく一時変数として利用しています。
オイラー法、修正オイラー法、解析解の結果を比較するコード例を以下に示します。
オイラー法（青線）の誤差が大きくなるのに対して、修正オイラー法（緑線）と解析解（赤線）はほぼ一致していることがわかります。 そのまま描画すると見分けがつかないので、ここでは修正オイラー法の線を太くしています。
このように、精度の次数が増えると誤差は大きく減少します。 単振動のような単純な系では修正オイラー法でも十分な精度となることがありますが、 複雑な系では実用的ではなく、実際の現場では通常四次以上の精度を持つ方法がとられます。