自由落下をオイラー法とルンゲクッタ法で解く

以前こちらで雑に考察した自由落下。

かれこれ3年以上前ということで感慨深いものがありましたが、

中途半端に終わっているので決着をつけようと思います。

 

1. 自由落下の運動方程式

2. 解析的手法

3. 解析的手法（跳ね返り）

4. オイラー法

5. オイラー法（跳ね返り）

6. ルンゲクッタ法

7. ルンゲクッタ法（跳ね返り）

8. C言語による実装（跳ね返り無し）

9. Javaによる実装（跳ね返りあり）

 

真空中を自由落下する物体には重力のみがはたらくので、次の運動方程式が成立します。

$$\displaystyle \frac{d^2}{dt^2}y = -g \ \ \ \ (1)$$

 

この式をオイラー法、ルンゲクッタ法で位置 $y$ について解き、解析的手法との比較を行います。

 

式(1)を1回積分、2回積分することで、次式が得られます。

$$\displaystyle v = v_0 – gt \ \ \ \ (2)$$
$$\displaystyle y = y_0 + v_0t – \frac{1}{2}gt^2 \ \ \ \ (3)$$

ここで $v$ は速度、 $v_0$ は初速度、 $y_0$ は初期位置を表します。

 

式(2)、(3)によってボールの位置が解析的に明らかになりました。

少し複雑になりますが、プログラム上で数値解法と比較するため、

時間ステップ $H$ で進めた際の、解析的な跳ね返りについて考えます。

一先ず落下のみ考えたいという方は、跳ね返りの部分は読み飛ばして下さい。

 

解析的な跳ね返りは次のような考えで実装しました。

位置がマイナスになったとき、次の処理を行います。

1. 初期位置 $y_0$ 、初速度 $v_0$ から、地面に接する（位置が $ 0$ になる）時刻 $t_c$ を求める。
2. 時刻 $t_c$ における速度を求め、跳ね返り係数を乗じて初速度 $v_0$ の値を更新する。
3. 初期位置 $y_0$ の値を $ 0$ に更新する。
4. 地面接触までの累計経過時刻 $t_n$ (初期値は $ 0$ )に、 $t_c$ を加算する。
5. 現在時刻 $t$ から $t_n$ を減算した時間（ $t – t_n$ ＜ $ H$ ）だけ進んだ位置を求め、$y$ を更新する。
6. 次回ステップ以降の計算においても、跳ね返りの瞬間時刻を起点とするため、 $t_n$ を減算した時間を用いて位置を算出する。

若干ややこしいですが、この考えによって、数値計算と解析的手法の比較が可能になります。

地面に接する時刻 $t_c$ は、式(3)の左辺に $ 0$ を代入して $t$ について解くことで、次のように求まります。

$$\displaystyle t_c = v_0 + \frac{\sqrt{v_0^2+2gy_0}}{g}$$

 

次はオイラー法による手法をまとめていきます。

式(1)の二階常微分方程式は、速度 $v$ を用いて次のように一階常微分方程式に分解できます。

 

$$\displaystyle\frac{d}{dt}y = v \ \ \ \ (4)$$

$$\displaystyle\frac{d}{dt}v = -g \ \ \ \ (5)$$

 

オイラー法では、微分方程式を差分方程式に近似することで、数値計算可能な形にします。

式(4)、(5)は、次のような差分方程式に近似されます。

 

$$\displaystyle y_{n+1} = y_n + v_n \times H \ \ \ \ (6)$$

$$\displaystyle v_{n+1} = v_n – g \times H \ \ \ \ (7)$$

 

ここで $H$ は時間ステップを表します。

式(6)、(7)の $n$ の値を増やしていくことで、数値的に自由落下の位置を計算できます。

実際のプログラム上ではメモリ節約のため、直前の $y_n$ 、 $v_n$だけ記憶させておき、

値の更新をすることによって変数を再利用するのが一般的です（多分）。

