$$
\def\quad{\;\;\;\;}
$$
$\to$ B series book
$\to$ Homepage
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- Solve for $b[1]$, $b[2]$, $b[3]$, $b[4]$ from (4), (5), (6), (7), and substitute into the remaining equations,
- Solve for $a[3,2]$, $a[4,2]$, $a[4,3]$, from (8), (9), (10) and substitute into (11),
- Simplify (11) and deduce that $c[4]=1$.
$\to$ B series book
$\to$ Homepage