【发布时间】：2021-11-27 16:33:00
【问题描述】：
我使用 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;
}
正如here 所提到的，我不需要 jacobian 作为显式求解器，这就是我为 jac 函数传递 NULL 指针的原因。 当我运行上面的代码时，我在所有时间点都得到了 -nan 值。
为了交叉检查，我用 python 编写了具有相同功能和相同参数的程序，使用 scipy.integrate.odeint 解决。它运行并提供了一个合理的答案。 按照我的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++ 代码中做错了什么？
【问题讨论】：
-
据我所知，您将驱动程序用作步长相当大的固定步长求解器。 python odeint 使用自适应步长，必要时更小，即使输出位于间隔较宽的点。
-
感谢您的回复。我还尝试为驱动程序使用更小的步长。我仍然得到 nan 作为回报（现在不是 -nan ）。但我仍然对它如何不起作用感到困惑。
-
@MuhammadMohsinKhan 这个问题与python无关。工作演示中使用的那个 python 并不意味着你使用那个标签。请阅读How to Ask 并查看tour
-
@eyllanesc 感谢您的编辑。我现在正在阅读这两个链接。