Deriving new explicit Runge--Kutta pairs using Test21

Jim Verner

Outline - Initiated 2022/03/20
Initiating the parameters required for deriving Runge--Kutta Pairs.
Specifying sets of nodes to generate rational coefficient methods.
Determining that two coefficients can be selected arbitrarily.

Our first B-series program, RKTest21.maple was a code written by John Butcher October 10, 2021, to determine the order of a Runge--Kutta method, and adapted it to solve the order conditions directly in subsets. In 1978 and 1979, the author showed how to derive pairs of Runge--Kutta methods of successive orders, and families of methods of all orders up to the design order by using the same derivative evaluations. Here, RKTest21 is adapted to derive pairs of Runge--Kutta methods of these types.

There are several ways to denote a pair of methods that use the same derivative evaluations in each step: here we shall refer to an $(s,q-p)$ procedure as one for which two solution approximations, one of order $q < p$, and another of order $p$, such that either approximation can be propagated, and the difference of the results from the two methods can be used to estimate the local truncation.

In the 1960s, Erwin Fehlberg showed that (6,4-5), (8,5-6), (10,6-7), (13,7-8) and (17,8-9) RK pairs could be derived using the same derivative evaluations. These were designed to be used so that the difference of the two approximations at each step point could be used to estimate the local truncation error. This error estimate would be used in turn to change the stepsize in order to obtain a solution as efficiently as possible. Unfortunately, for $q>=5$, the two sets of weights in the formulas computed the same value for any component of a differential system that was a function of the independent variable alone: as a result, such an error estimate may have allowed an increased stepsize that was not justified. Larry Shampine showed that this deficiency could be countered by using special quadrature formulas to avoid inappropriate stepsize reduction.

Methods of orders (5,3-4) and (6,4-5) found and studied independently by R.E. Merson, R.F. Scranton, F. Ceschino and R. England in the 1960s had different sets of weights for the two error formulas. Yet, higher order pairs of this type were unknown. In the early 1970s, Tom Hull, Wayne Enright and their colleagues undertook a project to compare various known methods for solving initial value problems. Their report DETEST reviewed a considerable variety of methods, and as a result, those authors solicited formulas similar to those obtained by Fehlberg that provided different sets of weights for the two methods. I believe there were three responses to their request: John Butcher had found a (9,5(6)) pair for which only eight derivative evaluations were required if the method of order 5 were propagated. This author derived an (8,5-6) pair, and this was subsequently coded and distributed in the IMSL Subroutine library as DVERK. A third formula, for which I no longer remember the author, was also sent to the Toronto group.

Subsequently, this author spent a year in New Zealand deriving parametric formulas for families of these as well as (10,6-7), (13,7-8) and (16,8-9) pairs with some assistance from John Butcher. The parameters for these methods were in terms of the nodes, and when $p>7,$ one more coefficient, and for each pair, there were some constraints imposed on the nodes to achieve the desired orders. Here suitable choices of nodes for these methods are referred to as 'confined'. These families formed the basis of a general theory for their derivation that was published in 1978. In 1979, the author showed that for each Runge--Kutta pair, sets of weights could be selected to give families of methods of all orders up to order p. Subsequently, some related formulas were studied by John Dormand and Peter Prince. Recently [2010], the author specified definitions for both robust and efficient formulas for RK pairs, and provided both exact and approximate values of coefficients of some particular pairs on the website jverner@math.sfu.ca (click on "Jim Verner's Refuge for Runge--Kutta pairs"). A number of these particular procedures have been coded for treating initial value problems in several software packages.

The derivation of both single Runge--Kutta methods and pairs of successive orders requires the solution of systems of algebraic equations. For high order explicit methods, these systems are a challenge to solve, and systems have been solved by making assumptions about the coefficients proposed by John for this purpose. Such 'simplifying conditions' have been particularly helpful in obtaining solutions to many systems of explicit Runge--Kutta order equations. For example, known explicit minimal-stage single methods of orders 8 (G.J. Cooper and J.H. Verner and A.R. Curtis) and methods of order 10 by A.R. Curtis, and then E. Hairer were found by making suitable choices of simplifying conditions.

