ode.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 *                                                                           *
 *                       / /   / /    __    \ \  / /                         *
 *                      / /__ / /   / _ \    \ \/ /                          *
 *                     /  ___  /   | |/_/    / /\ \                          *
 *                    / /   / /    \_\      / /  \ \                         *
 *                                                                           *
 *                         Jakub Benda (c) 2014                              *
 *                     Charles University in Prague                          *
 *                                                                           *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef HEX_ODE
#define HEX_ODE

#include <gsl/gsl_odeiv2.h>

#define DIVERGENCE_THRESHOLD    100

#define ABORT_ON_OVERFLOW       0
#define RETURN_ON_OVERFLOW      1
#define NORMALIZE_ON_OVERFLOW   2

#include "misc.h"

int solve2 (
    double xg[], int N, double h, double yg[][2],
    int (*derivs) (double, const double [2], double [2], void*),
    void * data, int flag
) {
    // initialize auxiliary variables
    size_t ntrial = 100000000;  // 10⁷
    double x0 = xg[0], x1 = xg[N - 2];
    double x = x0, xnext;
    int status = GSL_SUCCESS, first_status = GSL_SUCCESS;
    size_t nsteps = 0;
    xg[0] = x0;
    
    // setup the stepper
    const gsl_odeiv2_step_type * step_type = gsl_odeiv2_step_rkck;
    gsl_odeiv2_step * step = gsl_odeiv2_step_alloc(step_type, 2);
    gsl_odeiv2_control * control = gsl_odeiv2_control_y_new(1e-6, 0.0);
    gsl_odeiv2_evolve * evolve = gsl_odeiv2_evolve_alloc(2);
    
    // setup the system
    gsl_odeiv2_system sys;
    sys.function = derivs;
    sys.jacobian = nullptr;
    sys.dimension = 2;
    sys.params = data;
    
    // for all steps
    for(int i = 1; i < N - 1 and status == GSL_SUCCESS; i++)
    {
        xnext = x0 + (x1 - x0) * ((double)i)/((double)(N - 2));
        double (& y_row) [2] = yg[i];
        y_row[0] = yg[i-1][0];
        y_row[1] = yg[i-1][1];
        
        // Step until we reach the next grid point
        bool done = false;
        while (not done and status == GSL_SUCCESS)
        {
            status = gsl_odeiv2_evolve_apply (
                evolve,
                control,
                step,
                &sys,
                &x, xnext,
                &h,
                y_row
            );
            
            // check finiteness
            if (not std::isfinite(y_row[0]))
                throw exception("[solve2] Infinite result (%g) for i = %d", y_row[0], i);
            
            // advance number of steps
            nsteps++;
            if (status != GSL_SUCCESS and first_status != GSL_SUCCESS)
                first_status = status;
            if (nsteps > ntrial)
                std::cerr << "[solve2] Too many steps required (ntrial=" << ntrial << ", x=" << x << ").";
            
            // avoid overflow (normalize)
            if (fabs(y_row[0]) > DIVERGENCE_THRESHOLD)
            {
                if (flag == NORMALIZE_ON_OVERFLOW)
                {
                    double inverse_norm = 1./fabs(y_row[0]);
                    for (int j = 0; j <= i; j++)
                    {
                        yg[j][0] *= inverse_norm;
                        yg[j][1] *= inverse_norm;
                    }
                }
                if (flag == RETURN_ON_OVERFLOW)
                {
                    return i - 1;
                }
                if (flag == ABORT_ON_OVERFLOW)
                {
                    throw exception ("[solve2] Overflow error.");
                }
            }
            
            // decide if we have reached the grid point
            if (x1 > x0)
            {
                if (x >= xnext)
                    done = true;
            }
            else
            {
                if (x <= xnext)
                    done = true;
            }
            
        }
    }
    
    gsl_odeiv2_evolve_free(evolve);
    gsl_odeiv2_control_free(control);
    gsl_odeiv2_step_free(step);
    
    return N;
}

#endif