[Top][All Lists]

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

[Bug-gsl] Documentation bug - Ordinary Differential Equations

From: zkoza
Subject: [Bug-gsl] Documentation bug - Ordinary Differential Equations
Date: Wed, 19 May 2004 15:16:37 +0200

There seems to be a bug in the description of stepping functions in the ODE 
part of GSL, see 

For example, the manual says that the function gsl_odeiv_step_rk2 is 

"Embedded 2nd order Runge-Kutta with 3rd order error estimate.".

Well, the error estimate in this routine is computed 
through the standard mid-point Rung-Kutta, which is of order 2
(e.i. its error estimate is O(h^3). At the same time the value 
of the system state (= the solution) at t+h is here calculated via a third-
order formula (of accuracy O(h^4)). 
So I think the correct description of the function should read 

"Embedded 3rd order Runge-Kutta with 2nd order error estimate.".

The same remark refers also to all other embedded Runge-Kutta formulaes 
implemented in GSL:  rkf45, rkck, and rk8pd.

In the example program at 
at the last part 
(after the words "It is also possible to work with a non-adaptive 
integrator..." in the main text)
there are three suspected lines of code:
   double dfdy[4], dydt_in[2], dydt_out[2];

  /* initialise dydt_in */
  GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in);

Now, the array dfdy is never used in the program, actually it CANNOT be used 
in this kind of programming and hence should never appear in an example code.
Alas, as it turns out, actually this array is used -- exactly once -- in the 
above-mentioned line
GSL_ODEIV_JA_EVAL(&sys, t, y, dfdy, dydt_in);

So the array dfdy is being defined, initialized and... forgotten! 
I suggest replacing these lines with the following:
  double dydt_in[2], dydt_out[2];

  /* initialise dydt_in using function and params stored in sys */
  GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);

Now the array dfdy is gone for good. 
And there's a bit longer description of why we chose 
to call the user-supplied function 'func' through 
a starnge-looking macro rather than to call it explicitly.

Zbigniew Koza

reply via email to

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