大学の後輩が研究として、物理シミュレーションして、VRで遊べるアプリケーションを作っていました。空気抵抗を考える必要があり、複雑な微分方程式を解かなくてはならないとのことで、大学時代にカオス理論を専門にしていた自分にRunge-Kutta法を教えてほしいと相談されたので、プログラムを作ってみました。
目次
Wikipediaを参照すると以下のように書かれています。
数値解析においてルンゲ＝クッタ法（英: Runge–Kutta method）とは、初期値問題に対して近似解を与える常微分方程式の数値解法に対する総称である。この技法は1900年頃に数学者カール・ルンゲとマルティン・クッタによって発展を見た。
Wikipedia: ルンゲ=クッタ法
端的に言うと初期値と常微分方程式を与えることで、解となる関数の数値を求める方法です。(解となる関数を得られるわけではない)
言語はC++です。
#ifndef RungeKutta_hpp
#define RungeKutta_hpp
#include <stdio.h>
class RungeKutta {
float * result;
int num;
void calculateNextStep(float (*f[])(float t, float *x), float t0, float *x, float tn, int num);
public:
RungeKutta(){};
RungeKutta(float (*f[])(float t, float *x), float t0, float tn, float *x, int div, int num);
float getValue(int step, int index);
};
#endif /* RungeKutta_hpp */#include "RungeKutta.hpp"
float RungeKutta::getValue(int step, int index){
return result[step*num + index];
}
RungeKutta::RungeKutta(float (*f[])(float t, float *x), float t0, float tn, float *x, int div, int _num){
num = _num;
result = new float[div*num];
float h = (tn - t0)/div;
for(int i=0; i<div; ++i){
calculateNextStep(f, t0, x, t0+h, num);
for(int j=0; j<num; ++j){
result[i*num+j] = x[j];
t0 += h;
}
}
}
void RungeKutta::calculateNextStep(float (*f[])(float t, float *x), float t0, float *x, float tn, int num){
float k1[num], k2[num], k3[num], k4[num], tmp[num];
float h = (tn - t0);
float t = t0;
for(int j=0; j<num; j++){
k1[j] = (*f[j])(t, x);
tmp[j] = x[j] + h*k1[j]/2;
k2[j] = (*f[j])(t+h/2, tmp);
tmp[j] = x[j] + h*k2[j]/2;
k3[j] = (*f[j])(t+h/2, tmp);
tmp[j] = x[j] + h*k3[j];
k4[j] = (*f[j])(t+h, tmp);
x[j] += (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])*h/6;
}
}
減衰振動の以下の式をシミュレーション
以下のようにyを定義します。
すると、減衰振動の方程式はこれを代入することで以下のようになります。
このように2つの一回微分方程式に分解することができました。
この二つの関数を_f0, _f1として定義して、Runge-Kuttaにかけます。
#define LOOP 10000
float _f0(float t, float * x){
return x[1];
}
float _f1(float t, float * x){
return - x[0] - 0.2 * x[1];
}
float (*f[2])(float, float*);
f[0] = _f0;
f[1] = _f1;
float initialValues[] = {30.0, 0.0};
RungeKutta rk = RungeKutta(f, 0, 30, initialValues, LOOP, 2);
for(int i=0; i<LOOP; ++i){
cout << rk.getValue(i, 0) << endl;
}変数名が変わっているので分かりにくいですが、_f0, _f1の中のx[0]がx、x[1]がyに相等します。
RungeKuttaクラスのコンストラクタの引数の意味を以下にまとめておきます。
今回のサンプルコードでは2つの初期位置30cm、初速0cm/sとしています。連立常微分方程式を0〜30秒の時間を10000分割(0.003秒の粒度)で、計算するという意味になります。
rk.getValueの引数の意味は以下の通りです。
を表します。
今回はx(x[0])の値が欲しいので、第2引数を0としています。
もし、第2引数を1にすると、y(x[1])の値が取れるので、速度が計算できます。
openFrameworksでグラフを描画してやると以下のようになります。
ちゃんと減衰振動していることが確認できました。
http://hooktail.org/computer/index.php?Runge-Kutta%CB%A1