StepperDopr5
クラスが提供されている
(クラスの中身についてはNR3, 17.2.3の項を参照)。
これをドライバルーチンOdeint
と組み合わせることで
「ステップ幅の適応制御」付きでODEを解くことが出来る。
以下では
4次のルンゲ・クッタ法で解いたものと同じvdP方程式(12)を、
5次
ドルマンド・プリンス法によって解くことにする(
とする)。
#include "../code/nr3.h" #include "../code/stepper.h" #include "../code/stepperdopr5.h" #include "../code/odeint.h" class RHSvdP { // 「vdP方程式の右辺」という意味(名前は何でも良い) public: Doub mu; RHSvdP(Doub muu) : mu(muu) {} // 初期化リストを用いたコンストラクタ void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx) { dydx[0]= y[1]; dydx[1]= -y[0] + mu * y[1] * ( 1 - y[0]*y[0]); } }; Int main(void) { const int N=2; const Doub atol=1.0e-3, // 絶対許容誤差 rtol=atol, // 相対許容誤差 h1=0.01, // 初期ステップ hmin=0.0, x1=0.0, //積分範囲[x1, x2] x2=50.0; VecDoub ystart(N); ystart[0] = 1.0; // 初期値 ystart[1] = 0.0; Output out(-1); // out(-1)にした場合は、実際に使用したステップで出力 // out(20)とした場合は[x1, x2]の間を20等分した位置での値を出力 RHSvdP d(10.0); // μ=10.0 でvdP方程式の右辺を定義 Odeint<StepperDopr5<RHSvdP> > ode(ystart, x1, x2, atol, rtol, h1, hmin, out, d); ode.integrate(); //ここで積分している 出力: cout << "x\ty[0]\ty[1]\n"; for(Int i=0;i<out.count;i++){ cout << out.xsave[i] << "\t"; cout << out.ysave[0][i] << "\t"; cout << out.ysave[1][i] << endl; } return 0; }
これによって得られる結果を図2に示した。 図1と異なり、ステップ幅が一定でないことに注意。
|
なお、NR3では8次のドルマンド・プリンス法(最終点における誤差の見積もりにおいてのみ5次と3次の方法を使用している)のステッパルーチンStepperDopr853
クラスも提供されており[4]、上記プログラムにおいてインクルードするヘッダファイルをstepperdopr853.h
に、ステッパルーチン名をStepperDopr853
にそれぞれ変更すればそのまま実行出来る(実際上はこちらを利用すると良い)。
※GSL [5] でも同様のルーチンが存在するので試されたい