# Accelerated Runge-Kutta Methods

Authors: Firdaus  E. Udwadia; Artin Farahani
Journal: Discrete Dynamics in Nature and Society
 (2008)

---

## Abstract

Standard Runge-Kutta methods are explicit, one-step, and generally constant step-size numerical
integrators for the solution of initial value problems. Such integration schemes of orders 3, 4, and 5
require 3, 4, and 6 function evaluations per time step of integration, respectively. In this paper, we
propose a set of simple, explicit, and constant step-size Accerelated-Runge-Kutta methods that are two-step in nature. For orders 3, 4, and 5, they require only 2, 3, and 5 function evaluations per time step,
respectively. Therefore, they are more computationally efficient at achieving the same order of local
accuracy. We present here the derivation and optimization of these accelerated integration methods.
We include the proof of convergence and stability under certain conditions as well as stability regions
for finite step sizes. Several numerical examples are provided to illustrate the accuracy, stability, and
efficiency of the proposed methods in comparison with standard Runge-Kutta methods.


---

## Body

## 1. Introduction

Most of the
ordinary differential equations (ODEs) that model systems in nature and society
are nonlinear, and as a result are often extremely difficult, or sometimes
impossible, to solve analytically with presently available mathematical methods.
Therefore, numerical methods are often very useful for understanding the
behavior of their solutions. A very important class of ODEs is the initial
value problem (IVP):(1.1)dy˜(t)dt=f˜(t,y˜(t)),t0≤t≤tN,y˜(t0)=y˜0,where y˜, y˜0,
and f˜ are (m−1)-dimensional vectors, and t denotes the independent variable, time.A constant step numerical integrator approximates the
solutiony˜(t1) of the IVP defined by (1.1) at point t1=t0+h by y˜1.
This is referred to as “taking a step.” Similarly, at step n+1,
the numerical integrator approximates the solution y˜(tn+1) at point tn+1=tn+h by y˜n+1.
By taking steps successively, the approximate solutions y˜2,…,y˜N at points t2,…,tN are generated. The accuracy of a numerical
integrator is primarily determined by the order of the local error it
generates in taking a step. The local error at the end of step n+1 is the difference u˜(tn+1)−y˜n+1,
where u˜(t) is the local solution that satisfies the local
IVP:(1.2)du˜(t)dt=f˜(t,u˜(t)),u˜(tn)=y˜n.Standard Runge-Kutta (RK) methods are a
class ofone-step numerical integrators. That is, the approximate
solution y˜n+1 is calculated using y˜n and the function f˜.
An RK method that has a local error of O(hp+1) is said to be of order p and is denoted by RKp.
It has been established that RK3, RK4, and RK5 require 3, 4, and 6 function
evaluations per time step of integration, respectively [1–3].In this paper, we propose a new and simple set of
numerical integrators named the accelerated-Runge-Kutta (ARK) methods. We
derive these methods for the autonomous version of (1.1) given by(1.3)dy(t)dt=f(y(t)),t0≤t≤tN,y(t0)=y0,where y, y0,
and f are m-dimensional vectors, and t denotes the independent variable, time. We
assume that f is a (sufficiently smooth) Lipschitz
continuous function, and hence the solution of the IVP given by (1.3) is unique.Accelerated Runge-Kutta methods are inspired by
standard Runge-Kutta methods but aretwo-step in nature. That is, the
approximate solution yn+1 is calculated using yn, yn−1,
and the function f.
ARKp denotes an ARK method whose local error is of O(hp+1).
We will see in the next section that ARK3, ARK4, and ARK5 require 2, 3, and 5
function evaluations per time step of integration, respectively. Since function
evaluations are often the most computationally expensive part of numerically
solving differential equations (see Section 4 for further details), ARK methods
are expected to be more computationally efficient than standard RK methods.Various authors have attempted to increase the
efficiency of RK methods. As a result, numerous variations of two-step explicit
RK methods exist today. Byrne and Lambert [4] proposed 3rd- and 4th-order two-step Pseudo RK
methods. Byrne [5]
later proposed a set of methods that minimize a conservative bound on the
error, and Gruttke [6]
proposed a 5th-order Pseudo RK method. Costabile [7] introduced the Pseudo RK
methods of II species which involve many parameters to improve local accuracy.
Jackiewicz et al. [8]
considered the two-step Runge-Kutta (TSRK) methods up to order 4. Jackiewicz
and Tracogna [9]
studied a more general class of these TSRK mothods and derived order conditions
up to order 5 using Albrecht's method. Recently, Wu [10] has proposed a set of
two-step methods which make use of the derivative of the differential equation.The accelerated-Runge-Kutta (ARK) methods presented
here are along the lines first proposed in [11] and are simple and straightforward. They could be
considered special cases of the more general TSRK methods. The simplicity of
the ARK methods' construction not only reduces computational overhead cost, but
also leads to a simpler and more elegant set of order equations that can be
directly analyzed more precisely. The ARK methods also benefit from a more
precise and effective optimization technique which transforms them into
higher-order methods for some problems. In addition, we have presented here a
complete analysis of all aspects of the ARK methods, including their
motivation, derivation, accuracy, speedup, and stability.In Section2, we explain the central idea of ARK
methods followed by a description of how to derive and optimize ARK3, ARK4, ARK4-4, (a
more accurate ARK4), and ARK5 in Sections 2.1,
2.2, 2.3, and
2.4, respectively.
In Section 3, we use seven standard initial value problems to compare the
accuracy of ARK and RK methods. In Section 4, we present the computational
speedup of ARK methods over the RK methods using some of the same standard
problems. The hallmark of our ARK method is its simplicity. Tracogna and
Welfert [12], based on methods developed in [9], use a more sophisticated approach
using B-series to develop a 2-step TSRK algorithm for numerical integration. In
Section 5, we compare the numerical performance of our method with that of [12]
using the above-mentioned standard initial value problems.
Section 6.1 deals
with stability and convergence of ARK methods at small step sizes, and
Section 6.2
compares the stability region of ARK and RK methods at finite step sizes. We
end with our conclusions and findings, and include in the appendices, the
description of the standard initial value problems used, and the Maple code
that derives the ARK3 method.


## 2. Central Idea of Accelerated Runge-Kutta Methods

The idea behind the ARK methods proposed herein is perhaps
best conveyed by looking at RK3 and ARK3. An example of the RK3 method solving
the IVP given by (1.1) has the form(2.1)y˜n+1=y˜n+16k˜1+23k˜2+16k˜3for0≤n≤N−1,
where(2.2)k˜1=hf˜(tn,y˜n),k˜2=hf˜(tn+12h,y˜n+12k˜1),k˜3=hf˜(tn+h,y˜n−k˜1+2k˜2).The approximate solution y˜n+1 at tn+1 is calculated based on the approximate
solution y˜n at tn along with 3 function evaluations k˜1, k˜2,
and k˜3 in [tn,tn+1].
Hence, RK3 is a one-step method with a computational cost of 3 function
evaluations per time step.An example of theaccelerated-RK3 (ARK3) method
solving the IVP given by (1.3) has the form(2.3)yn+1=yn+12k1+12k−1+(k2−k−2)for1≤n≤N−1,where(2.4)k1=hf(yn),k−1=hf(yn−1),k2=hf(yn+512k1),k−2=hf(yn−1+512k−1).The approximate solution yn+1 at time tn+1 is calculated based on the approximate
solutions yn and yn−1 at times tn and tn−1, respectively, along with 2 function
evaluations k1 and k2 in [tn,tn+1] and 2 reused function evaluations k−1 and k−2 in [tn−1,tn].
In each step of integration from point tn to point tn+1,
only k1 and k2 are calculated, while k−1 and k−2 are reused from the previous step's k1 and k2.
Figure 1 illustrates this idea. The reuse of previously calculated data reduces
the number of function evaluations at each time step of integration from 3 for
RK3 to 2 for ARK3, making ARK3 more computationally efficient than RK3. This
increased efficiency is shared by ARK4 and ARK5 when compared to RK4 and RK5,
respectively, and is the central idea behind the ARK methods.Figure 1
The central idea of the ARK3 method is
illustrated here. Approximating the solutionyn+1 at tn+1 requires 2 function evaluations k1 and k2 (shown in red) and two reused values k−1 and k−2 (shown in blue) from the previous steps k1 and k2.
This process repeats for the next step (approximating the solution at tn+2).It is important to note that at the first step, there
is no previous step. Therefore, ARK methods cannot be self-starting. A one-step
method must supply the approximate solutiony1 at the end of the first step t1.
The one-step method must be of sufficient order to ensure that the difference y1−y(t1) is of order p or higher. In this case, the ARKp method will be convergent of order p (see Section 6.1, Theorem 6.4). For
example, the
ARK3 method can be started by the RK3 or RK4 methods. This extra computation
occurs only in the first step of integration; for all the subsequent steps, ARK
methods bear less computational cost than their RK counterparts.The general form of the ARK methods presented in this
paper is(2.5)yn+1=c0yn−c−0yn−1+c1k1−c−1k−1+∑i=2vci(ki−k−i)for1≤n≤N−1,where(2.6)k1=hf(yn),k−1=hf(yn−1),k2=hf(yn+a1k1),k−2=hf(yn−1+a1k−1),k3=hf(yn+a2k2),k−3=hf(yn−1+a2k−2),k4=hf(yn+a3k3),k−4=hf(yn−1+a3k−3),k5=hf(yn+a4k4),k−5=hf(yn−1+a4k−4).Here, v denotes the number of function evaluations
performed at each time step and increases with the order of local accuracy of
the ARK method. In each step of integration, only k1,k2,… are evaluated, while k−1,k−2,… are reused from the previous step.To determine the appropriate values of the parametersc's and a's, a three-part process is performed. (The
first and second parts of this process are shown in detail for ARK3 in
[11].)  First, the ARK expression (2.5) is expanded using the Taylor
series expansion. Second, after some algebraic manipulation, this expansion is
equated to the local solution u(tn+1) at tn+1 given by the Taylor series expansion (see
(1.2))(2.7)u(tn+1)=u(tn+h)=u(tn)+hu′(tn)+h22!u′′(tn)+h33!u′′′(tn)+⋯=yn+hyn′+h22!yn′′+h33!yn′′′+⋯.This results in the system of
nonlinear algebraic order equations
(2.8a)–(2.13a). Third, we
attempt to solve as many order equations as possible because the highest power
of h for which all of the order equations are
satisfied is the order of the resulting ARK method. The above process requires
a great deal of algebraic and numeric calculations which were mainly performed
using Maple (see Appendix B).In the following, we
will show what order of accuracy an ARK method can achieve with a given number
of function evaluationsv.
For v=1,
a second-order method ARK2 is possible. It is, however, omitted here because
such a low order of accuracy is of limited use in accurately solving
differential equations; we will look at cases v=2,…,5.


### 2.1. ARK3

