next up previous
: ブリルシュ・ストア法 : ステップ幅の適応制御 : 結果の出力

計算例

NR3においては、 5次ドルマンド・プリンス法を実装したステッパルーチンStepperDopr5クラスが提供されている (クラスの中身についてはNR3, 17.2.3の項を参照)。 これをドライバルーチンOdeintと組み合わせることで 「ステップ幅の適応制御」付きでODEを解くことが出来る。 以下では 4次のルンゲ・クッタ法で解いたものと同じvdP方程式(12)を、 5次 ドルマンド・プリンス法によって解くことにする( $ \mu =10,0, y(0)=1.0, z(0)=y^\prime (0)=0.0$ とする)。

#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と異なり、ステップ幅が一定でないことに注意。

図: vdP方程式(14)の5次ドルマンド・プリンス法による数値計算結果. $ \mu =10,0, y(0)=1.0, z(0)=y^\prime (0)=0.0$ としたときの結果. (A)数値解$ y[0](x)$ . (B) $ x=[17,20]$ の拡大図. ステップ幅が変動していることが分かる. (C) ステップ幅の変化.
(A) (B)
\includegraphics[width=8cm]{vdP_stepDopr5.eps} \includegraphics[width=8cm]{vdP_stepDopr5_inset.eps}
(C)
\includegraphics[width=8cm]{vdP_stepDopr5_stepsize.eps}  

なお、NR3では8次のドルマンド・プリンス法(最終点における誤差の見積もりにおいてのみ5次と3次の方法を使用している)のステッパルーチンStepperDopr853クラスも提供されており[4]、上記プログラムにおいてインクルードするヘッダファイルをstepperdopr853.hに、ステッパルーチン名をStepperDopr853にそれぞれ変更すればそのまま実行出来る(実際上はこちらを利用すると良い)。

※GSL [5] でも同様のルーチンが存在するので試されたい



ykagawa 平成20年7月29日