経緯
後輩が研究として、物理シミュレーションして、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でグラフを描画してやるとこんな感じ。
ちゃんと減衰振動してる!
コメントを残す