C++
openFrameworks
数学
ルンゲクッタ法

Runge-Kutta法で微分方程式を解いたらちょっと感動した話

More than 3 years have passed since last update.


Runge-Kutta法で微分方程式を解いたらちょっと感動した話


経緯

後輩が研究として、物理シミュレーションして、VRで遊べるアプリケーションを作っていた。空気抵抗とか結構考えないといけないみたいで、複雑な微分方程式を解かなきゃならないっぽい。んで、大学時代にカオスとかやってた自分に「Runge-Kutta法」教えて!みたいに来たので、プログラムを作ってみた。


Runge-Kutta法ってなんじゃい!

http://hooktail.org/computer/index.php?Runge-Kutta%CB%A1

↑を参照。プログラムもここを大いに参考にした。


プログラム

言語はC++です。

RungeKutta.hpp

#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 */

RungeKutta.cpp

#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つの1階微分方程式に分解できるので、それらを_f0, _f1として定義して、RungeKuttaにかける。

#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) + ofGetHeight()/2 << endl;
}


結果

openFrameworksでグラフを描画してやるとこんな感じ。

スクリーンショット 2015-12-15 12.45.13.png

ちゃんと減衰振動してる!


拙ブログからの転載