Forv=2,
from (2.5), we have the 6 parameters c0,c−0,c1,c−1,c2, and a1.
The derived order equations up to O(h4) are
(2.8a)O(h0)c0−c−0=1,(2.8b)O(h1)c−0+c1−c−1=1,(2.8c)O(h2)−12c−0+c−1+c2=12,(2.8d)O(h3)−112c−0+c2a1=512,(2.8e)O(h4)c2a12=13,(2.8f)O(h4)0=16.
Clearly, (2.8f) cannot be satisfied, so it is impossible to achieve a
4th-order method with 2 function evaluations per time step. It is, however,
possible to achieve a 3rd-order method ARK3 with 2 function evaluations per
time step by satisfying (2.8a)–(2.8d). We solve these equations in terms of c0,c−0,c1,
and c−1 and let the remaining two parameters c2 and a1 become free parameters.We use sets of free parameters to optimize the
accuracy of an ARK method by minimizing its local error. The local error is represented
by higher order equations (order 4 and above in case of ARK3). The closer these
order equations are to being satisfied, the smaller the local error is.
Therefore, we try to minimize the difference between the LHS and RHS of these
order equations. This can be done in several different ways. One popular
technique is to minimize a general bound on the error [5]. A second technique is to
minimize a norm of the error [13]. A third technique is to satisfy as many individual
order equations as possible (hinted to by [4]). In general, there is no way of knowing a priori
which technique yields the most accurate ARK method since the error ultimately
depends on the problem at hand. We have chosen the last technique for two
reasons. One, it is straightforward to implement exactly. Two, for some classes
of problems, the differential term multiplying the unsolved higher order
equations may vanish, which could effectively increases the order of the
method. For instance, the differential term multiplying (2.8f) may vanish, and
provided (2.8e) is already satisfied, the resulting ARK method (with 2 function
evaluations) will become a 4th-order method.Using the two free parametersc2 and a1,
we can solve two higher-order equations: (2.8e) and any one of the two 5th-order
equations,
(2.9a)O(h5)1120c−0+c2a13=31120,(2.9b)O(h5)c−0=31.
However, we cannot choose (2.9b) because to prove stability, we will need −1≤c−0<1 (see condition (6.16)). Therefore, we choose
(2.8e) and (2.9a), which leads to Set 2 of Table 1. We also present a set in which c−0=0 because doing so can increase the stability
(see Section 6) and speed (see Section 4) of an ARK method. Setting c−0=0 and solving (2.8a)–(2.8e) leads to Set 3 of Table 1. Set 1 of
Table 1 is the only set of parameters here that is not found by
solving all possible higher order equations. It is found by setting c−0=0,
solving (2.8a)–(2.8d) for c1, c−1,
and c2,
and searching for a value of a1 that leads to a highly accurate ARK method in
solving some standard IVPs (see Section 3 and Appendix A).Table 1
Optimized ARK3 parameter
sets.
Set 1Set 2Set 3c01−441−119+411c−00−541−79+410c11216(641−1)3(9+41)24748c−1−124(341−13)3(9+41)2−148c214003(9+41)22548a15129+412045


### 2.2. ARK4

Forv=3,
we have the 8 parameters c0,c−0,c1,c−1,c2,c3,a1,
and a2 along with order equations up to O(h5) which are
(2.10a)O(h0)c0−c−0=1,(2.10b)O(h1)c−0+c1−c−1=1,(2.10c)O(h2)−12c−0+c−1+c2+c3=12,(2.10d)O(h3)−112c−0+c2a1+c3a2=512,(2.10e)O(h4)c2a12+c3a22=13,(2.10f)O(h4)c3a1a2=16,(2.10g)O(h5)1120c−0+c2a13+c3a23=31120,(2.10h)O(h5)1240c−0+c3a1a22=31240,(2.10i)O(h5)1360c−0+c3a12a2=31360,(2.10j)O(h5)c−0=31.Since there are 10 order equations up toO(h5),
we cannot achieve a 5th-order method. However, we can achieve a 4th-order
method ARK4 by solving the first 6 order equations (2.10a)–(2.10f). Therefore, ARK4
requires 3 function evaluations per time step. We have chosen to solve these
equations for the c's while keeping a1 and a2 as free parameters because all order equations
are linear in terms of the c's which makes it straightforward to solve for
them.We usea1 and a2 to optimize ARK4's accuracy by solving two of
the O(h5) equations. We choose (2.10g) and
(2.10h) because
they were used in deriving the other O(h5) equations. Set 2 and Set 3 of
Table 2 are the
two solutions. Set 1 of
Table 2 is a set in which c−0=0 and is produced by solving (2.10b) to
(2.10g). It is
a set of approximate values and naturally contains some round-off error. The
effects of this round-off error (in parameter values) can be significant if the
error produced by the ARK methods is small enough to be of the same order. Therefore,
to measure the accuracy of the ARK methods correctly (see Section 3), we have
calculated the parameter values using 25 digits of computational precision in
Maple. Note, however, that such high precision is not required for the ARK
methods to be effective.Table 2
Optimized ARK4 parameter
sets.
Set 1Set 2Set 3c01.0−441−119+41−441−119+41c−00.0−541−79+41−541−79+41c11.017627673204495246749635163641−1(9+41)2163641−1(9+41)2c−10.0176276732044952467496350843341−13(9+41)243341−13(9+41)2c2−0.133003777809752528077129302003(9+41)2c30.61537610460525728132749424003(9+41)22003(9+41)2a10.35888611391988193765959429+41409+4120a20.75466023484835962323552579+41209+4120


### 2.3. ARK4-4

Forv=4,
we have the 10 parameters c0,c−0,c1,c−1,c2,c3,c4,a1,a2, and a3.
Order equations up to O(h5) are
(2.11a)O(h0)c0−c−0=1,(2.11b)O(h1)c−0+c1−c−1=1,(2.11c)O(h2)−12c−0+c−1+c2+c3+c4=12,(2.11d)O(h3)−112c−0+c2a1+c3a2+c4a3=512,(2.11e)O(h4)c2a12+c3a22+c4a32=13,(2.11f)O(h4)c3a1a2+c4a2a3=16,(2.11g)O(h5)1120c−0+c2a13+c3a23+c4a33=31120,(2.11h)O(h5)1240c−0+c3a1a22+c4a2a32=31240,(2.11i)O(h5)1360c−0+c3a12a2+c4a22a3=31360,(2.11j)O(h5)1720c−0+c4a1a2a3=31720.
Since there are exactly 10 order equations up to O(h5),
it would seem that we can achieve a 5th-order method. However, the only two
solutions of the above order equations have c−0=1 and c−0=25,
both of which produce unstable methods (see condition (6.16)).Therefore, we setc−0=0 which leaves us with the 8 parameters c1, c−1, c2, c3, c4, a1, a2,
and a3.
We solve (2.11b)–(2.11f) for the c's and solve three of the four O(h5) order equations using the free parameters a1, a2,
and a3.
The result is a 4th-order method called ARK4-4 that is more accurate than ARK4
for most of the IVPs that are investigated herein. The “-4” indicates
that this 4th-order method requires 4 function evaluations, unlike ARK4 which
requires only 3.Set 1 of Table3 satisfies the higher order equations
(2.11g),
(2.11h), and
(2.11j). Set 3 of Table 3 satisfies the higher order equations
(2.11g),
(2.11h), and
(2.11i). Set 2 of Table 3 satisfies the higher order equations
(2.11g),
(2.11j) and
(2.12b). To understand where
(2.12b) comes from, we first need to
explain a unique property of the O(h5) order equations.Table 3
Optimized ARK4-4 parameter
sets.
Set 1Set 2Set 3c01.01.01.0c−00.00.00.0c11.0228319288392032115814110.95999836297405233577612921.038087495003156301209584c−10.02283192883920321158141016−0.040001637025947664223868920.03808749500315630120958582c2−0.045158301883180231641969730.2483344505743049392964305−0.1206952296752875905594747c3−0.08618700613581317473462200−0.44002905880512272992927910.4307688535040614391640197c40.60851337917979019479518550.73169624525676545485671520.1518388811680698501858681a10.24641898480453520276639880.21280761842314480370072750.2340555618293773386595766a20.37942760708511201070162690.38075868967914793913977410.7532489015566390666145791a30.75675617797074070285366690.72620858035488573173473520.7932084970935761571360267TheO(h5) order equations  are fundamentally different
from lower order equations. If the order equations
for ARK4-4 are derived with a two-dimensional coupled ODE system, we get the
four O(h5) order equations presented above, that is,
(2.11g)–(2.11j). However, if the derivation is done with a one-dimensional ODE or
with a two-dimensional uncoupled ODE system, we get the three O(h5) equations:
(2.12a)O(h5)1120c−0+c2a13+c3a23+c4a33=31120,(2.12b)O(h5)190c−0+c3a1a2(a1+2a2)+c4a2a3(a2+2a3)=3190,(2.12c)O(h5)1720c−0+c4a1a2a3=31720.
Hence,  coupling of the ODE system affects the O(h5) order equations but not the O(h0),…,O(h4) order equations. We see that the difference
between the two sets of order equations is that (2.11h) and
(2.11i) are replaced by
(2.12b) = 2 × (2.11h) + (2.11i). Therefore, a parameter
set that satisfies (2.11h) and
(2.11i) also automatically satisfies
(2.12b), but not
vice versa. To ensure that the order equations are valid for all ODE systems,
we have considered a two-dimensional coupled ODE system in our derivations.Although we have not derived the order equations based
on a 3- or higher-dimensional coupled ODE system, we hypothesize that the order
equations would be the same as (2.11g)–(2.11j). We base this hypothesis on two
observations. One, unlike moving from a one-dimensional system to a
two-dimensional system where the coupling property emerges, no property emerges
when we move from a two-dimensional system to a three-dimensional system, so
there is no perceivable reason for the order equations to change. Two, using a
two-dimensional coupled ODE system, we have derived the same standard RK order
equations that others have reported [13].With that explanation, it might seem counterintuitive
to use a parameter set that satisfies (2.12b). After all, this equation only
applies for uncoupled systems. However, the thinking is that for uncoupled
systems, the resulting ARK4-4 becomes a 5th-order method. While for coupled
systems, the resulting method yields a very accurate 4th-order method because
it satisfies two of the four O(h5) equations as well as a weighted average of the
other two order equations. This is confirmed by our experiments in solving
standard problems (see Section 3) since the accuracy of Set 2 is comparable to that
of Set 1 and Set 3.


### 2.4. ARK5

