Re: GSL 1.4: BUG #8 [rk8pd.c]

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Luc Maisonobe]

Hello,

The remaining problem is still the same as one year ago: the comments
say the method is an 8(9) method (i.e. extrapolation beeing of order 8
and error estimation of order 9), but it is really an 8(7) method.
Despite I think nobody is interested anymore, I promised to check the
coefficients of the method, so here are my results.

The maximal error in rk8pd.c is about 3.92e-17.
The maximal error in rksuite.f is about 4.49e-29.

These values were computed by picking two coefficients from rksuite
(c(9) and a(8,4)) as the remaining free parameters of the model,
transforming them into rational numbers by continued fractions method
and using them to solve exactly all the equations of conditions and all
the additional equations from the article. The exact result is available
in http://www.spaceroots.org/software/rkcheck/DormandPrinceRKsuite.xml

                                                      Luc

Re: GSL 1.4: BUG #8 [rk8pd.c]

[M Joonas Pihlaja]

On Sat, 8 Nov 2003, Luc Maisonobe wrote:

[snip]
> Despite I think nobody is interested anymore, I promised to check the
> coefficients of the method, so here are my results.

I'm interested, at least, so thank you for sharing.

> The maximal error in rk8pd.c is about 3.92e-17.
> The maximal error in rksuite.f is about 4.49e-29.

What do you mean by maximal error here -- of the main estimator
itself or the second embeded estimator?  And how do you compute
these -- from the principal error function or some other method?
I'm very curious, because you say the results for the
coefficients are exact after fixing two free parameters.

Best Regards,

Joonas Pihlaja

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Luc Maisonobe]

Hello Joonas,


> > The maximal error in rk8pd.c is about 3.92e-17.
> > The maximal error in rksuite.f is about 4.49e-29.
> 
> What do you mean by maximal error here -- of the main estimator
> itself or the second embeded estimator?

It is simply the absolute value of the difference between the
coefficients I obtained and the coefficients in the source files.


>  And how do you compute
> these -- from the principal error function or some other method?
> I'm very curious, because you say the results for the
> coefficients are exact after fixing two free parameters.

The Dormand and Prince article roughly explains the method. You first
have a (somewhat huge) set of equations of conditions, the first ones
are very simple like 

\sum_{i=1}^{i=s} b_{i} = 1

The last ones are much more complex like

\sum_{i=6}^{i=s}\left(b_{i}
  \sum_{j=5}^{j=i-1}{\left(a_{i,j}
    \sum_{k=4}^{k=j-1}{\left(a_{j,k}
      \sum_{l=3}^{l=k-1}{\left(a_{k,l}
        \sum_{m=2}^{m=l-1}{\left(a_{l,m}
          \sum_{p=1}^{p=m-1}{\left(a_{m,p}
                               c_{p}
                             \right)}
                           \right)}
                         \right)}
                       \right)}
                     \right)}
                \right) = \frac{1}{5040}

For one order 8 method and one order 7 embedded method, there are 284
equations. to this you can add the traditional homogeneity equation :

  \sum_{j=1}^{j=s}a_{i,j} = c_{i}

In order to help solving these equations, Dormand and Prince added their
own set like for example:

  \sum_{j=1}^{j=i-1} a_{i,j} c_{i}^k = \frac{c_{i}^{k+1}}{k+1}

for the first (k,i) pairs. This greatly reduces the number of equations
of conditions since most of them become equivalent to other ones.

These equations bind most of the variables and the global model has only
10 degrees of freedom left (Dormand and Prince chose c_2, c_3, c_6, c_7,
c_8, c_10, c_11, \hat{b}_13, b_12 and a_{8,4}). All other variables can
be expressed from these ones.

I followed their method, but already had available the values of most of
the free variables as provided by their paper (c_8 = 93/200, c_10 =
13/20 ...). The only variables that were not obvious were c_11 and
a_{8,4}, so I had to chose these ones.

