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