From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 25590 invoked by alias); 8 Nov 2003 14:20:47 -0000 Mailing-List: contact gsl-discuss-help@sources.redhat.com; run by ezmlm Precedence: bulk List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sources.redhat.com Received: (qmail 25582 invoked from network); 8 Nov 2003 14:20:46 -0000 Received: from unknown (HELO ioskeha.hittite.isp.9tel.net) (62.62.156.27) by sources.redhat.com with SMTP; 8 Nov 2003 14:20:46 -0000 Received: from osts (unknown [62.39.214.79]) by ioskeha.hittite.isp.9tel.net (Postfix) with ESMTP id A412517B4C2; Sat, 8 Nov 2003 15:20:44 +0100 (CET) Subject: Re: GSL 1.4: BUG #8 [rk8pd.c] From: Luc Maisonobe To: M Joonas Pihlaja Cc: gsl-discuss@sources.redhat.com In-Reply-To: References: <1068274810.1218.19.camel@lehrin.local.spaceroots.org> Content-Type: text/plain Message-Id: <1068301243.2285.52.camel@lehrin.local.spaceroots.org> Mime-Version: 1.0 Date: Sat, 08 Nov 2003 14:20:00 -0000 Content-Transfer-Encoding: 7bit X-SW-Source: 2003-q4/txt/msg00054.txt.bz2 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