Recently, when John was developing an improved version of an algorithm for testing the order of Runge--Kutta methods, he realized that the order conditions might be used more directly for finding their solution. His new code, RKTest21.maple has been implemented to do this. To determine how this code might operate, I decided to use it to determine whether it might be used to derive some explicit Runge--Kutta pairs. At the outset, I thought it might not be completely successful, but expected that some known pairs could be obtained, and as well this might provide further insight in how to improve RKTest21.maple

The experiment was successful, and for given sets of confined nodes the following code when executed for selected values of $s$ and $p$, coefficients of pairs of methods will be computed and displayed. As well, the maximum local truncation error coefficient will also be displayed (although this has to be divided by suitable functions of each tree since the order equations are all normalized.)

Some modifications to RKTest21.maple are needed to derive pairs. In that program, subsets of equations are solved initially with MAPLE automatically selecting sufficient parameters from among those available to solve each subset when this is possible. For pairs of moderate orders, when nodes are rational requiring small numbers of digits for exact representation, solutions are obtained. For pairs of low to moderate orders, and exact rational values of confined nodes, pairs are obtained by MAPLE. For methods of orders 7 and 8 with 13 stages, in some cases, the object 'soln' consists of a several solution 'sets', each of which is a valid solution. Initially to accommodate this, each element was extracted, and this eventually led to a complete solution. What is remarkable, is that MAPLE was able to find parametric solutions for each particular set of rational confined nodes requiring few digits that was selected for testing. Later it was found that only more restricted sets of nodes yielded solutions.

In the article "Pairs with error estimates" [1978] by the author, it was shown that to achieve the orders for the numbers of stages specified, and to obtain different sets of weights for the two approximations, the nodes must be restricted to satisfy some 'nodal conditions' that determined sets of suitably confined nodes. For each of many sets of confined nodes for $p<8$, only a single solution pair was found using RKTest21.maple . However for (13,7-8) derivations, and several confined sets of nodes, some partial solutions provided several individual solution sets, and some of these led to different solutions. In particular, by isolating each solution set, and allowing MAPLE to proceed, we found that some solutions found could be represented with two essentially arbitrary coefficients a[7,6] and a[8,7] for each set of nodes. If these two parameters are suppressed within Test21.maple, then for each specified set of confined nodes, MAPLE algebraically generates a two-parameter family of solutions in terms of these two parameters.

Hence, this more general approach to solving the Runge--Kutta order conditions has already shown its value in leading to the discovery of new solutions. For such (13,7-8) families, Test21.maple seems not to be able to compute values of the coefficients for sets of confined nodes that require large numbers of digits for exact rational representation. Even so, it might be possible to obtain finite digit approximations for members of such larger families using fsolve in MAPLE to solve the order equations approximately; in this respect, N-digit approximations (N chosen) to the coefficients for methods with smaller local truncation errors might lead to more efficient production solvers for initial value problems. This discovery also suggests another problem of possible interest: it might be possible to develop a modified approach that would allow the direct derivation of these parametric formulas for all of the coefficients in terms of the confined nodes and the two entries a[7,6] and a[8,7] (or possibly two different arbitrarily chosen entries of the matrix A).

Runge-Kutta pairs of formulas of successive orders exist. Such pairs may be utilized to develop software for obtaining approximate solutions of ordinary differential equations with error control. Historically, to derive a pair of specific orders, the number of stages is minimized, and for this, usually constraints on the nodes have been imposed as shown below. Additionally, we attempt to have all nodes lie within the interval [0,1] to avoid derivative evaluations near singularities that might lie ahead of the computation points. Again, this imposes constraints on nodes. For example, to obtain (8,5(6)) pairs, select any value of c[3] in (0,1) other than 1/3 or 2/5 , c[4] is a specified value in [0,1], and other nodes can be selected arbitrarily in [0,1].

For p < 6, s=p+1 is sufficient to obtain either methods of order p, or (s,(p-1)p) RK pairs. For p > 5, s>p, and we can obtain (s,p) methods, or (s+1,(p-1)p) RK pairs. Here, we select s to initiate computation of coefficients of (s+1,(p-1)p) RK pairs.

restart;
with(LinearAlgebra);
Put  #  before any statement that is to be suppressed.

solvindex := 0;
maxsolno := 1;
s := 6; p := 5;
s := 7; p := 6;
s := 9; p := 7;
s := 12; p := 8;

