\documentstyle[psfig,leqno,seceqn]{article}
%\documentstyle{article}
\oddsidemargin=-0.10cm
\textwidth=15cm
%\textwidth=+1cm
%\textheight=+1cm
\newtheorem{teo}{Theorem}         

\title{Parallel Iterated Methods based on \\
 Multistep Runge-Kutta Methods of Radau type\thanks{draft3.0}}
\date{ }
\author{K. Burrage\thanks{Department of Mathematics, University
of Queensland, Australia} \and H. Suhartanto\thanks{Department of Mathematics,
University of Queensland, Australia}}
       
\begin{document}
\pssilent
\psfigurepath{/home/ciamp/hs/rs}
\maketitle
\begin{abstract}
This paper investigates iterated Multistep Runge-Kutta methods of Radau type
as a class of explicit methods suitable for parallel implementation.  
Using the idea of van der Houwen and Sommeijer \cite{hs90}, the method is designed
in such a way that the right-hand side evaluations can be computed in
parallel. We use  stepsize control and variable order based on
iterated approximation of the solution. A code is developed and its performance is
compared with codes based on iterated Runge-Kutta methods of Gauss type and various Dormand
and Prince pairs \cite{hvol1}.
The accuracy of some of our methods are  comparable with the PIRK10 methods
of  
van der Houwen and Sommeijer
\cite{hs90}, but require less processors. In addition at very stringent
tolerances these new methods are competitive with RK78 pairs in a sequential
implementation.
\end{abstract}

\section{Introduction}
The invention of parallel computers influences  the development of methods
for solving initial value problems for systems of ordinary differential
equations
\begin{equation}
\label{e1}
y'(x) = f(y(x)), \; y(x_{0})=y_{0}, \; f: R^{m} \rightarrow R^{m}.
\end{equation}
Implicit methods  such as Implicit Runge-Kutta (IRK) and Multistep Runge-Kutta
(MRK) methods
which were not previously of much interest due to the increased number
of coupled nonlinear algebraic equations, unless very special structures
were imposed on the coefficient matrix, has now begun to attract more people.
However, it was Jackson and N\o rsett \cite{r10jn88} and van der Houwen and Sommeijer \cite{hs88}  who introduced the 
iteration of IRK to take advantage of parallelism. Since
then various issues have been debated in \cite{r9in90,r11lie87} and
\cite{r12ns89}. Despite these debates, the results of van der Houwen and Sommeijer \cite{hs90}   
showed that codes based on parallel iteration are comparable in accuracy and
more efficient in terms of function evaluation than standard sequential codes
such as DOPRI8 which implements high order formulas of Prince and Dormand \cite{hvol1}. 
Van der Houwen and Sommeijer's work
also showed that it is possible to construct high order methods with a minimum
number of stage evaluations. Their contribution motivates us to construct
 higher order methods  than their methods with a  similar number of internal stages.
This is possible since our iterated methods will be based on
MRK methods which can be 
order of $2s+r-1$, see for example Burrage \cite{bu4}. 
The method is characterized by

\begin{eqnarray}
Y_{i} & = & \sum_{j=1}^{r} a_{ij}y_{n+1-j} + h \sum_{j=1}^{s} b_{ij} f(Y_{j}), \quad 
i=1,\ldots, s, \nonumber \\
y_{n+1} &  =  & \sum_{j=1}^{r} \alpha_{j} y_{n+1-j} + h\sum_{j=1}^{s} \gamma_{j} 
f(Y_{j}). \label{e2}
\end{eqnarray}

The method is an example of more general class of methods called a multivalue (general linear) methods. 
The $y_{n+1-i}, \, (i=1,\ldots,r)$ contain all the information from step to step
and are updated at the end
of a step by computing $f(Y_{1}), \ldots, f(Y_{s})$, which represent
approximations to the derivatives of the solution at $s$ internal points, and
taking appropriate linear combination. In section 2 we will review the concepts
of MRK methods and,  based on the idea of van der Houwen and
Sommeijer \cite{hs90}, we will develop iterated methods based on
MRK in section 3. The last two sections will be devoted
to implementation issues and numerical results.

\section{Multistep Runge-Kutta methods}
The general class of MRK methods have been studied
by Butcher \cite{but6,but7}, Cooper \cite{cop8}, Hairer and Wanner \cite{hw9}, Burrage
\cite{bu4,bu5} and Burrage and Moss \cite{bumos3}. In particular, Burrage 
\cite{bu4} has studied
the order conditions of these methods, and has shown that one can always
construct methods of order $2s+r-1$ based on collocation approximation.
By defining the following relationships: 

\begin{eqnarray*}
C(p): & & A \lambda^{k} = c^{k} - k B c^{k-1}, \quad k=0,\ldots,p,\\
B(p): & & \alpha^{T} \lambda^{k} = r^{k} -  k \gamma^{T} c^{k-1}, \quad
 k=0,\ldots, p, \\ 
D(p): & & k \gamma^{T} C^{k-1} B = r^{k}\gamma^{T} - \gamma^{T}C^{k},
\quad k=1,\ldots,p,
\end{eqnarray*}
where $A$ and $B$ are $s \times r$ and $s \times s$ matrices whose $(i,j)$
elements are $a_{ij}$ and $b_{ij}$ respectively, $\lambda = (0,-1,\ldots,
-r+1)^{T}, c=(c_{1},\ldots, c_{s})^{T}, \alpha = (\alpha_{1}, \ldots, \alpha_{r})^{T},
\gamma=(\gamma_{1},\ldots,\gamma_{s})^{T}, C=\mbox{diag}(c_{1},\ldots, c_{s})$
and multiplication of vectors is done componentwise , Burrage \cite{bu4} 
has proved

\begin{teo}.
A multistep Runge-Kutta method will be of order $w$ if
$B(w), C(\eta), D(\mu)$ hold where $\mu$ and $\eta$ are nonnegative integers
such that $\eta \geq r-1, \quad w \leq \mu + \eta + 1, \quad w \leq 2 \eta+2$.
\end{teo}

\begin{teo}.
A multistep Runge-Kuta method will be of order $s+r-1+t$ if
$C(s+r-1)$, and $B(s+r-1+t)$ hold, for $t=1,\ldots,s$.
\label{teo2}
\end{teo}

Hairer and Wanner \cite{hnw80} view the methods, characterized by
Theorem 2,  as multistep collocation
methods, assuming that one has $s$ real numbers $c_{1}, \ldots, c_{s}$
and $r$ solution values $y_{n}, \ldots, y_{n+1-r}$. Thus reverting to the 
non-autonomous version of (\ref{e1}) and defining  the 
corresponding  polynomial $u(x)$ of degree $s+r-1$ by

\begin{eqnarray}
\label{e323a}
u(x_{j}) & = & y_{j}, \quad j=n, \ldots, n+1-r  \\
\label{e323b}
u'(x_{n}+c_{i}h) & = & f(x_{n}+c_{i}h, u(x_{n}+c_{i}h)), \quad i=1, \ldots, s,
\end{eqnarray}
the numerical solution is then given by 
\begin{equation}
\label{e323c}
y_{n+1} = u(x_{n+1}).
\end{equation}

Now suppose that the derivatives $u'(x_{n}+c_{i}h)$ are known, then 
(\ref{e323a}) and (\ref{e323b}) form a Hermite interpolation problem
with insufficient data since the function values at $x_{n}+c_{i}h$ are
not available. In order to overcome this problem the dimensionless coordinate
$t=(x-x_{n})/h, \, x=x_{n}+th$, nodes $t_{1}=-r+1, \ldots,
t_{r-1}=-1, t_{r}=0$ are defined. In addition, the  polynomial 
$\varphi_{i}(t),  (i=1,\ldots, r)$ of degree $s+r-1$ is defined by

