[Help-gsl] ODE solvers and overdamped motion

[Jonny Taylor]

I am having a few problems with the GSL ODE solvers, and I was hoping  
somebody might be able to save me a lot of puzzling over this. I  
suspect that I am using the wrong tool for the job, or misusing the  
tool, but it's not obvious to me what I need to do.

My problem involves a multi-particle system, with complicated inter- 
particle forces, but the issue I am having can be seen in a simple  
overdamped spring calculation (the derivatives function simply returns  
dydt[i] = -y[i] * 9.0). I use the RK45 stepper, and error control  
gsl_odeiv_control_y_new(1e-2, 0). My motivation for this choice of  
error control is that I the motion I am interested in - and indeed the  
length scale of the interactions in the real system - is on the scale  
of a wavelength (normalized to 1). GSL is correctly respecting that  
error bound, but is not meeting my naive expectations of how I'd hoped  
it might behave.

The most obvious example of this is when the particle approaches its  
equilibrium point (zero in this simple case). Yes, I have specified an  
absolute error tolerance of 1e-2, but I would nevertheless have hoped  
that it might get closer to zero over time. This doesn't happen (it  
bounces around apparently randomly within +-1e-2). Again, I realise I  
cannot complain about this - it is following my specification.  
However, looking at the points at which it samples the velocity during  
each timestep (whether or not it is close to equilibrium) it looks  
like it is consistently overestimating the motion and then correcting  
itself. This makes me wonder if there is something I can do, or a  
different stepper I can use, that will improve the performance. I  
wonder if fundamentally the problem is that the motion involved is  
exponential decay, and so any sort of polynomial fit will always tend  
to over-estimate it.

The two things I would like to improve, if possible, in the  
performance are:
1. Greater inherent accuracy in the estimation, allowing larger  
timesteps (and therefore fewer samples) for a given error tolerance.
2. If possible, motion that gets closer to the "real" equilibrium  
point as time goes on

Can anybody suggest anything I could try to improve this?
Thanks in advance
Jonny