This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Stepping to negative direction in ODE
- To: gsl-discuss at sources dot redhat dot com
- Subject: Stepping to negative direction in ODE
- From: Timo Laitinen <timolai at green dot srl dot utu dot fi>
- Date: Wed, 3 Oct 2001 13:22:36 +0300 (EEST)
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