\documentclass[twoside]{article}
\usepackage{amssymb} % used for R in Real numbers
\usepackage{epsf} % used for two PostScript figures
\pagestyle{myheadings}
\markboth{\hfil Fifth-order Runge-Kutta \hfil}%
{\hfil David Goeken \& Olin Johnson \hfil}
\begin{document}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc 15th Annual Conference of Applied Mathematics, Univ. of Central Oklahoma,}
\hfill\break
Electronic Journal of Differential Equations, Conference~02, 1999, pp. 1--9.
\newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu (login: ftp)}
\vspace{\bigskipamount} \\
Fifth-order Runge-Kutta with higher order derivative approximations
\thanks{ {\em 1991 Mathematics Subject Classifications:} 65L06.
\hfil\break\indent
{\em Key words and phrases:} multistep Runge-Kutta, third-order method,
\hfil\break\indent
fourth-order method, fifth-order method, higher order derivatives.
\hfil\break\indent
\copyright 1999 Southwest Texas State University and University of
North Texas. \hfil\break\indent
Published November 23, 1999.} }
\date{}
\author{David Goeken \& Olin Johnson}
\maketitle
\begin{abstract}
Given $y'=f(y)$, standard Runge-Kutta
methods perform multiple evaluations of $f(y)$ in each integration
sub-interval as required for a given accuracy. Evaluations of
$y''=f_{y}f$ or higher derivatives are not considered due
to the assumption that the calculations involved in these functions
exceed those of $f$. However, $y''$ can be approximated to
sufficient accuracy from past and current evaluations of $f$ to
achieve a higher order of accuracy than is available through current functional
evaluations alone.
In July of 1998 at the ANODE (Auckland Numerical Ordinary Differential
Equations) Workshop, we introduced
a new class of Runge-Kutta methods based on this observation (Goeken~1999).
We presented a third-order method which requires
only two evaluations of $f$ and a fourth-order method which requires three.
This paper reviews these two methods and gives the general solution
to the equations generated by the fifth-order methods of this new class.
Interestingly, these fifth-order methods require only four functional
evaluations per step whereas standard Runge-Kutta methods require six.
\end{abstract}
\section{Third-order method}
We consider initial value problems expressed in autonomous form.
Starting with the non-autonomous form,
we assume that $f(x,y)$ is a continuous function with
domain $D$ in ${\mathbb R}^{n+1}$ where $x\in{\mathbb R}$, $y\in{\mathbb R}^n$ and $(x,y)\in D$
We assume that
$$
\|f(x,y_1)-f(x,y_2)\|_2\le L\|y_1-y_2\|_2
$$
for all $(x,y_1)$, $(x,y_2)\in D$; thus the problem
$$ \displaylines{
y'=f(x,y) \cr
y(x_{0})=y_{0}$ with $(x_{0},y_{0})\in D\cr}
$$
has a unique solution.
In autonomous form, $y$ and $f$ have $n+1$ components with $y_{n+1}=x$
and $f_{n+1}(y)=1$. The initial value problem is then written
$$\displaylines{
y'=f(y)\cr
y(x_{0})=y_{0}$ where $(y_{0})_{n+1}=x_{0}\,.\cr}
$$
Most efforts to increase the order of the Runge-Kutta methods have been
accomplished by increasing the number of Taylor's series terms used and
thus the number of functional evaluations required (Butcher~1987)
(Gear~1971). The use of higher order derivative terms has been proposed for
stiff problems (Rosenbrock~1963) (Enright~1974). Our method
adds higher order derivative terms to the Runge-Kutta $k_{i}$ terms ($i>1$)
to achieve a higher order of accuracy. For example, our new third-order
method, GJ3, for autonomous systems, lets
$$y_{n+1}=y_{n}+b_1k_1+b_2k_2$$
and $k_1=hf(y_{n})$.
However, we introduce additional terms by assigning
$$k_2=hf(y_{n}+a_{21}k_1+ha_{22}f_{y}(y_{n})k_1)$$
Using Taylor's series expansion techniques, the above is
uniquely satisfied to $O(h^3)$ as follows
\begin{eqnarray*}
k_1&=&hf(y_{n}) \\
k_2&=&hf(y_{n}+\frac{2}{3}k_1+\frac{2}{9}hf_{y}(y_{n})k_1) \\
y_{n+1}&=&y_{n}+{1\over 4}k_1+{3\over 4}k_2
\end{eqnarray*}
\section{Utilizing $f_{y}$}
The previous section developed a two-stage, third-order method;
however, it introduced a term with $f_{y}$.
The result is the addition of a higher derivative term to the
standard Runge-Kutta method.
The following describes three methods to utilize $f_{y}$.
\paragraph{Method 1:} If one knows or can generate $f_{y}$, and
if the evaluation of $f_{y}$ is cheaper than the evaluation of $f$,
then savings can be realized. For example, with a linear system of
equations, $y'=Ay$, $f_{y}$ is known and constant.
\paragraph{Method 2:} Since $y''=f'=f_{y}f$ for autonomous equations,
and since $k_1=hf$, $k_2$ can be replaced with
\begin{eqnarray*}
k_2&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}hf_{y}k_1) \\
&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}hf_{y}hf) \\
&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}h^2f_{y}f) \\
&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}h^2f')
\end{eqnarray*}
or
$$k_2=hf(y_{n}+{2\over 3}k_1+{2\over 9}h^2y'')\,.$$
Again, savings can be realized if one can formulate $y''$ (or $f'$)
and if it is cheaper to evaluate than $f$.
\paragraph{Method 3:} Building onto Method 2, one can approximate $y''$ (or $f'$)
by using the current and previous evaluations of $f$. For our third-order
method, this approximation must be of $O(h)$.
Since an $O(h)$ approximation of $f'$ is given by $f'=(f_{n}-f_{n-1})/h$,
one can approximate $k_2$ as follows
\begin{eqnarray*}
k_2&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}h^2f') \\
&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}h^2(f_{n}-f_{n-1})/h) \\
&=&hf(y_{n}+{2\over 3}k_1+{2\over 9}h(f_{n}-f_{n-1}))
\end{eqnarray*}
Since $f_{n}$ is calculated in the current step in the evaluation of
$k_1$, one only has to store the previous value, $f_{n-1}$.
In effect, the use of previous values for the approximation has
created a multistep Runge-Kutta method.
\section{Fourth-order method}
Similarly, our fourth-order method, GJ4, for autonomous systems, lets
$$y_{n+1}=y_{n}+b_1k_1+b_2k_2+b_3k_3$$
and
\begin{eqnarray*}
k_1&=&hf(y_{n}) \\
k_2&=&hf(y_{n}+a_{21}k_1+ha_{22}f_{y}(y_{n})k_1) \\
k_3&=&hf(y_{n}+a_{31}k_1+a_{32}k_2+ha_{33}f_{y}(y_{n})k_1
+ha_{34}f_{y}(y_{n})k_2)
\end{eqnarray*}
The Taylor's series expansion of these higher order methods
is tedious and error prone. We used modern symbolic math packages to
expand and then to solve the resulting systems of nonlinear equations
that were generated.
In this work, we used the symbolic math packages Reduce~(Reduce 1999),
PARI/GP~(PARI/GP 1999), and Octave~(Eaton 1997).
PARI/GP was used to generate the Taylor's series expansion of the
above, resulting in the following system of equations
$$\displaylines{
b_1+b_2+b_3=1 \cr
b_2a_{21}+b_3[a_{31}+a_{32}]=1/2 \cr
b_2a_{21}^2+b_3[a_{31}+a_{32}]^2=1/3 \cr
b_2a_{21}^3+b_3[a_{31}+a_{32}]^3=1/4 \cr
b_2a_{22}+b_3[a_{21}a_{32}+a_{33}+a_{34}]=1/6 \cr
b_3[a_{21}a_{34}+a_{22}a_{32}]=1/24 \cr
b_2a_{21}a_{22}+b_3[{a_{21}a_{32}(\frac{1}{2}a_{21}
+a_{31}+a_{32})+(a_{31}+a_{32})(a_{33}+a_{34})}]=1/6 \cr }
$$
However, in order to utilize Methods 2 and 3 of Section 2, we must restrict
the solution with $a_{34}=0$. The general solution to the above system
of equations (with $a_{34}=0$) has been found with Reduce
and example solutions are shown in Table 1.
\begin{table}[ht]
\caption{Example of fourth-order autonomous solutions}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
$b_1$&$b_2$&$b_3$&$a_{21}$&$a_{22}$&$a_{31}$&$a_{32}$&$a_{33}$\\
\hline
$1/6$&$1/6$&$2/3$&1&\hspace{1em}$1/2$&\hspace{1em}$3/8$&$\hspace{1ex}1/8$&\hspace{1em}0\\
$1/6$&$2/3$&$1/6$&$1/2$&\hspace{1em}$1/8$&\hspace{1em}-{1}&\hspace{1em}2&$-{1/2}$\\
$1/6$&$2/3$&$1/6$&$1/2$&$-{1/8}$&\hspace{1em}3&\hspace{1ex}-{2}&\hspace{1em}$5/2$\\
$1/10$&$1/2$&$2/5$&$1/3$&$\hspace{1em}1/18$&$-{25/24}$&$15/8$&$\hspace{1em}-{5/18}$\\
$1/10$&$1/2$&$2/5$&$1/3$&$-{1/6}$&$\hspace{1em}35/24$&$-{5/8}$&$\hspace{1em}5/6$\\
\hline
\end{tabular}
\end{center}
\end{table}
\section{Fifth-order method }
In July of 1998, the authors presented (Goeken~1999) this
numerical integration
technique at a meeting attended by Dr.\ John~Butcher.
Using his tree-based approach (Butcher~1987),
Dr.\ Butcher suggested a fifth-order
method. Since the meeting, his technique has been
verified using Taylor's series expansion techniques
to determine the general solution for our fifth-order methods.
Our fifth-order method, GJ5, for autonomous systems, lets
$$y_{n+1}=y_{n}+b_1k_1+b_2k_2+b_3k_3+b_4k_4$$
and
\begin{eqnarray*}
k_1&=&hf(y_{n}) \\
k_2&=&hf(y_{n}+a_{21}k_1+ha_{22}f_{y}(y_{n})k_1) \\
k_3&=&hf(y_{n}+a_{31}k_1+a_{32}k_2+ha_{33}f_{y}(y_{n})k_1) \\
k_4&=&hf(y_{n}+a_{41}k_1+a_{42}k_2+a_{43}k_3+ha_{44}f_{y}(y_{n})k_1)
\end{eqnarray*}
PARI/GP was used to generate the Taylor's series expansion of the
above, resulting in the following system of equations
$$ \displaylines{
b_1+b_2+b_3+b_4=1 \cr
b_2a_{21}+b_3[a_{31}+a_{32}]+b_4[a_{41}+a_{42}+a_{43}]=1/2 \cr
b_2a_{21}^2+b_3[a_{31}+a_{32}]^2+b_4[a_{41}+a_{42}+a_{43}]^2=1/3 \cr
b_2a_{22}+b_3[a_{21}a_{32}+a_{33}]+b_4[a_{21}a_{42}+a_{43}(a_{31}+a_{32})
+a_{44}]=1/6 \cr
b_2a_{21}^3+b_3[a_{31}+a_{32}]^3+b_4[a_{41}+a_{42}+a_{43}]^3=1/4 \cr
b_2a_{21}a_{22}+
b_3[{1\over 2}a_{21}^2a_{32}+(a_{31}+a_{32})(a_{21}a_{32}+a_{33})]+
{1\over 2}b_4[a_{21}^2a_{42} \cr
+a_{43}(a_{31}+a_{32})^2+
2(a_{41}+a_{42}+a_{43})(a_{21}a_{42}+(a_{31}+a_{32})a_{43}+a_{44})]=1/6 \cr
b_3a_{22}a_{32}+b_4[a_{21}a_{32}a_{43}+a_{22}a_{42}+a_{33}a_{43}]=1/24 \cr
b_2a_{21}^{4}+b_3[a_{31}+a_{32}]^{4}+b_4[a_{41}+a_{42}+a_{43}]^{4}=1/5 \cr
3b_2a_{21}^2a_{22}+
b_3[a_{21}^3a_{32}+3(a_{31}+a_{32})^2(a_{21}a_{32}+a_{33})]+
b_4[a_{21}^3a_{42} \cr
+(a_{31}+a_{32})^3a_{43}+
3(a_{41}+a_{42}+a_{43})^2(a_{21}a_{42}+(a_{31}+a_{32})a_{43}+a_{44})]=7/20 \cr
b_3a_{21}^2a_{32}(a_{31}+a_{32})+
b_4[(a_{41}+a_{42}+a_{43})(a_{21}^2a_{42}+(a_{31}+a_{32})^2a_{43})]=1/15 \cr
{1\over 2}b_2a_{22}^2+
b_3[a_{21}a_{32}({1\over 2}a_{21}a_{32}+a_{22}+a_{33})+
a_{22}a_{32}(a_{31}+a_{32})+{1\over 2}a_{33}^2] \cr
+b_4[{1\over 2}a_{21}^2(a_{32}a_{43}+a_{42}^2)+
(a_{31}+a_{32})(a_{21}(a_{32}a_{43}+a_{42}a_{43})+a_{43}(a_{33}+a_{44})\cr
+{1\over 2}(a_{31}+a_{32})a_{43}^2)+
a_{21}a_{42}(a_{22}+a_{44})+(a_{21}a_{32}a_{43}+a_{22}a_{42}\cr
+a_{33}a_{43})(a_{41}+a_{42}+a_{43})+
{1\over 2}a_{44}^2]=11/120 \cr
b_4a_{22}a_{32}a_{43}=1/120 \cr}
$$
The solution presented by Dr.\ Butcher and verified using the above system
of equations is
\begin{eqnarray*}
k_1&=&hf(y_{n}) \\
k_2&=&hf(y_{n}+{1\over 3}k_1+{1\over 18}hf_{y}k_1) \\
k_3&=&hf(y_{n}-{152\over 125}k_1+{252\over 125}k_2-{44\over 125}hf_{y}k_1) \\
k_4&=&hf(y_{n}+{19\over 2}k_1-{72\over 7}k_2+{25\over 14}k_3+{5\over 2}hf
_{y}k_1) \\
y_{n+1}&=&y_{n}+{5\over 48}k_1+{27\over 56}k_2+{125\over 336}k_3+{1\over 24
}k_4 \\
\end{eqnarray*}
Additional solutions to the above nonlinear system of equations have
been found using Octave. Three additional solutions are shown in Table 2.
\begin{table}[ht]
\caption{Example of fifth-order autonomous solutions}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
$b_1$&1/24&5/54&1/14 \\
\hline
$b_2$&125/336&250/567&32/81 \\
\hline
$b_3$&27/56&32/81&250/567 \\
\hline
$b_4$&5/48&1/14&5/54 \\
\hline
$a_{21}$&1/5&3/10&1/4 \\
\hline
$a_{22}$&1/50&9/200&1/32 \\
\hline
$a_{31}$&-52/27&-9/8&-329/250 \\
\hline
$a_{32}$&70/27&15/8&252/125 \\
\hline
$a_{33}$&-8/27&-9/32&-259/1000 \\
\hline
$a_{41}$&43/5&17/3&209/35 \\
\hline
$a_{42}$&-64/7&-490/81&-32/5 \\
\hline
$a_{43}$&54/35&112/81&10/7 \\
\hline
$a_{44}$&13/10&23/18&11/10 \\
\hline
\end{tabular}
\end{center}
\end{table}
\section{Numerical results}
To demonstrate that the new methods are of the order claimed, several
equations have been solved using the new third-, fourth-, and fifth-order
method on scalar autonomous equations, systems of autonomous equations,
and scalar non-autonomous equations.
Our previous paper (Goeken~1999) demonstrated the third- and fourth-order
methods utilizing $f_{y}$ and the approximation to $f'$.
Here we will concentrate on our new fifth-order method.
For scalar autonomous examples, we use the equations shown in Table 3,
with initial condition $y(0)=1$. These equations were solved using
a standard fifth-order Runge-Kutta method along with our fifth-order
method, GJ5, using $f_{y}$ directly.
Relative error was plotted against the step size
and is shown in Figures 1 and 2. Results are comparable to standard
fifth-order Runge-Kutta solution, thus demonstrating our claim.
The new method requires four functional evaluations of $f$ and
one of $f_{y}$ per step or four functional evaluations of $f$ and
three historical values of $f$;
whereas, the standard fifth-order Runge-Kutta method
requires six functional evaluations of $f$.
\begin{table}
\caption{Test problems}
\begin{center}
\begin{tabular}{|l|l|c|}
\hline
Function & Solution &$y(0)$\\
\hline
$y'=-y$ & $y=e^{-t}$ & 1\\
$y'={y\over{4}}-{y^2\over{80}}$ &$y={20\over{1+19e^{-t\over{4}}}}$ & 1\\
\hline
\end{tabular}
\end{center}
\end{table}
\section{Conclusions}
New third-, fourth-, and fifth-order
numerical integration techniques inspired by the Runge-Kutta method
have been presented.
The new methods exploit the use of higher order derivatives, specificly
$f_{y}$.
In particular, a technique utilizing an approximation to $y''$ has been
presented resulting in a multistep Runge-Kutta method.
Table 4 compares the
computational effort required for standard Runge-Kutta methods
with our methods.
Table 4 shows the cases where the proposed methods
are more efficient than the standard Runge-Kutta methods. Specifically,
the proposed methods are more efficient for cases where
\begin{itemize}
\item $f_{y}$ or $y''$ is cheaper to evaluate than $f$,
\item the use of historical values of $f$ is cheaper then evaluating $f$, and
\item for the fifth-order case, the number of total functional evaluations
can be reduced from 6 to 4 when using an approximation of $f'$.
\end{itemize}
\begin{table}[ht]
\caption{Number of evaluations comparison}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
&Standard&
\multicolumn{2}{c|}{Proposed Method}&
\multicolumn{2}{c|}{Proposed Method} \\
&Runge-Kutta&
\multicolumn{2}{c|}{Exact $f_y$}&
\multicolumn{2}{c|}{Approximating $f'$} \\
\cline{2-6}
&Num. $f$&Num. $f$&Num. $f_y$&Num. $f$&Num. $f_{n-i}$\\
Order&Evals.&Evals&Evals.&Evals.&Values\\
\hline
3&3&2&1&2&1\\
4&4&3&1&3&2\\
5&6&4&1&4&3\\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{figure*}
\epsfysize=72mm
\leavevmode
\epsfbox {fig1.ps}
\end{figure*}
\begin{figure*}
\epsfysize=72mm
\leavevmode
\epsfbox {fig2.ps}
\end{figure*}
\begin{thebibliography}{0}
\bibitem{jb} Butcher, J. C. (1987) {\it The Numerical Analysis of Ordinary Differential
Equations Runge-Kutta and General Linear Methods}, John Wiley \& Sons
Ltd., New York
\bibitem{je} Eaton, J. W. (1997) {\it GNU Octave A high-level interactive language
for numerical computations}, \hfill\newline
http://www.we.fh-osnabrueck.de/labsim/octave/manual/octave.html
\bibitem{we} Enright, W. H. (1974) {\it Second Derivative Multistep Methods for Stiff
Ordinary Differential Equations}, SIAM J. Numer. Anal., {\bf 11}, 321-331
\bibitem{cg} Gear, C. W. (1971) {\it Numerical Initial Value Problems in Ordinary
Differential Equations}, Prentice-Hall, Englewood Cliffs, New Jersey
\bibitem{dg} Goeken, D. and Johnson, O. (1999) {\it Runge-Kutta with Higher
Order Derivative Approximations}, submitted to Applied Numerical Mathematics
\bibitem{pari} PARI/GP (1999) \hfill\newline
http://hasse.mathematik.tu-muenchen.de/ntsw/pari/Welcome
\bibitem{red} Reduce (1999) \hfill\newline
http://www.sub.uni-goettingen.de/ssgfi/math/infodata/000838.html
\bibitem{hr} Rosenbrock, H. H. (1963) {\it Some General Implicit Processes for the
Numerical Solution of Differential Equations}, Comp. J., {\bf 5}, 329-330
\end{thebibliography}
\bigbreak
\noindent{\sc David Goeken} \\
Department of Computer Science \\
The University of Houston \\
Houston, TX 77204-3475, USA \\
e-mail: dgoeken@cs.uh.edu \\
Now with the LinCom Corporation, Houston, Texas \medskip
\noindent {\sc Olin Johnson} \\
Department of Computer Science \\
The University of Houston \\
Houston, TX 77204-3475, USA \\
e-mail: johnson@cs.uh.edu
\end{document}