test := proc(A, b, c, p, q, phi) 
local n, k, r, l, fact, Psi, pt, qt, first, one, R, compt, 
i, ctest, ttest, err; global s, last; 
s := Dimension(c); 
one := Vector(s, 1); 
ctest := phi~((A . one) - c); 
compt := 0; 
for i to s do 
   if ctest[i] <> 0 then 
      compt := compt + 1; 
      err[compt] := ctest[i]; 
   fi; 
od: 
ttest := phi~((b . one) - 1); 
if ttest <> 0 then 
   compt := compt + 1; 
   err[compt] := ttest; 
fi:
first[1] := 1; last[1] := 1; fact[1] := 1; 
pt[1] := 1; qt[1] := 0; 
Psi[1] := one; 
n := 1; 
for k from 2 to p do 
   first[k] := n + 1; 
   for r to last[k - 1] do 
      for l from first[k - pt[r]] to last[k - pt[r]] do 
         if l = 1 or r <= R[l] then 
            n := n + 1; 
            pt[n] := k; 
            if l = 1 then 
               qt[n] := 1 + qt[r]; 
            else 
               qt[n] := qt[l] + qt[r]; 
            fi: 
            if q < qt[n] then 
               n := n - 1; 
            else 
               R[n] := r; 
               fact[n] := fact[l]*fact[r]*k/pt[l]; 
               Psi[n] := phi~((DiagonalMatrix(Psi[l])).(A.(Psi[r]))); 
               ttest := phi((b.(Psi[n])) - 1/fact[n]); 
               if ttest <> 0 then 
                  compt := compt + 1; 
                  err[compt] := ttest; 
            	fi: 
         	fi; 
      	fi;
   	od; 
	od: 
	last[k] := n; 
od; 
Vector(compt, i -> err[i]); 
end proc;

phi := x -> factor(simplify(expand(x)));
nentries := proc() _npassed; end proc;
N := 100;

tester := proc(p, q := p) 
global A, b, c, e, f; 
e := test(A, b, c, p, q, phi); 
f := e; 
Matrix(min(Dimension(f), N), 2, (i, j) -> if j = 1 then i; else f[i]; fi:); 
end proc;

show := proc(n := Dimension(f)) 
local i; 
global f; 
for i to min(Dimension(f), n) do 
   printf("%g: %a\n", i, f[i]); 
od; 
end proc;

filter := proc(x) 
local temp, i, n; 
global e, f; 
n := Dimension(e); 
if type(x, numeric) then 
	f := Vector([e[x]]); 
fi; 
if type(x, algebraic) then 
	print(`type(x) is algebraic`); 
	temp := convert(e, set); 
	for i to n do 
		if coeff(e[i], x, 1) = 0 then 
			temp := temp minus {e[i]}; 
		fi; 
	od; 
	f := Vector(convert(temp, list)); 
fi; 
if type(x, list) then 
	f := Vector(Dimension(Vector(x)), i -> e[x[i]]); 
fi; 
if type(x, range) then 
	temp := Vector(min(convert(x, list)[2], n), i -> e[i]); 
	f := Vector(convert(temp, list)); 
fi; 
if type(x, set) then 
	f := Vector(convert(x, list)); 
fi; 
Matrix(min(Dimension(f), N), 2, (i, j) -> if j = 1 then i; else f[i]; fi); 
end proc;

clean := proc() 
local temp, i, n; 
global soln; 
temp := soln; 
n := numelems(soln); 
for i to n do 
	if soln[i] then 
		temp := temp minus {soln[i]}; 
	fi; 
od; 
print(`Leaving clean, type of soln is `*whattype(temp)); 
print(`The set of cleaned solution components   `); 
print(`has a number of new values equal to `); 
soln := temp; 
end proc;

solver := proc(variables) 
local i, n, solno; 
global soln, f, maxsolno, solvindex; 
solvindex := solvindex + 1; 
lprint(` `); 
print(`==========================================`); 
print(`        Execute solver iteration           `, solvindex); 
print(`The maximum number of solutions is  `, maxsolno); 
soln := solve(convert(f, set), variables); 
print(`In solver, the type of  soln is  `, whattype(soln)); 
n := nentries(soln); 
lprint(`Number of solutions for this iteration is  `, n); 
if n = 0 then 
	lprint(`-------`); 
	print(`soln are   `, soln); 
	print(`It appears for maxsolno = `, maxsolno, `  there are NO SOLUTIONS.`); 
	return ; 