Forv=5,
we have the 12 parameters c0,c−0,c1,c−1,c2,c3,c4,c5,a1,a2,a3, and a4.
Order equations up to O(h6) are
(2.13a)O(h0)c0−c−0=1,(2.13b)O(h1)c−0+c1−c−1=1,(2.13c)O(h2)−12c−0+c−1+c2+c3+c4+c5=12,(2.13d)O(h3)−112c−0+c2a1+c3a2+c4a3+c5a4=512,(2.13e)O(h4)c2a12+c3a22+c4a32+c5a42=13,(2.13f)O(h4)c3a1a2+c4a2a3+c5a3a4=16,(2.13g)O(h5)1120c−0+c2a13+c3a23+c4a33+c5a43=31120,(2.13h)O(h5)1240c−0+c3a1a22+c4a2a32+c5a3a42=31240,(2.13i)O(h5)1360c−0+c3a12a2+c4a22a3+c5a32a4=31360,(2.13j)O(h5)1720c−0+c4a1a2a3+c5a2a3a4=31720,(2.13k)O(h6)c2a14+c3a24+c4a34+c5a44=15,(2.13l)O(h6)c3a1a23+c4a2a33+c5a3a43=110,(2.13m)O(h6)c3a13a2+c4a23a3+c5a33a4=120,(2.13n)O(h6)c3a12a22+c4a22a32+c5a32a42=115,(2.13o)O(h6)c4a1a22a3+c5a2a32a4=140,(2.13p)O(h6)c4a1a2a32+c5a2a3a42=160,(2.13q)O(h6)c4a12a2a3+c5a22a3a4=130,(2.13r)O(h6)c5a1a2a3a4=1120.Solving the first 10 equations, (2.13a)–(2.13j), gives us
a 5th-order method ARK5 which requires 5 function evaluations per time step. We
solve the first 8 equations for the c's and solve equations (2.13i)–(2.13l) for a1,…,a4.
Among the multiple solutions of this system, we have found Set 3 in
Table 4 to
be the most accurate in solving the standard problems given in Appendix A.
Setting c−0=0 results in 10 parameters and 9 equations of
order up to O(h5). We use the one free parameter to solve (2.13k)
which leads to numerous solutions. Among them, Set 1 and Set 2 in
Table 4 are
the most accurate in solving the standard initial value problems in Section 3.Table 4
Optimized ARK5 parameter
sets.
Set 1Set 2Set 3c01.01.01.871204587171582065174140c−00.00.00.8712045871715820651713061c11.0555621513716989365889960.84781861161579177688825250.2696466886663821637128020c−10.05556215137169893658900796−0.15218138838420822311175440.1408512758379642288874380c2−0.15507826549018113423494420.63424822240505828729250600.3158759465556997630808750c30.42592470856062909111684540.051958763825071413882297940.3212830748049407866018770c40.1103009310583581269934950−0.25919009955146520907640610.1591061035393050004573704c50.063290474499494979535563050.2251645017055437310133241−0.001514107152118746437838297a10.21634433210095616972608890.97101495143869389525856860.5094586945643958664798805a20.7355421089142943499801371−0.25561031463318690045865660.5161588401001171574027862a30.70463958528507163869393351.0945995422706924901951021.041695566100089398625120a40.93551217959468840143281400.43431677438762241454203282.134538676833492640695294


## 2.1. ARK3

Forv=2,
from (2.5), we have the 6 parameters c0,c−0,c1,c−1,c2, and a1.
The derived order equations up to O(h4) are
(2.8a)O(h0)c0−c−0=1,(2.8b)O(h1)c−0+c1−c−1=1,(2.8c)O(h2)−12c−0+c−1+c2=12,(2.8d)O(h3)−112c−0+c2a1=512,(2.8e)O(h4)c2a12=13,(2.8f)O(h4)0=16.
Clearly, (2.8f) cannot be satisfied, so it is impossible to achieve a
4th-order method with 2 function evaluations per time step. It is, however,
possible to achieve a 3rd-order method ARK3 with 2 function evaluations per
time step by satisfying (2.8a)–(2.8d). We solve these equations in terms of c0,c−0,c1,
and c−1 and let the remaining two parameters c2 and a1 become free parameters.We use sets of free parameters to optimize the
accuracy of an ARK method by minimizing its local error. The local error is represented
by higher order equations (order 4 and above in case of ARK3). The closer these
order equations are to being satisfied, the smaller the local error is.
Therefore, we try to minimize the difference between the LHS and RHS of these
order equations. This can be done in several different ways. One popular
technique is to minimize a general bound on the error [5]. A second technique is to
minimize a norm of the error [13]. A third technique is to satisfy as many individual
order equations as possible (hinted to by [4]). In general, there is no way of knowing a priori
which technique yields the most accurate ARK method since the error ultimately
depends on the problem at hand. We have chosen the last technique for two
reasons. One, it is straightforward to implement exactly. Two, for some classes
of problems, the differential term multiplying the unsolved higher order
equations may vanish, which could effectively increases the order of the
method. For instance, the differential term multiplying (2.8f) may vanish, and
provided (2.8e) is already satisfied, the resulting ARK method (with 2 function
evaluations) will become a 4th-order method.Using the two free parametersc2 and a1,
we can solve two higher-order equations: (2.8e) and any one of the two 5th-order
equations,
(2.9a)O(h5)1120c−0+c2a13=31120,(2.9b)O(h5)c−0=31.
However, we cannot choose (2.9b) because to prove stability, we will need −1≤c−0<1 (see condition (6.16)). Therefore, we choose
(2.8e) and (2.9a), which leads to Set 2 of Table 1. We also present a set in which c−0=0 because doing so can increase the stability
(see Section 6) and speed (see Section 4) of an ARK method. Setting c−0=0 and solving (2.8a)–(2.8e) leads to Set 3 of Table 1. Set 1 of
Table 1 is the only set of parameters here that is not found by
solving all possible higher order equations. It is found by setting c−0=0,
solving (2.8a)–(2.8d) for c1, c−1,
and c2,
and searching for a value of a1 that leads to a highly accurate ARK method in
solving some standard IVPs (see Section 3 and Appendix A).Table 1
Optimized ARK3 parameter
sets.
Set 1Set 2Set 3c01−441−119+411c−00−541−79+410c11216(641−1)3(9+41)24748c−1−124(341−13)3(9+41)2−148c214003(9+41)22548a15129+412045


## 2.2. ARK4

Forv=3,
we have the 8 parameters c0,c−0,c1,c−1,c2,c3,a1,
and a2 along with order equations up to O(h5) which are
(2.10a)O(h0)c0−c−0=1,(2.10b)O(h1)c−0+c1−c−1=1,(2.10c)O(h2)−12c−0+c−1+c2+c3=12,(2.10d)O(h3)−112c−0+c2a1+c3a2=512,(2.10e)O(h4)c2a12+c3a22=13,(2.10f)O(h4)c3a1a2=16,(2.10g)O(h5)1120c−0+c2a13+c3a23=31120,(2.10h)O(h5)1240c−0+c3a1a22=31240,(2.10i)O(h5)1360c−0+c3a12a2=31360,(2.10j)O(h5)c−0=31.Since there are 10 order equations up toO(h5),
we cannot achieve a 5th-order method. However, we can achieve a 4th-order
method ARK4 by solving the first 6 order equations (2.10a)–(2.10f). Therefore, ARK4
requires 3 function evaluations per time step. We have chosen to solve these
equations for the c's while keeping a1 and a2 as free parameters because all order equations
are linear in terms of the c's which makes it straightforward to solve for
them.We usea1 and a2 to optimize ARK4's accuracy by solving two of
the O(h5) equations. We choose (2.10g) and
(2.10h) because
they were used in deriving the other O(h5) equations. Set 2 and Set 3 of
Table 2 are the
two solutions. Set 1 of
Table 2 is a set in which c−0=0 and is produced by solving (2.10b) to
(2.10g). It is
a set of approximate values and naturally contains some round-off error. The
effects of this round-off error (in parameter values) can be significant if the
error produced by the ARK methods is small enough to be of the same order. Therefore,
to measure the accuracy of the ARK methods correctly (see Section 3), we have
calculated the parameter values using 25 digits of computational precision in
Maple. Note, however, that such high precision is not required for the ARK
methods to be effective.Table 2
Optimized ARK4 parameter
sets.
Set 1Set 2Set 3c01.0−441−119+41−441−119+41c−00.0−541−79+41−541−79+41c11.017627673204495246749635163641−1(9+41)2163641−1(9+41)2c−10.0176276732044952467496350843341−13(9+41)243341−13(9+41)2c2−0.133003777809752528077129302003(9+41)2c30.61537610460525728132749424003(9+41)22003(9+41)2a10.35888611391988193765959429+41409+4120a20.75466023484835962323552579+41209+4120


## 2.3. ARK4-4

Forv=4,
we have the 10 parameters c0,c−0,c1,c−1,c2,c3,c4,a1,a2, and a3.
Order equations up to O(h5) are
(2.11a)O(h0)c0−c−0=1,(2.11b)O(h1)c−0+c1−c−1=1,(2.11c)O(h2)−12c−0+c−1+c2+c3+c4=12,(2.11d)O(h3)−112c−0+c2a1+c3a2+c4a3=512,(2.11e)O(h4)c2a12+c3a22+c4a32=13,(2.11f)O(h4)c3a1a2+c4a2a3=16,(2.11g)O(h5)1120c−0+c2a13+c3a23+c4a33=31120,(2.11h)O(h5)1240c−0+c3a1a22+c4a2a32=31240,(2.11i)O(h5)1360c−0+c3a12a2+c4a22a3=31360,(2.11j)O(h5)1720c−0+c4a1a2a3=31720.
Since there are exactly 10 order equations up to O(h5),
it would seem that we can achieve a 5th-order method. However, the only two
solutions of the above order equations have c−0=1 and c−0=25,
both of which produce unstable methods (see condition (6.16)).Therefore, we setc−0=0 which leaves us with the 8 parameters c1, c−1, c2, c3, c4, a1, a2,
and a3.
We solve (2.11b)–(2.11f) for the c's and solve three of the four O(h5) order equations using the free parameters a1, a2,
and a3.
The result is a 4th-order method called ARK4-4 that is more accurate than ARK4
for most of the IVPs that are investigated herein. The “-4” indicates
that this 4th-order method requires 4 function evaluations, unlike ARK4 which
requires only 3.Set 1 of Table3 satisfies the higher order equations
(2.11g),
(2.11h), and
(2.11j). Set 3 of Table 3 satisfies the higher order equations
(2.11g),
(2.11h), and
(2.11i). Set 2 of Table 3 satisfies the higher order equations
(2.11g),
(2.11j) and
(2.12b). To understand where
(2.12b) comes from, we first need to
explain a unique property of the O(h5) order equations.Table 3
Optimized ARK4-4 parameter
sets.
Set 1Set 2Set 3c01.01.01.0c−00.00.00.0c11.0228319288392032115814110.95999836297405233577612921.038087495003156301209584c−10.02283192883920321158141016−0.040001637025947664223868920.03808749500315630120958582c2−0.045158301883180231641969730.2483344505743049392964305−0.1206952296752875905594747c3−0.08618700613581317473462200−0.44002905880512272992927910.4307688535040614391640197c40.60851337917979019479518550.73169624525676545485671520.1518388811680698501858681a10.24641898480453520276639880.21280761842314480370072750.2340555618293773386595766a20.37942760708511201070162690.38075868967914793913977410.7532489015566390666145791a30.75675617797074070285366690.72620858035488573173473520.7932084970935761571360267TheO(h5) order equations  are fundamentally different
from lower order equations. If the order equations
for ARK4-4 are derived with a two-dimensional coupled ODE system, we get the
four O(h5) order equations presented above, that is,
(2.11g)–(2.11j). However, if the derivation is done with a one-dimensional ODE or
with a two-dimensional uncoupled ODE system, we get the three O(h5) equations:
(2.12a)O(h5)1120c−0+c2a13+c3a23+c4a33=31120,(2.12b)O(h5)190c−0+c3a1a2(a1+2a2)+c4a2a3(a2+2a3)=3190,(2.12c)O(h5)1720c−0+c4a1a2a3=31720.
Hence,  coupling of the ODE system affects the O(h5) order equations but not the O(h0),…,O(h4) order equations. We see that the difference
between the two sets of order equations is that (2.11h) and
(2.11i) are replaced by
(2.12b) = 2 × (2.11h) + (2.11i). Therefore, a parameter
set that satisfies (2.11h) and
(2.11i) also automatically satisfies
(2.12b), but not
vice versa. To ensure that the order equations are valid for all ODE systems,
we have considered a two-dimensional coupled ODE system in our derivations.Although we have not derived the order equations based
on a 3- or higher-dimensional coupled ODE system, we hypothesize that the order
equations would be the same as (2.11g)–(2.11j). We base this hypothesis on two
observations. One, unlike moving from a one-dimensional system to a
two-dimensional system where the coupling property emerges, no property emerges
when we move from a two-dimensional system to a three-dimensional system, so
there is no perceivable reason for the order equations to change. Two, using a
two-dimensional coupled ODE system, we have derived the same standard RK order
equations that others have reported [13].With that explanation, it might seem counterintuitive
to use a parameter set that satisfies (2.12b). After all, this equation only
applies for uncoupled systems. However, the thinking is that for uncoupled
systems, the resulting ARK4-4 becomes a 5th-order method. While for coupled
systems, the resulting method yields a very accurate 4th-order method because
it satisfies two of the four O(h5) equations as well as a weighted average of the
other two order equations. This is confirmed by our experiments in solving
standard problems (see Section 3) since the accuracy of Set 2 is comparable to that
of Set 1 and Set 3.


