* GSL 1.4: BUG #8 [rk8pd.c] @ 2003-11-02 12:57 Jerome BENOIT 2003-11-03 19:42 ` Gerard Jungman 2003-11-05 15:26 ` Brian Gough 0 siblings, 2 replies; 10+ messages in thread From: Jerome BENOIT @ 2003-11-02 12:57 UTC (permalink / raw) To: gsl-discuss; +Cc: Luc.Maisonobe 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 :() ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-02 12:57 GSL 1.4: BUG #8 [rk8pd.c] Jerome BENOIT @ 2003-11-03 19:42 ` Gerard Jungman 2003-11-04 11:27 ` MAISONOBE Luc 2003-11-05 15:26 ` Brian Gough 1 sibling, 1 reply; 10+ messages in thread From: Gerard Jungman @ 2003-11-03 19:42 UTC (permalink / raw) To: jgmbenoit; +Cc: gsl-discuss, Luc.Maisonobe 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-03 19:42 ` Gerard Jungman @ 2003-11-04 11:27 ` MAISONOBE Luc 0 siblings, 0 replies; 10+ messages in thread From: MAISONOBE Luc @ 2003-11-04 11:27 UTC (permalink / raw) To: Gerard Jungman; +Cc: jgmbenoit, gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-02 12:57 GSL 1.4: BUG #8 [rk8pd.c] Jerome BENOIT 2003-11-03 19:42 ` Gerard Jungman @ 2003-11-05 15:26 ` Brian Gough 2003-11-05 15:32 ` Jerome BENOIT 1 sibling, 1 reply; 10+ messages in thread From: Brian Gough @ 2003-11-05 15:26 UTC (permalink / raw) To: jgmbenoit; +Cc: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-05 15:26 ` Brian Gough @ 2003-11-05 15:32 ` Jerome BENOIT 2003-11-05 16:52 ` Brian Gough 0 siblings, 1 reply; 10+ messages in thread From: Jerome BENOIT @ 2003-11-05 15:32 UTC (permalink / raw) To: Brian Gough; +Cc: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-05 15:32 ` Jerome BENOIT @ 2003-11-05 16:52 ` Brian Gough 0 siblings, 0 replies; 10+ messages in thread From: Brian Gough @ 2003-11-05 16:52 UTC (permalink / raw) To: jgmbenoit; +Cc: gsl-discuss 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. ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] @ 2003-11-08 12:51 Luc Maisonobe 2003-11-08 13:21 ` M Joonas Pihlaja 0 siblings, 1 reply; 10+ messages in thread From: Luc Maisonobe @ 2003-11-08 12:51 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-08 12:51 Luc Maisonobe @ 2003-11-08 13:21 ` M Joonas Pihlaja 2003-11-08 14:20 ` Luc Maisonobe 2003-11-08 14:21 ` Jerome BENOIT 0 siblings, 2 replies; 10+ messages in thread From: M Joonas Pihlaja @ 2003-11-08 13:21 UTC (permalink / raw) To: Luc Maisonobe; +Cc: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-08 13:21 ` M Joonas Pihlaja @ 2003-11-08 14:20 ` Luc Maisonobe 2003-11-08 14:21 ` Jerome BENOIT 1 sibling, 0 replies; 10+ messages in thread From: Luc Maisonobe @ 2003-11-08 14:20 UTC (permalink / raw) To: M Joonas Pihlaja; +Cc: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
* Re: GSL 1.4: BUG #8 [rk8pd.c] 2003-11-08 13:21 ` M Joonas Pihlaja 2003-11-08 14:20 ` Luc Maisonobe @ 2003-11-08 14:21 ` Jerome BENOIT 1 sibling, 0 replies; 10+ messages in thread From: Jerome BENOIT @ 2003-11-08 14:21 UTC (permalink / raw) To: gsl-discuss 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 ^ permalink raw reply [flat|nested] 10+ messages in thread
end of thread, other threads:[~2003-11-08 14:21 UTC | newest] Thread overview: 10+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2003-11-02 12:57 GSL 1.4: BUG #8 [rk8pd.c] Jerome BENOIT 2003-11-03 19:42 ` Gerard Jungman 2003-11-04 11:27 ` MAISONOBE Luc 2003-11-05 15:26 ` Brian Gough 2003-11-05 15:32 ` Jerome BENOIT 2003-11-05 16:52 ` Brian Gough 2003-11-08 12:51 Luc Maisonobe 2003-11-08 13:21 ` M Joonas Pihlaja 2003-11-08 14:20 ` Luc Maisonobe 2003-11-08 14:21 ` Jerome BENOIT
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).