GSL ODE求解器返回-nan但同参数同ODE在Python中可正常求解

核心问题1：求解器绝对误差配置错误 #

你调用gsl_odeiv2_driver_alloc_y_new初始化求解器驱动时，最后一个参数（绝对误差）设置为0.0，这会直接导致求解器的误差控制逻辑失效，步长调整异常，迭代过程中很容易出现y值下溢为0甚至负数的情况，进而触发pow函数的域错误，返回NaN。
修复方式：将绝对误差改为适配你y值量级的小正数，例如1e-12，修改后的初始化代码如下：

// 原代码最后一个参数为0.0，修改为1e-12
gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 1e-12);

核心问题2：幂函数计算无边界保护 #

你的导数表达式中指数(m-1)/m ≈ -2.33为负数，当迭代过程中y[0]/(k*pow(s,n))趋近于0时，pow计算会返回无穷大，进一步导致y值计算异常。
修复方式：在计算幂函数前增加底数的最小值保护，修改func函数代码如下：

int func (double t, const double y[], double f[], void *params)
{
  (void)(t); /* avoid unused parameter warning */

  struct param_type *my_params_pointer = (param_type *)params;
  double k  = my_params_pointer->k;
  double n  = my_params_pointer->n;
  double m  = my_params_pointer->m;
  double s  = my_params_pointer->s;

  double base = y[0]/(k*pow(s,n));
  // 增加底数保护，避免趋近于0导致的无穷大
  if (base < 1e-15) base = 1e-15;
  f[0] = m*k*pow(s,n)*pow(base, (m-1)/m);

  return GSL_SUCCESS;
}

其他非致命问题 #

1. 全局声明的int * jac;没有任何作用，属于冗余代码可以直接删除，不会影响求解逻辑。
2. 你在循环结束后调用gsl_vector_add(time, fun_val)将时间数组和函数值数组相加后写入文件，如果你需要单独存储时间和求解结果，这部分业务逻辑需要调整，不过这不是求解返回NaN的原因。

修改完成后重新编译运行即可得到和Pythonscipy.integrate.odeint一致的结果。

内容的提问来源于stack exchange，提问作者Muhammad Mohsin Khan