Re: Stepping to negative direction in ODE

Re: Stepping to negative direction in ODE

[Brian Gough]

Timo Laitinen writes:
 >  I found a problem with the step size evolution mechanism in the
 > functions
 >  gsl_odeiv_evolve_apply and std_control_hadjust
 >  The problem is that they cannot handle negative step sizes, ie
 > they cannot be used for solving an equation backwards in time (or
 > in space!!!). My solution for the problem is to change in
 >  evolve.c, line 147

Thanks, I'll make those changes.

 >  There is also another problem in function gsl_odeiv_evolve_apply.
 > If the initial h0 is larger than the remaining time dt, the
 > final_step flag is set, on line 150. BUT, if the step is now
 > decreased by gsl_odeiv_control_hadjust, line 182, and the new step
 > h0<dt, the final_step flag is NOT RESET!. Thus, you end up taking a
 > small step *h, and updating the time *t to t1, which would give
 > strange results.

Right, you may want to grab the latest version from CVS as I checked
in a fix for this recently (also some improvements to the BSIMP algorithm).
Here it is,

Index: evolve.c
===================================================================
RCS file: /cvs/gsl/gsl/ode-initval/evolve.c,v
retrieving revision 1.11
retrieving revision 1.12
diff -r1.11 -r1.12
21c21
<  * RCS:     $Id: evolve.c,v 1.11 2001/06/25 10:27:45 bjg Exp $
---
>  * RCS:     $Id: evolve.c,v 1.12 2001/09/29 20:34:09 bjg Exp $
150a151,154
>     }
>   else
>     {
>       final_step = 0;


best regards
Brian Gough

Stepping to negative direction in ODE

[Timo Laitinen]

Hi,

I found a problem with the step size evolution mechanism in the functions 

gsl_odeiv_evolve_apply
and
std_control_hadjust

The problem is that they cannot handle negative step sizes, ie they cannot
be used for solving an equation backwards in time (or in space!!!). My
solution for the problem is to change in 

evolve.c, line 147

   if (h0 > dt)
to
   if (GSL_SIGN(h0)*h0 > GSL_SIGN(dt)*dt)

and in cstd.c, lines 109 and 118,

     *h = GSL_MAX_DBL(0.2 * h_old, h_new);
to 
     *h = (double) GSL_SIGN(*h)*GSL_MAX_DBL(0.2 * h_old, h_new); 

and

     *h = GSL_MIN_DBL(5.0 * h_old, h_new);

     *h = (double) GSL_SIGN(*h)*GSL_MIN_DBL(5.0 * h_old, h_new); 

, respectively. There may be a better/faster way of doing this so think
first (offhand, i can't remember why I chose not to use fabs in the first
one...)

There is also another problem in function gsl_odeiv_evolve_apply.
If the initial h0 is larger than the remaining time dt, the final_step
flag is set, on line 150. BUT, if the step is now decreased by
gsl_odeiv_control_hadjust, line 182, and the new step h0<dt, the
final_step flag is NOT RESET!. Thus, you end up taking a small step *h,
and updating the time *t to t1, which would give strange results.

Of course one should take care of not starting with a step size larger
than the integration period, but accidents happen (as in my case :-) ),
and tracking the problem is a real pain...

Otherwise it seems to work fine, thank you for a great work, and for 
saving my time considerably :-)

----------------
Timo                      FABRICATI DIEM, PVNK.
email:timo.laitinen@utu.fi