GSLで常微分方程式の初値問題を解く

GSLで常微分方程式の初値問題を解く
GNU Scientific Library(GSL)は科学計算用のC言語クラスライブラリである.ここでは,GSLを用いて常微分方程式の初期値問題の計算を行う方法を紹介する.
n次元の常微分方程式群は一般に[frac{{dy_i(t)}{{dx}=f_i(t,y_1(t),ldots y_n(t)]i=1,ldots,nに対応する初期条件がy_i (t_0) = Y_{i0}
GSLは,Runge‐Kutta法やBulirsch‐Stoer法のような下位層の方法を含む常微分方程式の初期問題を解決するための多くの方法を提供し,適応ステップ長制御などの上位層のアルゴリズムも含む.
GSLにおける1つの構造体gsl_odeiv_システムは常微分方程式のグループを記述し、

typedef struct  
{
  int (* function) (double t, const double y[], double dydt[], void * params);
  int (* jacobian) (double t, const double y[], double * dfdy, double dfdt[], void * params);
  size_t dimension;
  void * params;
}
gsl_odeiv_system;

functionは,ユーザ定義関数,t,y[],paramsを入力パラメータとして指す関数ポインタであり,dydtはn次元のベクトル場f_を返す.i(t,y,params)、計算に成功した場合、関数はGSL_を返すべきであるSUCCESS;
JAcobianも関数ポインタで、ユーザー定義の関数、t,y[],paramsを入力パラメータとして、dfdtはpartial f_を返します.i/partial t,dfdyはJacobian行列J_を返す{ij},記憶方式はJ(i,j)=dfdy[i*dimension+j]であり,計算に成功した場合,関数はGSL_を返すべきである.SUCCESS;
いくつかの簡単なアルゴリズムに対してJacobian行列を必要としないで、この時jacobianを空(NULL)に指向することができます.しかし、多くの効率的なアルゴリズムはJacobian行列を使うので、Jacobian行列を提供することができる時できるだけ提供するべきです.
dimensionはシステムの次元数を記述します.
paramsは任意のポインタであり,具体的な役割は以下の例で説明する.
以下、GSLで提供される関数について3つの部分に分けて説明します.
1ステップ関数
ステップ関数はユーザが呼び出すことができる最下層の機能関数であり,他の上位関数もステップ関数を呼び出すことによって最後の計算作業を完了する.
必要なスペースの割り当てと解放に使用

gsl_odeiv_step * gsl_odeiv_step_alloc(const gsl_odeiv_step_type * T, size_t dim);
int  gsl_odeiv_step_reset(gsl_odeiv_step * s);
void gsl_odeiv_step_free(gsl_odeiv_step * s);

ここでgsl_odeiv_step_typeは次の値になります.

gsl_odeiv_step_rk2;    //   Runge-Kutta  
gsl_odeiv_step_rk4;    //     Runge-Kutta  
gsl_odeiv_step_rkf45;  //Runge-Kutta-Fehlberg  
gsl_odeiv_step_rkck;  //Runge-Kutta Cash-Karp
gsl_odeiv_step_rk8pd; //Runge-Kutta Prince-Dormand
gsl_odeiv_step_rk2imp; //    Runge-Kutta 
gsl_odeiv_step_rk2simp;
gsl_odeiv_step_rk4imp; //Implicit 4th order Runge-Kutta at Gaussian points
gsl_odeiv_step_bsimp; //    Bulirsch-Stoer method of Bader and Deuflha （   Jacobian   ）
gsl_odeiv_step_gear1; //M=1 implicit Gear method
gsl_odeiv_step_gear2; //M=2 implicit Gear method

使用するアルゴリズムの情報を返します

const char * gsl_odeiv_step_name(const gsl_odeiv_step *);
unsigned int gsl_odeiv_step_order(const gsl_odeiv_step * s);

具体的なステップ計算関数は次のとおりです.

int  gsl_odeiv_step_apply(gsl_odeiv_step *s, double t, double h, double y[], double yerr[], const double dydt_in[], double dydt_out[], const gsl_odeiv_system * dydt);

計算t->t+h,t+hの時刻の値はy[]に格納され,yerr[]はy[]の計算誤差の推定値である.dydt_In[]nullではなくdydt_In[]はdydt_が提供されていない場合、t時刻の時間に対する導関数であるべきである.in[]の場合、この関数の内部で対応する計算が行われます.dydt_out[]t+h時刻対時間の導関数を出力する.
2次非線形Van der Pol方程式を考慮すると、[x'(t)+mu x'(t)(x(t)^2−1)+x(t)=0]という方程式は、次のような常微分方程式群に簡略化できる.
y_1' = -y_0 +\mu y_1 (1 - {y_0}^2) y_0' = y_1
Jacobian行列は[J=left({begin{array}{*{20}c}0&1{-1-2muy_0 y_1}&{mu(1-{y_0}^2)}\end{array}right)]

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

int func (double t, const double y[], double f[], void *params)
{//        
 double mu = *(double *)params;
 f[0] = y[1];
 f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
 return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
{
 double mu = *(double *)params;
 gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2);
 gsl_matrix * m = &dfdy_mat.matrix;
 gsl_matrix_set (m, 0, 0, 0.0);
 gsl_matrix_set (m, 0, 1, 1.0);
 gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
 gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
 dfdt[0] = 0.0;
 dfdt[1] = 0.0;
 return GSL_SUCCESS;
}

int main (void)
{
 const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);

 double mu = 10;
 gsl_odeiv_system sys = {func, jac, 2, &mu};
 double t = 0.0, t1 = 100.0;
 double h = 1e-6; //     h
 double y[2] = { 1.0, 0.0 }; //  
 double y_err[2];
 double dydt_in[2], dydt_out[2];
 
 GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in); //       dy/dt
 while (t < t1)
 {
  int status = gsl_odeiv_step_apply (
  s, t, h, y, y_err, dydt_in, dydt_out, &sys);
  if (status != GSL_SUCCESS) break;
  dydt_in[0] = dydt_out[0];
  dydt_in[1] = dydt_out[1];
  t += h;
  printf ("%.5e %.5e %.5e
", t, y[0], y[1]);
 }
 gsl_odeiv_step_free (s);
 return 0;
}

プログラムは簡単で、これ以上説明する必要はありません.この方法の欠点は積分ステップ長が自分で決定する必要があり,ステップ長が大きすぎると計算結果に大きな誤差があり,ステップ長が小さすぎると効率が低下し,ステップ長が小さすぎると誤差がかえって増大する可能性があることである.
したがって、一般的には、より高度なアルゴリズム、すなわち適応ステップアルゴリズムを選択します.
2適応ステップ制御

gsl_odeiv_control * gsl_odeiv_control_standard_new(double eps_abs, double eps_rel, double a_y, double a_dydt);
gsl_odeiv_control * gsl_odeiv_control_y_new(double eps_abs, double eps_rel);
gsl_odeiv_control * gsl_odeiv_control_yp_new(double eps_abs, double eps_rel);
gsl_odeiv_control * gsl_odeiv_control_scaled_new(double eps_abs, double eps_rel, double a_y, double a_dydt, const double scale_abs[], size_t dim);

これらの関数は大同小異で,いずれも1つのステップ関数のステップ後の誤差を検査し,ユーザが設定した誤差レベルと比較し,ステップ長を調整するかどうかを決定する.具体的な各パラメータの意味はマニュアルを参照してください.

gsl_odeiv_control * gsl_odeiv_control_alloc(const gsl_odeiv_control_type * T);

新しいステップ制御関数にスペースを割り当てるには、新しいステップ制御関数をカスタマイズしない限り、一般的なユーザーはあまり使用しません.

int gsl_odeiv_control_init(gsl_odeiv_control * c, double eps_abs, double eps_rel, double a_y, double a_dydt);

ステップ制御関数の初期化

void gsl_odeiv_control_free(gsl_odeiv_control * c);

対応するスペースを解放

int gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s, const double y0[], const double yerr[], const double dydt[], double * h);

この関数は、適切なステップhを計算するために使用され、ステップ長が大きくなるとGSLを返す.ODEIV_HADJ_INC、ステップサイズはGSLに戻るODEIV_HADJ_DEC、ステップ長が調整されていない場合はGSL_に戻るODEIV_HADJ_NIL.

const char * gsl_odeiv_control_name(const gsl_odeiv_control * c);

ステップ制御関数の名前を返します.
以下、上述したVan der Pol方程式を例にとる.

int main (void)
{
 const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
 gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-5, 0.0);
 double mu = 10;
 gsl_odeiv_system sys = {func, jac, 2, &mu};
 double t = 0.0, t1 = 100.0;
 double h = 1e-2; //     h
 double *ph = &h;
 double y[2] = { 1.0, 0.0 }; //  
 double y_const[2];
 double y_err[2];
 double dydt_in[2], dydt_out[2];
 int status;
 
 GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in); //       dy/dt
 while (t < t1)
 {
  gsl_odeiv_step_apply (s, t, h, y_const, y_err, dydt_in, dydt_out, &sys);
  status = gsl_odeiv_control_hadjust (c, s, y_const, y_err, dydt_in, ph);
  gsl_odeiv_step_reset (s);
  
  status = gsl_odeiv_step_apply (s, t, h, y, y_err, dydt_in, dydt_out, &sys);
  
  y_const[0] = y[0];
  y_const[1] = y[1];
    
  if (status != GSL_SUCCESS) break;
  dydt_in[0] = dydt_out[0];
  dydt_in[1] = dydt_out[1];
  t += h;
  printf ("%.5e %.5e %.5e
", t, y[0], y[1]);
 }
 gsl_odeiv_control_free(c);
 gsl_odeiv_step_free (s);
 return 0;
}

3進化アルゴリズム
進化アルゴリズムは最上位層のアルゴリズムであり,ステップ関数とステップ長制御関数を結合している.ステップ長を自動的に調整できます.
必要なスペースの割り当てと解放に使用

gsl_odeiv_evolve * gsl_odeiv_evolve_alloc (size t dim);
int gsl_odeiv_evolve_reset (gsl odeiv evolve * e);
void gsl_odeiv_evolve_free (gsl odeiv evolve * e);

int gsl_odeiv_evolve_apply (gsl odeiv evolve * e, gsl odeiv control * con, gsl odeiv step * step, const gsl odeiv system * dydt, double * t, double t1, double * h, double y[]);

この関数は次のステップの値を自動的に計算し、計算中にhを調整する可能性があるが、t+h<=t 1を保証するため、複数回計算すると正確にt 1に達する.

int main (void)
{
 const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 2);
 gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);
 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);
 double mu = 10;
 gsl_odeiv_system sys = {func, jac, 2, &mu};
 double t = 0.0, t1 = 100.0;
 double h = 1e-6;
 double y[2] = { 1.0, 0.0 };
 while (t < t1)
 {
  int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
  if (status != GSL_SUCCESS) break;
  printf ("%.5e %.5e %.5e
", t, y[0], y[1]);
 }
 gsl_odeiv_evolve_free (e);
 gsl_odeiv_control_free (c);
 gsl_odeiv_step_free (s);
 return 0;
}

微分方程式の解は、特定の時刻の値(例えば、等時間間隔の一連の時刻の値)を必要とすることがよくあります.この場合、進化アルゴリズムを簡単に採用するだけではだめです.ステップ長を制御できないためです.次のコードは、解決策を示しています.

for (i = 1; i <= 100; i++)
{
 double ti = i * t1 / 100.0;
 while (t < ti)
 {
  gsl_odeiv_evolve_apply (e, c, s, &sys, &t, ti, &h, y);
 }
 printf ("%.5e %.5e %.5e
", t, y[0], y[1]);
}