## 2.4. ARK5

Forv=5,
we have the 12 parameters c0,c−0,c1,c−1,c2,c3,c4,c5,a1,a2,a3, and a4.
Order equations up to O(h6) are
(2.13a)O(h0)c0−c−0=1,(2.13b)O(h1)c−0+c1−c−1=1,(2.13c)O(h2)−12c−0+c−1+c2+c3+c4+c5=12,(2.13d)O(h3)−112c−0+c2a1+c3a2+c4a3+c5a4=512,(2.13e)O(h4)c2a12+c3a22+c4a32+c5a42=13,(2.13f)O(h4)c3a1a2+c4a2a3+c5a3a4=16,(2.13g)O(h5)1120c−0+c2a13+c3a23+c4a33+c5a43=31120,(2.13h)O(h5)1240c−0+c3a1a22+c4a2a32+c5a3a42=31240,(2.13i)O(h5)1360c−0+c3a12a2+c4a22a3+c5a32a4=31360,(2.13j)O(h5)1720c−0+c4a1a2a3+c5a2a3a4=31720,(2.13k)O(h6)c2a14+c3a24+c4a34+c5a44=15,(2.13l)O(h6)c3a1a23+c4a2a33+c5a3a43=110,(2.13m)O(h6)c3a13a2+c4a23a3+c5a33a4=120,(2.13n)O(h6)c3a12a22+c4a22a32+c5a32a42=115,(2.13o)O(h6)c4a1a22a3+c5a2a32a4=140,(2.13p)O(h6)c4a1a2a32+c5a2a3a42=160,(2.13q)O(h6)c4a12a2a3+c5a22a3a4=130,(2.13r)O(h6)c5a1a2a3a4=1120.Solving the first 10 equations, (2.13a)–(2.13j), gives us
a 5th-order method ARK5 which requires 5 function evaluations per time step. We
solve the first 8 equations for the c's and solve equations (2.13i)–(2.13l) for a1,…,a4.
Among the multiple solutions of this system, we have found Set 3 in
Table 4 to
be the most accurate in solving the standard problems given in Appendix A.
Setting c−0=0 results in 10 parameters and 9 equations of
order up to O(h5). We use the one free parameter to solve (2.13k)
which leads to numerous solutions. Among them, Set 1 and Set 2 in
Table 4 are
the most accurate in solving the standard initial value problems in Section 3.Table 4
Optimized ARK5 parameter
sets.
Set 1Set 2Set 3c01.01.01.871204587171582065174140c−00.00.00.8712045871715820651713061c11.0555621513716989365889960.84781861161579177688825250.2696466886663821637128020c−10.05556215137169893658900796−0.15218138838420822311175440.1408512758379642288874380c2−0.15507826549018113423494420.63424822240505828729250600.3158759465556997630808750c30.42592470856062909111684540.051958763825071413882297940.3212830748049407866018770c40.1103009310583581269934950−0.25919009955146520907640610.1591061035393050004573704c50.063290474499494979535563050.2251645017055437310133241−0.001514107152118746437838297a10.21634433210095616972608890.97101495143869389525856860.5094586945643958664798805a20.7355421089142943499801371−0.25561031463318690045865660.5161588401001171574027862a30.70463958528507163869393351.0945995422706924901951021.041695566100089398625120a40.93551217959468840143281400.43431677438762241454203282.134538676833492640695294


## 3. Accuracy

We consider
seven initial value problems (IVPs) to illustrate the accuracy of the ARK
methods. These problems are taken from examples used by previous researchers
[11, 14]. They are named IVP-1
through IVP-7, and are listed in Appendix A. Most of the problems model real-world
phenomena (suggested by [14]) and have been used by other researchers to test
numerical integrators. Specifically, IVP-3 models the Euler equations for a rigid body
without external forces, IVP-4 and IVP-5 model a one-body gravitational problem
in the plane, IVP-6 models a radioactive decay chain, and IVP-7 models the
gravitational problem for five planets around the sun (and four inner planets).
The first two example problems (IVP-1 and IVP-2) are one-dimensional ODEs—the
first, autonomous, the second, nonautonomous. The dimension of the rest of the
examples, which include linear and non-linear ODE systems, ranges from 3 to 30.The approximate solution to each example ODE problem
is computed for0≤t≤15 at the step sizes: 0.001, 0.0025, 0.005, 0.01,
0.025, 0.05, and 0.1 using ARK3, ARK4, ARK4-4, and ARK5. For comparison, we
also show the corresponding results using the RK2, RK3, RK4, and RK5 methods.
Each ARK method is started with an RK method of the same order, taking 10 steps
in the first time step at 1/10th the value of step size. For ARK methods, Set 1
parameter values of Tables 1,
2, 3, and
4 are used, and for RK methods, the
parameter values in
Table 5 are used.Table 5
RK parameters in Butcher's
notation [15]: (a) RK2,
(b) RK3, (c) RK4, and (d) RK5.The accuracy plots for problems IVP-1 through IVP-7
are given in Figures2,
3,
4,
5,
6,
7, and
8, respectively. Along the ordinate
is plotted an average of the 2-norm of the global error, which is computed by
taking the 2-norm of the difference between the solution y(tn) and approximate solution yn,
and then averaging this norm over the interval 10≤t≤15.
The exception is IVP-5 for which the average absolute value of global error  in all four
components of the solution are plotted separately. The results obtained when
using RK methods are shown by lines with circles along them, while the results
for ARK methods are shown by lines with squares along them. Lines with the same
color indicate methods that have the same number of function evaluations.Figure 2
Accuracy plot for IVP-1 using log-log scale. The RK
methods are shown by lines with circles. The ARK methods are shown by lines
with squares. Lines of the same color correspond to methods with the same
number of function evaluations per time step.Figure 3
Accuracy plot for IVP-2 (see caption of
Figure2).Figure 4
Accuracy plot for IVP-3 (see caption of
Figure2).Figure 5
Accuracy plot for IVP-4 (see caption of Figure2).Accuracy plots for IVP-5: (a)y1,
(b) y2,
(c) y3,
and (d) y4 (see caption of Figure 2).
(a)(b)(c)(d)Figure 7
Accuracy plot for
IVP-6 (see caption of Figure2).Figure 8
Accuracy plot for IVP-7 (see caption of
Figure2).Due to small step sizes, for example10−3, and the high order of local accuracy for some
of the RK and ARK methods, the global error for some problems reaches 10−24.
This is beyond MATLAB's precision capabilities. MATLAB has a minimum round-off
error of about eps≈2.3×10−16 (at magnitude 1) and a maximum of 17
significant digits of accuracy. Therefore, a proper analysis of the accuracy of
the ARK and RK methods cannot be done in MATLAB. In fact, we found that for
small step sizes, the higher order ARK and RK plot lines produced by MATLAB are
erroneous; they are jagged lines with slopes different than the order of the
method. That is why the PARI mathematics environment [16] which provides 28 digits of accuracy by default was
used for all calculations in this section. This ensures that the computational
round-off error is far smaller than the methods' global error. Note however,
that using high precision computation is not necessary for the ARK methods to
be effective. The ARK methods should retain their accuracy and efficiency advantages
over the RK methods at any precision.The slopes of the plot lines confirm the order of our
ARK methods; the slopes of ARK3, ARK4, ARK4-4, and ARK5 lines (except when the
time steps are large) are 3, 4, 4, and 5, respectively. The slope of the ARK5
lines for problems IVP-3 to IVP-7 (Figures4 to
8) is especially important. We
hypothesized in Section 2.3 that although our O(h5) order equations were derived based on a
two-dimensional ODE system, they would also apply to higher-dimensional
problems. The fact that problems IVP-3 to IVP-7 are 3-dimensional or higher and
their ARK5 accuracy plot lines have slope 5 (for small step sizes) confirms our
hypothesis.Comparing ARK and RK methods, the plots show that
generally speaking, the global error of the ARK methods is greater than that of
the corresponding RK methods by about one order of magnitude when comparing
methods of the same order, such as comparing ARK3 versus RK3. Although great
effort was made to find the parameter sets that yield the most accurate ARK
methods, it is possible that better sets can be found which would considerably
improve the accuracy of the ARK methods. Set 1 of ARK4 exemplifies this claim.
It is seen from the accuracy plots that for problems IVP-2, IVP-3, IVP-5, and
IVP-7, ARK4 has virtually the same accuracy as RK4.A more useful comparison of ARK and RK methods would
be to compare those that have the same number of function evaluations per time
step. Table6 summarizes the order and number of function evaluations for the
ARK and RK methods. In comparing ARK3 versus RK2, ARK4 versus RK3, and ARK4-4
versus RK4 in the accuracy plots, ARK methods prove to be more accurate. Because
ARK3 and ARK4 are higher-order methods than RK2 and RK3, respectively, their
gain in accuracy grows as the step size decreases. For example, for IVP-5 (see
Figure 6) at step size 0.1, ARK3 and ARK4 are less than one order of magnitude
more accurate than RK2 and RK3, respectively, but at step size 0.001, they are
three and four orders of magnitude more accurate than RK2 and RK3,
respectively. A similar gain of accuracy can occur even for ARK4-4 versus RK4
when for some problems (see Figures 2 and
7), ARK4-4 becomes a 5th-order
method.Table 6
Summary of
order and number of function evaluations per time step for RK methods and ARK
methods presented herein.
RK2ARK3RK3ARK4RK4ARK4-4RK5ARK5Order of methodO(h2)O(h3)O(h3)O(h4)O(h4)O(h4)O(h5)O(h5)No.
of function evaluations22334465