\begin{eqnarray}
\varphi_{i}(t_{j}) & = & \left\{ \begin{array}{lll}
                            0 & \mbox{if} & i \neq j \\
                            1 & \mbox{if} & i=j \end{array} 
                         \right. \; j=1,\ldots,r \label{e324}  \\
\varphi_{i}'(c_{j}) & = & 0, \quad j=1, \ldots, s \nonumber
\end{eqnarray}
and polynomial $\psi_{i}(t) \, (i=1,\ldots,s)$ is defined by
\begin{eqnarray}
\psi_{i}(t_{j}) & = &  0,  \quad  j=1,\ldots,r \nonumber  \\
\psi_{i}'(c_{j}) & = & \left\{ \begin{array}{lll}
1 & \mbox{if} & i=j \\
0 & \mbox{if} & i\neq j \end{array} \right. \; j=1,\ldots,s. \label{e325}
\end{eqnarray}
Now polynomial $u(x)$ can be written as
\begin{equation}
\label{e326}
u(x_{n}+th) = \sum_{j=1}^{r} \varphi_{j}(t)y_{n+1-j} + h \sum_{j=1}^{s}
\psi_{j}(t)u'(x_{n}+c_{j}h).
\end{equation}
In (\ref{e326}) if we put $t=c_{i}$, and writing $u(x_{n}+c_{i}h) = v_{i}$ then
inserting the collocation condition (\ref{e323b}) gives 

\begin{eqnarray}
v_{i} & =  & \sum_{j=1}^{r} \varphi_{j}(c_{i}) y_{n+1-j} +
h \sum_{j=1}^{s} \psi_{j}(c_{i}) f(x_{n}+c_{j}h,v_{j}), \quad i=1, \ldots, s
\label{e327a} \\
y_{n+1} & = & \sum_{j=1}^{r} \varphi_{j}(1) y_{n+1-j} + h 
\sum_{j=1}^{s} \psi_{j}(1) f(x_{n}+c_{j}h, v_{j}), \label{e327b} 
\end{eqnarray}
which is the MRK  formula 
(\ref{e2}) where   $a_{ij} = \varphi_{j}(c_{i}), b_{ij} =
\psi_{j}(c_{i}), \alpha_{j}=\varphi_{j}(1),$ and $\gamma_{j} =
\psi_{j}(1)$. 

The highest possible order of the method depends on the way one chooses 
the nodes $c_{i}$.  Burrage\cite{bu4}, Guillou \& Soule \cite{gs69} and Lie \&
N\o rsett \cite{ln89} constructed a method of maximal order $p=2s+r-1$. 
But these methods are not stiffly stable and are not suitable for
stiff problems. Since this paper represents the first of a series
of papers in which we hope to develop a type sensitive code based on
MRKs, it will be necessary to have a stiffly stable corrector for the stiff
part, and so we will also use
this corrector for the nonstiff part. Thus
we consider
methods of Radau type by choosing 
$c_{s}=1$ and try to determine the remaining nodes $c_{i}$ so that 
we obtain a method of order $p=2s+r-2$. 
The methods are characterized by the tableau form
\begin{equation}
\begin{tabular}{c}
c \\ \hline 
A \\ \hline
B
\end{tabular}.
\label{tab1}
\end{equation}
Hairer and Wanner \cite{hnw80}
using the idea of Krylov \cite{k59} and adapting to
quadrature problems have proved that one can obtain
a method of order $p=2s+r-2$ where the nodes are computed by
solving nonlinear equations $\chi'_{i}(c_{i})=0, \quad i=1,\ldots,s-1$, where
$\chi_{i}(t)$ is defined as
\begin{equation}
\label{e335}
\chi_{i}(t) = C\,\prod_{j=1}^{r}(t-t_{j})\,\prod_{j=1,j\neq i}^{s}
(t-c_{j})^{2},
\end{equation}
and where $C$ is determined by  $\chi_{i}(c_{i})=1$. The solution to the ensuing 
nonlinear equation has  $\left( \begin{array}{c} s+r-1 \\ r-1 \end{array}
\right)$ solutions, however we consider only 
the solutions $c_{i}$ with $0<c_{i}<1$. 
Appendix A defines the method parameters for the $(r,s)$ pairs $(2,2)$,
$(2,3)$, and $(3,3)$ as well as for what appears to be the three
most effective methods namely the $(4,2)$, $(5,3)$, and $(7,3)$
methods.

% move these tables into appendix A


Note that when $r=1$ we obtain the class of Radau Runge-Kutta 
methods of order $2s-1$ that are stiffly accurate, while when
$s=1$ we obtain the family of BDF methods. We also note that
the class of stiffly accurate methods of order $2s+r-2$
have been implemented by Schneider \cite{sch93} who shows that
this class of methods can be very effective at stringent
tolerances for stiff problems.

\section{Iterated Multistep Runge-Kutta Methods}
We start with the iterated method based on general MRK methods,
not necessarily of Radau type,  by using a Kronecker tensor
product form and following the notation of van der Houwen and Sommeijer \cite{hs90}.
Thus we write (\ref{e2}) as 
\begin{equation}
\label{mrk1}
y_{n+1} =  (\alpha^{T} \otimes I) y^{(n)} + h ( \gamma^{T} \otimes I)  F(Y_{n}),
\end{equation}
where 
\[
Y_{n}  =  (A \otimes I )y^{(n) } + h (B \otimes I) F(Y_{n}). 
\]
Here  $Y_{n}$ represents $s$ stages vectors
$(Y_{n,1}^{T}, \ldots, Y_{n,s}^{T})^{T}$, $F(Y_{n}) = (f(Y_{n,1})^{T}, \ldots, f(Y_{n,s})^{T})^{T}$,  $y^{(n)}$ represents $r$ previous solutions
$(y_{n}^{T}, \ldots, y_{n+1-r}^{T})^{T}$ and $I$ denotes the identity
matrix of order $m$. 
Since we are only concerned with solving non-stiff problems, the
implicitly defined equations above can be solved by some form of
fixed-point iteration without placing too severe a restriction on
the stepsize in terms of a  Lipschitz constant for $f$. 
Thus after $L$ corrections we obtain the method

\begin{eqnarray}    
\begin{array}{lll}
\begin{array}{lll}
Y_{n}^{(j)}  & = &  (A \otimes I) y^{(n)} + h ( B \otimes I) F(Y_{n}^{(j-1)}), 
\end{array}  & &  \quad j=1, \ldots, L  \label{it0} \\
y_{n+1}   =   (\alpha^{T} \otimes I)  y^{(n)}
+ h (\gamma^{T} \otimes I)  F(Y_{n}^{(L)}). & &  
\end{array} 
\end{eqnarray}
The computational time needed for one iteration of 
(\ref{it0}) is equivalent
to the time required to evaluate one right-hand side function on 
a sequential computer assuming that there are $s$ processors available and
the component of $F(Y_{n}^{(j)})$ is done in parallel.

We call a method providing $F(Y_{n}^{(0)})$  the predictor method and
(\ref{mrk1}) the corrector method, and the resulting parallel, iterated MRK
method will be called PIMRK.
Let $F(Y_{n}^{(0)})$ be an approximation to $F(Y_{n})$ satisfying the condition
\begin{equation}
\label{p1}
f(Y^{(0)}_{n,i}) = f(Y_{n,i}) + O(h^{q}), \quad i=1,\ldots,s,
\end{equation}
resulting in $y_{n+1}^{(0)} = y_{n+1} + O(h^{q+1})$, then the  predictor method satisfying
(\ref{p1}) is called a predictor method of order $q$.

Since we are using corrector methods of Radau type of order $p=2s+r-2$,
we have that  $\alpha_{j} = a_{sj}, \gamma_{j} = b_{sj}$ and 
so $y_{n+1} = (e_{s}^{T} \otimes I) Y^{(L)}$. 
Thus  assuming that we have a predictor of order $q-1$, we can generalize
the result of Jackson and  N\o rsett  \cite{r10jn88}
and Burrage \cite{bu7}  that the (global)
order of $y_{n+1}$ equals $p^{*} = \mbox{min}(p,q+L)$. 
Note that for  $F(Y^{(0)})= (e \otimes I) f(y_{n})$ one has $F(Y^{(0)}_{n}) =F(Y_{n})+O(h)$, i.e., 
$q=1$, and then we have

\begin{teo}.
Let ${A,B,\alpha,\gamma}$ define an $r$-step, $s$-stage MRK methods, and
denote $F(Y_{n}^{(i)})$ by $F_{n}^{(i)}$, then the PIMRK method defined by
\begin{eqnarray}
\begin{array}{lll}
F_{n}^{(0)} &    = &  (e \otimes I) f(y_{n}),  \\
Y_{n}^{(j)}  & = & (A \otimes I) y^{(n)} + h (B \otimes I) F_{n}^{(j-1)}, \quad  j=1, \ldots, L, \\
y_{n+1} &   = & (e_{s}^{T} \otimes I) Y_{n}^{(L)}, 
\end{array}
 \label{teo1} 
\end{eqnarray}
represents an $(L+1)s$-stage explicit MRK method of order
$p^{*} = min(p,L)$ requiring $L$ parallel stages. Here $e^{T} = (1,1,\ldots,1)$
and $ e_{s}^{T} = (0,0,\ldots,0,1).$
\end{teo}

\subsection{Predictor methods}
Unlike van der Houwen and Sommeijer \cite{hs90} who used predictor methods
by performing the auxiliary vector recursion we instead use predictor
methods which are based on past information. Note that Burrage \cite{bu6}
suggested predictor methods based only on past information on node points,
hence  the stability function is easily obtained. In addition,  
Jackiewicz and Tracogna \cite{jt88} studied a class of two-step Runge-Kutta
methods which depend on stage values at two consecutive steps. Their methods,
which do not belong in the class of MSRKs,
are quite efficient with respect to the number of function evaluations.
Our methods are more general than theirs for we are not limited
to two past stages values.

With this information we actually
could predict either $Y_{n}^{(0)}$ or  $F_{n}^{(0)}$. In the former case, we should
evaluate $F(Y_{n}^{(0)})$ then start the iteration, while in the latter case
we immediately start the iteration saving one function evaluation
on every step. We have conducted a number of experiments comparing these two
approaches. In terms of accuracy they are comparable, except for
predictors using the information on both the  node and internal points and then
the former case is more accurate than the latter. However, the latter
approach is more efficient in terms of the number of function evaluations
for predictors using various values of past information.  
Observing this phenomena,  we do not consider predicting $Y_{n}^{(0)}$.

The trivial predictor $P_{0}$ for $F_{n}^{(0)}$ is defined as
$(f(y_{n})^{T}, \ldots, f(y_{n})^{T})^{T}$  and gives an iterated method of 
order $p^{*} = min(p,L)$. 
We denote predictor $P_{a}$, where $a$ represents classes of predictor
using specific past information, with  $a$ taking the values $1,2$
, $3$ and $4$.  $P_{1}$ means the predictor uses derivatives on $r$ previous
node points,
$P_{2}$ indicates that the predictor uses both derivatives on $r$ previous 
node points and $s$ stage points in the previous step,
$P_{3}$ means the predictor uses only derivatives on $s$ stage points, and
$P_{4}$ indicates that the predictor uses derivatives on $s$ previous stage points
and $t_{n-1}$.

For example, we consider predictor $P_{1}$ which uses $r$ previous 
derivatives $f_{n}, \ldots, f_{n+1-r}$. We will determine the $s \times r$  matrix
$V_{1}$ such that
\begin{eqnarray}
F_{n}^{(0)} & = & (V_{1} \otimes I) f^{(n)},  \label{ep10}
\end{eqnarray}
where $f^{(n)}$ denotes $(f(y_{n})^{T}, \ldots, f(y_{n+1-r})^{T})^{T}$.
By expanding the solution and the derivatives at $t_{n}$ we obtain the order
conditions,
\begin{eqnarray}
p c^{p-1} & = & p V_{1}  \lambda^{p-1}, \label{ep2}
\end{eqnarray}
where $\lambda = (0, -1, \ldots, -r+1)^{T}$. It is obvious that
$F_{n}^{(0)}$ is of order $r-1$. Note that (\ref{ep2}) can be
simplified into $c^{p-1} = V_{1} \lambda^{p-1}$.
Since both $c$ and $\lambda$ are available, then $V_{1}$ can be obtained
by solving systems of linear equations  to give an $s \times r$ matrix.
For $P_{2}$, $P_{3}$ and $P_{4}$, similar systems of equations arise with  $\lambda$  defined as
$\lambda= (-1+c_{1}, \ldots, -1+c_{s-1}, 0, -1, \ldots, -r+1)^{T}$, 
$\lambda = (-1+c_{1}, \ldots, -1+c_{s-1}, 0)^{T}$, and 
$\lambda = ( -1+c_{1},\ldots, -1+c_{s-1},0,-1)^{T}$, respectively, where in 
each of the cases, 
$V_{i}, i=2,3,4$ is an $s \times (s+r-1), s \times s,$ and $s \times (s+1)$
matrix, 
respectively.

\subsection{Stability}
We consider now the linear stability of a PIMRK (based on the 
trivial predictor)  with respect to the test equation
\begin{equation}
\label{te}
y'(t)=\lambda y(t).
\end{equation}
It is obvious that the application of (\ref{teo1}) yields the  recursion
\begin{equation}
Y_{n}^{(L)} =( I +z B+z^{2}B^{2}+\ldots +z^{L-1}B^{L-1} ) A  y^{(n)} + z^{L} B^{L} e\, y_{n},
\end{equation}
where $z = \lambda h$. Let $m(z) = (I+zB + \ldots + z^{L-1}B^{L-1})A$ then
\begin{equation}
y_{n+1} = e_{s}^{T} m(z) y^{(n)} + e_{s}^{T}z^{L} B ^{L} e y_{n},
\end{equation}
and this takes the form of the recurrence relation
\begin{equation}
y_{n+1} - \alpha_{1}(z)y_{n} - \ldots - \alpha_{r}(z) y_{n+1-r} = 0,
\end{equation}
where 
\[
\alpha_{1}(z)=e_{s}^{T}z^{L}B^{L}\bar{e}+e_{s}^{T}m(z)\bar{e}, \quad
\alpha_{j}(z) = e_{s}^{T} m(z) \bar{e_{j}}, \quad j=2, \ldots, r.
\]
where $\bar{e}^{T} = (1,\ldots,1) \in R^{r}$ and $\bar{e_{j}}^{T}$
are the standard basis vectors for $R^{r}$. 
If we choose $L=p$, where $p$ is the order of the corrector then we obtain
a stability polynomial of degree $p$.
The stability regions of PIMRK methods
with pair values $(r,s)$ of $(2,2), (4,2)$
and the PIRK method of order 4 based on iterating the $s$ stage Gauss method after $L=3$ are given in Figure (\ref{sr}). The PIRK method appears to have
a slightly larger stability region than the PIMRK methods but there is little
difference.

\begin{figure}
\centerline{
\psfig{figure=gs2rs22p0.ps,height=3in,width=3in}}
\caption{Stability region for (2,2), (4,2), and PIRK, after $L=3$}
\label{sr}
\end{figure}

Now, we consider the effects of various non trivial predictors on  the
stability polynomial. Since we are using $s$ past derivatives values and
$r$ node points
for estimating $F_{n}^{(0)}$, we write
\begin{equation}
\label{ep}
F_{n}^{(0)} = ( W_{i}\otimes I) \, f(Z_{n}), \; \quad i=1,\ldots,4,
\end{equation} 
where $W_{i}$ depends on whether we are using
$P_{1}, P_{2}, P_{3}$ or $P_{4}$, and 

\[
Z_{n} = (Y_{n-1,1},\ldots, Y_{n-1,s},
y_{n-1},\ldots, y_{n+1-r})^{T} \in R^{s+r-1}.
\]
 
Note that $Y_{n-1,s}$ is $ y_{n}$ since the method is stiffly accurate.
 Applying (\ref{te}) we obtain
$F^{(0)} = \lambda W_{i}\, Z_{n}$. For the case of $P_{1}, P_{2}, P_{3}$
and $P_{4}$, respectively, $W_{i}$ is defined as
\[
\left[ \begin{array}{ll}
        0_{s \times (s-1)} & V_{1} \end{array} \right], \quad
\left[ \begin{array}{l}
         V_{2} \end{array} \right], \quad
\left[ \begin{array}{ll}
         V_{3} & 0_{s \times (r-1)}\end{array} \right], \quad
\left[ \begin{array}{ll}
         V_{4} & 0_{s \times (r-2)}\end{array} \right], \quad
\]
where $0_{i \times j}$ stands for zero matrix of dimension
 $i$ by $j$.
So after $L$ iterations, we rewrite the corrector as
\begin{equation}
Y_{n}^{(L)} = m(z) y^{(n)} + z^{L} B^{L} W_{i} Z_{n}, \quad i=1,\ldots,4.
\end{equation}
And this can be further rewritten as
\begin{equation}
Z_{n+1} = M_{i}(z) Z_{n}, \quad i=1,\ldots,4,
\label{l3.13}
\end{equation}
where 
\[
M_{i}(z) = 
\left[ \begin{array}{ll}
        0_{s \times (s-1)} & m(z) \\
        0_{(r-1) \times (s-1)} & J \end{array}
\right] + z^{L}\left[ \begin{array}{l} B^{L} W_{i} \\ 0_{(r-1) \times (s+r-1)} \end{array} \right], \quad i=1,\ldots, 4,
\]
and $J$ is a matrix of dimension $r-1 \times r$ as follows

\[
J = \left[ \begin{array}{lllll}
             1 & 0 &  \cdots & 0 & 0 \\
             0 & 1 &  \cdots & 0 & 0 \\
             0 & 0 &  \ddots & 0 & 0 \\
             0 & 0 &  \cdots & 1 & 0 \end{array}
    \right].
\]
The characteristic polynomial associated with (\ref{l3.13}) is 
\begin{equation}
p_{i}(\zeta,z) = det(I_{r+s-1}*\zeta - M_{i}(z)), \quad i=1,\ldots,4,
\end{equation}
where $I_{r+s-1}$ is the identity matrix of order $r+s-1$. We have 
computed 
the stability regions for $(r,s)=(2,2)$ after $L=3$ iterations using 
predictors
$P_{0}, P_{1}, P_{2}, P_{3} $ and $P_{4}$,  see figure (\ref{f2}). 
Predictors $P_{0}$ and $P_{1}$ give the largest stability regions in this case.
Note that the two separate lines in rightmost plots represent boundaries
from the Schur's criterion and can be ignored.

\begin{figure}
\centerline{\hbox{
\psfig{figure=p22p0p1.ps,height=2.5in,width=2.5in,angle=270}
\psfig{figure=p22p2.ps,height=2.5in,width=2.5in,angle=270}}}
\centerline{\hbox{
\psfig{figure=p22p3.ps,height=2.5in,width=2.5in,angle=270}
\psfig{figure=p22p4.ps,height=2.5in,width=2.5in,angle=270}
}}
\caption{(2,2) with $P_{0}, P_{1}, P_{2}, P_{3}$ and $P_{4}$}
\label{f2}
\end{figure}      

\subsection{The rate of convergence of the inner iterations}
Following the definition of the rate convergence of the inner iterations
by Cong \& Mitsui \cite{cong1}, the convergence bound of a PMIRK method
is defined in a similar way as for methods proposed in \cite{cong2,cong3,hcong}.
We use a test equation $y'(x) = \lambda y(t)$, where $\lambda$ runs 
through the eigenvalues of the Jacobian matrix $\partial f/ \partial y$,
then we obtain the iteration error equation

\begin{equation}
Y^{(j)}_{n} - Y_{n} = z B [ Y^{(j-1)}_{n} - Y_{n}], \quad z = \lambda h, \quad
j =1, \ldots, L.
\label{congh1}
\end{equation}

So the spectral radius $\rho(B)$ of the matrix $B$ dominates the rate 
of convergence. We call $\rho(B)$ the convergence factor of the PIMRK method.
Requiring that $\rho(z B) < 1$, gives us the convergence condition in the inner
iterations

\begin{equation}
| z | < \frac{1}{\rho(B)} \quad \mbox{or} \quad h < \frac{1}{\rho(B)\rho(\partial f / \partial y)}.
\end{equation}

The convergence factors for particulars PIMRK methods are discussed in 
the next two section. 

\section{Implementation considerations}
In this section we describe a simple strategy for implementing PIMRK with 
a variable stepsize in order to control the local truncation error. 
$r-1$ starting procedures based on IRK of Radau type of order $2s^{*}-1$
are used to compute $y_{1}, \ldots, y_{r-1}$ where $s^{*}$ is determined so that it is at least
 the order of the main MRK method, i.e.,  $2s+r-2$. Hence $s^{*} =
\lceil \frac{2s+r-1}{2} \rceil$.

Our
strategy is similar to the one implemented by van der Houwen and Sommeijer \cite{hs90}
which is also implemented in DOPRI8. 
If the trivial predictor is used then we define $L=p$ so that $y_{n+1}$ is of $O(h^{p+1})$, since in each
iteration we have $Y_{s}^{(j)}$, which is the approximation of $y_{n+1}$ of
$O(h^{j+1})$, then at the end of iteration we have $Y_{s}^{(p)}$ of
$O(h^{p+1})$. As the result we define 
an estimate of local error $\epsilon$ from step $n$ to step $n+1$
as
\begin{equation}
\epsilon  = \| y_{n+1} - Y_{s}^{(p-1)} \|,
\label{est}
\end{equation}
for some norm $\| . \|$. It is obvious that $\epsilon = O(h^{p})$.
As a matter of fact, one could consider Richardson's strategy for
estimating local error, but our few experiments indicated that this
approach is more expensive than (\ref{est}).

In the fixed iteration implementation,
a step is accepted when $\epsilon \leq TOL$ and rejected otherwise. 
The new step size is chosen as
\begin{equation}
\label{step}
h_{\mbox{new}} = h \times \mbox{min} \{ 6, \mbox{max}\{ \frac{1}{3},
0.9(\frac{TOL}{\epsilon})^{\frac{1}{p}}\}\}.
\end{equation}
The constants $6$ and $\frac{1}{3}$ in the expression are used to prevent
an abrupt change in the stepsize and the safety factor $0.9$ is added
to increase the probability that the next step will be accepted. 

The PIMRK method also allows us to implement variable order techniques. This 
can be done by avoiding fixing the number of iterations and using two successive
approximation to the solution as the estimate of local error, i.e,
\[
\epsilon^{(j)} = \| Y_{s}^{(j)} - Y_{s}^{(j-1)} \|.
\]
If during iteration the condition $\epsilon^{j} \leq TOL$ is satisfied
for some $j=j_{0}$ then accept $Y_{s}^{(j_{0})} $ as the numerical
solution $y_{n+1}$.  The next step size will be determined as in (\ref{step})
but $\epsilon^{j_{0}}$ and $p^{*}=\mbox{min}\{p,q+j_{0}\}$ is used 
instead of $\epsilon$ and $p$. If the tolerance condition is not
satisfied after $M$ iterations, the step will be rejected, $h$
is redefined and the step is restarted. In the implementation
we define $M=2s+r-2$. In this way a variable order PIMRK
can be implemented.


\section{Numerical experiments}
In this section we present some numerical experiments which are done
in Matlab, on a Sun Sparc2 workstation, whose unit roundoff error is $2.22e^{-16}$. The coding in Matlab is convenient  since
one can perform vector and matrix operation very simply hence
requiring less effort compared with coding in a conventional programming language, such as Fortran.
The performance of parallel computation is estimated by observing
the number of function evaluations per processor.

\subsection{Convergence factors}
Table (\ref{tcf}) is the result of a direct numerical 
computation on the convergence
factor of PIMRK methods defined on section 4. Our convergence factors
are slightly larger than those of PIRK methods (a class of
iterated Runge-Kutta methods, see van der Houwen and Sommeijer \cite{hs90})  and PITRK (
a class of iterated two step Runge-Kutta methods, see Cong and
Mitsui \cite{cong1}). However our numerical
experiments indicated that these differences do not greatly affect the overall 
performance of the method.

% table goes here..

\begin{table}[p]
\caption{Convergence factors for various PIMRK, PIRK, PTIRK methods.}
\label{tcf}
\begin{center}
\begin{tabular}{lllllllll}
\hline
order   & p=4 & p= 5 & p= 6 & p= 6 & p= 8 & p= 9 & p= 10 & p= 12 \\
\hline
PIMRK   &(2,2)&(3,2) &(4,2) &(2,3) & (4,3)& (3,4)& (4,4) & (6,4) \\
        &.357 &.283  &.330  & .276 & .219 & .188 & .175  & .159 \\
PIRK    &.289 &      &.     & .215 & .165 &      & .137  &      \\
PITRK   &.193 &      &      & .136 & .106 &      & .085  &  \\
\hline
\end{tabular}
\end{center}
\end{table}

One way of reducing the convergence factors of PIMRK methods is by reducing
the order of the methods. For example, for $r=3, s=2$  we have a method
of order $5$, but reducing it to order $4$ we will have one free parameter
$c_{1}$. The parameter is computed in such a way that all eigenvalues
of a matrix $\hat{B}$ of the resulting method are minimum. 
In fact by assuming order $s+r-1$ (with $c_{s}=1$) there are $s-1$ free
parameters $c_{1}, \ldots, c_{s-1}$.
We performed this approach on methods
$(r,s) = (3,2)$ and $(4,2)$ which are of order $5$ and $6$. 
The convergence factors of these methods after reducing to order $4$ and $5$
are respectively $.257$ and $.273$. They are indeed smaller than those
convergence factors of the original methods, see Table (\ref{tcf}). 

We performed numerical experiments by 
using the original methods $(3,2)$ and $(4,2)$
and their companion methods of order one lower with variable stepsize and 
variable
order strategies on several test problems such as problem
A1, A2, A3, Fehlberg from DTEST \cite{dtest}, and a constant problem TC
$ y_{1}' = 0.5, y_{2}'=100$, 
 with various values of
tolerances. Tables (\ref{z1},\ref{z2},\ref{z3},\ref{z4},\ref{z5}) show the 
number of (sequential) function evaluations required using the predictor
$P_{3}$. 
In general, the results show that there is no significant improvement obtained 
by the companion methods. Several cases in which the companion methods are
 more efficient  are in Table \ref{z1} for method $(3,2)$ at log tolerance
$-4$, Table \ref{z2} for method $(3,2)$ at log tolerance
$-4, -6, -8$ and for method $(4,2)$ at log tolerance
$-8, -10, -12$, Table (\ref{z3}) for method $(3,2)$
at log tolerance  $-4, -6, -12$ and for method $(4,2)$ at log tolerance
$-4, -6$, Table (\ref{z4}) for method $(3,2)$ at
log tolerance $-4, -6$ and for method $(4,2)$ at 
log tolerance $-4, -6$, and $-8$, and Table (\ref{z5}) for method
$(3,2)$ at log tolerance $-4,-8$ and for method $(4,2)$
at log tolerance $-8$. The star $*$ symbol at Table
(\ref{z5}) indicate that the integration is  aborted due to
too many steps required. Note that in table (\ref{z4}) since the original
method is higher order than the companion and the problem solved is 
a constant ODE, its costs is independent on tolerances. Note also that
in table (\ref{z5}), we integrate the problem on a longer interval than
in table (\ref{t4}), hence the costs on the table differ significantly.
The results also indicate that
the convergence factors do not affect too much on the performance
of the method, for example, method $(3,4)$ of order $9$
with convergence factor $.188$ is less efficient than method
$(6,3)$ of order $10$ with convergence factor $.249$. The stability
of the method is more dominant than the convergence factor. 


\subsection{Comparison among $(r,s)$}
We have solved Van der Pol's equation
\begin{eqnarray*}
y_{1}' & = & y_{1}(1-y_{2}^{2})-y_{2}, \\
y_{2}' & = & y_{1},
\end{eqnarray*}
where the interval of integration is $[0 ,20]$ and $y_{0} = (0 \;, .25)^{T}$. 
Several $(r,s)$ pairs for
$r$-steps $s$-stages PIMRK are implemented with the trivial predictor.
 
Table (\ref{t2}) is the result
of a fixed number of iterations per step, while Table (\ref{t3}) is the result of the variable order implementation.  Various values of tolerances are used.
The character * in the tables indicates
that the program is aborted due to too many steps required.

For the sequential cost, we see in Table (\ref{t2}) that higher order 
methods $(5,3)$ and $(7,3)$
with order $9$ and $11$, respectively are suitable for stringent tolerances,
but at less stringent tolerance the lower order methods such as $(4,2)$ are preferred. 
Meanwhile  in 
parallel costing, the high order methods such as
$(4,4)$ and $ (6,4)$ are preferred
for all values of tolerances. 


We observe from Table (\ref{t3}) that the variable order implementation of PIMRK works
as expected, it is more efficient than the fixed number iteration implementation.
The results also shows that higher order PIMRK such as $(6,4),
(7,3)$ and $(7,4)$ are suitable both for
sequential and parallel implementation with fair to stringent tolerances. 

\subsection{Comparing predictors}
We considered the accucary and the
efficiency of the methods. We solved problem A1, A2, A3, and A4 of
\cite{hefs72}. For accuracy comparison, fixed stepsize $h$ and fixed iteration methods are used.
In the
test we use $r$-step $s$-stages PIMRKs with values of $(r,s)$ as $(4,2), (5,3)$
and $(7,3)$. 
 In Tables (\ref{a1}), (\ref{a2}), (\ref{a3}) and (\ref{a4})  we listed values 
TP,
NFS,
NFP, 
and D which, respectively, denote the type of predictor, the number of function evaluations
for sequential machines, the number of function evaluations for
parallel machines, and the number of similar digits to the exact solution at
the end point of integration (absolute accuracy). 
Comparing experiments at fixed stepsize $h$, we observed 
that all type of predictors have comparable accuracy,
except  $P_{2}$ which is the least accurate amongst the other type
predictors. 


Next we implement variable order methods with stepsize selection strategy
described earlier.
On the test,  we compare the performance of low order method $(4,2)$, higher
order methods $(5,3)$ and $(7,3)$,  Gauss3, Gauss5 (van der Houwen  and
Sommeijer \cite{hs90}) and RK78 (Dormand and Prince \cite{hvol1}).
The results for 
the Fehlberg problem are given in Table (\ref{t4}),
the results for the motion of rigid body 
without external forces, problem B5 \cite{hefs72}, are presented
in Table (\ref{t6}),
and the results to the Orbit problem, problem D2 \cite{hefs72}
are presented 
in Table (\ref{t7}).
These tables also show that the low order methods 
are suitable for fair tolerances while higher order methods
are suitable for stringent tolerances. 
In a sequential environment the $(7,3)$ methods are as effective as
the RK78 pair at very stringent tolerances, while for a parallel implementation
a number of our methods are comparable with Gauss5 methods. The advantage
being that they are require fewer processors than the $5$ processors required
for the implementation of Gauss5.
The results also
indicate that $P_{3}$ is the most robust predictor.

\section{Conclusion}
We have designed parallel iterated methods based on multistep Runge-Kutta
of  Radau type.
The methods are implemented, and both sequential and parallel costs are
examined. The sequential cost indicates that the maximum performance
are attained by methods $(5,3)$ and $(7,3)$. The performance of
methods with orders higher than these are less efficient. 
In general a method $(r_{2},s_{2})$
is less efficient than $(r_{1},s_{1})$ where $s_{2} > s_{1}$  even though
they are the same order, this is due to the fact that $(r_{2},s_{2})$
requires more function evaluation than $(r_{1},s_{1})$,
see for example $(4,2)$ and $(2,3)$, $(2,4)$ and $(4,3)$,
and $(3,4)$ and $(5,3)$ which are respectively, of orders
$6$, $8$, and $9$.  As we expected
the higher methods are preferred for stringent tolerances, in particular
those with predictor $P_{3}$ and $P_{4}$ have sequential
performances comparable to RK78. In the Table (\ref{t4}), at tolerance
$10^{-14}$ for instance,
their sequential performances are even better than those of RK78. 

We have also seen  (from Tables  (\ref{t2}) and (\ref{t3}), for example)
that the variable order implementation is both more efficient and 
robust than the corresponding fixed order implementation at all
levels of tolerance and suggests that the variable order implementation
is the method of choice. In this case the $(6,4)$ and $(7,4)$
variable order methods (which have maximum order $12$ and $13$ respectively)
appear to be the most suitable candidate both sequentially and in parallel. 

In the case of parallel cost estimates, we observed that even though $s_{2} > s_{1}$,
method $(r_{2},s_{2})$ is mostly better than $(r_{1},s_{1})$. Unlike the performance
in sequential cost, the higher order methods gain more efficiency than lower
methods, this motivates us to continue our works on real parallel implementation
of the methods so that we could find suitable methods $(r,s)$ for
a multiprocessor environment, in which, for example the 
the computation of the right-hand side is dominant and costly.
This  will be considered in a later paper.

\begin{thebibliography}{99}

\bibitem{bumos3} Burrage, K. \& Moss, P. M., (1980),
{\em Simplifying assumptions for the order of partioned multivalue
methods}, BIT {\bf 20}, 452--465.

\bibitem{bu5} Burrage, K. (1985),
{\em Order and stability of explicit multivalue methods},
Appl. Numer. Math., {\bf 1}, 363--379.

\bibitem{bu4} Burrage, K. (1988),
{\em Order properties of multivalue methods}, IMA J. Numer. Anal. {\bf 8},
43--69.


\bibitem{bu6} Burrage, K. (1993),
{\em The Search for the Holy Grail, or: Predictor-Corrector methods
for solving ODEIVPs}, App. Num. Math. {\bf 11}, 125--141.

\bibitem{bu7} Burrage, K. (1995),
{\em Parallel and sequential methods for Ordinary Differential Equations},
Oxford University Press, New York.

\bibitem{but6} Butcher, J.C. (1966),
{\em On the convergence of numerical solutions to ordinary differential
equations}, Math. Comp., {\bf 20}, 1--10.

\bibitem{but7} Butcher, J. C. (1973),
{\em The order of numerical methods for ordinary differential equations},
Math. Comp., {\bf 27}, 793--806.

\bibitem{cong1} Cong, Nguyen huu \&  Mitsui, Taketomo (1996), 
{\em A class of explicit parallel two-step Runge-Kutta methods},
in preparation.

\bibitem{cong2} Cong, Nguyen huu (1994), 
{\em Parallel iteration of symmetric Runge-Kutta methods for non stiff
initial value problems}, J. Comput. Appl. Math., {\bf 51}, 117--125.

\bibitem{cong3} Cong, Nguyen huu (1995),
{\em Explicit parallel two-step Runge-Kutta-Nystrom methods},
to appear in Comput. Math. Appl.

\bibitem{cop8} Cooper, G. J. (1978),
{\em The order of convergence of general linear methods for
ordinary differential equations}, SIAM J. Numer. Anal, {\bf 15}, 643--661.

\bibitem{dtest} Enright, W.H \& Pryce, J.D. (1987),
{\em Two Fortran Packages for Assessing Initial Value Methods}, 
ACM Transactions on Mathematical Software,
{\bf 13}, 1--27. 

\bibitem{gs69} Guillou, A. \& Soul\'{e}, J.L. (1969),
{\em La r\'{e}solution num\'{e}rique des probl\`{e}mes diff\'{e}rentiels
aux conditions initiales par des m\'{e}thodes de collocation},
R.I.R.O., {\bf R-3}, 17--44.

\bibitem{hw9} Hairer, E. \& Wanner, G. (1974),
{\em On the Butcher group and general multivalue methods},
Computing, {\bf 13}, 1--15.

\bibitem{hvol1} Hairer, E. \& Wanner, G. (1987),
{\em Solving Ordinary Differential Equations I:  Nonstiff Problems},
Springer Series in Comp. Math.,
Springer-Verlag, Berlin.

\bibitem{hnw80} Hairer, E., N\o rsett, S. P.,  \& Wanner, G. (1991),
{\em Solving Ordinary Differential Equations II: Stiff and
Differential-Algebraic Problems}, Springer Series in Comp. Math.,
Springer-Verlag, Berlin.

\bibitem{hs88} Houwen, P. J. van der \& Sommeijer, B. P. (1988),
{\em Variable step integration of high order Runge-Kutta methods on 
parallel computers}, Rep. NM-R8817, CWI, Amsterdam, The Netherlands.

\bibitem{hs90} Houwen, P.J. van der \& Sommeijer, B.P. (1990),
{\em Parallel iteration of high-order Runge-Kutta methods with
stepsize control}, J. Comput. Appl. Math. {\bf 29}, 111--127.


\bibitem{hcong} Houwen, P. J. van der \& Cong, Nguyen huu (1993),
{\em Parallel block predictor-corrector methods of Runge-Kutta type},
Appl. Numer. Math., {\bf 13}, 109--123.


\bibitem{hefs72} Hull, T.E., Enright, W.E., Fellen, B.M.
\& Sedgwick, A.E. (1972), {\em  Comparing numerical methods for
ordinary differential equations}, SIAM J. Numer. Anal. {\bf 9},
603--637.

\bibitem{r9in90} Iserles, A. \& N\o rsett, S.P. (1990),
{\em On the theory of parallel Runge-Kutta methods}, IMA J. Numer. Anal.
{\bf 10}, 463--488.

\bibitem{jt88} Jackiewicz, Z. \& Tracogna, S. (1988),
{\em A general class of two-step Runge-Kutta methods for ordinary
differential equations},
SIAM J. Numer. Anal. {\bf 1}, 1--38.

\bibitem{r10jn88} Jackson, K.R. \& N\o rsett, S.P. (1995),
{\em The potential for parallelism in Runge-Kutta methods. Part I:
RK formula in standard forms }, SIAM J. Numer. Anal. {\bf 32}, 49--82.

\bibitem{k59} Krylov, V.I. (1959),
{\em Priblizhennoe Vyschisslenie Integralov}, Goz. Izd. Fiz.-Mat, Lit.,
Moscow. English translation: Approximate calculation of 
integrals. Macmillan, New York, 1962.

\bibitem{r11lie87} Lie, I. (1987),
{\em Some aspects of parallel Runge-Kutta methods},
Report No. 3/87, University of Trondheim, Division Numerical
Mathematics, Norway.

\bibitem{ln89} Lie, I. \& N\o rsett, S.P. (1989),
{\em Superconvergence for multistep collocation},
Math. of Comput., {\bf 52}, 65--79.

\bibitem{r12ns89} N\o rsett, S.P. \& Simonsen, H.H. (1989),
{\em Aspects of parallel Runge-Kutta methods}, in A. Bellen (ed):
Workshop on Numerical Methods for Ordinary Differential Equations,
L`Aquila, 1987, Lecture Notes in Mathematics, Vol. 1386, Springer-Verlag,
Berlin, 103--117.

\bibitem{sch93} Schneider, S. (1993), 
{\em Numerical experiments with a multistep Radau method},
BIT {\bf 33}, 332--350.

\end{thebibliography}

\appendix
\section{Appendix: method parameters}
\subsection{MRK with $s=2, r=2$}
\begin{flushleft}
\begin{tabular}{cc}
\hline
0.39038820  & 1.00000000 \\
\hline 
 1.04671555 & -0.04671555 \\
1.02010510 & -0.02010510 \\
\hline
0.40044075 & -0.05676810\\
 0.77072386 & 0.20917105\\
\hline
\end{tabular}
\end{flushleft}

\subsection{MRK with $s=3,r=2$}
\begin{flushleft}
\begin{tabular}{ccc}
\hline
0.17789172 & 0.67323526  & 1.00000000 \\
\hline
1.00622582  &  -0.00622582 & \\  
 0.99832216 &  0.00167784 & \\  
1.00092398  &  -0.00092398 &   \\
\hline
 0.20550986 & -0.04964153 & 0.01579757\\
 0.43921211 & 0.26945763 & -0.03375664\\
 0.41315606 & 0.48645093 & 0.09946902\\
\hline
\end{tabular}
\end{flushleft}

\subsection{MRK with $s=3,r=3$}
\begin{flushleft}
\begin{tabular}{ccc}
\hline
0.19216964 & 0.68931797 & 1.00000000 \\
\hline
0.20964870 & -0.04220518 & 0.01246369\\
0.46587564 & 0.25655438 & -0.02988950\\
0.43415567 & 0.47073597 & 0.09305311\\
\hline
 1.013199 &  -0.01413601 & 0.00093679 \\
0.99626059 & 0.00425626 & -0.00051685 \\
1.00211147 &  -0.00216768 & 0.00005621\\
\hline
\end{tabular}
\end{flushleft}

\newpage
\subsection{MRK with $s=2,r=4$}
\begin{flushleft}
\begin{tabular}{cccc}
\hline
0.44784375 &  1.00000000 & & \\
\hline
1.14343426 & -0.18153852 &  0.04369589 & -0.00559162 \\
1.06570627 & -0.07647305 & 0.01189055 &  -0.00112376 \\
\hline
0.37688667 &  -0.03996454 & &  \\
0.76866717 & 0.17526959 & & \\
\hline
\end{tabular}
\end{flushleft}

\subsection{MRK with $s=3,r=5$}
\begin{flushleft}
\begin{tabular}{ccccc}
\hline
0.21139546 &  0.70879842 &  1.00000000 & & \\
\hline
1.02761716 & -0.03323146 & 0.00669661 & -0.00119344 & 0.00011114 \\
0.99181431 & 0.01152114 & -0.00411899 & 0.00087188 & -0.00008834 \\
1.00488332 & -0.00526501 & 0.00041916 & -0.00003976 & 0.00000229 \\
\hline
0.00011114 & -0.00119344 & 0.00669661\\ 
-0.00008834 & 0.00087188 & -0.00411899\\ 
0.00000229 & -0.00003976 & 0.00041916\\  
\hline
\end{tabular}
\end{flushleft}

\subsection{MRK with $s=3,r=7$}
\begin{flushleft}
\begin{tabular}{ccccccc}
\hline
0.22489755 & 0.72107268 &  1.00000000 & & & & \\
\hline
1.04191312 & -0.05559195 & 0.01884687 & -0.00674189 & 0.00188717 & -0.00034292 & 0.00002959 \\
0.98734001 & 0.02115880 & -0.01244334 & 0.00526817 & -0.00160232 & 0.00030596 & -0.00002729  \\
1.00793593 & -0.00894651 & 0.00120534 & -0.00023009 & 0.00003998 & -0.00000494 & 0.00000031\\ 
\hline
0.21523002 & -0.03021097 & 0.00776626 & & & & \\
0.52080767 & 0.23077049 & -0.02347204 & & & & \\
0.47508152 & 0.43681502 & 0.08101438 & & & & \\
\hline
\end{tabular}
\end{flushleft}

\begin{table}[p]
\caption{The number of sequential function evaluations for problem A1}
\label{z1}
\begin{center}
\begin{tabular}{|l|l|lllll|}
\hline
Method & Order & \multicolumn{5}{|c|}{-log(Tol)} \\
       &       & 4   & 6   & 8   & 10   & 12   \\ \hline
(3,2)  & 5     & 155 & 289 & 749 & 1921 & 4843 \\ \cline{2-7}
       & 4     & 127 & 395 & 1141& 4085 & 13959\\ \hline
(4,2)  & 6     & 138 & 258 & 490 & 1046 & 2254 \\ \cline{2-7}
       & 5     & 192 & 312 & 708 & 1766 & 4394 \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{The number of sequential function evaluations for problem A2}
\label{z2}
\begin{center}
\begin{tabular}{|l|l|lllll|}
\hline
Method & Order & \multicolumn{5}{|c|}{-log(Tol)} \\
       &       & 4    & 6   & 8   & 10   & 12   \\ \hline
(3,2)  & 5     & 93   & 193 & 275 & 493 & 1145 \\ \cline{2-7}
       & 4     & 89   & 179  & 247  & 505  & 1045 \\ \hline
(4,2)  & 6     & 104  & 170  & 282  & 530  & 958 \\ \cline{2-7}
       & 5     & 120  & 194  & 256  & 472  & 868  \\ \hline
\end{tabular}
\end{center}
\end{table}
 
\begin{table}[p]
\caption{The number of sequential function evaluations for problem A3}
\label{z3}
\begin{center}
\begin{tabular}{|l|l|lllll|}
\hline
Method & Order & \multicolumn{5}{|c|}{-log(Tol)} \\
       &       & 4    & 6   & 8   & 10   & 12   \\ \hline
(3,2)  & 5     & 367  & 679 & 1125  & 2165  & 4667 \\ \cline{2-7}
       &   4   & 343 & 639  & 1143  & 2255  &  4543 \\ \hline
(4,2)  & 6     & 348  & 634  & 1072  & 1974  & 3720 \\ \cline{2-7}
       & 5     & 340 & 624  & 1196  & 2212  & 4398  \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{The number of sequential function evaluations for problem TC}
\label{z4}
\begin{center}
\begin{tabular}{|l|l|lllll|}
\hline
Method & Order & \multicolumn{5}{|c|}{-log(Tol)} \\
       &       & 4    & 6   & 8   & 10   & 12   \\ \hline
(3,2)  & 5     & 19 & 19  &19  & 19  & 19  \\ \cline{2-7}
       & 4     & 17 & 17  & 37  & 19 & 19  \\ \hline
(4,2)  & 6     & 26  & 26  & 26  & 26  & 26  \\ \cline{2-7}
       & 5     & 24 & 24  & 24  & 52  & 64  \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{The number of sequential function evaluations for Fehlberg's problem} 
\label{z5}
\begin{center}
\begin{tabular}{|l|l|lllll|}
\hline
Method & Order & \multicolumn{5}{|c|}{-log(Tol)} \\
       &       & 4    & 6   & 8   & 10   & 12   \\ \hline
(3,2)  & 5     & 8112 & 15971  & 32233  & *  &  *\\ \cline{2-7}
       & 4     & 7940  & 17351  & 28843  & *  & * \\ \hline
(4,2)  & 6      & 7320 & 15976  & 30136  & & \\ \cline{2-7}
       & 5     & 8550 & 17074  & 28266  & *  &  * \\ \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS and NFP by fixed iteration and trivial predictor for Van der Pol's problem}
\label{t2}
\begin{center}
\begin{tabular}{|l|r|rr|rr|rr|rr|rr|rr|} \hline \hline
\multicolumn{2}{|c|}{-log(TOL)}  & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{6} & \multicolumn{2}{c}{8}
& \multicolumn{2}{c}{10} & \multicolumn{2}{c}{12}
& \multicolumn{2}{c|}{14} \\  \hline
Method & P & NFS&  NFP & NFS & NFP& NFS & NFP & NFS & NFP & NFS & NFP 
& NFS & NFP \\ \hline
(22) &4 &1118 & 558 &2942 &1470 &8516 &4255 &* & * &* & *&*&*\\
(32) &5 &1051 & 523 &2227 &1111 &4855 &2423 & 11611 &5799 &* & *&*&*\\
(42) &6 & 930 & 443 &1636 & 802 &2806 &1387 &5080 &2518 & 10200 &5078&21514&10729\\
(23) &6 &1556 & 518 &2846 & 948 &5141 &1713 & 10010 &3334 & 20795 &6929&44159&14715\\
(33) &7 &1257 & 417 &2013 & 669 &3327 &1107 &5547 &1845 & 10065 &3351&18567&6183\\
(24) &8 &2254 & 563 &3318 & 829 &5334 &1333 &8582 &2145 & 13818 &3454&23490&5870\\
(43) &8 &1256 & 400 &1802 & 582 &2497 & 819 &4051 &1337 &6634 &2198&11021&3655\\
(34) &9 &1523 & 379 &2227 & 555 &3699 & 923 &5299 &1323 &7859 &1963&12059&3011\\
(53) &9 &1165 & 365 &1829 & 581 &2629 & 853 &3933 &1293 &5709 &1885&8653&2861\\
(63) & 10 &1392 & 408 &2175 & 669 &2823 & 885 &4005 &1289 &5868 &1910&8433&2765\\
(44) & 10 &1516 & 367 &2236 & 547 &3064 & 754 &4180 &1033 &6016 &1492&8464&2104\\
(73) & 11 &1477 & 427 &2197 & 667 &2887 & 887 &3757 &1187 &5137 &1657&7207&2347\\
(64) & 12 &* & * &2714 & 638 &3418 & 814 &4254 &1023 &5266 &1276&6982&1705\\
(74) & 13 &* & * &2383 & 547 &3391 & 799 &4495 &1075 &5359 &1291&7087&1723\\
GS3&6&1305&475&2519&919&4576&1686&8646&3216&17626&6586&36981&13851\\
GS5&10&1784&380&2515&535&3522&750&5216&1112&7646&1634&10775&2315\\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS and NFP by variable order and trivial predictor  for Van der Pol's problem}
\label{t3}
\begin{center}
\begin{tabular}{|l|r|rr|rr|rr|rr|rr|rr|} \hline \hline
\multicolumn{2}{|c|}{-log(TOL)}  & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{6} & \multicolumn{2}{c}{8}
& \multicolumn{2}{c}{10} & \multicolumn{2}{c}{12} & \multicolumn{2}{c|}{14} \\  \hline
Method & P & NFS&  NFP & NFS & NFP& NFS & NFP & NFS & NFP & NFS & NFP &NFS & NFP\\ \hline
(22)&4&1109&554&2940&1469&8516&4255&*&*&*&*&*&*\\
(32)&5&951&474&2225&1110&4855&2423&11611&5799&*&*&*&*\\
(42)&6&802&394&1562&771&2768&1371&5078&2517&10200&5078&21514&10729\\
(23)&6&1250&417&2733&911&5131&1710&10010&3334&20792&6928&44156&14714\\
(33)&7&1044&348&1979&659&3280&1092&5544&1844&10062&3350&18567&6183\\
(24)&8&1721&431&3038&760&5155&1289&8405&2101&13746&3436&23430&5855\\
(43)&8&934&308&1627&537&2599&859&4084&1350&6607&2189&10997&3647\\
(34)&9&1153&289&2023&506&3597&899&5137&1283&7711&1926&11983&2992\\
(53)&9&905&297&1543&507&2490&820&3778&1244&5724&1890&8662&2864\\
(63)&10&933&300&1596&516&2496&811&3756&1221&5742&1878&8190&2689\\
(44)&10&1102&274&1748&434&2914&724&3774&936&5684&1412&8426&2096\\
(73)&11&940&300&1519&487&2377&767&3409&1099&4897&1589&7036&2296\\
(64)&12&1019&248&1630&397&2569&628&3115&757&4898&1199&6569&1613\\
(74)&13&1101&267&1615&391&2493&606&3293&797&4715&1148&6413&1568\\
GS3 &6&1016 &372 &1941 &705 &3577 & 1297 &5868 & 2136 & 10537 & 3853 & 19282 & 7078\\
GS5 &10&1171 &255 &2021 &441 &3170 &694 &4103 &883 &6102 & 1310 &8818 & 1882\\
RK78 &8& 702 & 702&1027 &1027&1482 &1482&2288 &2288&3614 &3614&5876 &5876 \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS, NFP and D for problem A1, * denotes a negative value}
\label{a1}
\begin{center}
\begin{tabular}{|r|l|rrr|rrr|}
\hline \hline
    &  & \multicolumn{3}{c|}{$h^{-1}$=2} & \multicolumn{3}{c|}{$h^{-1}$=4} \\
$(r,s)$ &  TP &  NFS   & NFP    &  D   &   NFS &  NFP &  D \\ \hline
$(42)$ & 0  &    446 &    207 &  10.4 &    846 &    407 &  11.4\\ 
&  1  &    224 &     96 &  10.1 &    384 &    176 &  11.7\\ 
&  2  &    158 &     63 & * &    238 &    103 & *\\ 
&  3  &    374 &    171 &  10.4 &    694 &    331 &  11.4\\ 
&  4  &    302 &    135 &  10.4 &    542 &    255 &  11.4\\  \hline
$(53)$ & 0  &   1029 &    325 &  12.9 &   1989 &    645 &  14.5\\ 
&  1  &    597 &    181 &  12.6 &   1077 &    341 &  13.8\\ 
&  2  &    399 &    115 &  * &    639 &    195 &  10.3\\ 
&  3  &    819 &    255 &  12.9 &   1539 &    495 &  14.5\\ 
&  4  &    714 &    220 &  12.7 &   1314 &    420 &  14.4\\  \hline
$(73)$ & 0  &   1387 &    407 &  15.9 &   2587 &    807 &  17.1\\ 
&  1  &    775 &    203 &  12.4 &   1255 &    363 &  14.0\\ 
&  2  &    595 &    143 &  *  &    835 &    223 &   4.0\\ 
&  3  &   1189 &    341 &  15.9 &   2149 &    661 &  17.1\\ 
&  4  &   1090 &    308 &  14.7 &   1930 &    588 &  17.8\\  \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS, NFP and D for problem A2, * indicates a negative value}
\label{a2}
\begin{center}
\begin{tabular}{|r|l|rrr|rrr|}
\hline \hline
    &  & \multicolumn{3}{c|}{$h^{-1}$=2} & \multicolumn{3}{c|}{$h^{-1}$=4} \\
$(r,s)$ &  TP &  NFS   & NFP    &  D   &   NFS &  NFP &  D \\ \hline
$(42)$ & 0  &    446 &    207 &   5.3 &    846 &    407 &   6.0\\ 
&  1  &    224 &     96 &   4.4 &    384 &    176 &   5.1\\ 
&  2  &    158 &     63 &   * &    238 &    103 &   3.6\\ 
&  3  &    374 &    171 &   5.3 &    694 &    331 &   6.0\\ 
&  4  &    302 &    135 &   5.3 &    542 &    255 &   6.0\\  \hline
$(53)$ & 0  &   1029 &    325 &   8.4 &   1989 &    645 &   9.3\\ 
&  1  &    597 &    181 &   7.3 &   1077 &    341 &   8.1\\ 
&  2  &    399 &    115 &   4.2 &    639 &    195 &   4.6\\ 
&  3  &    819 &    255 &   8.4 &   1539 &    495 &   9.3\\ 
&  4  &    714 &    220 &   8.4 &   1314 &    420 &   9.2\\  \hline
$(73)$ & 0  &   1387 &    407 &  10.2 &   2587 &    807 &  10.9\\ 
&  1  &    775 &    203 &   7.7 &   1255 &    363 &   8.3\\ 
&  2  &    595 &    143 &   5.2 &    835 &    223 &   4.6\\ 
&  3  &   1189 &    341 &  10.2 &   2149 &    661 &  10.9\\ 
&  4  &   1090 &    308 &  10.2 &   1930 &    588 &  10.9\\  \hline
\end{tabular}
\end{center}
\end{table}
 
\begin{table}[p]
\caption{Values of NFS, NFP and D for problem A3, * indicates a negative value}
\label{a3}
\begin{center}
\begin{tabular}{|r|l|rrr|rrr|}
\hline \hline
    &  & \multicolumn{3}{c|}{$h^{-1}$=2} & \multicolumn{3}{c|}{$h^{-1}$=4} \\
$(r,s)$ &  TP &  NFS   & NFP    &  D   &   NFS &  NFP &  D \\ \hline
$(42)$ & 0  &    446 &    207 &   2.0 &    846 &    407 &   3.2\\ 
&  1  &    224 &     96 &   1.5 &    384 &    176 &   3.4\\ 
&  2  &    158 &     63 &  * &    238 &    103 &  *\\ 
&  3  &    374 &    171 &   2.0 &    694 &    331 &   3.1\\ 
&  4  &    302 &    135 &   2.0 &    542 &    255 &   3.1\\  \hline
$(53)$ &  0  &   1029 &    325 &   4.4 &   1989 &    645 &   5.6\\ 
&  1  &    597 &    181 &   3.7 &   1077 &    341 &   5.5\\ 
&  2  &    399 &    115 &   0.3 &    639 &    195 &   3.4\\ 
&  3  &    819 &    255 &   4.4 &   1539 &    495 &   5.6\\ 
&  4  &    714 &    220 &   4.8 &   1314 &    420 &   5.6\\  \hline
$(73)$ &  0  &   1387 &    407 &   4.9 &   2587 &    807 &   7.0\\ 
&  1  &    775 &    203 &   3.3 &   1255 &    363 &   6.5\\ 
&  2  &    595 &    143 &  * &    835 &    223 &   3.4\\ 
&  3  &   1189 &    341 &   4.9 &   2149 &    661 &   7.0\\ 
&  4  &   1090 &    308 &   4.9 &   1930 &    588 &   7.0\\  \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS, NFP and D for problem A4.}
\label{a4}
\begin{center}
\begin{tabular}{|r|l|rrr|rrr|}
\hline \hline
    &  & \multicolumn{3}{c|}{$h^{-1}$=2} & \multicolumn{3}{c|}{$h^{-1}$=4} \\
$(r,s)$ &  TP &  NFS   & NFP    &  D   &   NFS &  NFP &  D \\ \hline
$(42)$ & 0  &    446 &    207 &   5.9 &    846 &    407 &   6.7\\ 
&  1  &    224 &     96 &   3.5 &    384 &    176 &   4.4\\ 
&  2  &    158 &     63 &   1.0 &    238 &    103 &   1.8\\ 
&  3  &    374 &    171 &   5.9 &    694 &    331 &   6.7\\ 
&  4  &    302 &    135 &   5.8 &    542 &    255 &   6.7\\  \hline
$(53)$ & 0  &   1029 &    325 &   9.2 &   1989 &    645 &  10.5\\ 
&  1  &    597 &    181 &   6.8 &   1077 &    341 &   8.2\\ 
&  2  &    399 &    115 &   2.6 &    639 &    195 &   3.4\\ 
&  3  &    819 &    255 &   9.2 &   1539 &    495 &  10.5\\ 
&  4  &    714 &    220 &   9.0 &   1314 &    420 &  11.1\\  \hline
$(73)$ & 0  &   1387 &    407 &  10.4 &   2587 &    807 &  12.3\\ 
&  1  &    775 &    203 &   6.8 &   1255 &    363 &   8.2\\ 
&  2  &    595 &    143 &   2.5 &    835 &    223 &   3.2\\ 
&  3  &   1189 &    341 &  10.4 &   2149 &    661 &  12.3\\ 
&  4  &   1090 &    308 &  10.4 &   1930 &    588 &  12.3\\  \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[p]
\caption{Values of NFS and NFP for Fehlberg's problem.}
\label{t4}
\begin{tabular}{|l|l|rr|rr|rr|rr|rr|rr|} \hline \hline
\multicolumn{2}{|c|}{-log(TOL)}  & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{6} & \multicolumn{2}{c}{8}
& \multicolumn{2}{c}{10} & \multicolumn{2}{c}{12} & \multicolumn{2}{c|}{14}\\  \hline
Method & TP & NFS&  NFP & NFS & NFP& NFS & NFP & NFS & NFP & NFS & NFP &NFS & NFP\\ \hline
42 &0&654 &324 &1108 &549 &2008 &999 &3776 &1881 &7816 &3899&16728&8355\\
 &1&550 &272 &1054 &522 &2054 &1022 &4642 &2314 & 10152 &5067&21622&10802\\
 &2&582 &288 &1090 &540 &2228 &1109 &4466 &2226 & 10130 &5056&21612&10797\\
 &3&560 &277 &1056 &523 &1900 &945 &3722 &1854 &8398 &4190&20066&10024\\
 &4&538 &266 &978 &484 &1872 &931 &4232 &2109 &9406 &4694&19674&9828\\ \hline
53 &0&746 &246 &1094 &360 &1778 &588 &2762 &914 &4055 &1343&6201&2059\\
 &1&866 &286 &1142 &376 &1730 &572 &2483 &821 &4109 &1361&6216&2064\\
 &2&680 &224 &947 &311 &1637 &541 &2627 &869 &3785 &1253&6324&2100\\
 &3&665 &219 &965 &317 &1448 &478 &2348 &776 &3545 &1173&6024&2000\\
 &4&677 &223 &983 &323 &1532 &506 &2354 &778 &3740 &1238&6066&2014\\ \hline
73 &0&832 &261 &1135 &358 &1573 &504 &2407 &777 &3577 &1162&5053&1657\\
 &1&874 &275 &1357 &432 &1789 &576 &2434 &786 &3520 &1143&4972&1630\\
 &2&778 &243 &1153 &364 &1567 &502 &2353 &759 &3448 &1119&4936&1618\\
 &3&691 &214 &1069 &336 &1459 &466 &2059 &661 &2866 &925&4360&1426\\
 &4&745 &232 &1000 &313 &1381 &440 &1942 &622 &3076 &995&4396&1438\\ \hline
GS3&0&854 & 310 &1648 & 592 &2635 & 955 &4454 &1626 &7902 &2904&14905&5487\\
 &3&754 & 252 &1309 & 437 &2503 & 835 &4429 &1477 &8440 &2814&18100&6034\\
 &4&632 & 232 &1213 & 439 &2356 & 850 &4651 &1665 &9368 &3342&18178&6478\\ \hline
GS5&0&1021 & 221 &1325 & 285 &2081 & 445 &3181 & 677 &4719 &1003&6357&1357\\
 &3&771 & 155 &1151 & 231 &1646 & 330 &2406 & 482 &3856 & 772&5776&1156\\
 &4&762 & 174 &1152 & 252 &1692 & 368 &2423 & 519 &4051 & 863&5378&1150\\ \hline
RK78 & &559 &559&741 &741&1118 &1118&1716 &1716&2782 &2782&4615&4615 \\
\hline
\end{tabular}
\end{table}

\begin{table}[p]
\caption{Values of NFS and NFP for Rigid body problem.}
\label{t6}
\begin{tabular}{|l|l|rr|rr|rr|rr|rr|rr|} \hline \hline
\multicolumn{2}{|c|}{-log(TOL)}  & \multicolumn{2}{c}{4} & \multicolumn{2}{c}{6} & \multicolumn{2}{c}{8}
& \multicolumn{2}{c}{10} & \multicolumn{2}{c}{12}& \multicolumn{2}{c|}{14} \\  \hline
Method & TP & NFS&  NFP & NFS & NFP& NFS & NFP & NFS & NFP & NFS & NFP &NFS & NFP\\ \hline
42 &0&238 &112 &446 &213 &784 &379 &1466 &717 &2930 &1443&6010&2983\\
 &1&216 &101 &408 &194 &812 &393 &1780 &874 &3694 &1825&7656&3806\\
 &2&226 &106 &436 &208 &854 &414 &1734 &851 &3586 &1771&7736&3846\\
 &3&214 &100 &376 &178 &720 &347 &1454 &711 &3374 &1665&7166&3561\\
 &4&208 & 97 &398 &189 &810 &392 &1668 &818 &3254 &1605&7242&3599\\ \hline
53 &0&293 & 93 &502 &160 &759 &243 &1127 &363 &1609 &521&2562&836\\
 &1&341 &109 &469 &149 &714 &228 &995 &319 &1612 &522&2391&779\\
 &2&323 &103 &427 &135 &597 &189 &911 &291 &1372 &442&2652&866\\
 &3&260 & 82 &430 &136 &594 &188 &875 &279 &1453 &469&2082&676\\
 &4&284 & 90 &433 &137 &600 &190 &935 &299 &1420 &458&2199&715\\ \hline
73 &0&325 & 95 &523 &155 &715 &213 &1015 &307 &1441 &443&2050&640\\
 &1&430 &130 &616 &186 &832 &252 &1099 &335 &1516 &468&2065&645\\
 &2&361 &107 &589 &177 &802 &242 &1162 &356 &1435 &441&2044&638\\
 &3&313 & 91 &481 &141 &694 &206 &892 &266 &1210 &366&1708&526\\
 &4&313 & 91 &472 &138 &655 &193 &976 &294 &1312 &400&1855&575\\ \hline
GS3&0&300 & 112 &621 & 225 &989 & 361 &1655 & 603 &2949 &1079&5349&1965\\
 &3&280 &94 &442 & 148 &745 & 249 &1591 & 531 &3217 &1073&6361&2121\\
 &4&250 &98 &524 & 192 &889 & 321 &1816 & 648 &3193 &1141&5761&2065\\ \hline
GS5&0&406 &94 &573 & 129 &913 & 205 &1369 & 301 &1830 & 398&2573&553\\
 &3&361 &73 &596 & 120 &696 & 140 &1126 & 226 &1521 & 305&2026&406\\
 &4&355 &87 &675 & 159 &922 & 210 &1137 & 253 &1609 & 349&2195&471\\ \hline
RK78 & &182 &182&234 &234&364 &364&611 &611&949 &949&1508&1508 \\
\hline
\end{tabular}
\end{table}

\begin{table}[p]
\caption{Values of NFS and NFP for Orbit problem.}
\label{t7}
\begin{tabular}{|l|l|rr|rr|rr|rr|rr|rr|} \hline \hline
\multicolumn{2}{|c|}{-log(TOL)}& \multicolumn{2}{c}{4} & \multicolumn{2}{c}{6} & \multicolumn{2}{c}{8}
& \multicolumn{2}{c}{10} & \multicolumn{2}{c}{12}& \multicolumn{2}{c|}{14} \\\hline
Method & TP & NFS&NFP & NFS & NFP& NFS & NFP & NFS & NFP & NFS & NFP &NFS & NFP\\ \hline
42 &0& 672 & 326 &1024 & 499 &1850 & 903 &3988 &1972 &8594 &4269&18448&9190\\
 &1& 586 & 283 &1100 & 537 &2228 &1092 &5210 &2583 & 11400 &5672&24106&12019\\
 &2& 558 & 269 &1146 & 560 &2508 &1232 &5062 &2509 & 10914 &5429&23562&11747\\
 &3& 630 & 305 &1138 & 556 &2236 &1096 &4410 &2183 &9428 &4686&23750&11841\\
 &4& 520 & 250 &1042 & 508 &2058 &1007 &4710 &2333 & 10040 &4992&22586&11259\\ \hline
53&0&679&219&1419&463&1828&594&2593&841&4300&1410&7157&2357\\
&1&775&251&1050&340&1663&539&2689&873&4246&1392&7355&2423\\
&2&637&205&1005&325&1591&515&2845&925&4435&1455&7562&2492\\
&3&622&200&1008&326&1531&495&2491&807&3874&1268&7241&2385\\
&4&652&210&1035&335&1537&497&2527&819&4000&1310&7037&2317\\\hline
73&0&724&222&1519&481&1774&554&2656&836&3691&1175&5305&1703\\
&1&871&271&1240&388&1801&563&2563&805&3841&1225&5647&1817\\
&2&757&233&1234&386&1672&520&2668&840&3943&1259&5620&1808\\
&3&685&209&1015&313&1552&480&2254&702&3499&1111&4849&1551\\
&4&652&198&1027&317&1570&486&2323&725&3385&1073&5011&1605\\\hline
GS3&0&848&314&1655&601&2375&875&4579&1687&8815&3247&17003&6263\\
&3&640&214&1240&414&2434&812&5176&1726&10807&3603&19549&6517\\
&4&696&260&1179&431&2636&948&5406&1932&10823&3857&21053&7499\\\hline
GS5&0&1081&245&1707&379&2505&541&3588&764&4881&1041&6120&1320\\
&3&851&171&1316&264&2066&414&2791&559&3941&789&5951&1191\\
&4&941&225&1528&352&2110&462&2876&620&4168&896&6308&1348\\\hline
RK78&&468&468&741&741&1105&1105&1638&1638&2587&2587&4511&4511\\
\hline
\end{tabular}
\end{table}



\end{document}