経緯

後輩が研究として、物理シミュレーションして、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つの一回微分方程式に分解できるので、それらを_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

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