## 4. Speedup

As mentioned
before, the ARK3, ARK4, and ARK5 methods require 2, 3, and 5 function
evaluations per time step of integration, respectively, while the RK3, RK4, and
RK5 methods require 3, 4, and 6 function evaluations per time step,
respectively. It is commonly believed that function evaluations constitute the
main computational cost of RK-type methods. Hence, we can expect that compared
to the RK3, RK4, and RK5 methods, the ARK3, ARK4, and ARK5 methods have
speedups of approximately(3−2)/3≈33%, (4−3)/4=25%,
and (6−5)/6≈17%, respectively.To verify these speedups experimentally, we have
implemented the methods in MATLAB. We have measured the time it takes to solve
IVP-2 and IVP-4 for0≤t≤1000 (with such a large time span needed to get a reliable
time measurement) at the 7 different step sizes used in Section 3. (We have
used the MATLAB utility TIC and TOC for time measurements). We repeated this
integration 100 times and took the average of integration times before
calculating the speedup (tRK−tARK)/tRK.
Figures 9 and
10 display the results. Speedup depends on the problem being
solved, and among IVP-1 through IVP-7, IVP-2 leads to
the lowest speedup, while IVP-1 leads to the highest. Speedup also depends on
the parameter set used because it affects the computational overhead. The ARK
implementations here use values of Set 1 in Tables 1,
2, and
4, and the RK
implementations use parameter values in Table 5. The plots show a linear least
square fit to the points. We expect this line to be fairly flat because the
number of steps for the ARK and RK methods are only slightly different (caused
by the startup of the ARK methods at the first time step) and because speedup
does not depend on the step size.Figure 9
Speedup of ARK methods compared with RK methods for
IVP-2.Figure 10
Speedup of ARK methods compared with RK methods for
IVP-4.The speedups seen in Figures9 and
10 are lower than
expected. Moreover, we have found the time measurement and consequently the
speedup to be machine dependent. Through measurements of the computation times
at different points in our computer codes, we have found two main causes for
the less-than-expected speedup of our ARK methods. First, the ARK methods
require additional multiplication and addition operations when evaluating the
expression for yn+1 than the RK methods require. This accounts for
about half of the loss in the expected speedup. Second, the ARK methods require
a number of assignment operations to set yn−1,k−1,… to yn,k1,… in every step of integration which the RK
methods do not require. Somewhat surprisingly, this takes a significant amount
of computation time in the MATLAB environment and is responsible for about
the other half of the loss in the expected speedup.


## 5. Comparison with other Two-Step Methods

We now compare
the ARK4 method (with parameters of Set 1 of Table2) to the more sophisticated
TSRK4 method presented in [12]. Both are fourth-order methods with 3 function
evaluations per time step. The accuracy of the two methods is calculated in the
same way as in Section 3 but implemented in MATLAB. Because of MATLAB's lower
precision, we have not calculated the accuracy at the two smallest step sizes.
The average of the 2-norm of the global error for problems IVP-1, IVP-2, and
IVP-4 is presented in Figures
11,
12, and
13. The TSRK4 method is more than one
order of magnitude more accurate in solving IVP-1; the ARK4 method is less than
one order of magnitude more accurate than the TSRK4 method in solving IVP-2;
and the TSRK4 method is slightly more accurate than ARK4 in solving IVP-4. The
TSRK4 method is also more accurate in solving the other problems listed in
Appendix A. This is because the TSRK4 method is more general and possesses more
parameters, which have been effectively optimized. In
Figure 14, we have
plotted the relative change in computation time when solving IVP-1, IVP-2, and
IVP-4 (similar to Section 4). The lower computation time (see Figure 14)
required by the ARK4 method is due to its lower computational overhead, which
is produced by its simpler form as given in (2.5).Figure 11
Comparison of the accuracy of
the ARK4  (blue) and TSRK4 (red) methods in solving IVP-1.Figure 12
Comparison of the accuracy of the ARK4 (blue) and TSRK4
(red) methods in solving IVP-2.Figure 13
Comparison of the accuracy of the ARK4 (blue) and TSRK4 (red) methods in solving IVP-4.Figure 14
Relative change in computation time of the TSRK4
method compared with the ARK4 method in solving IVP-1 (blue), IVP-2 (red), and IVP-4 (black),
showing that TSRK4 requires greater computation time.


## 6. Stability

We have
already established that the error of an ARKp method over a single time step (local error)
has the order O(hp+1) (see Section 1). In this section, we will look
at how this error accumulates over many time steps (global error) by looking at
two cases, small step sizes (approaching
zero, D-stability) and large step
sizes. Similar results are provided in [12] in the form of (a) a bound on the global error for
small step sizes and (b) a comparable stability region for finite step sizes.


### 6.1. Small Step Sizes