また、式(6)、(7)では、自由落下（ $v_0 = 0$ ）の際には1番初めのステップを進めても変化がありません。

それに対して、解析解でははじめから少しずつ落下します。ここで時間 $H$ のズレが発生します。

これを解消するために、まずは $v_{n+1}$ の値を更新し、最新の $v_{n+1}$ の値を、 $y_{n+1}$ を求める際に利用します。

数式は以下のようになります。

 

$$\displaystyle v_{n+1} = v_n – g \times H$$

$$\displaystyle y_{n+1} = y_n + v_{n+1} \times H \ \ \ \ (8)$$

 

私は昔、位置がマイナスになったら位置を $ 0$ とし、

速度に $-e$ （反発係数）を乗じるだけ、という方法をとっていましたが、

さすがにお粗末なので、今回はちゃんとした跳ね返りを考えます。

と言っても実装上は解析的手法の跳ね返りよりも簡単で、

位置がマイナスになったとき、次のような手順を行います。

1. 直前の位置 $y_n$ と速度 $v_n$ を用いて、地面に接する（位置が $ 0$ になる）時刻 $t_c$ を求める。
2. 時刻 $t_c$ を用いて、地面に衝突した瞬間の速度 $v_c$ を計算する。
3. 位置を $ 0$ 、跳ね返り直後の速度を $-v_c\times e$とする。
4. 次のステップまでの残り時間 $H-t_c$ だけ時間を進め、 $y_{n+1}$ 、 $v_{n+1}$ を算出する。

以上となります。

（時間を巻き戻して）地面に接する瞬間の時刻、次のステップまでに跳ね返っている位置を求める点で、

解析的手法と同様ではありますが、時間の起点を考慮しなくて良い分考えやすくなっています。

地面に接する時刻 $t_c$ は、式(8)の左辺に $ 0$ を代入して $t$ について解くことで、次のように求まります。

$$\displaystyle t_c = \frac{y_n}{v_{n+1}}$$

 

そしてルンゲクッタ法。

ルンゲクッタ法は複数勾配の重み付け平均を利用することで、高精度な近似を行う手法です。

今回は一般的な4次のものを利用します。

詳細は割愛して、今回のケースで利用する式をどんどん書いていきます。

 

二階微分方程式のルンゲクッタ法ということで混乱しそうになるかもしれませんが、

分解した各一階微分方程式の勾配が、計算過程で相互に影響し合うように書きます。

 

$$\displaystyle y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 +2 k_3 + k_4) \ \ \ \ (9)$$

$$\displaystyle v_{n+1} = v_n + \frac{1}{6}(l_1 + 2l_2 +2 l_3 + l_4) \ \ \ \ (10)$$

 

と定義したとき、各勾配は次のように計算していきます。

 

$$\displaystyle k_1 = H\times v_n$$

$$\displaystyle l_1 = H\times (-g) $$

$$\displaystyle k_2 = H\times \left( v_n +\frac{ l_1}{2} \right)$$

$$\displaystyle l_2 = H\times (-g) $$

$$\displaystyle k_3 = H\times \left( v_n + \frac{l_2}{2} \right)$$

$$\displaystyle l_3 = H\times (-g) $$

$$\displaystyle k_4 = H\times (v_n + l_3) $$

$$\displaystyle l_4 = H\times (-g) $$

 

これら式を用いることで、ルンゲクッタ法による自由落下を計算することができます。

 

ルンゲクッタ法による跳ね返りは、オイラー法の時と同様の流れになります。

しかし、地面に接する時刻 $t_c$ を求めるのは若干ややこしくなります。

勾配の式を整理しながら式(9)に代入していくと、次のようになります。

$$\displaystyle y_{n+1} = \frac{v_n H + 4\left( v_n – \frac{gH}{2} \right) + H(v_n – gH)}{6}$$

$$\displaystyle y_{n+1} = y_n + v_n H – \frac{1}{2}gH^2 \ \ \ \ (11)$$

