[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Bug in odeiv integrator
From: |
Taneli Kalvas |
Subject: |
[Bug-gsl] Bug in odeiv integrator |
Date: |
Wed, 15 Apr 2009 13:53:31 +0300 |
User-agent: |
Thunderbird 2.0.0.21 (X11/20090318) |
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