Despite the article lists c_11 as a free parameter and hence c_9 appears
only has a dependant parameter, I think from observing the float values
in rksuite.f that c_11 was computed from c_9. This is because c_9 as
exactly the same value (5490023248/9719169821) as in the paper whereas
c_11 is slightly different from the one published (difference is
2.33e-18). They are tightly related to each other by the following
simple relation which can be derived from the equations of conditions
and the other parameters:

 c_11 = (10009221 c_9 - 8842703) / (11614350 c_9 - 10009221)

So I opted for c_9 = 5490023248 / 9719169821

The only remaining parameter was a_{8,4}, I simply used a continued
fraction approximation to retrieve it from the rksuite floating value.

All other values have been computed from these ones.

Changing the free parameters lead to other 8(7) methods. The choice of
the free parameters was done by Dormand and Prince in order to minimize
the truncation error (i.e. the first equations of conditions that are
*not* verified). Having chosen (almost) the same parameters as in
rksuite mean I did not redo this minimiation step. In fact, I tried to
do it for a{8,4}, but it occurred to me that the first few largest
residuals of the 9th order equations are symetrical (i.e. the largest
negative is exactly the opposite of the largest positive, and also the
second largest and third largest) and do not depend on a{8,4}.

If you want, I can provide you some intermediate results, like the
solving steps in Maxima syntax (Maxima is the free version of good old
Macsyma computer algebra system) and the one degree of freedom model
where only a{8,4} is still unbound.

Dormand and Prince did not provide exact values because (quoting the
paper) :

"it proved too expensive computationally to retain exact rational
parameters as in the previous case."

However, what was "too expensive computationally" in 1980 is now easy
for any home PC (I did all that on my personal PC running GNU/Linux and
only freeware tools).

                                                     regards,
                                                       Luc

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Jerome BENOIT]

M Joonas Pihlaja wrote:
> On Sat, 8 Nov 2003, Luc Maisonobe wrote:
> 
> [snip]
> 
>>Despite I think nobody is interested anymore, I promised to check the
>>coefficients of the method, so here are my results.
> 
> 
> I'm interested, at least, so thank you for sharing.

I am still interested because I currently using it:
the rk4 series did not pass my test whereas the rk8pd
passed it. Note that my test is rather specific to what
I am doing: the test just checks the conservstion of
one quatity (namely the energy).
Nevertheless I was confused by the BUG report,
hence my email to gsl-discuss.
So any precision is welcome.

> 
> 
>>The maximal error in rk8pd.c is about 3.92e-17.
>>The maximal error in rksuite.f is about 4.49e-29.
> 
> 
> What do you mean by maximal error here -- of the main estimator
> itself or the second embeded estimator?  And how do you compute
> these -- from the principal error function or some other method?
> I'm very curious, because you say the results for the
> coefficients are exact after fixing two free parameters.
> 
> Best Regards,
> 
> Joonas Pihlaja
> 

Thanks,
Jerome BENOIT

GSL 1.4: BUG #8 [rk8pd.c]

[Jerome BENOIT]

Hello GSL,

what is the current status of the BUG #8 ?

Is there any available patch which the coefficients ?

Thanks in advance,
Jerome BENOIT

PS:
note that the Paper

High order embedded Runge-Kutta formulae
P. J. Princea and J. R. Dormanda
J. Comp. Appl. Math., 7 (1981), pp 67-76.

is now available at Science Direct:

http://dx.doi.org/10.1016/0771-050X(81)90010-3

(but unfortunately I cannot get it :()

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Gerard Jungman]

On Sun, 2003-11-02 at 05:57, Jerome BENOIT wrote:

