Задача с решением задачи Коши с использованием GSL - PullRequest
Новая
       23

Задача с решением задачи Коши с использованием GSL

0 голосов
/ 10 июня 2019

Я пытаюсь решить проблему Коши, используя GSL.У меня есть примитивная функция с одним параметром.Я думаю, что проблема может быть в моих параметрах, но ошибка 

gsl: driver.c:356: ERROR: integration limits and/or step direction not consistent
Default GSL error handler invoked.

выдает в любое время.

Вот код. 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
// GSL lib includes
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>


struct Dots {
    double par;
    double x;
    double y;
};

int ode_func (double x, const double y[], double f[], void *params)
{
    double mu = *(int *)params;
    f[0] = (x + 2 * y[0]) / (1 + mu * mu);
    return GSL_SUCCESS;
}

void calc_cauchy_problem(double x_start, double x_end, double y_start,
                         int count, int param1, int param2) {
    int dim = 1;
    double x = x_start;
    double y[1] = {y_start};
    int param = param1;
    int status = 0;

    for (param = param1; param <= param2; param++) {
        gsl_odeiv2_system sys = {ode_func, NULL, dim, &param};
        gsl_odeiv2_driver * d =
                gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rkf45, 1e-6, 1e-6, 0.0);

        for (int i = 1; i <= count; i++) {
            double xi = x_start + i * (x_end - x_start) / count;

            int status = gsl_odeiv2_driver_apply(d, &x, xi, y);
            if (status != GSL_SUCCESS)
            {
                printf ("error, return value=%d\n", status);
                break;
            }
        }
        gsl_odeiv2_driver_free (d);
    }
}

int main() {
    double start_time = omp_get_wtime();
    double x_start = 0;
    double x_end = 2;
    double y_start = 0;
    const int count = 50;
    int param1 = 0;
    int param2 = 10;
    calc_cauchy_problem(x_start, x_end, y_start, count, param1, param2);
    printf("Elapsed time = %f\n", omp_get_wtime() - start_time);
    return 0;
}

В настоящее время нет манипуляцийс этими данными.Я обнаружил, что если у меня есть только один цикл цикла for (param = param1; param <= param2; param++) {, он работает нормально.Проблема появляется во втором цикле цикла в этой строке int status = gsl_odeiv2_driver_apply(d, &x, xi, y);

...