タケウチCarlの地平線

都内在住31歳の独身サラリーマンが、日々木工や読書、散歩などを楽しみつつ、いつか脱サラして小屋暮らしや旅暮らしをすることを夢見るブログ

ルンゲ=クッタ法①~ルンゲ=クッタ法を使うと物体の運動をコンピュータ上でシミュレーションすることができます

ルンゲ=クッタ法と微分方程式

ルンゲ=クッタ法は微分方程式の初期値問題に対して近似解を与える一連の方法です。


さて、物体の運動はニュートン方程式F=maによって近似的に表されます。ここでFは力、mは質量、aは加速度であり、質量mの物体に力が働いた際にどのような加速度を受けるのか、反対に加速度があるときにどのような力が働いているのか、といった「力と加速度の関係」を表す方程式です。


そして、加速度aとは位置xの時間tによる2階微分\frac{dx^2}{d^2t})ですが、このように、微分された関数による方程式を一般に微分方程式と呼んでいます。


ということは、ルンゲ=クッタ法を使えばニュートン方程式の近似解を求めることができるのですから、物体の運動をコンピュータ上でシミュレーションすることができます。これは工作をする上で何かに役立ちそうです、というか本当のところを言うと、地震が起こったときに本棚がどのように動くのかをシミュレーションしたいがためにルンゲ=クッタ法を調べたのでした。


微分方程式の解を数値的に計算する方法はたくさんありますが、ルンゲ=クッタ法は計算公式がシンプルでプログラムに実装しやすく精度もよいため、一般的に使われているとのことです。

ルンゲ=クッタ法の計算公式*1

ここでは微分方程式にyの関数やyの時間微分しか出てこない場合に、ルンゲ=クッタ法の計算公式を書いてみます。要は「一変数の場合」ということですね~。


実際にシミュレーションするときには、縦×横の2次元空間や、縦×横×高さの3次元空間の中の物体の運動をシミュレーションすることが多そうなので、これだけでは実用に足らないかもしれません、多変数の場合は次回にします。


さて、初期条件と時刻tでの位置yの時間変化がf(t,y)で与えらている微分方程式を考えます。数式で書くと
\frac{dy}{dt}=f(t,y),y(t_0)=y_0
ですねー。


ルンゲ=クッタ法の計算公式では時間の刻み幅をhとして、ある時刻の位置y_nから、次の時刻の位置y_{n+1}を近似的に計算するのですが、それは
y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)
と計算されます。


ここで、k_1k_4は時刻t_nとそのときの位置y_nからfによって計算される値であり、
k_1=f(t_n,y_n)
k_2=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1)
k_3=f(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2)
k_4=f(t_n+h,y_n+k_3)
です。

補足説明

  • k_1は時刻t_n位置y_nでの変化率です。
  • k_2は次の時間\frac{h}{2}の変化率をk_1で近似した場合の、時刻t_n+\frac{h}{2}での変化率です。
  • k_3は次の時間\frac{h}{2}の変化率をk_2で近似した場合の、時刻t_n+\frac{h}{2}での変化率です。
  • k_4は次の時間hの変化率をk_3で近似した場合の、時刻t_n+hでの変化率です。

これら4種類の変化率の加重平均(\frac{k_1+2k_2+2k_3+k_4}{6})をとることで、時間hの「正しい変化率」を近似計算しているわけです。

ルンゲ=クッタ法は精度がいい

ルンゲ=クッタ法から計算される物体の位置y_nはあくまで近似解ですので、厳密な物体の位置y_nからは徐々にずれてきます。しかしながら、両者をテイラー展開して比較すると4次の項まで一致していることが確かめられます。

(続く)

*1:古典的な4次の