fi; 
solno := min(maxsolno, n); 
if 1 < n then 
	lprint(`Here we compute soln #`, solno, ` out of `, n, `  solutions.`); 
fi; 
if whattype(soln) = exprseq then 
   soln := soln[solno]; 
fi; 
if soln <> {} then 
   clean(); 
fi; 
n := numelems(soln); 
end proc;

fix := proc() 
global A, b, c, soln; 
b := phi~(subs(soln, b)); 
c := phi~(subs(soln, c)); 
A := phi~(subs(soln, A)); 
end proc;

A := Matrix(s, s, (i, j) -> if j < i then a[i, j]; else 0; fi);
AA := Matrix(s, s, (i, j) -> a[i, j]);
Aij := r -> convert(SubMatrix(AA, r[1] .. r[2], r[3] .. r[4]), set);
all := Aij([1, s, 1, s]);
row := i -> Aij([i, i, 1, s]);
col := j -> Aij([1, s, j, j]);

For s > 8, some variables have been assumed to have values equal to zero in order to simplify the solution by the MAPLE solve operator. It is possible that more solutions exist than those possible here. Observe, for s=12, we have found new parametric solutions, when this restriction was not totally imposed.

if s >= 9 then 
	all := all minus {A[7, 2], A[8, 2], A[9, 2], a[4, 2], a[5, 2], a[6, 2]}; 
	for i from 4 to 9 do 
		a[i, 2] := 0; 
	od; 
fi;
if s >= 12  then
    for i from 6 to 12 do 
    	 all := all minus {a[i, 2], a[i, 3]}; 
    	 a[i, 2] := 0; 
    	 a[i, 3] := 0; 
    od;
    all := all minus {a[7, 6]};
    a[7, 6] := a76;
    all := all minus {a[8, 7]};
    a[8, 7] := a87;
fi;

if s = 5 then
  	c := Vector([0, 1/3, 1/4, 2/3, 3/4]);
   b := Vector[row]([b1, b2, b3, b4, b5]);
elif s = 6 then
   c := Vector([0, 3/10, 2/5, 1, 39/40, 1/40]);
   c4 := c[3]/(10*c[3]^2 - 8*c[3] + 2);
   if c4 <> c[4] then
      lprint(`This choice of nodes will fail to yield a method `);
      lprint(`of order 6 for the proposed   (8,5-6)  pair.`);
      lprint(`We have changed   c[4] = `, c[4], ` to `, c4);
      lprint(` in order to get such a method.`);
      c[4] := c4;
   fi;
   b := Vector[row]([b1, 0, b3, b4, b5, b6]);
elif s = 7 then
   c := Vector([0, 1/6, 4/15, 2/3, 3/4, 1/15, 1]);
   c4 := c[3]/(15*c[3]^2 - 10*c[3] + 2);
   if c4 <> c[4] then
      lprint(`This choice of nodes will fail to yield a method `);
      lprint(`of order 6 for the proposed   (8,5-6)  pair.`);
      lprint(`We have changed   c[4] = `, c[4], ` to `, c4);
      lprint(` in order to get such a method.`);
      c[4] := c4;
   fi;
   b := Vector[row]([b1, 0, b3, b4, b5, b6, b7]);
elif s = 9 then
   c := Vector([0, 1/18, 1/6, 1/4, 4/5, 8/9, 3/4, 1/8, 1]);
   c3 := (2*c[4])/3;
   sm := c[4] + c[5];
   pr := c[4]*c[5];
   c6 := (sm - pr*(12 - 7*sm))/(3 - sm*(12 - 14*sm + 70*pr) + pr*(24 + 105*pr));
   if c3 <> c[3] or c[6] <> c6 then
      lprint(`This choice of nodes will fail to yield a method`);
      lprint(`of order 7 for the proposed   (10,6-7)  pair.`);
      lprint(`We have changed   c[3] = `, c[3], ` to `, c3);
      lprint(` and c[6] = `, c[6], ` to `, c6);
      c[3] := c3;
      c[6] := c6;
   fi;
   b := Vector[row]([b1, 0, b3, b4, b5, b6, b7, b8, b9]);
   for i from 4 to 9 do
      a[i, 2] := 0;
   od;
