[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Bug-gsl] rk4

From: Brian Gough
Subject: Re: [Bug-gsl] rk4
Date: Tue, 03 Jul 2007 10:21:01 +0100
User-agent: Wanderlust/2.14.0 (Africa) Emacs/22.1 Mule/5.0 (SAKAKI)

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

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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]