用 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 来描述常微分方程组，

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 可以为以下值：

返回一些所用算法的信息

具体的步进计算函数如下

计算 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)
\]