The ARK method
given by expression (2.5) yields an approximate solution {yn}(2≤n≤N) to the initial value problem (1.3) and can be
written as(6.1)[yn+1yn]=[c0I−c−0II0][ynyn−1]+[c1k1−c−1k−1+∑i=2vci(ki−k−i)0]for 1≤n≤N−1,
where I is an m by m identity matrix. We can write (6.1)
as(6.2)y¯n+1=Qy¯n+ϕ¯(y¯n,f,h)for1≤n≤N−1,where(6.3)y¯n+1=[yn+1yn],y¯n=[ynyn−1],Q is the 2m by 2m block matrix given by(6.4)Q=[c0I−c−0II0],and ϕ¯ is the 2m by 1 vector defined as(6.5)ϕ¯=[ϕ0]=[c1k1−c−1k−1+∑i=2vci(ki−k−i)0].With a slightly different
initial condition z0=y0+Δ0 at time t0,
a slightly different approximate solution z1=y1+Δ1 at time t1,
and a slightly different recipe ϕ(z¯n,f,h)+hδn(1≤n≤N−1),
the ARK method yields a “perturbed” approximate solution {zn}(2≤n≤N),
and can be written as(6.6)[zn+1zn]=[c0I−c−0II0][znzn−1]+[c1k1−c−1k−1+∑i=2vci(ki−k−i)0]+h[δn0]for 1≤n≤N−1, or can be written as(6.7)z¯n+1=Qz¯n+ϕ¯(z¯n,f,h)+hδ¯nfor1≤n≤N−1,where(6.8)δ¯n=[δn0].To prove stability, we will show that under some
conditions,∥zn−yn∥∞ is bounded. To prove convergence, we will show
that ∥y(tn)−yn∥∞→0 as h→0,
where y(tn) is the exact solution of the IVP defined by
(1.3). In order to do this, we will assume that the ARK method is stable and
that as h tends to zero, ∥y(t1)−y1∥∞ tends to zero, where t0 is the initial time and t1=t0+h.
The stability analysis comprises of two lemmas and two theorems. In Lemma 6.1, we
will consider the eigenvalue problem for the matrix Q.
In Lemma 6.2, we will show that the function ϕ¯ is Lipschitz continuous. In Theorem 6.3, we will
establish stability of the ARK methods by looking at the difference of (6.2) and
(6.7). In Theorem 6.4, we will prove convergence for ARK methods using the result
of Theorem 6.3. This stability analysis is inspired by Shampine [17].Lemma 6.1.
The eigenvalue problem for the matrixQ given in (6.4)  can be written as(6.9)QW=WD.If −1≤c−0<1,
the matrix W is nonsingular and satisfies (6.10)∥W−1QW∥=1,(6.11)∥W∥=2,(6.12)w:=∥W−1∥=2|c−0−1|, where here and throughout the
stability analysis, all norms denote the infinity norm ∥⋅∥∞.Proof.
Taking(6.13)W=[Ic−0III],D=[I00c−0I],where we have made use of the O(h0) order equation, c0=1+c−0,
we observe that (6.9) is satisfied. Noting that(6.14)W−1=1c−0−1[−Ic−0II−I],we require that c−0≠1 for W−1 to exist. Recalling that ∥⋅∥≡∥⋅∥∞,
we can write(6.15)∥W−1QW∥=∥D∥=max(absoluterowsum)=max(1,|c−0|).With the
condition(6.16)−1≤c−0<1,we can see from (6.15) that ∥W−1QW∥=1.
From condition (6.16) and expression (6.13), we find that(6.17)∥W∥=max(1+|c−0|,2)=2,and from expression (6.14), we
obtain(6.18)w:=∥W−1∥=1|c−0−1|max(1+|c−0|,2)=2|c−0−1|.Thus, Lemma 6.1 is proven.Lemma 6.2.
Iff is a Lipschitz continuous function so
that(6.19)∥f(zn)−f(yn)∥≤L∥zn−yn∥foraconstantL,andanyzn,yn(1≤n≤N),then, so is ϕ¯,
that is,(6.20)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≤hL^∥z¯n−y¯n∥,foraconstantL^,andanyz¯n,y¯n(1≤n≤N).Proof.
First, we will
show by induction that(6.21)∥ki(zn)−ki(yn)∥≤hγi∥zn−yn∥foraconstantγi(1≤i≤v),where we recall that ∥⋅∥≡∥⋅∥∞, and that v is the number of function evaluations of the
ARK method. Also, here and throughout this proof, 1≤n≤N.
For i=1, we have(6.22)∥k1(zn)−k1(yn)∥=∥hf(zn)−hf(yn)∥≤hL∥zn−yn∥=hγ1∥zn−yn∥,where γ1=L.
Assuming inequality (6.21) is true for i=j,
we have(6.23)∥kj+1(zn)−kj+1(yn)∥=∥hf(zn+ajkj(zn))−hf(yn+ajkj(yn))∥≤hL∥zn−yn+aj(kj(zn)−kj(yn))∥≤hL[∥zn−yn∥+|aj|hγj∥zn−yn∥]=hL[1+hγj|aj|]∥zn−yn∥=hγj+1∥zn−yn∥,where γj+1=L[1+hγj|aj|].
Hence, inequality (6.21) holds for j+1 and, therefore, it holds for 1≤i≤v.
Similarly, we can show that(6.24)∥k−i(zn−1)−k−i(yn−1)∥≤hγi∥zn−1−yn−1∥for1≤i≤v.
Functionϕ can be written as(6.25)ϕ=c1k1−c−1k−1+∑i=2vci(ki−k−i)=∑i=1vciki−∑i=1vc−ik−i,wherec−i=cifor2≤i≤v=ϕ1−ϕ−1.Using inequality (6.21), we can
write(6.26)∥ϕ1(zn)−ϕ1(yn)∥=∥∑i=1vciki(zn)−∑i=1vciki(yn)∥≤∑i=1v|ci|∥ki(zn)−ki(yn)∥≤h[∑i=1v(|ci|γi)]∥zn−yn∥.Similarly, using inequality
(6.24), we can write(6.27)∥ϕ−1(zn−1)−ϕ−1(yn−1)∥≤h[∑i=1v(|c−i|γi)]∥zn−1−yn−1∥.Since ∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≡∥ϕ¯(z¯n)−ϕ¯(y¯n)∥∞=max(∥ϕ(z¯n)−ϕ(y¯n)∥,0),
we have(6.28)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥=∥ϕ(z¯n)−ϕ(y¯n)∥=∥ϕ1(zn)−ϕ−1(zn−1)−(ϕ1(yn)−ϕ−1(yn−1))∥≤∥ϕ1(zn)−ϕ1(yn)∥+∥ϕ−1(zn−1)−ϕ−1(yn−1)∥≤h[∑i=1v(|ci|γi)]∥zn−yn∥+h[∑i=1v(|c−i|γi)]∥zn−1−yn−1∥.Because ∥z¯n−y¯n∥=max(∥zn−yn∥,∥zn−1−yn−1∥),
we have(6.29)∥zn−yn∥≤∥z¯n−y¯n∥,∥zn−1−yn−1∥≤∥z¯n−y¯n∥.Therefore,(6.30)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≤h[∑i=1v|ci|γi+∑i=1v|c−i|γi]∥z¯n−y¯n∥=hL^∥z¯n−y¯n∥,where L^=∑i=1v(|ci|+|c−i|)γi.
Thus, Lemma 6.2 is proven.Theorem 6.3.
Suppose that an ARK method is used to solve the IVP
given in (1.3), and that f is a Lipschitz continuous function throughout
the integration span. If −1≤c−0<1,
the ARKp method is stable for all sufficiently small
step sizes, and
(6.31)∥zn−yn∥≤Emax(maxr=0,1∥zr−yr∥,max1≤m<n∥δm∥)for2≤n≤N,where(6.32)E=4we2L^w(tn−t1)max(tn−t1,1).Proof.
Subtracting (6.7) from (6.2), we have(6.33)z¯n+1−y¯n+1=Q(z¯n−y¯n)+ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)+hδ¯nfor1≤n≤N−1.We multiply this equation by W−1 to get(6.34)W−1(z¯n+1−y¯n+1)=W−1Q(z¯n−y¯n)+W−1[ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)]+hW−1δ¯n.Let(6.35)dn=W−1(z¯n−y¯n)sothatWdn=(z¯n−y¯n)for1≤n≤N,and we get(6.36)dn+1=W−1QWdn+W−1[ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)]+hW−1δ¯n.
Taking the infinity norm and using the Lipschitz
condition onϕ¯ (Lemma 6.2), we can write(6.37)∥dn+1∥≤∥W−1QW∥∥dn∥+hL^∥W−1∥∥z¯n−y¯n∥+h∥W−1∥∥δ¯n∥for1≤n≤N−1.Considering expressions (6.10) and
(6.12), we have(6.38)∥dn+1∥≤∥dn∥+hL^w∥z¯n−y¯n∥+hw∥δ¯n∥.Also, using expressions (6.35) and
(6.11), we have(6.39)∥z¯n−y¯n∥=∥Wdn∥≤∥W∥∥dn∥=2∥dn∥,so inequality (6.38)
becomes(6.40)∥dn+1∥≤(1+2hL^w)∥dn∥+hw∥δ¯n∥for1≤n≤N−1.
Relation (6.40) establishes a bound on ∥dn+1∥ based on ∥dn∥.
This is not very useful because we do not know how large ∥dn∥ is. Instead, we want to find a bound based on ∥d1∥.
We will now show by induction that this bound can be written as(6.41)∥dn∥≤e2L^w(tn−t1)∥d1∥+hw∑m=1n−1e2L^w(tn−tm+1)∥δ¯m∥for2≤n≤N.
In the following, we will use the
inequality(6.42)1≤1+t≤etfort≥0,which follows from the
expression et=1+t+(1/2)t2+⋯.
For n=2,
the right hand side of inequality (6.41) becomes(6.43)e2L^w(t2−t1)∥d1∥+hwe2L^w(t2−t2)∥δ¯1∥=e2L^wh∥d1∥+hw∥δ¯1∥.Inequality (6.40) can be written
as(6.44)∥dn∥≤(1+2hL^w)∥dn−1∥+hw∥δ¯n−1∥for2≤n≤N,and for n=2 as(6.45)∥d2∥≤(1+2hL^w)∥d1∥+hw∥δ¯1∥≤e2L^wh∥d1∥+hw∥δ¯1∥.Comparing expressions (6.45) and
(6.43), we can see that inequality (6.41) holds for n=2.
Assuming it holds for n=k,
we have from inequality (6.40) that(6.46)∥dk+1∥≤(1+2hL^w)∥dk∥+hw∥δ¯k∥for1≤k≤N−1≤e2hL^w(e2L^w(tk−t1)∥d1∥+hw∑m=1k−1e2L^w(tk−tm+1)∥δ¯m∥)+hw∥δ¯k∥.Since for any r, h+(tk−tr)=tk+1−tr,
we have that(6.47)∥dk+1∥≤e2L^w(tk+1−t1)∥d1∥+hw∑m=1k−1e2L^w(tk+1−tm+1)∥δ¯m∥+hwe2L^w(tk+1−tk+1)∥δ¯k∥=e2L^w(tk+1−t1)∥d1∥+hw∑m=1ke2L^w(tk+1−tm+1)∥δ¯m∥for1≤k≤N−1.Therefore, inequality (6.41) holds
for k+1 and by induction, it holds for 2≤n≤N.
We can simplify inequality (6.41) by
writing(6.48)∥dn∥≤e2L^w(tn−t1)∥d1∥+hw∑m=1n−1e2L^w(tn−tm+1)∥δ¯m∥for2≤n≤N≤e2L^w(tn−t1)∥d1∥+hwe2L^w(tn−t2)(n−1)max1≤m<n∥δ¯m∥≤e2L^w(tn−t1)∥d1∥+hwe2L^w(tn−t1)(n−1)max1≤m<n∥δ¯m∥=e2L^w(tn−t1)∥W−1(z¯1−y¯1)∥+we2L^w(tn−t1)(tn−t1)max1≤m<n∥δ¯m∥≤e2L^w(tn−t1)w[∥z¯1−y¯1∥+(tn−t1)max1≤m<n∥δ¯m∥]≤e2L^w(tn−t1)w[max(tn−t1,1)∥z¯1−y¯1∥+max(tn−t1,1)max1≤m<n∥δ¯m∥]≤2we2L^w(tn−t1)max(tn−t1,1)max(∥z¯1−y¯1∥,max1≤m<n∥δ¯m∥).Using the inequality (6.39), we
have(6.49)∥z¯n−y¯n∥≤4we2L^w(tn−t1)max(tn−t1,1)max(∥z¯1−y¯1∥,max1≤m<n∥δ¯m∥).Since ∥zn−yn∥≤∥z¯n−y¯n∥,
and ∥δ¯m∥=∥δm∥,
we have(6.50)∥zn−yn∥≤Emax(maxr=0,1∥zr−yr∥,max1≤m<n∥δm∥)for2≤n≤N,where(6.51)E=4we2L^w(tn−t1)max(tn−t1,1).Hence, stability of ARK methods
is established. Recall that the ARK methods require to be provided with the
initial condition y0 at the beginning of the first time step t0 and the approximate solution y1 at the end of the first time step t1.
The maxr=0,1∥zr−yr∥ term in expression (6.31) refers to the
perturbations in these values.Theorem 6.4.
Consider an ARKp method that is used to solve the IVP given by
(1.3), where f is a sufficiently smooth Lipschitz continuous
function throughout the integration span. If −1≤c−0<1 and the approximate solution y1 at time t1 is accurate to O(hq),
the ARKp method is convergent to order min(p,q),
and (6.52)∥y(tn)−yn∥≤Emax(O(hq),O(hp))for2≤n≤N,where(6.53)E=4we2L^w(tn−t1)max(tn−t1,1).In particular, ∥y(tn)−yn∥→0 as h→0.Proof.
Consider again
the approximate solution{yn}(2≤n≤N) and the perturbed solution {zn}(2≤n≤N) of the ARKp method given by (6.54)yn+1=c0yn−c−0yn−1+ϕ(yn,yn−1,f,h)for1≤n≤N−1,(6.55)zn+1=c0zn−c−0zn−1+ϕ(zn,zn−1,f,h)+hδnfor1≤n≤N−1,where z0=y0+Δ0, z1=y1+Δ1,
and zn's recipe is perturbed by hδn.
The truncation error hδn* is defined as the difference between the exact
and approximate solutions at the end of the current step if the method is
provided with the exact solution at the beginning of the current and previous
steps. Therefore, the exact solution y(tn) satisfies(6.56)y(tn+1)=c0y(tn)−c−0y(tn−1)+ϕ(y(tn),y(tn−1),f,h)+hδn*.We will now show that {y(tn)}(2≤n≤N) is a perturbed solution of the ARK method. For n=1, we have (6.57)y2=c0y1−c−0y0+ϕ(y1,y0,f,h),(6.58)z2=c0z1−c−0z0+ϕ(z1,z0,f,h)+hδ1,(6.59)y(t2)=c0y(t1)−c−0y(t0)+ϕ(y(t1),y(t0),f,h)+hδ1*.We set(6.60)z0=y(t0)=y0,z1=y(t1)=y1+Δ1*,δ1=δ1*,and from (6.58)
and (6.59), we have z2=y(t2).
Similarly, we set δn=δn* for 1≤n≤N−1,
and from (6.55) and (6.56), we get zn=y(tn) for 2≤n≤N.
Since we have made the same assumptions as Theorem6.3,
relations (6.31) and (6.32) hold, and we have(6.61)∥y(tn)−yn∥≤Emax(∥y(t1)−y1∥,max1≤m<n∥δm*∥)for2≤n≤N.For ARKp, δm*=O(hp).
The approximate solution y1 at time t1 which is required by the ARKp method is supplied by a one-step method such
as RKp.
If y1 is accurate of order q,
that is, ∥y(t1)−y1∥=O(hq),
the above inequality shows that the ARKp method is convergent of order min(p,q),
that is,(6.62)∥y(tn)−yn∥≤Emax(O(hq),O(hp))for2≤n≤N.In particular, ∥y(tn)−yn∥→0 as h→0,
and Theorem 6.4 is proven.


### 6.2. Large Step Sizes