elif s = 12 then
   maxsolno := 3;           # Choose this to be the upper bound on the number of solutions
   c := Vector(12, [0, 2/27, 1/9, 1/6, 5/12, 1/2, 5/6, 1/6, 2/3, 1/4, 1/3, 1]);
   a76 := 125/54;
   a87 := 13/900;
   c[4] := c[6]*(3*c[6] - 4*c[5])/(4*c[6] - 6*c[5]);
   c[3] := 2/3*c[4];
   s1 := c[6] + c[7] + c[8];
   s2 := c[6]*(c[7] + c[8]) + c[7]*c[8];
   s3 := c[6]*c[7]*c[8];
   tu := 20 - 28*s1 + 42*s2 - 70*s3;
   tv := 28 - 42*s1 + 70*s2 - 140*s3;
   tw := 15 - 24*s1 + 42*s2 - 84*s3;
   tx := 24 - 42*s1 + 84*s2 - 210*s3;
   ty := 4 - 7*s1 + 14*s2 - 35*s3;
   tz := 3 - 6*s1 + 14*s2 - 42*s3;
   c[9] := (tu*tz - tw*ty)/(tv*tz - tx*ty);
   print(`Nodes  c  = `, Vector[row](c));
   b := Vector[row]([b1, 0, 0, b4, b5, b6, b7, b8, b9, b10, b11, b12]);
   for i from 6 to 12 do
      a[i, 2] := 0;
      a[i, 3] := 0;
   od;
   a[4, 2] := 0;
   a[5, 2] := 0;
fi;

tester(p, 0);
filter(1 .. s - 1);
solver(col(1));
fix();
if s = 5 then
    vars := {b1, b3, b4, b5};
    b2 := 0;
elif s = 6 then
    vars := {b1, b3, b4, b5, b6};
    b2 := 0; bb2 := 0; b[2] := 0; bb[2] := 0;
    bb6 := 0; bb[6] := 0;
elif s = 7 then
    vars := {b1, b3, b4, b5, b6, b7};
    b2 := 0; bb2 := 0; b[2] := 0; bb[2] := 0;
    bb6 := 0; bb[6] := 0;
    bb7 := 0; bb[7] := 0;
elif s = 9 then
    vars := {b1, b4, b5, b6, b7, b8, b9};
    b2 := 0; bb2 := 0; b[2] := 0; bb[2] := 0;
    b3 := 0; bb3 := 0; b[3] := 0; bb[3] := 0;
    bb8 := 0; bb[8] := 0;
    bb9 := 0; bb[9] := 0;
elif s = 12 then
    vars := {b1, b10, b11, b12, b6, b7, b8, b9};
    b2 := 0; bb2 := 0; b[2] := 0; bb[2] := 0;   
    b3 := 0; bb3 := 0; b[3] := 0; bb[3] := 0;   
    b4 := 0; bb3 := 0; b[4] := 0; bb[4] := 0;
    b5 := 0; bb5 := 0; b[5] := 0; bb[5] := 0;
    bb11 := 0; bb[11] := 0;
    bb12 := 0; bb[12] := 0;
fi;

tester(p, 1);
solver(vars);
fix();
print(`Weights = `, Vector(b));

tester(p, 2);
solver(all);
fix();

tester(p, 3);
solver(all);
fix();

if 3 < p then tester(p, 4); end if;
solver(all);
fix();

if 4 < p then tester(p, 5); end if;
solver(all);
fix();

if 5 < p then tester(p, 6); end if;
solver(all);
fix();

if 6 < p then tester(7); end if;
solver(all);
fix();

if 7 < p then tester(8); end if;
solver(all);
fix();

######################################################
#
#                  The ERROR ESTIMATING METHOD
#
######################################################

# First, find the quadrature method on   ce:=c   minus    c[s-1]

# Then solve        SUM    beAe[j] - be(1-ce[j])  = 0   for   j=1..s-2.

