Re: [Bug-gsl] rk4

[Brian Gough]

At Mon, 2 Jul 2007 15:30:23 -0400 (EDT),
Brian Powell wrote:
> Thanks in advance for your help.  I've attached a tar.gz file with a
> sample code and makefile.  Simply run "make test".  The code solves a
> system of ODEs and outputs the solution in the file "out.dat"  The final
> column is the time variable.
> 
> I'm interested in seeing if you can recreate the problem.  You should get
> different evolutions under gsl-1.5 and gsl-1.8.

Thanks for the program.  I can see there is a difference (below).

I believe it's because the RK4 stepper error estimates were improved
in version 1.6.  In version 1.5 the error was overestimated, which
meant the overall solution was actually more accurate than it should
be for the same tolerance, because smaller steps were taken by the
solver.

If the tolerance is reduced to gsl_odeiv_control_y_new(1e-8,1e-8) the
two solutions converge reasonably at t=0 (strict example below).

If you use a stricter tolerance do you still see the same problem of
the forward and backward evolution not matching up?

-- 
Brian Gough
(GSL Maintainer)

Network Theory Ltd,
Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/


==> base/out-1.5.dat <==
5.498446e+00 1.097511e-04 1.558111e-01 -1.352890e+00 2.331965e-01 -6.074385e-02 
-2.438367e-03 3.734914e-03 0.000000

==> base/out-1.6.dat <==
5.472987e+00 1.147616e-04 1.199053e-01 -1.081635e+00 1.844766e-01 -4.511818e-02 
-1.017724e-03 2.115210e-03 0.000000

==> strict/out-1.5.dat <==
5.498446e+00 1.097511e-04 1.558111e-01 -1.352890e+00 2.331965e-01 -6.074385e-02 
-2.438368e-03 3.734914e-03 0.000000

==> strict/out-1.6.dat <==
5.497661e+00 1.099080e-04 1.545413e-01 -1.343615e+00 2.315072e-01 -6.018281e-02 
-2.378950e-03 3.669031e-03 0.000000