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