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
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
