This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: bug list summary


Brian Gough wrote:

> BUG#8 -- inexact coefficients in rk8pd.c

I finally got the reference article from Dormand and Prince (many thanks
again Joonas). In fact, the article itself gives these false
approximations and explicitly says "they were computed using a precision
of about 24 significant digits and the rationals presented in table 2
are continued fractions accurate to 18 significant figures".

I still think there is an error in the first and fourth elements of the
b7 array, but this error was in the original article. I also think the
real coefficients given in rksuite are more accurate, and since they
thank Dormand and Prince for their help, I think they got the original
coefficients from the authors. I don't know if all 30 significant
figures given in rksuite are good or if only the first 24 ones are good.

I will check this as time permit. For the moment, I propose the
following patch, which only add some comments (corected order of the
method, article reference) and correct two coefficients. I have just
written it now to answer this mail, not tried it (and even no compiled
it, so beware of typos ...)

                                                         Luc
--- gsl-1.2/ode-initval/rk8pd.c~	Mon Nov 19 22:39:32 2001
+++ gsl-1.2/ode-initval/rk8pd.c	Wed Oct  9 09:49:26 2002
@@ -17,7 +17,12 @@
  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  */
 
-/* Runge-Kutta 8(9), Prince-Dormand */
+/* Runge-Kutta 8(7), Prince-Dormand
+ * The original article is :
+ *  High order embedded Runge-Kutta formulae
+ *  P.J. Prince and J.R. Dormand
+ *  Journal of Computational and Applied Mathematics, volume 7, no. 1, 1981
+ */
 
 /* Author:  G. Jungman
  */
@@ -80,14 +85,21 @@
 static const double b4[] = { 1.0 / 32.0, 0.0, 3.0 / 32.0 };
 static const double b5[] = { 5.0 / 16.0, 0.0, -75.0 / 64.0, 75.0 / 64.0 };
 static const double b6[] = { 3.0 / 80.0, 0.0, 0.0, 3.0 / 16.0, 3.0 / 20.0 };
+
+/* the b7[0] and b7[3] coefficients were corrected by Luc Maisonobe
+ * since the approximation given in the original article seemed to be
+ * wrong, because the sum b7[0]+b7[1]+ ... +b7[5] was not equal to
+ * ah[5] as it should be.
+ */
 static const double b7[] = {
-  29443841.0 / 614563906.0,
+  215595617.0 / 4500000000.0, /* was 29443841/614563906 in the article */
   0.0,
   0.0,
-  77736538.0 / 692538347.0,
+  202047683.0 / 1800000000.0, /* was 77736538/692538347 in the article */
   -28693883.0 / 1125000000.0,
   23124283.0 / 1800000000.0
 };
+
 static const double b8[] = {
   16016141.0 / 946692911.0,
   0.0,

Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]