[Bug-gsl] Bug in odeiv integrator

[Taneli Kalvas]

Hi!

The behaviour of gsl_odeiv_evolve_apply() differs from documentation. 
The documentation says:

"If the user-supplied functions defined in the system dydt return a 
status other than GSL_SUCCESS the step will be aborted. In this case, t 
and y will be restored to their pre-step values and the error code from 
the user-supplied function will be returned."

The y values are restored, but the t value is not. This happens in case, 
where the control function has first rejected a step because of too 
large error and when taking a new step, the user function returns an 
error. See the following print for a simple example:

address@hidden:~/src/gsl_oveiv_bug$ gsl-config --version
1.12
address@hidden:~/src/gsl_oveiv_bug$ gcc --version
gcc (GCC) 4.2.4 (Ubuntu 4.2.4-1ubuntu3)
Copyright (C) 2007 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

address@hidden:~/src/gsl_oveiv_bug$ cat test.c
#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>


#define ERR 201


int dfunc( double t, const double *x, double *dxdt, void *data )
{
    static calls = 0;
    printf( "Evaluate at t=%lf\n", t );

    if( t < 1.0 )
    dxdt[0] = t;
    else
    dxdt[0] = 1.0;

    calls++;
    if( calls == 12 )
    return( ERR );
    return( GSL_SUCCESS );
}


int main( void )
{
    gsl_odeiv_system      system;
    gsl_odeiv_step       *step;
    gsl_odeiv_control    *control;
    gsl_odeiv_evolve     *evolve;

    system.jacobian  = NULL;
    system.params    = NULL;
    system.function  = dfunc;
    system.dimension = 1;

    step    = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, 
system.dimension );
    control = gsl_odeiv_control_standard_new( 1.0e-6, 1.0e-6, 1.0, 1.0 );
    evolve  = gsl_odeiv_evolve_alloc( system.dimension );

    double t    = 0.0;
    double x    = 0.0;
    double dt   = 1.2;
    double maxt = 2.0;
    while( t < 2.0 ) {
    int retval = gsl_odeiv_evolve_apply( evolve, control, step, &system,
                         &t, maxt, &dt, &x );
    if( retval == ERR )
        printf( "Breaking at error: t=%12.6lf x=%12.6lf\n", t, x );
    else
        printf( "Accept: t=%12.6lf x=%12.6lf\n", t, x );
    }


    gsl_odeiv_evolve_free( evolve );
    gsl_odeiv_control_free( control );
    gsl_odeiv_step_free( step );

    return( 0 );
}
address@hidden:~/src/gsl_oveiv_bug$ gcc -lgsl -lgslcblas test.c -o test
address@hidden:~/src/gsl_oveiv_bug$ ./test
Evaluate at t=0.000000
Evaluate at t=0.240000
Evaluate at t=0.360000
Evaluate at t=0.720000
Evaluate at t=1.200000
Evaluate at t=1.050000
Evaluate at t=1.200000
Evaluate at t=0.056886
Evaluate at t=0.085328
Evaluate at t=0.170657
Evaluate at t=0.284428
Evaluate at t=0.248874
Breaking at error: t=    1.200000 x=    0.000000
Evaluate at t=1.200000
Evaluate at t=1.256886
Evaluate at t=1.285328
Evaluate at t=1.370657
Evaluate at t=1.484428
Evaluate at t=1.448874
Evaluate at t=1.484428
Accept: t=    1.484428 x=    0.284428
Evaluate at t=1.484428
Evaluate at t=1.587542
Evaluate at t=1.639100
Evaluate at t=1.793771
Evaluate at t=2.000000
Evaluate at t=1.935553
Evaluate at t=2.000000
Accept: t=    2.000000 x=    0.800000
address@hidden:~/src/gsl_oveiv_bug$

It is easy to get around this problem by resetting the t in case of user 
function errors, so this is easy to fix, but I hope the problem would be 
fixed internally in gsl_odeiv_evolve_apply() to prevent further headaches.

Thanks and best regards,
Taneli Kalvas

--
Taneli Kalvas
M.Sc., Researcher
Physics Department, room FL114
P.O. Box 35 (YFL)
40014 University of Jyväskylä, Finland
Phone: +358-44-314-1602
Fax:   +358-14-260-2351
Email: address@hidden