こんにちは、Rafkaです。
今回は、今まで紹介してきた数値解法の誤差比較を、斜方投射を題材に行っていきます。
問題設定
斜方投射は以下の図のように、よく高校で習うような一般的な形でモデル化しました。
この時の位置に関する厳密解は初期位置が \(\boldsymbol{p}(0) = [p_x(0), p_y(0)]^\text{T} = [0, 0]^\text{T}\) ならば以下のようになりますね。
$$\begin{align}p_x(t) &= v_0t\cos(\theta)\\p_y(t) &= v_0t\sin(\theta)-\frac{1}{2}gt^2\end{align}$$
各数値解法による値と、この厳密解の値とを比較することで、どの程度誤差が発生しているのかを見ていきます。
オイラー法等で特には微分方程式で表現されたモデルが必要なので、厳密解が分かっているこの問題ですが、今回はあえて、ニュートンの運動方程式から出発して微分方程式で表してみます。
ニュートンの運動方程式より、以下の関係が成り立ちます。
$$m\boldsymbol{a}(t) = \boldsymbol{F}(t) = [0, -mg]^\text{T}$$
ここで、\({}^\text{T}\) は行列の転置を表しています。上で時間 \(t\) での位置を \(\boldsymbol{p}(t)\) と表したので、加速度 \(\boldsymbol{a}(t)\) について以下の式が成り立ちます。
$$ \boldsymbol{a}(t) = \frac{d^2\boldsymbol{p}(t)}{dt^2} = \frac{d\boldsymbol{v}(t)}{dt}$$
ですので、上で示した関係式は以下のように変形できます。
$$\frac{d^2\boldsymbol{p}(t)}{dt^2} = [0, -g]^\text{T}$$
ちゃっかり \(m\) で割ってます。
このままでは２階微分の方程式になってしまうので、\(\frac{d\boldsymbol{p}(t)}{dt} = \boldsymbol{v}(t)\) であることを利用して、以下のように１階微分の方程式に変形します。
$$\left\{\begin{align}\frac{d\boldsymbol{v}(t)}{dt} &= [0, -g]^\text{T}\\\frac{d\boldsymbol{p}(t)}{dt} &= \boldsymbol{v}(t) = [v_x(t), v_y(t)]^\text{T}\end{align}\right.$$
更にここで、\(\boldsymbol{x}(t) = [\boldsymbol{v}(t), \boldsymbol{p}(t)]^\text{T} = [v_x(t), v_y(t), p_x(t),p_y(t)]^\text{T}\) とすることで、以下のように一纏めにしてしまいます。
$$\frac{d\boldsymbol{x}(t)}{dt} = [0, -g, v_x(t), v_y(t)]^\text{T}$$
これで微分方程式を用いた斜方投射のモデル化が完了しました。
プログラム
前回の最後にプログラムを一新するかもと書きましたが、一新しました。ついでに本記事シリーズで使用するプログラムを纏めてGitHubに上げました。
一応決定版の数値解法クラスのソースコードを以下に示します。
変更点は以下の通りです。
- 各解法毎に実装していたループの共通部分を一つのメソッドに分離しました。
- その分離したメソッドの内部で、指定された解法に応じて更新メソッドを選択するようにしました。
- １ステップの更新量を計算するメソッドはPrivateだったものを、ユーザが１ステップ毎に独自の処理を挟めるようにPublicにアクセスできるようにしました。
この解法プログラムを用いて、斜方投射のシミュレーションを行うために作成したプログラムは以下の通りです。
１０行目から２４行目に定義されているクラスが、斜方投射の数理モデルを表すクラスです。NumericalModelWithODEクラスを継承し、１６行目で先程示した微分方程式をFunctionメソッドに実装しました。また、今回は厳密解との比較を行うので、時間 \(t\) での厳密解を返すメソッドを ExactSolutions として実装しました。１３行目から１５行目はこのクラスのコンストラクタの定義ですが、初期値の設定を行っており、コンストラクタで値を渡さなければ、\(v_0 = 5.0 \text{m}/\text{s},\;\theta = 45^\circ\) として初期値を設定するようにしています。
今回は複数の解法を組み合わせてプロットするので、各解法の値を同じファイルに出力できるように、１ステップの更新量を計算するメソッドを使用しました。実際にシミュレーションを行いながらファイルに出力を行っているのは、３６行目から４６行目のソースコードになります。
実行結果
実行した結果をグラフにプロットしました。今回のグラフは”再生”ボタンを押すことで、各時間までの軌跡を描きながらアニメーションする機能や、サイドバーを動かして指定した時間までの軌跡を描画することができるので、是非遊んでみてください。
結果を見てみると、オイラー法の結果は厳密解である赤線から大きく外れてしまっている事が分かります。一方で、ホイン法やルンゲ＝クッタ法の結果は厳密解と殆ど同じ軌跡を辿っており、これらの手法の精度の高さを確認することができました。ホイン法とルンゲ＝クッタ法の精度の差を見ることが難しい結果となってしまいましたが、これは時間刻みをより荒くした時に顕著に確認できるようになると思われます。
まとめ
今回は、今まで紹介してきた数値計算手法の誤差の比較を、斜方投射を題材に行いました。
オイラー法は単純な手法なだけあって誤差が大きくなってしまうことや、ホイン法やルンゲ＝クッタ法はオイラー法に比べて誤差を抑えられることが確認できました。
今回の題材は、簡単な物理で厳密解もあっさりと分かってしまうものだったので、次回は数値計算でなければ解くのが難しい現象を題材にシミュレーションを行い、その後、脳の数理モデルの話に移りたいと思います。
次回 → 数値シミュレーション －二重振り子－