|
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 |
----------------------------------------------------------------------
|