シミュレーションについて軽くおさらい
微分方程式を立てて解く
世の中の様々な物理現象は、基本的に微分方程式を立てて、それを解くことで明らかにすることができます。微分方程式はニュートンの運動方程式とかナビエ・ストークス方程式などがあり、式さえ立てれば、後はそれを解くだけです\(^o^)/。
2つの解き方がある
解くための方法として、解析解と数値解の2つあります。解析解は微分方程式を数学的操作のみを使って求めます。ゆえに、厳密な結果(誤差0)を導き出せますな。一方、数値解は、微分方程式を四則演算の関係で表す差分方程式に変換して、数値的に求めます。
数値解ってな~に?
数値解(オイラー法)
数値解の求め方ってなんなの?って思われるかもしれないので説明しておきます。例としてニュートンの運動方程式F = ma を挙げます。ここでは、数値解として変位xを求めてみましょう。以下の流れで求まりますね。
a = F/m → v = v + aΔt → x = x + vΔt ※ ” = ” は右辺を左辺に代入
これをオイラー法と言います。 v = v + aΔtでは、加速度a をΔtの間だけ一定とみなしているので、速度vは時刻tに関して直線的に増減します。ここで解析解との誤差が生じることがわかりますね。x = x + vΔt についても同様です。
計算誤差
複雑な物理現象になると、微分方程式も複雑なものになり解析解を求めることが困難になります。このとき、数値解を使わざるを得ないが、数値解には誤差がつきもの(T_T)。
より精度の高い数値解(ルンゲクッタ法)
数値解における計算誤差を減らしてできるだけ精度の高い計算をしようとするのが、ルンゲクッタ法です。以下のアルゴリズムで求めます。
V(t, x, v) = v , A(t ,x, v) = F(t, x, v)/m
v(1) = V(t, x, v)
a(1) = A(t, x, v)
v(2) = V(t+Δt/2, x+v(1)Δt/2, v+a(1)Δt/2)
a(2) = A(t+Δt/2, x+v(1)Δt/2, v+a(1)Δt/2)
v(3) = V(t+Δt/2, x+v(2)Δt/2, v+a(2)Δt/2)
a(3) = A(t+Δt/2, x+v(2)Δt/2, v+a(2)Δt/2)
v(4) = V(t, x+v(3)Δt, v+a(3)Δt)
a(4) = A(t, x+v(3)Δt, v+a(3)Δt)
x = x + Δt/6 (v(1)+2v(2)+2v(3)+v(4))
v = v + Δt/6 (a(1)+2a(2)+2a(3)+a(4))
いざ、シミュレーション!!
ルンゲクッタでシミュレーションしてみました。Pygameというモジュールを使いました。
う~ん、うまく行かないなぁ。物体がたくさんあるとそれだけ誤差が蓄積されますね。
物体と物体をつなぐ糸の長さはすべて等しいという体で運動方程式を立てたのですが、なんかの伸び縮みしてますね… 。力技で物体間の糸の長さを一定に保つように、PD制御で張力にフィードバックさせたのですが、うまく行かなかったです(TдT)。まぁ、コレは今後の課題ということで(笑)
おまけ