if s = 6 then
   ct := c;
   ce := c;
   bt := b;
   be := Vector(s, i -> if i = 2 then 0; elif i < s then bb[i]; else 0; fi);
else
   ct := Vector(s + 1, i -> if i <= s then c[i]; else 1; fi);
   ce := Vector(s + 1, i -> if i < s - 1 then c[i]; 
                            elif i = s + 1 then 1; else 0; fi);
   bt := Vector(s + 1, i -> if i <= s then b[i]; else 0; fi);
   be := Vector(s + 1, i -> if i < s - 1 then bb[i];
                            elif i = s + 1 then bb[s + 1]; else 0; fi);
fi;

vars := {bb[1]};
eqns := {};
if s = 6 then
   terr[1] := add(bb[i], i = 1 .. 5) - 1;
   eqns := {terr[1], op(eqns)};
   for k from 2 to 4 do
      terr[k] := add(bb[i]*ce[i]^(k - 1), i = 2 .. 5) - 1/k;
      eqns := {terr[k], op(eqns)};
      vars := {bb[k + 1], op(vars)};
   od;
else
   for k to p - 1 do 
      terr[k] := add(bb[i]*ce[i]^(k - 1), i = s - p + 2 .. s - 2) 
                + bb[s + 1]*ce[s + 1]^(k - 1) - 1/k; 
      if k = 1 then 
         terr[1] := terr[1] + bb[1]; 
         vars := {bb[s + 1], op(vars)}; 
      fi; 
      eqns := {terr[k], op(eqns)}; 
      if k < p - 2 then 
         vars := {bb[s - p + k + 1], op(vars)}; 
      fi; 
   od;
fi;

bsol := solve(eqns, vars);
assign(bsol);

if s = 6 then
   Ae := A;
else
   Ae := Matrix(s + 1, s + 1, (i, j) -> if i <= s and j < i then A[i, j]; 
                                        else 0; fi);
   for j from s - 1 by -1 to 1 do
      Ae[s + 1, j] := (bb[j]*(1 - ce[j]) - add(bb[i]*A[i, j], i = j + 1 .. s - 2))/bb[s + 1];
   od;
fi;
 
#   NOW DISPLAY COEFFICIENTS of an ERROR ESTIMATING 
#             RUNGE-KUTTA PROCEDURE

print(Vector(ct), Matrix(Ae));
if s = 12 then
   print(c[11], `         `, Row(Ae, 11));
   print(c[12], `                     `, Row(Ae, 12));
   print(ce[13], `                `, Row(Ae, 13));
fi;
print(`            `, Vector[row](bt));
print(`            `, Vector[row](be));

# We now verify the order of each method by showing that there are no 
# errors up to order  p for the high order method:
test(Ae, bt, ct, p, p, phi);

# and no errors up to order p-1 for the low order (embedded) method:
test(Ae, be, ct, p - 1, p - 1, phi);

# Further evidence is given by showing error values in the order conditions 
# of order p  for the method  of order   p-1   that are NOT satisfied:

if 11 < s then k := 5; else k := 1; fi;

for j to k do
    maxerr := 0;
    if 11 < s then
        a76 := 13./(900);
        a87 := (3.);
    fi;
    terrors := test(Ae, be, ct, p, p, phi);
    n := numelems(terrors);
    for i to n do
      maxerr := max(maxerr, abs(terrors(i)));
    od;
    lprint(`maxerr(`, j, `) is  `, evalf[5](maxerr));
od;
lprint(`Number of error terms is   `, n);
for j to k do 
	maxerr := 0; 
	if 11 < s then a76 := 1./(1000*j); a87 := (-1.)/(2000*j); fi; 
	terrors := test(Ae, bt, ct, p + 1, p + 1, phi); 
	n := numelems(terrors); 
	for i to n do 
		maxerr := max(maxerr, abs(terrors(i))); 
	od; 
	lprint(`maxerr(`, j, `) is  `, evalf[5](maxerr)); 
od;
lprint(`Number of error terms is   `, n);

print(`For this set of error terms, the maximum number of non-zero \rerror conditions should be `, last[p], ` - `, last[p - 1], ` = `, last[p] - last[p - 1]);
print(`BUT, for some methods, some errors of order `, p, ` are zero; for the current pair, the number of non-zero terms is shown above.`);