To produce the
stability region of ARK methods for finite step sizes, we look at the
approximate solution of the linear one-dimensional problemy′=λy,
where λ is a complex number. After some simplification
using the order equations, the ARK method's approximate solution of this problem
becomes(6.63)[yn+1yn]=[1+c−0+12(3−c−0)(λh)+S−c−0−12(1+c−0)(λh)−S10][ynyn−1],where for ARK3,(6.64)S=112(5+c−0)(λh)2,for ARK4,(6.65)S=112(5+c−0)(λh)2+16(λh)3,for ARK4-4,(6.66)S=112(5+c−0)(λh)2+16(λh)3+c4a1a2a3(λh)4,and for ARK5,(6.67)S=112(5+c−0)(λh)2+16(λh)3+1720(31−c−0)(λh)4+c5a1a2a3a4(λh)5.For stability, we require that the modulus of the eigenvalues
of the matrix in (6.63) be smaller than 1. Figures
15, 16,
17, and 18 show
the stability regions of ARK and RK methods produced in this way. The curves
represent the points at which the computed maximum modulus of the eigenvalues
equals 1. The region of stability is the area enclosed by the curves and the
real axis. It is clear from (6.63) that the stability region of ARK3 and ARK4
depends only on c−0,
but the stability of ARK4-4 and ARK5 depends also on the products c4a1a2a3 and c5a1a2a3a4, respectively.Figure 15
Stability region of ARK3 sets of Table1 and RK3 set
of Table 5.Figure 16
Stability region of ARK4
sets of Table2 and RK4 set of Table
5.Figure 17
Stability region of ARK4-4 sets of Table3 and RK4 set
of Table 5.Figure 18
Stability region of ARK5
sets of Table4 and RK5 set
of Table 5.


## 6.1. Small Step Sizes

The ARK method
given by expression (2.5) yields an approximate solution {yn}(2≤n≤N) to the initial value problem (1.3) and can be
written as(6.1)[yn+1yn]=[c0I−c−0II0][ynyn−1]+[c1k1−c−1k−1+∑i=2vci(ki−k−i)0]for 1≤n≤N−1,
where I is an m by m identity matrix. We can write (6.1)
as(6.2)y¯n+1=Qy¯n+ϕ¯(y¯n,f,h)for1≤n≤N−1,where(6.3)y¯n+1=[yn+1yn],y¯n=[ynyn−1],Q is the 2m by 2m block matrix given by(6.4)Q=[c0I−c−0II0],and ϕ¯ is the 2m by 1 vector defined as(6.5)ϕ¯=[ϕ0]=[c1k1−c−1k−1+∑i=2vci(ki−k−i)0].With a slightly different
initial condition z0=y0+Δ0 at time t0,
a slightly different approximate solution z1=y1+Δ1 at time t1,
and a slightly different recipe ϕ(z¯n,f,h)+hδn(1≤n≤N−1),
the ARK method yields a “perturbed” approximate solution {zn}(2≤n≤N),
and can be written as(6.6)[zn+1zn]=[c0I−c−0II0][znzn−1]+[c1k1−c−1k−1+∑i=2vci(ki−k−i)0]+h[δn0]for 1≤n≤N−1, or can be written as(6.7)z¯n+1=Qz¯n+ϕ¯(z¯n,f,h)+hδ¯nfor1≤n≤N−1,where(6.8)δ¯n=[δn0].To prove stability, we will show that under some
conditions,∥zn−yn∥∞ is bounded. To prove convergence, we will show
that ∥y(tn)−yn∥∞→0 as h→0,
where y(tn) is the exact solution of the IVP defined by
(1.3). In order to do this, we will assume that the ARK method is stable and
that as h tends to zero, ∥y(t1)−y1∥∞ tends to zero, where t0 is the initial time and t1=t0+h.
The stability analysis comprises of two lemmas and two theorems. In Lemma 6.1, we
will consider the eigenvalue problem for the matrix Q.
In Lemma 6.2, we will show that the function ϕ¯ is Lipschitz continuous. In Theorem 6.3, we will
establish stability of the ARK methods by looking at the difference of (6.2) and
(6.7). In Theorem 6.4, we will prove convergence for ARK methods using the result
of Theorem 6.3. This stability analysis is inspired by Shampine [17].Lemma 6.1.
The eigenvalue problem for the matrixQ given in (6.4)  can be written as(6.9)QW=WD.If −1≤c−0<1,
the matrix W is nonsingular and satisfies (6.10)∥W−1QW∥=1,(6.11)∥W∥=2,(6.12)w:=∥W−1∥=2|c−0−1|, where here and throughout the
stability analysis, all norms denote the infinity norm ∥⋅∥∞.Proof.
Taking(6.13)W=[Ic−0III],D=[I00c−0I],where we have made use of the O(h0) order equation, c0=1+c−0,
we observe that (6.9) is satisfied. Noting that(6.14)W−1=1c−0−1[−Ic−0II−I],we require that c−0≠1 for W−1 to exist. Recalling that ∥⋅∥≡∥⋅∥∞,
we can write(6.15)∥W−1QW∥=∥D∥=max(absoluterowsum)=max(1,|c−0|).With the
condition(6.16)−1≤c−0<1,we can see from (6.15) that ∥W−1QW∥=1.
From condition (6.16) and expression (6.13), we find that(6.17)∥W∥=max(1+|c−0|,2)=2,and from expression (6.14), we
obtain(6.18)w:=∥W−1∥=1|c−0−1|max(1+|c−0|,2)=2|c−0−1|.Thus, Lemma 6.1 is proven.Lemma 6.2.
Iff is a Lipschitz continuous function so
that(6.19)∥f(zn)−f(yn)∥≤L∥zn−yn∥foraconstantL,andanyzn,yn(1≤n≤N),then, so is ϕ¯,
that is,(6.20)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≤hL^∥z¯n−y¯n∥,foraconstantL^,andanyz¯n,y¯n(1≤n≤N).Proof.
First, we will
show by induction that(6.21)∥ki(zn)−ki(yn)∥≤hγi∥zn−yn∥foraconstantγi(1≤i≤v),where we recall that ∥⋅∥≡∥⋅∥∞, and that v is the number of function evaluations of the
ARK method. Also, here and throughout this proof, 1≤n≤N.
For i=1, we have(6.22)∥k1(zn)−k1(yn)∥=∥hf(zn)−hf(yn)∥≤hL∥zn−yn∥=hγ1∥zn−yn∥,where γ1=L.
Assuming inequality (6.21) is true for i=j,
we have(6.23)∥kj+1(zn)−kj+1(yn)∥=∥hf(zn+ajkj(zn))−hf(yn+ajkj(yn))∥≤hL∥zn−yn+aj(kj(zn)−kj(yn))∥≤hL[∥zn−yn∥+|aj|hγj∥zn−yn∥]=hL[1+hγj|aj|]∥zn−yn∥=hγj+1∥zn−yn∥,where γj+1=L[1+hγj|aj|].
Hence, inequality (6.21) holds for j+1 and, therefore, it holds for 1≤i≤v.
Similarly, we can show that(6.24)∥k−i(zn−1)−k−i(yn−1)∥≤hγi∥zn−1−yn−1∥for1≤i≤v.
Functionϕ can be written as(6.25)ϕ=c1k1−c−1k−1+∑i=2vci(ki−k−i)=∑i=1vciki−∑i=1vc−ik−i,wherec−i=cifor2≤i≤v=ϕ1−ϕ−1.Using inequality (6.21), we can
write(6.26)∥ϕ1(zn)−ϕ1(yn)∥=∥∑i=1vciki(zn)−∑i=1vciki(yn)∥≤∑i=1v|ci|∥ki(zn)−ki(yn)∥≤h[∑i=1v(|ci|γi)]∥zn−yn∥.Similarly, using inequality
(6.24), we can write(6.27)∥ϕ−1(zn−1)−ϕ−1(yn−1)∥≤h[∑i=1v(|c−i|γi)]∥zn−1−yn−1∥.Since ∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≡∥ϕ¯(z¯n)−ϕ¯(y¯n)∥∞=max(∥ϕ(z¯n)−ϕ(y¯n)∥,0),
we have(6.28)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥=∥ϕ(z¯n)−ϕ(y¯n)∥=∥ϕ1(zn)−ϕ−1(zn−1)−(ϕ1(yn)−ϕ−1(yn−1))∥≤∥ϕ1(zn)−ϕ1(yn)∥+∥ϕ−1(zn−1)−ϕ−1(yn−1)∥≤h[∑i=1v(|ci|γi)]∥zn−yn∥+h[∑i=1v(|c−i|γi)]∥zn−1−yn−1∥.Because ∥z¯n−y¯n∥=max(∥zn−yn∥,∥zn−1−yn−1∥),
we have(6.29)∥zn−yn∥≤∥z¯n−y¯n∥,∥zn−1−yn−1∥≤∥z¯n−y¯n∥.Therefore,(6.30)∥ϕ¯(z¯n)−ϕ¯(y¯n)∥≤h[∑i=1v|ci|γi+∑i=1v|c−i|γi]∥z¯n−y¯n∥=hL^∥z¯n−y¯n∥,where L^=∑i=1v(|ci|+|c−i|)γi.
Thus, Lemma 6.2 is proven.Theorem 6.3.
Suppose that an ARK method is used to solve the IVP
given in (1.3), and that f is a Lipschitz continuous function throughout
the integration span. If −1≤c−0<1,
the ARKp method is stable for all sufficiently small
step sizes, and
(6.31)∥zn−yn∥≤Emax(maxr=0,1∥zr−yr∥,max1≤m<n∥δm∥)for2≤n≤N,where(6.32)E=4we2L^w(tn−t1)max(tn−t1,1).Proof.
Subtracting (6.7) from (6.2), we have(6.33)z¯n+1−y¯n+1=Q(z¯n−y¯n)+ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)+hδ¯nfor1≤n≤N−1.We multiply this equation by W−1 to get(6.34)W−1(z¯n+1−y¯n+1)=W−1Q(z¯n−y¯n)+W−1[ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)]+hW−1δ¯n.Let(6.35)dn=W−1(z¯n−y¯n)sothatWdn=(z¯n−y¯n)for1≤n≤N,and we get(6.36)dn+1=W−1QWdn+W−1[ϕ¯(z¯n,f,h)−ϕ¯(y¯n,f,h)]+hW−1δ¯n.
Taking the infinity norm and using the Lipschitz
condition onϕ¯ (Lemma 6.2), we can write(6.37)∥dn+1∥≤∥W−1QW∥∥dn∥+hL^∥W−1∥∥z¯n−y¯n∥+h∥W−1∥∥δ¯n∥for1≤n≤N−1.Considering expressions (6.10) and
(6.12), we have(6.38)∥dn+1∥≤∥dn∥+hL^w∥z¯n−y¯n∥+hw∥δ¯n∥.Also, using expressions (6.35) and
(6.11), we have(6.39)∥z¯n−y¯n∥=∥Wdn∥≤∥W∥∥dn∥=2∥dn∥,so inequality (6.38)
becomes(6.40)∥dn+1∥≤(1+2hL^w)∥dn∥+hw∥δ¯n∥for1≤n≤N−1.
Relation (6.40) establishes a bound on ∥dn+1∥ based on ∥dn∥.
This is not very useful because we do not know how large ∥dn∥ is. Instead, we want to find a bound based on ∥d1∥.
We will now show by induction that this bound can be written as(6.41)∥dn∥≤e2L^w(tn−t1)∥d1∥+hw∑m=1n−1e2L^w(tn−tm+1)∥δ¯m∥for2≤n≤N.
In the following, we will use the
inequality(6.42)1≤1+t≤etfort≥0,which follows from the
expression et=1+t+(1/2)t2+⋯.
For n=2,
the right hand side of inequality (6.41) becomes(6.43)e2L^w(t2−t1)∥d1∥+hwe2L^w(t2−t2)∥δ¯1∥=e2L^wh∥d1∥+hw∥δ¯1∥.Inequality (6.40) can be written
as(6.44)∥dn∥≤(1+2hL^w)∥dn−1∥+hw∥δ¯n−1∥for2≤n≤N,and for n=2 as(6.45)∥d2∥≤(1+2hL^w)∥d1∥+hw∥δ¯1∥≤e2L^wh∥d1∥+hw∥δ¯1∥.Comparing expressions (6.45) and
(6.43), we can see that inequality (6.41) holds for n=2.
Assuming it holds for n=k,
we have from inequality (6.40) that(6.46)∥dk+1∥≤(1+2hL^w)∥dk∥+hw∥δ¯k∥for1≤k≤N−1≤e2hL^w(e2L^w(tk−t1)∥d1∥+hw∑m=1k−1e2L^w(tk−tm+1)∥δ¯m∥)+hw∥δ¯k∥.Since for any r, h+(tk−tr)=tk+1−tr,
we have that(6.47)∥dk+1∥≤e2L^w(tk+1−t1)∥d1∥+hw∑m=1k−1e2L^w(tk+1−tm+1)∥δ¯m∥+hwe2L^w(tk+1−tk+1)∥δ¯k∥=e2L^w(tk+1−t1)∥d1∥+hw∑m=1ke2L^w(tk+1−tm+1)∥δ¯m∥for1≤k≤N−1.Therefore, inequality (6.41) holds
for k+1 and by induction, it holds for 2≤n≤N.
We can simplify inequality (6.41) by
writing(6.48)∥dn∥≤e2L^w(tn−t1)∥d1∥+hw∑m=1n−1e2L^w(tn−tm+1)∥δ¯m∥for2≤n≤N≤e2L^w(tn−t1)∥d1∥+hwe2L^w(tn−t2)(n−1)max1≤m<n∥δ¯m∥≤e2L^w(tn−t1)∥d1∥+hwe2L^w(tn−t1)(n−1)max1≤m<n∥δ¯m∥=e2L^w(tn−t1)∥W−1(z¯1−y¯1)∥+we2L^w(tn−t1)(tn−t1)max1≤m<n∥δ¯m∥≤e2L^w(tn−t1)w[∥z¯1−y¯1∥+(tn−t1)max1≤m<n∥δ¯m∥]≤e2L^w(tn−t1)w[max(tn−t1,1)∥z¯1−y¯1∥+max(tn−t1,1)max1≤m<n∥δ¯m∥]≤2we2L^w(tn−t1)max(tn−t1,1)max(∥z¯1−y¯1∥,max1≤m<n∥δ¯m∥).Using the inequality (6.39), we
have(6.49)∥z¯n−y¯n∥≤4we2L^w(tn−t1)max(tn−t1,1)max(∥z¯1−y¯1∥,max1≤m<n∥δ¯m∥).Since ∥zn−yn∥≤∥z¯n−y¯n∥,
and ∥δ¯m∥=∥δm∥,
we have(6.50)∥zn−yn∥≤Emax(maxr=0,1∥zr−yr∥,max1≤m<n∥δm∥)for2≤n≤N,where(6.51)E=4we2L^w(tn−t1)max(tn−t1,1).Hence, stability of ARK methods
is established. Recall that the ARK methods require to be provided with the
initial condition y0 at the beginning of the first time step t0 and the approximate solution y1 at the end of the first time step t1.
The maxr=0,1∥zr−yr∥ term in expression (6.31) refers to the
perturbations in these values.Theorem 6.4.
Consider an ARKp method that is used to solve the IVP given by
(1.3), where f is a sufficiently smooth Lipschitz continuous
function throughout the integration span. If −1≤c−0<1 and the approximate solution y1 at time t1 is accurate to O(hq),
the ARKp method is convergent to order min(p,q),
and (6.52)∥y(tn)−yn∥≤Emax(O(hq),O(hp))for2≤n≤N,where(6.53)E=4we2L^w(tn−t1)max(tn−t1,1).In particular, ∥y(tn)−yn∥→0 as h→0.Proof.
Consider again
the approximate solution{yn}(2≤n≤N) and the perturbed solution {zn}(2≤n≤N) of the ARKp method given by (6.54)yn+1=c0yn−c−0yn−1+ϕ(yn,yn−1,f,h)for1≤n≤N−1,(6.55)zn+1=c0zn−c−0zn−1+ϕ(zn,zn−1,f,h)+hδnfor1≤n≤N−1,where z0=y0+Δ0, z1=y1+Δ1,
and zn's recipe is perturbed by hδn.
The truncation error hδn* is defined as the difference between the exact
and approximate solutions at the end of the current step if the method is
provided with the exact solution at the beginning of the current and previous
steps. Therefore, the exact solution y(tn) satisfies(6.56)y(tn+1)=c0y(tn)−c−0y(tn−1)+ϕ(y(tn),y(tn−1),f,h)+hδn*.We will now show that {y(tn)}(2≤n≤N) is a perturbed solution of the ARK method. For n=1, we have (6.57)y2=c0y1−c−0y0+ϕ(y1,y0,f,h),(6.58)z2=c0z1−c−0z0+ϕ(z1,z0,f,h)+hδ1,(6.59)y(t2)=c0y(t1)−c−0y(t0)+ϕ(y(t1),y(t0),f,h)+hδ1*.We set(6.60)z0=y(t0)=y0,z1=y(t1)=y1+Δ1*,δ1=δ1*,and from (6.58)
and (6.59), we have z2=y(t2).
Similarly, we set δn=δn* for 1≤n≤N−1,
and from (6.55) and (6.56), we get zn=y(tn) for 2≤n≤N.
Since we have made the same assumptions as Theorem6.3,
relations (6.31) and (6.32) hold, and we have(6.61)∥y(tn)−yn∥≤Emax(∥y(t1)−y1∥,max1≤m<n∥δm*∥)for2≤n≤N.For ARKp, δm*=O(hp).
The approximate solution y1 at time t1 which is required by the ARKp method is supplied by a one-step method such
as RKp.
If y1 is accurate of order q,
that is, ∥y(t1)−y1∥=O(hq),
the above inequality shows that the ARKp method is convergent of order min(p,q),
that is,(6.62)∥y(tn)−yn∥≤Emax(O(hq),O(hp))for2≤n≤N.In particular, ∥y(tn)−yn∥→0 as h→0,
and Theorem 6.4 is proven.


