用 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 中用一个结构体 gsl_odeiv_system 来描述常微分方程组，

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 中提供的函数。

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[] 应为 t 时刻的对时间的导数，如果没有提供 dydt_in[] 的话，这个函数的内部会作相应的计算。dydt_out[] 输出t + h 时刻对时间的导数。

下面是一个简单的例子：
考虑二阶的非线性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 - 2\mu y_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\n", 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);

这几个函数大同小异，都是检查一个步进函数步进后的误差，与用户设定的误差水平进行比较，然后决定是否要调整步长。具体各个参数的含义请看手册。

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\n", 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 <= t1，因此多次计算后会精确的到达t1。

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\n", 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\n", t, y[0], y[1]);
}