GSL 自适应积分给出无效指针错误

Question

我正在尝试使用内置的 RK2imp 方法来解决 GSL sample examples 中给出的范德波尔问题的略微修改版本，以使用自适应求解器进行集成。但是，我不断收到此“错误：无效指针”。 我在这里附上我的代码

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

#define ABTOL 1e-18

#define A 2.00861986087484313650940188
#define mu 1.0
#define TEND 6.6632868593231301896996820305

double eps_abs=ABTOL, eps_rel=0.0, tend=TEND;

int func (double t, const double y[], double f[], void *params)
{
    (void)(t); /* avoid unused parameter warning */
    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)
{
    (void)(t); /* avoid unused parameter warning */
    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_odeiv2_step_type * T = gsl_odeiv2_step_rk2imp;

    gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, 2);
    gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
    gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (2);

    gsl_odeiv2_system sys = {func, jac, 2, NULL};

    double t = 0.0, t1 = TEND;
    double h = 1e-6;
    double y[2] = { A, 0.0 };

    while (t < t1)
    {
        int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &h, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error: %s\n", gsl_strerror (status));
            break;
        }

        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

    gsl_odeiv2_evolve_free (e);
    gsl_odeiv2_control_free (c);
    gsl_odeiv2_step_free (s);
    return 0;
}
 

这是我得到的输出：

 error: invalid pointer

我知道我可能犯了一个非常简单（而且很愚蠢）的错误，但由于某种原因我看不到它。我真的很感激任何帮助！谢谢。

Answer 1

据我了解，不恰当的 stepping function 用于数值积分方法 (Runge-Kutta)。这里最好使用显式方法，甚至可以使用gsl_odeiv2_step_rk2步进函数，但解会收敛很长时间。 或者你可以使用隐式集成方法，但是你需要使用step driver：

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

#define ABTOL 1e-18

#define A 2.00861986087484313650940188
#define mu 1.0
#define TEND 6.6632868593231301896996820305

double eps_abs=ABTOL, eps_rel=0.0, tend=TEND;

int func (double t, const double y[], double f[], void *params)
{
    (void)(t); /* avoid unused parameter warning */
    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)
{
    (void)(t); /* avoid unused parameter warning */
    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_odeiv2_step_type * T = gsl_odeiv2_step_msadams;
    // with "gsl_odeiv2_step_rk2imp" stepping function the solution
    // converges VERY SLOWLY and may be "failure" error at the end

    gsl_odeiv2_step * s = gsl_odeiv2_step_alloc (T, 2);
    gsl_odeiv2_control * c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
    gsl_odeiv2_evolve * e = gsl_odeiv2_evolve_alloc (2);

    gsl_odeiv2_system sys = {func, jac, 2, NULL};

    gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(&sys, T, 1e-6, 1e-6, 1e-6 );
    gsl_odeiv2_step_set_driver(s, d);

    double t = 0.0, t1 = TEND;
    double h = 1e-6;
    double y[2] = { A, 0.0 };

    while (t < t1)
    {
        int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &h, y);

        if (status != GSL_SUCCESS)
        {
            printf ("error: %s\n", gsl_strerror (status));
            break;
        }

        printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

    gsl_odeiv2_evolve_free (e);
    gsl_odeiv2_control_free (c);
    gsl_odeiv2_step_free (s);
    return 0;
}