Runge-Kutta order test 2021

$$ \def\quad{\;\;\;\;} $$

Runge-Kutta order test 2021

In Our first B-series program Jim Verner and I described how we worked on an APL code to test the order of a method given by the three arrays $A$, $b^T$, $c$. Another way of verifying the order of a given method, is to modify Easy Octave algorithm 1 so that order conditions are evaluated recursively, as trees up to the required order are being generated. In addition to the order conditions, it is convenient that the inner consistency conditions are also checked. The Maple code given here, which incorporates a preliminary version of the procedure test, verifies the order of the method

\[ \begin{array}{c|ccc} 0 & \\ \frac13 & \frac13\\ \frac23 & 0 & \frac23\\ \hline & \frac14& 0 & \frac34 \end{array} \]

restart:
with(LinearAlgebra):

test := proc(A,b,c,p)
local i,s,n,k,r,l,fact,Psi,pt,first,last,one,err,R,ctest;
   s := Dimension(c);
   one := Vector(s,1);     
   ctest := A.one-c;
   for i from 1 to s do
      if ctest[i] <> 0 then
         print(ctest[i])
      fi;
   od;
   err := b.one-1;
   if err <> 0 then
      print(err) 
   fi;	
   first[1] := 1;
   last[1] := 1;
   fact[1] := 1;
   pt[1] := 1;
   Psi[1] := one;
   n := 1;
   for k from 2 to p do
      first[k] := n + 1:
      for r from 1 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;
               R[n] := r;
               fact[n] := fact[l]*fact[r]*k/(pt[l]);
               Psi[n] := DiagonalMatrix(Psi[l]).(A.Psi[r]);
               err := b.Psi[n]-1/fact[n];
               if err <> 0 then
                  print(err)
               fi
            fi
         od
      od;
      last[k] := n
   od;
end:

A := Matrix([[0,0,0],[1/3,0,0],[0,2/3,0]]):
c := Vector([0,1/3,2/3]):
b := Vector[row]([1/4,0,3/4]):

test(A, b, c, 3):
'test complete';

The absence of any output indicates that all conditions are satisfied exactly.

A second experiment using

A := Matrix(4,4,(i,j)->if i>j then a[i,j] else 0 fi):
bb := Vector[row](4,i->b[i]):
cc := Vector(4,i->if i>1 then c[i] else 0 fi):

test(A, bb, cc, 4):
gives the conditions that a 4 stage method has order 4.

\begin{align} &a[2,1]-c[2], \tag{1}\\ &a[3,1]+a[3,2]-c[3], \tag{2}\\ &a[4,1]+a[4,2]+a[4,3]-c[4], \tag{3}\\ &b[1]+b[2]+b[3]+b[4]-1,\\ &b[2]*a[2,1]+b[3]*\big(a[3,1]+a[3,2]\big)\\ &\quad+b[4]*\big(a[4,1]+a[4,2]+a[4,3]\big)-1/2,\\ &b[2]*a[2,1]^2+b[3]*\big(a[3,1]+a[3,2]\big)^2\\ &\quad+b[4]*\big(a[4,1]+a[4,2]+a[4,3]\big)^2-1/3\\ &b[3]*a[3,2]*a[2,1]+b[4]*a[4,2]*a[2,1]\\ &\quad +b[4]*a[4,3]*\big(a[3,1]+a[3,2]\big)-1/6\\ &b[2]*a[2,1]^3+b[3]*\big(a[3,1]+a[3,2]\big)^3\\ &\quad +b[4]*\big(a[4,1]+a[4,2]+a[4,3]\big)^3-1/4\\ &b[3]*a[3,2]*a[2,1]*\big(a[3,1]+a[3,2]\big)\\ &\quad+b[4]*a[4,2]*a[2,1]*(a[4,1)+a[4,2]+a[4,3])\\ &\quad+b[4]*a[4,3]*\big(a[3,1]+a[3,2]\big)\\ &\quad\quad*\big(a[4,1]+a[4,2]+a[4,3]\big)-1/8\\ &b[3]*a[3,2]*a[2,1]^2+b[4]*a[4,2]*a[2,1]^2\\ &\quad +b[4]*a[4,3]*\big(a[3,1]+a[3,2]\big)^2 -1/12\\ &b[4]*a[4,3]*a[3,2]*a[2,1]-1/24 \end{align}

This is a system of 11 equations in 13 unknowns. Although it can be solved using Maple, we will use the approach used by Kutta. First however, we solve (1), (2), (3) for a[2, 1], a[3,1], a[4,1] and substitute into the remaining equations, reordered for convenience

\begin{align} &b[1]+b[2]+b[3]+b[4]-1 \tag{4}\\ &b[2]*c[2]+b[3]*c[3]+b[4]*c[4]-1/2 \tag{5}\\ &b[2]*c[2]^2+b[3]*c[3]^2+b[4]*c[4]^2-1/3 \tag{6}\\ &b[2]*c[2]^3+b[3]*c[3]^3+b[4]*c[4]^3-1/4\tag{7}\\ &b[3]*a[3,2]*c[2]+b[4]*a[4,2]*c[2]\\ &\quad+b[4]*a[4,3]*c[3]-1/6\tag{8}\\ &b[3]*a[3,2]*c[2]*c[3] +b[4]*a[4,2]*c[2]*c[4]\\ &\quad+b[4]*a[4,3]*c[3]*c[4]-1/8\tag{9}\\ &b[3]*a[3,2]*c[2]^2+b[4]*a[4,2]*c[2]^2\\ &\quad+a[4,3]*c[3]^2-1/12\tag{10}\\ &b[4]*a[4,3]*a[3,2]*c[2] -1/24 \tag{11} \end{align}

The traditional Kutta approach is to carry out three steps
  1. Solve for $b[1]$, $b[2]$, $b[3]$, $b[4]$ from (4), (5), (6), (7), and substitute into the remaining equations,
  2. Solve for $a[3,2]$, $a[4,2]$, $a[4,3]$, from (8), (9), (10) and substitute into (11),
  3. Simplify (11) and deduce that $c[4]=1$.
This process, and similar manipulations for higher order and implicit methods can be made more convenient by modifying the procedure test in a number of ways. These are described in detail, with many examples of method derivation and analysis in a document to be written
$\to$ B series book
$\to$ Homepage