\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}