用 GSL 解常微分方程初值问题
来源:互联网 发布:阿里云飞天系统 编辑:程序博客网 时间:2024/04/27 14:16
用 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-Karpgsl_odeiv_step_rk8pd; //Runge-Kutta Prince-Dormandgsl_odeiv_step_rk2imp; //隐式的 Runge-Kutta gsl_odeiv_step_rk2simp;gsl_odeiv_step_rk4imp; //Implicit 4th order Runge-Kutta at Gaussian pointsgsl_odeiv_step_bsimp; //隐式的 Bulirsch-Stoer method of Bader and Deuflha (需要 Jacobian 矩阵)gsl_odeiv_step_gear1; //M=1 implicit Gear methodgsl_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]);}
- 用 GSL 解常微分方程初值问题
- 用 GSL 解常微分方程初值问题
- 求解常微分方程初值问题之Runge_Kutta法
- 求解常微分方程初值问题之Runge_Kutta_Fehlberg法
- 简单的常微分方程(组)初值问题
- 四阶龙格-库塔法求解常微分方程的初值问题
- 求解常微分方程初值问题之多变量Runge_Kutta_Gill法
- 求解常微分方程初值问题之多步Euler预报-校正法
- 求解常微分方程初值问题之Milne预报-校正法
- 四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序
- GSL微分方程
- Matlab解常微分方程
- 求解常微分方程初值问题之改进Euler法:预报-校正公式
- Eular方法解常微分方程
- 常微分方程数值解上机
- Python3 SciPy解常微分方程 用Matplotlib演示
- 常微分方程
- 常微分方程
- 栈的应用举例:数制转换,表达式求值
- doxygen filter
- 男士必须收藏:男士健身方案
- jQuery操作数组的方法
- linux 计算时间 time
- 用 GSL 解常微分方程初值问题
- 使用GIMP使图片背景透明化
- JAVA动态绑定详解
- 推荐算法学习(1)
- Android中 检查网络连接状态的变化,无网络时跳转到设置界面
- c语言对时间的处理函数和计时的实现
- spring mvc annotation 的例子
- p2p打洞原理
- window7 上装SecureCRT,连接小鸡