Blame ode-initval/TODO

Packit 67cb25
# -*- org -*-
Packit 67cb25
#+CATEGORY: ode-initval
Packit 67cb25
Packit 67cb25
* Add a higher level interface which accepts a 
Packit 67cb25
Packit 67cb25
start point, 
Packit 67cb25
end point, 
Packit 67cb25
result array (size N,  y0, y1, y2 ... ,y(N-1))
Packit 67cb25
desired relative and absolute errors epsrel and epsabs
Packit 67cb25
Packit 67cb25
it should have its own workspace which is a wrapper around the
Packit 67cb25
existing workspaces
Packit 67cb25
Packit 67cb25
* Implement other stepping methods from well-known packages such as
Packit 67cb25
RKSUITE, VODE, DASSL, etc
Packit 67cb25
Packit 67cb25
* Roundoff error needs to be taken into account to prevent the
Packit 67cb25
step-size being made arbitrarily small
Packit 67cb25
Packit 67cb25
* The entry below has been downgraded from a bug.  We use the
Packit 67cb25
coefficients given in the original paper by Prince and Dormand, and it
Packit 67cb25
is true that these are inexact (the values in the paper are said to be
Packit 67cb25
accurate 18 figures).  If somebody publishes exact versions we will
Packit 67cb25
use them, but at present it is better to stick with the published
Packit 67cb25
versions of the coefficients them use our own.
Packit 67cb25
----------------------------------------------------------------------
Packit 67cb25
BUG#8 -- inexact coefficients in rk8pd.c
Packit 67cb25
Packit 67cb25
From: Luc Maisonobe <Luc.Maisonobe@c-s.fr>
Packit 67cb25
To: gsl-discuss@sources.redhat.com
Packit 67cb25
Subject: further thoughts about Dormand-Prince 8 (RK8PD)
Packit 67cb25
Date: Wed, 14 Aug 2002 10:50:49 +0200
Packit 67cb25
Packit 67cb25
I was looking for some references concerning Runge-Kutta methods when I
Packit 67cb25
noticed GSL had an high order one. I also found a question in the list
Packit 67cb25
archive (April 2002) about the references of this method which is
Packit 67cb25
implemented in rk8pd.c. It was said the coefficients were taken from the
Packit 67cb25
"Numerical Algorithms with C" book by Engeln-Mullges and Uhlig.
Packit 67cb25
Packit 67cb25
I have checked the coefficients somewhat with a little java tool I have
Packit 67cb25
developped (see http://www.spaceroots.org/archive.htm#RKCheckSoftware)
Packit 67cb25
and found they were not exact. I think this method is really the method
Packit 67cb25
that is already in rksuite (http://www.netlib.org/ode/rksuite/) were the
Packit 67cb25
coefficients are given as real values with 30 decimal digits. The
Packit 67cb25
coefficients have probably been approximated as fractions later on.
Packit 67cb25
However, these approximations are not perfect, they are good only for
Packit 67cb25
the first 16 or 18 digits depending on the coefficient.
Packit 67cb25
Packit 67cb25
This has no consequence for practical purposes since they are stored in
Packit 67cb25
double variables, but give a false impression of beeing exact
Packit 67cb25
expressions. Well, there are even some coefficients that should really
Packit 67cb25
be rational numbers but for which wrong numerators and denominators are
Packit 67cb25
given. As an example, the first and fourth elements of the b7 array are
Packit 67cb25
given as 29443841.0 / 614563906.0 and 77736538.0 / 692538347, hence the
Packit 67cb25
sum off all elements of the b7 array (which should theoretically be
Packit 67cb25
equal to ah[5]) only approximate this. For these two coefficients, this
Packit 67cb25
could have been avoided using  215595617.0 / 4500000000.0 and
Packit 67cb25
202047683.0 / 1800000000.0, which also looks more coherent with the
Packit 67cb25
other coefficients.
Packit 67cb25
Packit 67cb25
The rksuite comments say this method is described in this paper :
Packit 67cb25
Packit 67cb25
   High Order Embedded Runge-Kutta Formulae
Packit 67cb25
   P.J. Prince and J.R. Dormand
Packit 67cb25
   J. Comp. Appl. Math.,7, pp. 67-75, 1981
Packit 67cb25
Packit 67cb25
It also says the method is an 8(7) method (i.e. the coefficients set
Packit 67cb25
used to advance integration is order 8 and error estimation is order 7).
Packit 67cb25
If I use my tool to check the order, I am forced to check the order
Packit 67cb25
conditions numerically with a tolerance since I do not have an exact
Packit 67cb25
expression of the coefficients. Since even if some conditions are not
Packit 67cb25
mathematically met, the residuals are small and could be below the
Packit 67cb25
tolerance. There are tolerance values for which such numerical test
Packit 67cb25
dedeuce the method is of order 9, as is said in GSL. However, I am not
Packit 67cb25
convinced, there are to few parameters for the large number of order
Packit 67cb25
conditions needed at order 9.
Packit 67cb25
Packit 67cb25
I would suggest to correct the coefficients in rk8pd.c (just put the
Packit 67cb25
literal constants of rksuite) and to add the reference to the article.
Packit 67cb25
Packit 67cb25
----------------------------------------------------------------------