> note that the Paper
> 
> High order embedded Runge-Kutta formulae
> P. J. Princea and J. R. Dormanda
> J. Comp. Appl. Math., 7 (1981), pp 67-76.
> 
> is now available at Science Direct:
> 
> http://dx.doi.org/10.1016/0771-050X(81)90010-3
> 
> (but unfortunately I cannot get it :()

We have a subscription, but I don't seem to have
access to electronic versions before 1995. So
I will wander over to the library and see if
it looks helpful. As you pointed out, the easiest
thing to do might be to grab the numbers out of
rksuite. I wouldn't trust anything printed in the
journal, unless somebody out there can vouch for
the accuracy. Has anybody listening here
checked the numbers in the journal?

Thanks.

-- 
Gerard Jungman <jungman@lanl.gov>
Los Alamos National Laboratory

Re: GSL 1.4: BUG #8 [rk8pd.c]

[MAISONOBE Luc]

Gerard Jungman a écrit :

> We have a subscription, but I don't seem to have
> access to electronic versions before 1995. So
> I will wander over to the library and see if
> it looks helpful. As you pointed out, the easiest
> thing to do might be to grab the numbers out of
> rksuite. I wouldn't trust anything printed in the
> journal, unless somebody out there can vouch for
> the accuracy. Has anybody listening here
> checked the numbers in the journal?

Since my last message on the subject 
(http://sources.redhat.com/ml/gsl-discuss/2002-q4/msg00020.html) I have 
check some of the coefficients.

As I said in that message, the article clearly says these numbers are 
approximations. The values from rksuite.f in Netlib are more accurate 
(Dormand and Prince seem to have provided better values to Brankin, 
Gladwell and Shampine).

I have tried to check these values using my rkcheck tool 
(http://www.spaceroots.org/software/rkcheck/index.html), and corrected 
some of them. Since Benoit's message yesterday, I have restarted this 
work and corrected more terms. What I mean here is *only* to replace 
float numbers with exact rational numbers, in order to check exactly the 
order. This is only of theoretical interest, for someone as picky as me 
;-) the values provided in the article are accurate to about 17 digits 
(just have a look to the first coefficients I corrected in my message 
last year) and the values in rksuite are accurate to about 30 digits. My 
  rational number are *not* better than this ones, they only form a 
globally consistent set, using some values as free parameters (for 
example 5490023248/9719169821 which is exactly the same in the article 
and in rksuite) and computing the other ones by solving the conditions 
equations and the additional equations from the article). The rational 
numbers I have computed so far (all the vector elements and some matrix 
elements) sometime involve really huge integers (more than 130 digits).

Again, I don't think my work is of *practical* interest. It will only be 
an external validation (if I succeed ...). I will post the xml file with 
the exact coefficients when it will be ready (for the moment, the one 
you can find at my site only contains the coefficients I had last year, 
i.e. the float values from rksuite and the two corrected rationals from 
my patch).

I am not subscribed to gsl list anymore, so this message will probably 
bounce. If either Gerard or Benoit think it could interest the list, 
please forward it by yourself.

Apart for Benoit :

I have brought the article to work this morning, if you want a copy, 
send me directly a postal adress or a fax number, and I will be happy to 
send you one.

                                                         Luc

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Brian Gough]

Jerome BENOIT writes:
 > Hello GSL,
 > 
 > what is the current status of the BUG #8 ?

Although it is listed as a "bug" the original coefficients are
accurate to 16 digits, which should be good enough for most practical
purposes.

-- 
Brian Gough

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Jerome BENOIT]

Thank you very much for your replies.

Brian Gough wrote:
> Jerome BENOIT writes:
>  > Hello GSL,
>  > 
>  > what is the current status of the BUG #8 ?
> 
> Although it is listed as a "bug" the original coefficients are
> accurate to 16 digits, which should be good enough for most practical
> purposes.
> 

May be this should be added in the BUG file to avoid confussion.

Thanks,
Jerome

Re: GSL 1.4: BUG #8 [rk8pd.c]

[Brian Gough]

Jerome BENOIT writes:
 > Brian Gough wrote:
 > > 
 > > Although it is listed as a "bug" the original coefficients are
 > > accurate to 16 digits, which should be good enough for most practical
 > > purposes.
 > > 
 > 
 > May be this should be added in the BUG file to avoid confussion.
 > 

It is there already -- the comments in the file give the full
description.