GSL ODE 求解器 returns -nan 尽管在 python 中求解具有相同参数的相同 ODE

Question

我使用 python 使用 scipy.integrate.odeint 求解 ODE。目前，我正在做一个小项目，我在 C++ 中使用 gsl 来解决 ODE。我正在尝试求解 ODE，但求解器为每个时间点返回 -nan。以下是我的代码：

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

struct param_type {
    double k;
    double n;
    double m;
    double s;
};


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;

  f[0] = m*k*pow(s,n)*pow((y[0]/(k*pow(s,n))),(m-1)/m);

  return GSL_SUCCESS;
}


int * jac;


int main ()
{
  struct param_type mu = {1e-07, 1.5, 0.3, 250};
  gsl_odeiv2_system sys = {func, NULL, 1, &mu};

  gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
  int i;
  double t = 0.0, t1 = 10.0;
  double step_size = 0.1;
  double y[1] = { 1e-06 };

  gsl_vector *time = gsl_vector_alloc ((t1 / step_size) + 1);
  gsl_vector *fun_val = gsl_vector_alloc ((t1 / step_size) + 1);

  for (i = 1; i <= t1/step_size; i++)
    {
      double ti = i * t1 / (t1 / step_size);
      int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

      if (status != GSL_SUCCESS)
        {
          printf ("error, return value=%d\n", status);
          break;
        }

      printf ("%.5e %.5e\n", t, y[0]);
      gsl_vector_set (time, i, t);
      gsl_vector_set (fun_val, i, y[0]);
    }

  gsl_vector_add(time, fun_val);

  {
     FILE * f = fopen ("test.dat", "w");
     gsl_vector_fprintf (f, time, "%.5g");
     fclose (f);
  }

  gsl_odeiv2_driver_free (d);
  gsl_vector_free (time);
  gsl_vector_free (fun_val);

  return 0;
}

如前所述，我不需要 jacobian 作为显式求解器，这就是我为 jac 函数传递 NULL 指针的原因。当我运行上面的代码时，我在所有时间点都得到 -nan 值。

为了交叉验证，我在python中编写了具有相同功能和相同参数的程序，使用scipy.integrate.odeint求解。它运行s 并提供了一个似是而非的答案。按照我的 python 代码：

import numpy as np
from scipy.integrate import odeint

def nb(y, t, *args):
    k = args[0]
    n = args[1]
    m = args[2]
    s = args[3]
    return m*k*s**n*(y/(k*s**n))**((m-1)/m)

t = np.linspace(0,10,int(10/0.1))
y0 = 1e-06

k = 1e-07
n = 1.5
m = 0.3
s = 250

res = odeint(nb, y0, t, args=(k,n,m,s)).flatten()
print(res)

谁能帮我弄清楚，我在使用 GSL 求解 ODE 的 C++ 代码中做错了什么？

Answer 1

你的问题在这里：

f[0] = m*k*pow(s,n)*pow((y[0]/(k*pow(s,n))),(m-1)/m);

随着求解器的进行，它可能想要对 y[0] 的负值进行采样。在 Python 中这没有问题，在 C++ 中它产生 NAN。

要处理这个问题，您可以模仿 Python 的行为：

  auto sign = (y[0] < 0) ? -1.0 : 1.0;
  f[0] = sign*m*k*pow(s,n)*pow((std::abs(y[0])/(k*pow(s,n))),(m-1)/m);

甚至将 sign 有效设置为 1:

  f[0] = m*k*pow(s,n)*pow((std::abs(y[0])/(k*pow(s,n))),(m-1)/m);

毕竟，将负值提高到非整数次方是错误的，除非考虑复数，而事实并非如此。

请注意 y[0] 使用 std::abs 保护。

GSL ODE 求解器 returns -nan 尽管在 python 中求解具有相同参数的相同 ODE

GSL ODE solver returns -nan although same ODE with same parameters is being solved in python

c++

ode

gsl