Адаптивная интеграция GSL выдает ошибку Invalid Pointer - PullRequest
Новая
       103

Адаптивная интеграция GSL выдает ошибку Invalid Pointer

0 голосов
/ 14 июля 2020

Я пытаюсь использовать встроенный метод RK2imp для решения слегка измененной версии задачи Ван дер Пола, приведенной в примерах GSL , для использования адаптивных решателей для интеграции. Однако я продолжаю получать эту «ошибку: неверный указатель». Я прикрепляю сюда свой код 

#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

Я знаю, что могу совершить очень простую (и глупую) ошибку, но для некоторых причина, по которой я не могу этого видеть. Буду очень признателен за любую помощь! Спасибо.

1 Ответ

1 голос
/ 15 июля 2020

Насколько я понимаю, неподходящая пошаговая функция для метода численного интегрирования (Рунге-Кутта). Здесь лучше использовать явные методы, можно даже использовать пошаговую функцию gsl_odeiv2_step_rk2, но решение будет сходиться очень долго. ИЛИ вы можете использовать неявные методы интеграции, но вам необходимо использовать драйвер шага : 

#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;
}
...