Re: [Maxima-discuss] rkf45 share

[Michel T. <ta...@lp...>]

Le 18/01/2017 à 06:57, Robert Dodier a écrit :
> On 2017-01-17, Dimiter Prodanov <dim...@gm...> wrote:
>
>> I think there was a numerical issue in the package of Panagiotis
>> Papasotiriou.
>
> I'll be happy to fix the coefficients if you can be more specific about
> what is the problem.
>
>> a2:float([ 1/4,  1/4]), /* 2^2 */
>> a3:float([ 3/8, 3/32, 9/32]), /* 2^5 */
>> a4:float([ 12/13, 1932/2197, -7200/2197, 7296/2197]),  /* 13^3 */
>> a5:float([ 1, 439/216, -8, 3680/513, -845/4104]),
>> a6:float([ 1/2, -8/27, 2, -3544/2565, 1859/4104, -11/40]),
>> b1:float([ 25/216, 1408/2565, 2197/4101, -1/5]),
>> b2:float([ 16/135, 6656/12825, 28561/56430, -9/50, 2/55])
>
> b1 and b2 are supposed to sum to 1, right? Looks like b2 is OK but b1
> doesn't sum to 1. I think it should have 2197/4104 instead of 2197/4101.

 From
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
it is indeed 4104 which is correct.

>
> I don't know how to check the a's, maybe you can suggest something.

Otherwise the a and the b are like in Wikipedia.
>
> Also, I don't understand how these coefficients fit into the existing
> rkf45.mac. I think it's best for this purpose if you'll provide a patch
> that shows what lines need to be changed.
>

Well all those rationals have been computed as double floats in 
rkf45.mac which is a good thing so that everything is computed as double 
float by contagion. Otherwise if the function is given in terms of 
rationals maxima would compute evertything to rationals and it would be 
terrible.

Here is the Wikipedia table computed to double precision:

First lines a

0

0.25	0.25

0.375	0.09375	 0.28125

0.9230769230769231 0.8793809740555303 −3.277196176604461
3.320892125625853 <--- in fact this is 54 if computed with bfloat

1.0 2.032407407407407 −8 7.173489278752436 <-- here bfloat gives 37
−0.2058966861598441

0.5 −0.2962962962962963	2.0 −1.381676413255361 0.4529727095516569
−0.275

These are all the numbers i see in rkf45.mac between
/* Compute solution at xi+h */  and
/* Estimate local truncation error */

Then lines b

First b line
0.1185185185185185 0 0.5189863547758284 0.5061314903420167 −0.18	
0.03636363636363636
I don't know where these numbers appear, the last one is in truc_error.

second b line
0.1157407407407407 0 0.5489278752436647 0.5353313840155945 <-- 46 bfloat
-0.2	0
These are the numbers appearing in /* Step accepted, do it */
In particular this is computed with the correct 2197/4104.

At first sight coefficients are correct in rkf45.mac. It is difficult to 
compare to netlib, because the computations seem organized quite 
differently here. There is a more readable C++ version doing as in 
netlib here:
https://people.sc.fsu.edu/~jburkardt/cpp_src/rkf45/rkf45.cpp
after
//
//    Output, float S[NEQN], the estimate of the solution at T+H.
//
{
   float ch;
   int i;

....

One can see it is not directly comparable.

-- 
Michel Talon