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: GSL 1.4: BUG #8 [rk8pd.c]


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



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