## 6.2. Large Step Sizes

To produce the
stability region of ARK methods for finite step sizes, we look at the
approximate solution of the linear one-dimensional problemy′=λy,
where λ is a complex number. After some simplification
using the order equations, the ARK method's approximate solution of this problem
becomes(6.63)[yn+1yn]=[1+c−0+12(3−c−0)(λh)+S−c−0−12(1+c−0)(λh)−S10][ynyn−1],where for ARK3,(6.64)S=112(5+c−0)(λh)2,for ARK4,(6.65)S=112(5+c−0)(λh)2+16(λh)3,for ARK4-4,(6.66)S=112(5+c−0)(λh)2+16(λh)3+c4a1a2a3(λh)4,and for ARK5,(6.67)S=112(5+c−0)(λh)2+16(λh)3+1720(31−c−0)(λh)4+c5a1a2a3a4(λh)5.For stability, we require that the modulus of the eigenvalues
of the matrix in (6.63) be smaller than 1. Figures
15, 16,
17, and 18 show
the stability regions of ARK and RK methods produced in this way. The curves
represent the points at which the computed maximum modulus of the eigenvalues
equals 1. The region of stability is the area enclosed by the curves and the
real axis. It is clear from (6.63) that the stability region of ARK3 and ARK4
depends only on c−0,
but the stability of ARK4-4 and ARK5 depends also on the products c4a1a2a3 and c5a1a2a3a4, respectively.Figure 15
Stability region of ARK3 sets of Table1 and RK3 set
of Table 5.Figure 16
Stability region of ARK4
sets of Table2 and RK4 set of Table
5.Figure 17
Stability region of ARK4-4 sets of Table3 and RK4 set
of Table 5.Figure 18
Stability region of ARK5
sets of Table4 and RK5 set
of Table 5.


## 7. Conclusion

In this paper,
we have developed a simple set of constant step size, explicit, accelerated
Runge-Kutta (ARK) methods for the numerical integration of initial value problems.
They rely on the general form posited in (2.5). The simplicity of this form
facilitates their detailed treatment, an aspect missing in previous work on
similar multistep methods. Specifically, we have provided a study of their
motivation, derivation, optimization, accuracy, speedup, stability, and
convergence.Our main conclusions are as follows.(1)
The ARK methods developed in this paper are
simple to analyze and implement and they lead to an elegant set of order
equations. Three accuracy-optimized parameter sets for ARK3, ARK4, ARK4-4, and
ARK5 have been provided. The accuracy of the methods obtained from these sets
is confirmed by numerically solving a standard set of seven initial value
problems.(2)
The ARK methods are more computationally
efficient than RK methods at providing the same order of local accuracy because
they reuse function evaluations from previous steps of integration. ARK3, ARK4,
and ARK5 require one less function evaluation per time step of integration than
RK3, RK4, and RK5, respectively.(3)
Numerical examples show that the
accuracy-optimized ARK3, ARK4, ARK4-4, and ARK5 methods presented here are
superior to RK methods that computationally cost the same in terms of the
number of function evaluations. Also, the ARK methods are comparable to RK
methods that have the same order of local accuracy.(4)
ARK methods require to be initially started by
a one-step method of equal or higher order such as an RK method; however, this
extra computation occurs only at the first step of integration and is
insignificant compared to the computational savings of ARK methods over the
subsequent steps.(5)
In solving the standard initial value problems
considered here, the ARK3, ARK4, and ARK5 methods exhibit minimum speedups of
19%, 17%, and 15% compared to RK3, RK4, and RK5 methods, respectively, on our
computer. However, the theoretical speedups of the ARK methods are
approximately 33%, 25%, and 17%, respectively, based on the smaller number of
function evaluations required. This reduction in speedup is found to be caused
by the higher computational overheads in the ARK methods, when compared to the
RK methods.(6)
Convergence and stability for small step sizes
are proven for ARK methods with some conditions. The stability regions of ARK
methods for large step sizes are shown to be generally smaller than, but
comparable to, those of RK methods.



---

*Source: 790619-2008-09-16.xml*