式(11)の左辺を $ 0$ とすると、その時の $H$ が $t_c$ であるため、この連立方程式を解くと

$$\displaystyle t_c = \frac{v_n + \sqrt{v_n^2 – 2gy_n}}{g}$$

となります。

・・・ここまで計算して、やってしまった感が出てきました。

4次のルンゲクッタの式を変形した結果、 $t_c$ が解析解のものと一致しています。

自由落下の問題において、4次のルンゲクッタによる数値解は、解析解に等しいということに帰着してしまいそうです。

そもそも式(11)の段階でそのことが半分証明されていました。もう証明してしまいましょう。

勾配の式を整理しながら式(10)に代入していくと、次のようになります。

$$\displaystyle v_{n+1} = v_n – gH \ \ \ \ (12)$$

式(12)は公差 $-gH$ の等差数列なので、初項を $v_0$ とすると

$$ v_n = v_0 – gH(n-1) $$

となります。

また、式(11)の $n$ の値をずらしながら方程式を立て、加算していくことで次式が得られます。

$$\displaystyle y_{n+1} = y_0 + \sum_{i=0}^n v_n H – \frac{1}{2}gH^2n$$

$$\displaystyle y_{n+1} = y_0 + \sum_{i=0}^n \{v_0-gH(n-1) \} H – \frac{1}{2}gH^2n$$

$$\displaystyle y_{n+1} = y_0 + (v_0 + gH) Hn – \sum_{i=0}^n gH^2n – \frac{1}{2}gH^2n$$

$$\displaystyle y_{n+1} = y_0 + v_0 Hn + gH^2n – gH^2\frac{(n+1)n}{2} – \frac{1}{2}gH^2n$$

$$\displaystyle y_{n+1} = y_0 + v_0 Hn – \frac{1}{2}g(Hn)^2$$

ここで、時刻 $t$ は、ステップ $H$ を $n$ 回進めた時間 $Hn$ に等しいため、

$$\displaystyle y_{n+1} = y_0 + v_0 t – \frac{1}{2}gt^2$$

以上により、自由落下の解析解と、ルンゲクッタによる数値解が等しいことが証明されました。

証明されてしまいましたが、見なかったことにして実装していきます。
 

跳ね返りを考えない自由落下を、C言語で実装すると次のようになります。

計算結果を「freefall.csv」に書き出すプログラムとなっています。

自由落下だけということもあり、割と短いソースコードになりました。

※コードをコピーする際は、タブ文字がスペースに置換されない単純ビュー（画面ダブルクリック）にすることをお勧めします。

書き出されたcsvをExcelで開き、グラフ化してみます。

解析解とオイラー法では僅かに誤差が観察されました。

ルンゲクッタ法とは線分が完全に一致してしまっています。

次に、解析解との差分を計算し、誤差曲線を描いてみました。

オイラー法では誤差は線形に増加していくのがわかります。

同じスケールでは、ルンゲクッタ法の誤差は全くわかりません。

拡大すると次のようになりました。

10のマイナス13乗のスケールでこつこつと誤差が増えていきます。

これは計算順序の違いによって発生したdouble型の精度の問題と考えられます。

 

次は、Javaで跳ね返りありのアニメーション付きで実装してみました。

最近NetBeansのGUIエディタに頼ってばっかりだったので、今回はEclipseでせっせとSwingコンポーネント書きました。

シミュレーション結果は次のようになりました。

尚、反発係数は0.9、時間ステップは0.05、0.5としています。

時間ステップを増加させると、Euler法では誤差が大きくなっていくことがわかります。

一方で、Runge-Kutta法では全くといっていいほど誤差が見られません。

数学的に同一の処理をしているので当然と言えば当然なのですが、Runge-Kutta法の実装には細かいミス多発で結構手こずりました。

また時間があったら、もっと複雑な微分方程式を解いてみて、性能検証してみようと思います。疲れた。