remarkRemark hypothesisHypothesis claimClaim Walking into the complex plane to
‘order’ better time integratorsJ. D. George, S. Y, Jung, and N. M. Mangan
ex_supplement

# Walking into the complex plane to ‘order’ better time integrators

Jithin D. George Department of Engineering Sciences and Applied Mathematics,
Northwestern University, Evanston, IL.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Samuel
Y. Jung 11footnotemark: 1 Niall M. Mangan11footnotemark: 1

###### Abstract

Most numerical methods for time integration use real time steps. Complex time
steps provide an additional degree of freedom, as we can select the magnitude
of the step in both the real and imaginary directions. By time stepping along
specific paths in the complex plane, integrators can gain higher orders of
accuracy or achieve expanded stability regions. Complex time stepping also
allows us to break the Runge-Kutta order barrier, enabling 5th order accuracy
using only five function evaluations. We show how to derive these paths for
explicit and implicit methods, discuss computational costs and storage
benefits, and demonstrate clear advantages for complex-valued systems like the
Schrodinger equation.

###### keywords:

complex time steps, complex time integrators, high order methods, absolute
stability, order barrier

30-08, 65E05, 65L05, 65L04, 65M12, 65M20

##  1 Introduction

The real line is a tiny sliver in the whole of the complex plane. Wandering
into the complex plane has often been employed to solve real problems. Two
direct examples are the solution of complicated real-valued integrals via
contour integrals [1] and the usage of the intuition gained from poles in the
complex plane in control-theory applications [2]. The complex plane provides
wonderful insights to real problems, allowing mathematicians to impact
practical applications. One of the motivations for this paper was a simple
numerical application highlighted in a 2018 SIAM News article [3] where Nick
Higham describes how taking a numerical complex derivative neatly gets rid of
round-off error. However, most methods for numerically solving a differential
equation from time $a$ to time $b$, construct a path utilizing only the real
line. We show that we can use the complex plane to our advantage when it comes
to numerical time integration by taking time steps in the complex plane.

This paper is not the first to talk about complex time steps. Early work on
time stepping in the complex plane focused on systematically avoiding
singularities during numerical integration of singular differential equations
[4] and as techniques to achieve superconvergence for Runge-Kutta methods [5].
A use of complex time steps most similar to our work was in [6], where
Chambers explores operator splitting in symplectic integrators using complex
substeps. Taking real substeps often results in negative substeps and substeps
larger than the original step causing the symplectic integrator go back and
forth along the direction of integration. By using complex substeps, Chambers
showed how to avoid negative and unstable substeps, significantly reducing the
truncation error. Chambers’s work inspired others to consider complex time
steps as part of operator splitting methods for various parabolic and
nonlinear evolution problems [7, 8, 9, 10, 11]. Concurrent with our work, a
very recent paper [12] in this area uses complex coefficients for high order
composition methods. Another recent paper develops time integration schemes
using nodes in the complex plane based on Cauchy’s integral formula and using
the complex nodes to facilitate parallel BDF formulation [13]. Rather than
focusing on a specific application or numerical scheme as has been done
previously, we ask fundamental questions about the utility of numerical time-
stepping in the complex plane, including: ‘What paths lead to the least
error?’ and ‘Are there general strategies and trade offs for complex time-
stepping?’. Here we show that stepping into the complex plane provides an
extra dimension which can be exploited to improve accuracy, design expanded
stability regions and even break through the order barrier [14, 15]. We
compare against existing Runge-Kutta methods and demonstrate improvement for
particular proof of concept cases. Finally we analyze increased operation
counts due to complex operations, potential storage benefits, extensions to
stiff solvers, and strong stability preserving properties. The code to
reproduce the numerical experiments in this paper is available on GitHub
111https://github.com/Dirivian/complex_time_integrators .

##  2 Taking complex time steps with the Euler method

The idea of a complex time step can be unusual and unintuitive. To make it
accessible, we first illustrate simple paths with complex time steps and
demonstrate that applying the forward Euler method along the complex paths
improves accuracy compared to forward Euler with real steps. The paths we
choose in the complex plane must intersect the real line at the time points
where we desire a solution. Our primary question is ‘what is the optimal
n-step path in the complex plane over which to integrate between two real time
points?’

To illustrate the simplest possible case, suppose we want to solve the
equation $\dot{y}=\lambda y$ with the initial condition $y_{0}$ at time $t=0$.
To find the solution $y(\Delta t)$, we can either take two real time steps of
length $\frac{\Delta t}{2}$ or take a complex time step of length
$\frac{\Delta t}{2}+i\frac{\Delta t}{2}$ followed by a complex time step of
length $\frac{\Delta t}{2}-i\frac{\Delta t}{2}$ returning back to the real
line. These two paths are shown in Fig. 1 (top). We numerically integrate the
differential equations $\dot{y}=y$ along these paths as shown in Fig. 1
(bottom). By eye, the complex time path clearly does much better than the real
time path in this linear case. The improvement with a complex path can hold
for nonlinear ODEs, second order ODEs, non-autonomous ODEs (Fig. 2) and
nonlinear PDEs, such as the viscous Burgers equation (Fig. 3). In section 3,
we describe the generalizations needed to solve nonlinear equations, but first
we provide theoretical insight into the simple linear case.

Figure 1: Taking complex time steps that return to the real line results in
increased accuracy. Figure 2: Comparison of real-valued classic Euler method
and 2-step complex-valued Euler method. The number of real timepoints for the
classic scheme are roughly double, $2n-1$, the number for the 2-step complex
scheme so that the net number of timesteps and function evaluations are the
same for both schemes. Figure 3: The complex time steps reduce the error for
the nonlinear viscous Burgers equation.

Why does time stepping on the complex path improve accuracy of the numerical
solution so significantly? To answer that, we analyze the linear problem
$\dot{y}=y$ with initial condition $y(t_{0})=y_{0}$. Consider time stepping
from the initial condition $y_{0}$ at time $t_{0}$ to time $t+\Delta t$ using
2 real time steps. The discrete values of $y(t)$ are given by the following at
$t=t_{0},\;t_{0}+\frac{\Delta t}{2},$ and $t_{0}+\Delta t$:

$\displaystyle y(t_{0})$ $\displaystyle=y_{0}$ $\displaystyle
y(t_{0}+\frac{\Delta t}{2})$ $\displaystyle=y_{0}\bigg{(}1+\frac{\Delta
t}{2}\bigg{)}$ (1) $\displaystyle y(t_{0}+\Delta t)$
$\displaystyle=y_{0}\bigg{(}1+\frac{\Delta t}{2}\bigg{)}\bigg{(}1+\frac{\Delta
t}{2}\bigg{)}=y_{0}\bigg{(}1+\Delta t+\frac{\Delta t^{2}}{4}\bigg{)}.$

Now consider timestepping from $y_{0}$ at time $t_{0}$ to time $t+\Delta t$
using the two complex time steps.

$\displaystyle y(t_{0})$ $\displaystyle=y_{0}$ $\displaystyle
y\Big{(}t_{0}+\frac{\Delta t}{2}+i\frac{\Delta t}{2}\Big{)}$
$\displaystyle=y_{0}\bigg{(}1+\frac{\Delta t}{2}+i\frac{\Delta t}{2}\bigg{)}$
(2) $\displaystyle y(t_{0}+\Delta t)$
$\displaystyle=y_{0}\bigg{(}1+\frac{\Delta t}{2}+i\frac{\Delta
t}{2}\bigg{)}\bigg{(}1+\frac{\Delta t}{2}-i\frac{\Delta
t}{2}\bigg{)}=y_{0}\bigg{(}1+\Delta t+\frac{\Delta t^{2}}{2}\bigg{)}$

The exact solution to the problem is given by the Taylor expansion

(3) $\displaystyle y(t_{0}+\Delta t)=y_{0}e^{\Delta t}=y_{0}\Big{(}1+\Delta
t+\frac{\Delta t^{2}}{2}+\ldots\Big{)}.$

The secret behind the higher accuracy of the complex 2-step method is that the
polynomial obtained through the numerical approximation matches the Taylor
series of the exact solution to the second order term, making forward Euler
method gain second order accuracy. This is one of two unique 2-step paths in
the complex plane which has this property. The other unique path uses the
other permutation of the complex time steps $(\frac{\Delta t}{2}-i\frac{\Delta
t}{2},\frac{\Delta t}{2}+i\frac{\Delta t}{2})$ instead of $(\frac{\Delta
t}{2}+i\frac{\Delta t}{2},\frac{\Delta t}{2}-i\frac{\Delta t}{2})$. Both
permutations of these steps will lead to the same polynomial with coefficients
that match the Taylor expansion to 2nd order.

By stepping into the complex plane, we gain an extra degree of freedom in the
coefficients that allows us to design time-steps with improved accuracy. This
extra dimension can be used to construct higher order paths for linear and
nonlinear differential equations (Section 3), design improved stability
(Section 6) and even break through the order barrier for Runge-Kutta methods
(Section 5). The paths can be combined with other methods, such as the
implicit midpoint method (Runge-Kutta-2) as demonstrated in Section 8. In the
next section, we demonstrate how to extend the linear example to construct
paths with higher order accuracy.

##  3 Constructing a higher order time integrator

Is it possible then to construct 3-step path on which forward Euler would be
3rd order accurate? Would it be nth order accurate on an n-step complex path?
To build intuition for the simplest higher order construction, we continue
analyzing the linear problem in Subsection 3.1, and then prove that achieving
higher order requires complex time-stepping. In Subsection 3.2, we extend the
construction of complex time-steppers to the general nonlinear problem.

###  3.1 Higher order integrators for linear differential equations

In this subsection, we show how to construct a third-order time integrator for
the ode $\dot{y}=\lambda y$ by taking 3 forward Euler time steps of variable
sizes in the complex plane that return to the real line at every desired time
point.

If we take the three complex steps $w_{1}\Delta t,w_{2}\Delta t$ and
$w_{3}\Delta t$ to timestep from $y_{0}$ at $t_{0}$ to $t_{0}+\Delta t$, the
approximate numerical solution $\bar{y}$ at time $t+\Delta t$ is given by

(4) $\displaystyle\bar{y}(t_{0}+\Delta t)=y_{0}(1+w_{1}\lambda\Delta
t)(1+w_{2}\lambda\Delta t)(1+w_{3}\lambda\Delta t)$

where the $w_{i}$’s are complex variables. Collecting $\lambda\Delta t$ terms
results in

(5) $\displaystyle\bar{y}(t_{0}+\Delta t)$
$\displaystyle=y_{0}(1+(w_{1}+w_{2}+w_{3})\lambda\Delta
t+(w_{1}w_{2}+w_{2}w_{3}+w_{1}w_{3})\lambda^{2}\Delta
t^{2}+w_{1}w_{2}w_{3}\lambda^{3}\Delta t^{3}).$

The exact solution is given by

(6) $\displaystyle y(t_{0}+\Delta t)=y_{0}e^{\lambda\Delta
t}=y_{0}(1+\lambda\Delta t+\frac{1}{2}\lambda^{2}\Delta
t^{2}+\frac{1}{6}\lambda^{3}\Delta t^{3}+\ldots)$

For the numerical approximation to be third-order accurate, the coefficients
of $y_{0}(\lambda\Delta t)^{n}$ for $n=1,2,3$ in equations (5) and (6) need to
match. This results in the following order conditions:

$\displaystyle w_{1}+w_{2}+w_{3}$ $\displaystyle=1$ $\displaystyle
w_{1}w_{2}+w_{2}w_{3}+w_{1}w_{3}$ $\displaystyle=\frac{1}{2}$ (7)
$\displaystyle w_{1}w_{2}w_{3}$ $\displaystyle=\frac{1}{6}.$

One of the solutions to the above nonlinear system is

(8)
$\displaystyle(w_{1},w_{2},w_{3})=(0.186731+0.480774i,0.626538,0.186731-0.480774i).$

However, as the nonlinear system 3.1 is symmetric about $w_{1},w_{2}$ and
$w_{3}$, all 6 permutations of
$(0.186731+0.480774i,0.626538,0.186731-0.480774i)$ result in unique 3-step 3rd
order complex time stepping paths. In Fig 4, all six 3-step 3rd order paths
are shown along with the two 2-step 2nd order paths and the single step
forward Euler path.

For a general n-step nth order path, the order conditions would be

$\displaystyle\sum_{i}w_{i}$ $\displaystyle=1$ $\displaystyle\sum_{i\neq
j}w_{i}w_{j}$ $\displaystyle=\frac{1}{2}$ $\displaystyle\sum_{i\neq j\neq
k}w_{i}w_{j}w_{k}$ $\displaystyle=\frac{1}{6}$ $\displaystyle\vdots$ (9)
$\displaystyle\prod_{i=1}^{n}w_{i}$ $\displaystyle=\frac{1}{n!}.$

For higher orders, the solutions take on a complicated form and solving the
nonlinear system becomes a tough symbolic task. However, we only need to find
one solution for high-order approximations and such these solution can be
found numerically and be stored in a library for general use.

Figure 4: All possible first, second, and third-order complex Euler methods
for linear differential equations. Each order requires an additional step and
adds an additional constraint on the coefficients.

These results beg the theoretical question: is there any way to take steps of
varying size on the real line and achieve a higher-order time integrator? or
is it only possible to do so by going into the complex plane?

######  Theorem 3.1.

It is impossible to construct a time integrator with order greater than or
equal to 2 using real-valued Euler time steps of variable sizes.

######  Proof 3.2.

We prove the above theorem by showing that it is impossible to construct a
second-order time integrator using forward Euler steps of variable size on the
real line.

Suppose we use $n$ time steps of variable size given by
$w_{1},w_{2},\ldots,w_{n}$. The approximate solution to $\dot{y}=\lambda y$
given by our time integrator $\bar{y}$ would be

$\bar{y}(\Delta t)=(1+w_{1}\lambda\Delta t)(1+w_{2}\Delta\lambda
t)\ldots(1+w_{n}\Delta\lambda t)$

where we assume the first order condition is satisfied,

$w_{1}+w_{2}+\ldots+w_{n}=1.$

Gathering together different order $\lambda\Delta t$ terms,

$\displaystyle\bar{y}(\Delta t)$
$\displaystyle=1+(w_{1}+w_{2}+\ldots+w_{n})\lambda\Delta
t+(w_{1}w_{2}+w_{1}w_{3}+\ldots w_{n-1}w_{n})\lambda^{2}\Delta t^{2}+\ldots$
$\displaystyle=1+\lambda\Delta t+(w_{1}w_{2}+w_{1}w_{3}+\ldots
w_{n-1}w_{n})\lambda^{2}\Delta t^{2}+\ldots$

where we have substituted in our first order condition. For the time stepper
to be second order, we need the coefficient on $\lambda^{2}\Delta t^{2}$to be

$w_{1}w_{2}+w_{1}w_{3}+\ldots+w_{2}w_{3}+\ldots+w_{1}w_{n}+w_{2}w_{n}+\ldots+w_{n-1}w_{n}=\frac{1}{2}.$

Factoring out a $w_{n}$ results in

$w_{1}w_{2}+w_{1}w_{3}+\ldots+w_{2}w_{3}+\ldots+(w_{1}+w_{2}+\ldots+w_{n-1})w_{n}=\frac{1}{2}.$

Solving the first order condition for $w_{n}$ and incorporating produces

$w_{1}w_{2}+w_{1}w_{3}+\ldots+w_{2}w_{3}+\ldots+(w_{1}+w_{2}+\ldots+w_{n-1})(1-w_{1}-w_{2}\ldots-
w_{n-1})=\frac{1}{2},$

which rearranges to

$w_{1}+w_{2}+\ldots+w_{n-1}-w_{1}^{2}-w_{2}^{2}\ldots-
w_{n-1}^{2}-w_{1}w_{2}-w_{1}w_{3}\ldots-w_{n-1}w_{n-2}=\frac{1}{2}.$

Since all the $w$s are real, we can express them as

$w_{j}=k_{j}w_{1}\quad j=1,\ldots,n-1$

where all the $k$s are real. Then, we get

$(k_{1}+k_{2}+\ldots+k_{n-1})w_{1}-w_{1}^{2}(k_{1}^{2}+k_{2}^{2}\ldots+k_{n-1}^{2}+k_{1}k_{2}+k_{1}k_{3}\ldots+k_{n-1}k_{n-2})=\frac{1}{2}.$

Rearranging in terms of powers of $w_{1}$ one has

$w_{1}^{2}(k_{1}^{2}+k_{2}^{2}\ldots+k_{n-1}^{2}+k_{1}k_{2}+k_{1}k_{3}\ldots+k_{n-1}k_{n-2})-(k_{1}+k_{2}+\ldots+k_{n-1})w_{1}+\frac{1}{2}=0.$

This is quadratic equation to solve for $w_{1}$ in the form of
$ax^{2}+bx+c=0$. For $w_{1}$ to only have real solutions, we need a positive
discriminant $\Delta=b^{2}-4ac\geq 0$. In our case, the discriminant,
$\Delta$, is

$\displaystyle\Delta$
$\displaystyle=(k_{1}+k_{2}+\ldots+k_{n-1})^{2}-4(k_{1}^{2}+k_{2}^{2}\ldots+k_{n-1}^{2}+k_{1}k_{2}+k_{1}k_{3}\ldots+k_{n-1}k_{n-2})\frac{1}{2}$
$\displaystyle=-(k_{1}^{2}+k_{2}^{2}+\ldots+k_{n-1}^{2})$

Since the discriminant is always negative, $w_{1}$ cannot be a purely real
time step and must to be complex to satisfy our second order condition.
Therefore, not only can we achieve higher order integrators through complex
time stepping, but we must step into the complex plane to achieve higher
order.

###  3.2 Higher order integrators for general nonlinear non-autonomous
differential equations

In the previous subsection, we showed how to construct a 3-step 3rd order
complex path on which to time integrate using forward Euler steps. We can
experimentally demonstrate that the order of accuracy is obeyed by every
linear ordinary differential equation and partial differential equation that
we have tested it on. However, the derivation assumed a linear problem and
does not guarantee the same order of accuracy for nonlinear problems. Indeed,
when we tested the complex time-stepping schemes on nonlinear and non-
autonomous equation, the 3 step paths achieved only second order accuracy. To
achieve higher order accuracy for nonlinear, non-autonomous equations, we
generalize the procedure we described so far to obtain 3rd order accuracy.

A general nonlinear ode is given by

(10) $\displaystyle\dot{y}=f(t,y),\quad y(t_{0})=y_{0}.$

Let us assume the function $f$ is sufficiently analytic near the region of
integration i.e it does not have any singularities close to the path of
integration. At time $t=t_{0}+\Delta t$, we have

(11) $\displaystyle y(t_{0}+\Delta t)=y_{0}+\Delta
tf(t_{0},y_{0})+\frac{1}{2}\Delta t^{2}\dot{f}(t_{0},y_{0})+\frac{1}{3!}\Delta
t^{3}\ddot{f}(t_{0},y_{0})+\ldots$

where

(12)
$\displaystyle\dot{f}(t,y)=\frac{d}{dt}f(t,y)\and[\Bigg{]}{t=t_{0},y=y_{0}}=f{\left(t_{0},y_{0}\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}\and[\Bigg{]}{t=t_{0},y=y_{0}}+\frac{\partial}{\partial
t}f{\left(t,y\right)\and[\Bigg{]}{t=t_{0},y=y_{0}}}.$

Expanding equation 11, we get

$\displaystyle y(t_{0}+\Delta t)=$ $\displaystyle y_{0}+\Delta
tf(t_{0},y_{0})+\frac{1}{2}\Delta
t^{2}\bigg{(}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}+\frac{\partial}{\partial
t}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{0},y=y_{0}}$
$\displaystyle+\frac{1}{3!}\Delta
t^{3}\bigg{(}f^{2}{\left(t,y\right)}\frac{\partial^{2}}{\partial
y^{2}}f{\left(t,y\right)}+\left(\frac{\partial}{\partial
y}f{\left(t,y\right)}\right)^{2}f{\left(t,y\right)}$ (13)
$\displaystyle+2f{\left(t,y\right)}\frac{\partial^{2}}{\partial y\partial
t}f{\left(t,y\right)}+\frac{\partial}{\partial
t}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}+\frac{\partial^{2}}{\partial
t^{2}}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{0},y=y_{0}}$
$\displaystyle+\ldots$

Now, consider a complex time stepper that takes 3 complex time steps
$w_{1}\Delta t,w_{2}\Delta t,w_{3}\Delta t$ starting from $y_{0}$ to go to
$y_{1},y_{2},y_{3}$. These steps are are given by:

(14) $\displaystyle y_{0}$ $\displaystyle=y(t_{0})$ (15) $\displaystyle y_{1}$
$\displaystyle=y_{0}+w_{1}\Delta tf(t_{0},y_{0})$ (16) $\displaystyle y_{2}$
$\displaystyle=y_{1}+w_{2}\Delta tf(t_{0}+w_{1}\Delta t,y_{1})$ (17)
$\displaystyle y_{3}$ $\displaystyle=y_{2}+w_{3}\Delta tf(t_{0}+w_{1}\Delta
t+w_{2}\Delta t,y_{2}).$

One can expand $f(t_{n},y_{n})$ around $f(t_{0},y_{0})$ using the multivariate
Taylor series. Doing so,

$\displaystyle y_{3}=$ $\displaystyle y_{0}+\Delta
t(w_{1}+w_{2}+w_{3})f{\left(t_{0},y_{0}\right)}$ $\displaystyle+\Delta
t^{2}(w_{1}w_{2}+w_{1}w_{3}+w_{2}w_{3})\bigg{(}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}+\frac{\partial}{\partial
t}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{0},y=y_{0}}$
$\displaystyle+\Delta
t^{3}\bigg{(}\frac{1}{2}(w_{1}^{2}w_{2}+w_{1}^{2}w_{2}+2w_{1}w_{2}w_{3}+w_{2}^{2}w_{3})f^{2}{\left(t,y\right)}\frac{\partial^{2}}{\partial
y^{2}}f{\left(t,y\right)}$
$\displaystyle+w_{1}w_{2}w_{3}f{\left(t,y\right)}\left(\frac{\partial}{\partial
y}f{\left(t,y\right)}\right)^{2}+w_{1}w_{2}w_{3}\frac{\partial}{\partial
t}f{\left(t,y\right)}\frac{\partial}{\partial y}f{\left(t,y\right)}$
$\displaystyle+(w_{1}^{2}w_{2}+w_{1}^{2}w_{2}+2w_{1}w_{2}w_{3}+w_{2}^{2}w_{3})f{\left(t,y\right)}\frac{\partial^{2}}{\partial
y\partial t}f{\left(t,y\right)}$ (18)
$\displaystyle+\frac{1}{2}(w_{1}^{2}w_{2}+w_{1}^{2}w_{2}+2w_{1}w_{2}w_{3}+w_{2}^{2}w_{3})\frac{\partial^{2}}{\partial
t^{2}}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{0},y=y_{0}}+O\left(\Delta
t^{4}\right).$

For this three step complex integrator to have 3rd order accuracy, the
coefficients of the partial derivatives in Eq 3.2 should match that of Eq 3.2.
Equating coefficients, we get this system of nonlinear equations.

$\displaystyle w_{1}+w_{2}+w_{3}=1$ $\displaystyle
w_{1}w_{2}+w_{1}w_{3}+w_{2}w_{3}=\frac{1}{2}$ $\displaystyle
w_{1}w_{2}w_{3}=\frac{1}{6}$ (19) $\displaystyle
w_{1}^{2}w_{2}+w_{1}^{2}w_{3}+2w_{1}w_{2}w_{3}+w_{2}^{2}w_{3}=\frac{2}{3}$

Unfortunately, this system is overconstrained and does not have a solution,
implying that we may need more than 3 steps for 3rd order accuracy. However,
for systems with only real solutions we can take further advantage of the
degree of freedom in the imaginary part of the step and find many 3-step paths
with 3rd order accuracy. If we allow the truncation error at the order of
$\Delta t^{3}$ to be purely imaginary, we can effectively remove the error at
the end of the 3 step sequence by taking only the real part of the numerical
solution at the end of the complex path. Let $y$ be the true solution and
$\bar{y}$ be the approximate numerical solution. Then, with our setup, we
would have

(20) $\displaystyle y(t+\Delta t)=\bar{y}(t+\Delta t)+ik\Delta t^{3}+O(\Delta
t^{4})$

where $k$ is real. Taking the real part on both sides, we get

(21) $\displaystyle\textbf{Re}(y(t+\Delta t))=\textbf{Re}(\bar{y}(t+\Delta
t))+O(\Delta t^{4})$

If the true solution to our differential equation and its initial condition
are known to be real-valued at real time points, then we have $y(t+\Delta
t)=\textbf{Re}(y(t+\Delta t))$ and taking the real part of the solution at the
end of complex path results in 3rd order accurate numerical solution. This
formulation, allowing the error to be purely imaginary at $O(\Delta t^{3})$
has fewer constraints on the coefficients, which must now satisfy the
following relaxed conditions:

$\displaystyle w_{1}+w_{2}+w_{3}=1$ $\displaystyle
w_{1}w_{2}+w_{1}w_{3}+w_{2}w_{3}=\frac{1}{2}$
$\displaystyle\textbf{Re}(w_{1}w_{2}w_{3})=\frac{1}{6}$ (22)
$\displaystyle\textbf{Re}(w_{1}^{2}w_{2}+w_{1}^{2}w_{3}+2w_{1}w_{2}w_{3}+w_{2}^{2}w_{3})=\frac{2}{3}.$

There are many paths satisfying these equations. Of the 3rd order paths in Fig
4 that were found to satisfy the linear order conditions, the top and bottom
(cyan and brown) paths also satisfy the 3rd order nonlinear order conditions.
We ran convergence tests for various nonlinear problems listed in Section 3.3,
and demonstrated 3rd order accuracy with those paths (Fig 5).

In this subsection, we have derived paths in the complex plane that satisfy
the order conditions for scalar non-autonomous differential equations. A
natural follow-up question is whether these paths return the same order for
vector valued non-autonomous differential equations, expanding our result to
systems of differential equations and partial differential equations. Butcher
showed in [15] that the order conditions for scalar non-autonomous
differential equations are the same as those for vector valued differential
equations up to order 4. As the paths derived here are of order 3 or less,
they hold for systems of differential equations and partial differential
equations as seen in the next subsection. To construct a fifth order or higher
complex path, we would need extract the necessary order conditions using
Butcher’s rooted tree approach [15].

###  3.3 Numerical Results

We have tested our complex time-stepping schemes on a variety of differential
equations to show that they satisfy the order conditions. The equations are
listed in Table 1 along with a unique label corresponding to their legend on
the convergence plot in Fig.5.

Figure 5: Convergence of various ODEs and PDEs (Sec. 3.2) for 1-step (forward
Euler), 2-step and 3-step complex time steps. The path/method used is
displayed above the corresponding convergence plot. The orders of accuracy are
as expected. Table 1: List of differential equations 1\. Linear ODEs • The
Dahlquist test problem (‘linear’) $\displaystyle\dot{y}=\lambda y,\quad
y(0)=1$ • Simple harmonic oscillator (‘shm’) $\displaystyle\ddot{y}=-y,\quad
y(0)=1,\quad y^{\prime}(0)=0$ 1\. Real unequal eigenvalues of same sign
---
2\. Nonlinear ODEs • Nonlinear autonomous ode 1(‘square’)
$\displaystyle\dot{y}=-y^{2},\quad y(0)=1$ • Nonlinear autonomous ode 2
(‘exp’) $\displaystyle\dot{y}=-e^{y},\quad y(0)=1$ • Nonlinear sine(‘nlsin’)
$\displaystyle\dot{y}=4y\sin{(t)}^{3}\cos{(t)},\quad y(0)=1$ • The Van Der Pol
Oscillator (vdp”)
$\displaystyle\dot{y_{1}}=y_{2},\dot{y_{2}}=\mu(1-y_{1}^{2})y_{2}-y_{1},\quad(y_{1},y_{2})(0)=(2,0)$
2\. Real unequal eigenvalues of opposite sign
3\. Linear PDE : The advection equation (‘wave’) This is a linear PDE given by
(23) $\displaystyle u_{t}=cu_{x},\quad x\in(0,2\pi),\quad t\in(0,1)$ This PDE
is solved with periodic boundary conditions and the initial condition
$u(x,t=0)=e^{7(x-\pi)^{2}}$. The exact solution is
$u(x,t)=e^{7(x-t-\pi)^{2}}$. The parameter $c$ is set to 1. The Fourier
spectral method is used for the spatial differencing with 70 modes.
3\. Nonlinear PDE : Viscous Burgers equation (‘burgers’) This is a nonlinear
PDE given by (24) $\displaystyle u_{t}+uu_{x}=\nu u_{xx},\quad
x\in(0,2\pi),\quad t\in(0,2)$ This PDE is solved with periodic boundary
conditions and the initial condition $u(x,t=0)=2\mu\sin x/(1.5+\cos x)$. The
exact solution is $u(x,t)=2\mu e^{-\mu t}\sin x/(1.5+e^{-\mu t}\cos x)$. The
parameter $\nu$ is set to 0.1. The Fourier spectral method is used for the
spatial differencing with 70 modes.

The error for all ODEs and PDEs apart from the Van Der Pol oscillator is
calculated by comparing the numerical solution to the exact solution. For the
Van Der Pol oscillator, the numerical solution using complex paths is compared
to the numerical solution obtained by using the ‘RK44’ 4th order method from
Nodepy [16] with $\Delta t=10^{-6}$.

###  3.4 Higher order integrators tailored to specific problems

For specific problems, it is possible to derive complex time steppers that are
uniquely accurate with fewer steps. For example, we can derive a 2-step
complex Euler method with 3rd order accuracy for the following nonlinear
problem,

(25) $\displaystyle y^{\prime}=-y^{2}\quad y(t_{0})=y_{0}.$

The true solution at time $t_{0}+\Delta t$ can be approximated to third order
accuracy using Eq (3.2) to obtain

(26) $\displaystyle y(t_{0}\Delta t)=\frac{y_{0}}{1+y_{0}\Delta
t}=y_{0}-y_{0}^{2}\Delta t+y_{0}^{3}\Delta t^{2}-y_{0}^{4}\Delta
t^{3}+O(\Delta t^{4}).$

Note that one does not need the true solution to obtain Eq (26) as this
expansion can be obtained from the right hand side of the differential
equation.

Now, we apply the following two-step complex method to integrate from $y_{0}$
at time $t=t_{0}$ to $y_{2}$ at time $t=t_{0}+\Delta t$ using forward Euler
steps.

$\displaystyle y_{0}$ $\displaystyle=y_{0}$ $\displaystyle y_{1}$
$\displaystyle=y_{0}+w_{1}\Delta ty_{0}^{2}$ (27) $\displaystyle y_{2}$
$\displaystyle=y_{1}+w_{2}\Delta ty_{1}^{2}$

where $w_{1}\Delta t,w_{2}\Delta t$ are complex time steps. Writing $y_{2}$ in
terms of the initial condition, we get

(28) $\displaystyle y_{2}=y_{0}+(w_{1}+w_{2})y_{0}^{2}\Delta
t+2w_{1}w_{2}y_{0}^{3}\Delta t^{2}+w_{1}^{2}w_{2}y_{0}^{4}\Delta t^{3}.$

This 2 step complex Euler method would have 3rd order accuracy if the
coefficients in Eq. 26 and 28 match. However, there are no values
$w_{1},w_{2}$ that satisfy these conditions. So, we assume the case in which
the true solution to (25) is real at real time points (valid if the initial
condition is real). For this case, we can let the error be purely imaginary
beyond order 1 and use the trick of taking the real part of the numerically
calculated values every time we complete a complex path and return to the real
line. This relaxed constraint produces the following a nonlinear system of
equations that can be solved,

$\displaystyle w_{1}+w_{2}=1$
$\displaystyle\textbf{Re}(w_{1}w_{2})=\frac{1}{2}$ (29)
$\displaystyle\textbf{Re}(w_{1}^{2}w_{2})=1.$

One of the many solutions to this system is

(30) $\displaystyle w_{1}=1-\frac{1}{\sqrt{2}}i\quad
w_{2}=\frac{1}{\sqrt{2}}i.$

Applying this new 2-step time stepper to the nonlinear problem, it indeed has
3rd order accuracy (Fig 6). Similar 2-step paths are derived for
$\dot{y}=f(t)$ and $\dot{y}=-e^{y}$ and they also show expected 3rd order
accuracy in Fig. 6.

Figure 6: Two step complex method Euler methods that have third order accuracy
for two specific problems.

In the process of constructing this integrator, we needed to plug in the right
hand side of the differential equation into Eq. 3.2 to derive the expansion.
This may be an intensive symbolic task for more complicated problems. However,
the right complex steps could be of ‘learned’ for specific and difficult
problems through approaches similar to [17] where a machine learning approach
is used to obtain Runge-Kutta methods tailored to particular ODE families.

In this section, we have shown the process of constructing higher order paths
for explicit Runge-Kutta type integrators. Later, in Section 8, we extend it
to implicit methods. In the next section, we discuss the benefits and costs
when compared to equivalent classical Runge-Kutta methods.

##  4 Operation counts and low storage benefits

Now, we discuss a potential drawback of the complex timestepping schemes, the
increased computational cost for performing complex operations, and a
potential benefit, the low storage cost compared to the Runge-Kutta methods.
In Sec. 4.1, we demonstrate that low storage cost is the main benefit for
complex-valued systems.

Consider the forward Euler method applied to the differential equation
$\dot{y}=f(t,y)$, yielding the step $y_{i+1}=y_{i}+\Delta tf(t,y)$. If the
timestepping is done along the real line, the computational cost is one real
multiplication, one real addition, and the operation cost of one real function
evaluation. If the timestepping is done along a path in the complex plane, the
computational cost is one complex multiplication, one complex addition and the
operation cost of one complex function evaluation. A complex multiplication
costs 4M + 2A or 3M+5A (Gauss’s trick), where M and A are the costs of real
multiplication and addition. The storage also needs to be doubled to store a
complex-valued vector compared to the real-valued vector.

Consider taking a complex timestep $w\Delta t$ with the ordinary differential
equation $\dot{y}=\lambda y$. We can express it as the following system.

(39) $\displaystyle\left[\begin{array}[]{ccccccccccccc}y_{r}\\\
y_{i}\end{array}\right]^{k+1}=\left[\begin{array}[]{ccccccccccccc}y_{r}\\\
y_{i}\end{array}\right]^{k}+\left[\begin{array}[]{ccccccccccccc}(\lambda
w)_{r}\Delta t&-(\lambda w)_{i}\Delta t\\\ (\lambda w)_{i}\Delta t&(\lambda
w)_{r}\Delta t\end{array}\right]\left[\begin{array}[]{ccccccccccccc}y_{r}\\\
y_{i}\end{array}\right]^{k}$

where $w$ represents the complex time step and the subscripts $r$ and $i$
represent the real and imaginary parts. This results in 4M+4A cost of
operations as opposed to M+A operations for purely real time steps. Based on
Eq (39), the complex time stepping process shows promise for parallelization,
which could reduce computational time by up to half.

Similarly, if we integrate the system of N equations $\dot{y}=Ay$ with complex
time steps, it becomes a system of $2N$ equations. In that scenario, the
dominant cost is the matrix-vector multiplication ($Ay$) which takes $N^{2}$
operations with real time steps and $4N^{2}$ operations with complex time
steps. For a k-step Complex Euler method,there are k matrix-vector
multiplications and the operation count is $4kN^{2}$ compared to $kN^{2}$ for
a k stage Runge-Kutta Method.

For nonlinear systems, the dominant cost comes from function evaluations.
Functions can be designed to minimize the computational cost, but the real
equivalent of any complex-valued function can always be calculated using fewer
operations. So, the complex time stepper will always be more expensive than an
explicit Runge-Kutta method if the number of function evaluations is the same.
In the next section, we explore how complex time integrators can require less
function evaluations than the classical Runge-Kutta methods of the same order
by breaking the Runge-Kutta order barrier [15]. In that case, the complex
integrators may end being the less expensive option.

Another potential advantage, is that the storage requirement for complex
stepping is a vector of size 2N while general Runge Kutta scheme store the
intermediate stages. For example, a 5th order complex timestepping scheme
needs to storage a vector of size $2N$ at all times while a classically
implemented 5th order Runge-Kutta method requires memory for vector of size
$5N$. However, there are low-storage implementations of Runge-Kutta schemes
that use fewer registers [18, 19, 20]. Storage is not much of an issue for
typical odes and simple systems but it becomes quite relevant in various
practical settings like in computational fluid dynamics [19] where partial
differential equations like the Navier Stokes equations are solved in 3
dimensions at high resolution.

###  4.1 The linear Schrodinger equation : A clear computational advantage for
the 3-step complex Euler method

In this subsection, we demonstrate an example case where our complex methods
have a clear benefit over existing methods. We compare the 3-step complex
Euler method to a standard RK-3 method (Ralston’s method) for a specific
problem. We choose to examine the linear Schrodinger equation because it is a
partial differential equation with complex values, so the cost of complex
operations is the same in both methods. We show that methods are of the same
order, but because 3-step complex Euler has lower storage requirements, it is
advantageous.

The linear Schrodinger equation is given by

(40) $\displaystyle u_{t}=iu_{xx},\quad x\in(0,2\pi).$

We solve the equation from $t=0$ to $t=3$ on the domain $x\in(0,2\pi)$ with
periodic boundary conditions. The initial condition is given by
$u(t=0)=e^{ix}$ and the Fourier spectral method is used for the spatial
differencing. The exact solution is given by $u(t)=e^{i(x-t)}$.

Figure 7: The linear Schrodinger is solved using both Runge-Kutta 3 (Ralston’s
method) and the complex 3-step method. Both methods result in nearly identical
error ($10^{-9}$ using the infinity norm).  Figure 8: Runge-Kutta 3 (Ralston’s
method) and the complex 3-step method require the same the number of function
calls and result in nearly identical error when computing the linear
Schrodinger equation.

Figs 7 and 8 show that the two methods have near identical error and require
equal number of function evaluations when integrating the linear Schrodinger
equation. The 3-step complex method has one advantage over the Runge-Kutta-3
method: Since the calculating the next step only requires the numerical
solution at the current step, the required storage is the size of the initial
condition vector. The Runge-Kutta-3 method requires thrice the storage to
store the 3 stages during each function evaluation in its basic
implementation.

##  5 Breaking the order barrier for explicit Runge-Kutta methods with complex
coefficients

One can construct an $s$ stage explicit Runge-Kutta method with order $p$
where $s=p$ up to order $p=4$. But for order 5 and beyond, the number of
stages needs to be more than the order desired [15]. This is known as the
order barrier for Runge-Kutta methods [14].

Specifically for the 5th order Runge-Kutta-5 method, 6 stages are required and
hence 6 function evaluations. This is because for 5th order accuracy, there
are 17 order conditions and the 5 stages only provide 15 free parameters to
satisfy them.

Can the order barrier be broken with complex time stepping? If we take $n$
complex Euler steps and use the trick of taking the real part and keeping the
$>1$ order error purely imaginary, we have $2n-1$ free parameters ($n$ complex
variables that need to return to the real line). $2n-1$ free parameters are
insufficient even for 4th order accuracy because there are 8 order conditions
and only 7 free parameters for $n=4$.

Complex stepping can be applied to numerical methods apart from Forward Euler
as well. What if we allowed some coefficients in a 5-stage Runge Kutta step to
be complex variables and constrained the error to be purely imaginary? That
easily increases the number of free parameters from 15 to 17 or higher in
order to solve the 17 order conditions!

Another approach to break the $p=5$ order barrier is to take two complex
Runge-Kutta steps: one Runge-Kutta-2 step and then one Runge-Kutta-3 step. If
we require that the two complex Runge-Kutta steps come back to the real line,
the total number of free parameters is 10 (3 free parameters for Runge-
Kutta-2, 6 free parameters for Runge-Kutta-3 and the two time steps sizes with
the constraint that they add up to $\Delta t$ (purely real)). If all 10 of
those free parameters are complex valued and the $>1$ order error is purely
imaginary, we get 20 free parameters which are sufficient to solve the 17
order conditions with 5th order accuracy. Using this method, we can achieve
5th order accuracy with only 5 function evaluations, breaking the order
barrier. As these are complex function evaluations, the cost is higher than
the real equivalent, but there are problems and architectures for which this
is likely not an issue as discussed in Sec. 4.

As a proof of concept, we find a 5th order accurate complex Runge-
Kutta-2/Runge-Kutta-3 method for the special case of scalar autonomous ODEs.
We choose this particular case because it is more tractable than the 17 order
conditions for the general case, which require a more complicated derivation
involving a rooted tree approach [15]. Even with the scalar autonomous case,
there are 12 order conditions to satisfy (which cannot be done with only 10
free parameters of the purely real version of the Runge-Kutta-2/Runge-Kutta-3
method). Here, we show that the complex Runge-Kutta-2/Runge-Kutta-3 method
returns a 5th order accurate numerical solution for scalar autonomous ODEs
with 5 function evaluations, something that cannot be done with its real
equivalent.

Figure 9: The complex two step Runge-Kutta-2/Runge-Kutta-3 method with complex
Runge-Kutta coefficients has 5th order accuracy for autonomous scalar ODEs.

The following equations describe our two step complex method. First an Runge-
Kutta-2 step of size $(b_{11}+b_{12})\Delta t$ goes from $y_{n}$ to $y_{m}$
and then an Runge-Kutta-3 step of size $(b_{21}+b_{22}+b_{23})\Delta t$ to
goes from $y_{m}$ to $y_{n+1}$. The steps and stages are

(41) $\displaystyle y_{n}$ $\displaystyle=y_{n}$ (42) $\displaystyle k_{11}$
$\displaystyle=f(y_{n})$ (43) $\displaystyle k_{12}$
$\displaystyle=f(y_{n}+a_{121}\Delta tk_{11})$ (44) $\displaystyle y_{m}$
$\displaystyle=y_{n}+b_{11}\Delta tk_{11}+b_{12}\Delta tk_{132}$ (45)
$\displaystyle k_{21}$ $\displaystyle=f(y_{m})$ (46) $\displaystyle k_{22}$
$\displaystyle=f(y_{m}+a_{221}\Delta tk_{21})$ (47) $\displaystyle k_{23}$
$\displaystyle=f(y_{m}+a_{231}\Delta tk_{21}+a_{232}\Delta tk_{22})$ (48)
$\displaystyle y_{n+1}$ $\displaystyle=\textbf{Re}(y_{m}+b_{21}\Delta
tk_{21}+b_{22}\Delta tk_{22}+b_{23}\Delta tk_{23}).$

Ensuring the method returns to the real line at an effective time step of
$\Delta t$ introduces the constraint $b_{11}+b_{12}+b_{21}+b_{22}+b_{23}=1$.

With these 9 complex variables, we solve the 12 order conditions and the 1
constraint guaranteeing that the time-stepper returns to the real line. The
resulting scheme provides 5th order accuracy for scalar autonomous ODEs with
only 5 function evaluations (Fig. 9).

##  6 Expanded regions of stability with complex time-steps

In the previous sections, we designed paths in the complex plane to achieve
higher order. In this section, we leverage complex paths in an alternative way
and expand the time steppers’s region of absolute stability. First, we discuss
the linear stability regions of the complex integrators developed in Sec. 3.1.
Next, we show how these regions can be expanded by designing the complex steps
to optimize a desired stability polynomial. Expanded stability regions have
been studied in a similar manner previously in the context of Runge-Kutta
methods where the order of accuracy for the Runge-Kutta method is sacrificed
to obtain a larger stability region [21, 22, 23, 24, 25, 26, 27]. In most
previous works, a new Runge-Kutta integrator (different integration method)
was derived that had a stability function/polynomial corresponding to the
desired stability region. Here, we keep the same integration method and
traverse a path in the complex plane to obtain the desired stability region
(different integration domain). We illustrate this approach using the forward
Euler method. We also show that complex stepping may enable even larger
regions designed for specific problems.

Consider the linear equation $\dot{y}=\lambda y$ studied in Sec. 3.1. Most
integrators, including those discussed here, convert the linear equation into
the difference equation $y(t+\Delta t)=\Phi(\lambda\Delta t)y(t)$ where
$\Phi(z)$ is the stability function of the integrator [28]. The region
$z\in\mathbb{C}$ where $|\Phi(z)|\leq 1$ is the region of absolute stability.
For a complex $n$-step $n$th order method, the stability function is given by
the first $n+1$ terms in the Taylor series of $e^{\Delta t}$. For the 3 step
3rd order method, the stability function is

(49) $\displaystyle\Phi(z)=1+z+\frac{z^{2}}{2}+\frac{z^{3}}{3!}.$

The regions of absolute stability for the 1 step 1st order, 2 step 2nd order
and 3 step 3rd order methods (Fig 10) are the same as those for the n-stage
Runge-Kutta because they have the same stability functions [28]. Note that
$\Delta t$ in all these methods refers to the step size separating the
beginning and end points of the complex path and hence a real time step. Since
the 1-step, 2-step and 3-step methods require a different number of steps and
hence a different number of function evaluations to take a net step of size
$\Delta t$, a standard metric for comparison among methods is the effective
step size $\Delta t_{eff}=\Delta t/n$ where $n$ is the number of steps in the
complex integrator case or the number of stages in the case of Runge Kutta
methods [25].

Figure 10: The 1-step first order, 2-step 2nd order and 3-step 3rd order
complex Euler methods share the same stability regions as the 1,2 and 3 stage
Runge-Kutta methods (left). Corresponding complex time-stepping paths are
given on right.

We would like to expand the regions in Fig 10 by optimizing the stability
function. If we relax our requirement for a 3-step integrator to be 3rd order
accurate, and require only 2nd order accuracy, its stability function would be
given by

(50) $\displaystyle\Phi(z)=1+z+\frac{z^{2}}{2}+kz^{3}.$

The coefficient $k$ in the above expression is the same as the coefficient in
front of the $\lambda^{3}\Delta t^{3}$ in Eq (5): $w_{1}w_{2}w_{3}$. We can
choose the complex steps $w_{1},w_{2}$ and $w_{3}$ such that they satisfy the
first and second order conditions
($w_{1}+w_{2}+w_{3}=1,w_{1}w_{2}+w_{2}w_{3}+w_{1}w_{3}=1/2$)and make
$w_{1}w_{2}w_{3}$ equal to our desired value of $k$. The optimal $k$ value is
the one that would allow us to take the largest stable time step. For problems
with purely real and imaginary eigenvalues, the optimal $k$ values and the
associated optimal stability functions/polynomials have been long-studied [21,
24, 23]. In the past decade, Ketcheson and co-authors [26, 27, 25] have
explored optimal stability functions in the context of $\dot{y}=Ly$ where $L$
is a matrix with quite arbitrary spectra. In [25], Ketcheson and Ahmadia
describe an optimization algorithm to obtain the optimal stability polynomial
that allows the maximal time step depending on the spectra of $L$. This
algorithm (RKOpt [29]) is publicly available through Nodepy [16].

By choosing the value $k$ that results in the largest permissible time step
along the negative real eigenvalue line, we achieve a 3-step 2nd order method
(green) with expanded stability along the negative real axis compared to the
3-step 3rd order method (red Fig 11). Requiring only 1st order accuracy for
the 3-step method can nearly triple the size of the stability region along the
negative real axis (gold) and surprisingly requires only real time-steps (Fig
11).

Figure 11: By sacrificing accuracy, one can expand the stability region (left)
using complex time steps, as demonstrated by these optimized 3-step complex
Euler paths (right) with 3rd, 2nd, and 1st order accuracy.

To our knowledge, our work is the first to use complex time steps to create
stability polynomials, enabling us to expand the stability region. Next, we
will demonstrate a potential path for even larger stability regions via
complex stepping.

The use of complex time steps in an integrator allows to expand the optimal
stability polynomials beyond those with real coefficients. Consider the
challenge of finding a 3-step 2nd order integrator with an optimal stability
region that encompasses the largest possible timestep along the eigenvalue
$\lambda=-1-2i$ (red dots Fig 12). The red dots further away from the origin
correspond to larger $\Delta t$ values along the eigenvalue. We want to find
the optimal stability polynomial/function $\Phi(\lambda\Delta t)$ that permits
the largest $\Delta t$ for a 3 stage/step 2nd order method in this scenario.

The optimal stability polynomial returned by RKOpt[29] (magenta Fig. 12) is

(51) $\displaystyle\Phi(z)=1+z+\frac{1}{2}z^{2}+0.1134z^{3}.$

Through trial and error, we find that an even better stability polynomial
(brown Fig. 12) given by

(52) $\displaystyle\Phi(z)=1+z+\frac{1}{2}z^{2}+0.1134z^{3}-0.06iz^{3}$

The complex coefficient allows for a larger stable time step size for a
problem with this eigenvalue. A trade-off is that the optimal region is no-
longer symmetric. So such regions would be primarily useful for complex
problems, where complex eigenvalues appear outside conjugate pair.
Alternatively, better stability region designs may be possible using more
sophisticated optimization methodology.

Figure 12: Using complex time steps produces stability polynomials with
complex coefficients. The resulting stability regions (left) can expand
further in a desired direction than those using purely real coefficients.
Corresponding paths on right.

This section has been an exploration of the paths in the complex plane on
which an integrator has improved linear stability properties. In the next
section, we consider nonlinear stability properties particularly strong
stability properties. Later, in Section 8, we discuss the derivation of
implicit complex methods with properties like L-stability.

##  7 Remarks on Strong Stability Preserving (SSP) properties

Time integrators with SSP properties are desirable because they extend they
provide guarantees of stability for nonlinear problems and non-smooth
solutions. In this section, we highlight the need for new analysis for
deriving practical SSP conditions for complex timesteppers and we show
numerically that complex timestepping can result in larger SSP time steps.

Consider the following differential equation system,

$u_{t}=F(u).$

If one use the forward Euler method is used to integrate the above system, the
system would be said to satisfy a SSP property if

$|u_{n+1}|=|u_{n}+\Delta tF(u_{n})|\leq|u_{n}|\quad\text{for}\quad 0\leq\Delta
t\leq\Delta t_{FE}$

where the maximal SSP time step is given by

$|\Delta t_{FE}|=\min_{n}\frac{2}{\bigg{|}\frac{F(u_{n})}{u_{n}}\bigg{|}}.$

Much effort has focused on obtaining optimal SSP Runge-Kutta methods in
various contexts [30, 20]. The general argument used to show SSP is as
follows; if forward Euler is SSP for some $\Delta t\leq\Delta t_{FE}$, these
optimal SSP Runge-Kutta methods can be expressed as a convex combination of
the forward Euler steps and guaranteed to be SSP for some $\Delta t\leq
C\Delta t_{FE}$. This $C$ is known as the SSP coefficient.

Following this line of argument, consider our $k$ step general complex Euler
time stepper expressed as

$\displaystyle u^{0}$ $\displaystyle=u_{n}$ $\displaystyle u^{i}$
$\displaystyle=u^{i-1}+w_{i}\Delta tF(u^{i-1}),i=1,\ldots,k$ $\displaystyle
u_{n+1}$ $\displaystyle=u^{k}.$

Here, $u_{n},u_{n+1}$ represent the solution at $n$ and $n+1$ time points on
the real line. The $u^{i}$ represent the solution at the complex time points
in between and the $w_{i}\Delta t$ are the complex time steps. Unlike the
Runge-Kutta methods in their Shu-Osher form, finding an SSP coefficient using
convex combinations of forward Euler is less straightforward because of the
complex nature of $w_{i}$ and $F(u^{i-1})$.

A strict condition for an SSP complex time integrator can be derived if we
require

$|u^{i}|=|u^{i-1}+w_{i}\Delta tF(u^{i-1})|\leq|u^{i-1}|,\quad i=1,\ldots,k.$

This results in a very restrictive time step.

$\Delta
t\leq\min_{i}\left(-2\frac{\textbf{Re}\bigg{(}\frac{w_{i}F(u^{i-1})}{u^{i-1}}\bigg{)}}{|w_{i}|^{2}\bigg{|}\frac{F(u^{i-1})}{u^{i-1}}\bigg{|}^{2}}\right).$

Notably, the lower bound on the largest possible SSP timestep can be as small
as zero if $\textbf{Re}\bigg{(}\frac{w_{i}F(u^{i-1})}{u^{i-1}}\bigg{)}$
vanishes. This is too strict of a condition to obtain a practical SSP
timestep. We can obtain a larger timestep with a looser SSP condition if we
require that only the solutions at the real time points at the end of the path
in the complex plane satisfy the SSP property $|u_{n+1}|\leq|u_{n}|$. For the
two step second order complex time integrator which takes the complex steps
$w$ and $\bar{w}$ at each real time point, where $w,\bar{w}$ are
$\frac{1}{2}+\frac{1}{2}i$ and $\frac{1}{2}-\frac{1}{2}i$, , this translates
to the following condition.

$|u_{n}+w\Delta tF(u_{n})+\bar{w}F(u_{n}+w\Delta tF(u_{n}))|\leq|u_{n}|$

We can also take the $\bar{w}$ step first and then the $w$ step before
returning to the real line. Taking different paths with the same number of
steps and the same order results in different SSP conditions. However, we have
not yet been able to obtain a clean expression for an upper bound on the SSP
timestep using this approach.

Using numerical experiments, we find the complex timestepping methods could
allow for a larger SSP stable timestep. We demonstrate for the simple case,

(53) $\displaystyle\dot{y}=f(y)=-ye^{-y}.$

We consider three 2nd order timestepping methods.

  * •

The midpoint method:

(54) $\displaystyle u_{n+1}=u_{n}+\Delta t_{n}f(u_{n}+\Delta t_{n}/2f(u_{n}))$

  * •

The optimal SSP Runge-Kutta 2 method: This is the second order RK method with
the optimal SSP coefficient of 1. [30]

$\displaystyle u_{m}$ $\displaystyle=u_{n}+\Delta t_{n}f(u_{n})$ (55)
$\displaystyle u_{n+1}$
$\displaystyle=\frac{u_{n}}{2}+\frac{u_{m}}{2}+\frac{1}{2}\Delta
t_{n}f(u_{m})$

  * •

The complex 2-step method: We assume our initial condition is real. So, we
make sure to take the real part of the solution at the end of the two steps.

$\displaystyle u_{m}$ $\displaystyle=u_{n}+w_{1}\Delta t_{n}f(u_{n})$ (56)
$\displaystyle u_{n+1}$ $\displaystyle=\textbf{Re}(u_{m}+w_{2}\Delta
t_{n}f(u_{m}))$

We numerically compute the largest timestep $\Delta t_{n}$ for each method for
various $u_{n}$ values such that for $\Delta t\leq\Delta t_{n}$, the SSP
condition ($|u_{n+1}|\leq|u_{n}|$) is satisfied. Different methods permit the
largest time-step for different sizes of $u_{n}$. However, for most values,
the complex 2-step method enables comparable maximal time-steps and for large
values of $u_{n}$, the complex 2-step method appears to enable larger maximum
time-steps (Fig. 13). While we have not yet been able to prove SSP properties
for complex methods, these numerical experiments show it may have desirable
properties.

Figure 13: For various values of $u_{n}$, the complex 2-step method allows us
to take a larger SSP timestep.

##  8 Dealing with stiff systems: Implicit methods

When dealing with stiff systems, expanding the stability region of the
integrator is rarely sufficient and implicit methods with properties like
L-stability are required. In this section, we derive a higher order complex
path for an implicit method: The implicit midpoint method. First we study
linear ODEs (Sec. 8.1) and then we extend to general, nonlinear, non-
autonomous differential equations (Sec. 8.2).

###  8.1 Linear differential equations

Consider the linear equation $\dot{y}=\lambda y$. Applying the implicit
midpoint method to this equation produces the following difference equation,

(57) $\displaystyle y_{n+1}=y_{n}+\frac{\lambda\Delta t}{2}(y_{n}+y_{n+1})$

which can also be written as

(58) $\displaystyle y_{n+1}=y_{n}\frac{1+\lambda\Delta t/2}{1-\lambda\Delta
t/2}.$

Now, suppose we take two complex steps $w_{1}\Delta t$ and $w_{2}\Delta t$ to
go from $y_{n}$ to $y_{n+1}$. Our numerical approximation at $y_{n+1}$ is

(59) $\displaystyle y_{n+1}=y_{n}\frac{(1+\lambda w_{1}\Delta t/2)(1+\lambda
w_{2}\Delta t/2)}{(1-\lambda w_{1}\Delta t/2)(1-\lambda w_{2}\Delta t/2)}.$

If $y(t)=y_{n}$, the true solution at time $t+\Delta t$ would be
$y_{n}e^{\lambda\Delta t}$. This solution can be represented to $4$th order
accuracy by the following rational function expansion of the exponential
function, known as the (2,2) Pade approximate,

(60) $\displaystyle y_{t+\Delta t}=y_{n}\frac{1+\lambda\Delta
t/2+\lambda^{2}\Delta t^{2}/12}{1-\lambda\Delta t/2+\lambda^{2}\Delta
t^{2}/12}+O(\Delta t^{5}).$

Setting the right hand sides of (59) and (60) equal to each other, gives us
the 2-step path in the complex plane on which the implicit midpoint method has
$4$th order accuracy.

By solving

(61) $\displaystyle w_{1}+w_{2}=1\quad,\quad w_{1}w_{2}=\frac{1}{3},$

we find that the two 2-step 4th order paths are
$(\frac{1}{2}+i\frac{1}{2\sqrt{3}},\frac{1}{2}-i\frac{1}{2\sqrt{3}})$ and
$(\frac{1}{2}-i\frac{1}{2\sqrt{3}},\frac{1}{2}+i\frac{1}{2\sqrt{3}})$.

Equations (59) and (60), provide a clear connection between the optimal k-step
complex path for the implicit midpoint method and the (k,k) Pade approximant.
We can derive the optimal k-step complex path for the Backward Euler method in
a similar way, by setting it equal to the (0,k) Pade approximate. This results
in the same paths derived in Subsection 3.1 for the forward Euler method.

We test out our newfound implicit paths on the heat equation with Dirichlet
boundary conditions.

(62) $\displaystyle u_{t}=u_{xx},\quad x\in(0,1)$

We choose the initial condition $u(x,0)=\sin(\pi x)$, which has the exact
solution $u(x,t)=e^{-\pi^{2}t}\sin(\pi x)$. The spatial discretization is
performed using 4th order finite differencing, periodic boundary conditions
and a grid size of 10000 cells. As expected, the 2-step complex implicit
midpoint method has 3rd order accuracy and the complex 3-step backward Euler
method has 4th order accuracy (Fig. 14).

Figure 14: The 2-step implicit midpoint method has 4th order accuracy while
the 3-step Backward Euler has 3rd order accuracy for the heat equation.

###  8.2 General nonlinear non-autonomous differential equations

In the previous subsection, we found that the path
$(\frac{1}{2}+i\frac{1}{2\sqrt{3}},\frac{1}{2}-i\frac{1}{2\sqrt{3}})$ allows
the implicit midpoint method to be 4th order accurate when solving linear
differential equations. Does the expected 4th order accuracy of the 2-step
complex path for the implicit midpoint method derived in the previous
subsection hold for general nonlinear differential equations? Numerically, we
confirm this to be true, provided we use the trick of taking the real part of
the numerically calculated values at the end of the complex path. To prove 4th
order accuracy analytically, we analyze the midpoint-method following a
similar line of reasoning as our previous nonlinear analysis (Equations (11) -
(3.2)). Suppose that we take two implicit steps $w_{1}\Delta t$ and
$w_{2}\Delta t$ from $y_{n}$ at $t_{n}$ to $y_{n+1}$ at at $t+\Delta t$ with
$y_{m}$ being the intermediate solution at at $t+w_{1}\Delta t$.

(63) $\displaystyle y_{n}$ $\displaystyle=y(t_{n})$ (64) $\displaystyle y_{m}$
$\displaystyle=y_{n}+w_{1}\Delta tf(t_{n}+w_{1}\Delta t/2,(y_{n}+y_{m})/2)$
(65) $\displaystyle y_{n+1}$ $\displaystyle=y_{m}+w_{2}\Delta
tf(t_{n}+w_{1}\Delta t+w_{2}\Delta t/2,(y_{m}+y_{n+1})/2)$

Unlike the previous analysis, equations (64) and (65) are implicit equations
and cannot be directly Taylor expanded around $f(t,y_{n})$. To overcome this
difficulty and calculate the desired $w_{1}$ and $w_{2}$ step sizes, we
replace the two implicit equations with explicit versions that are accurate to
the original implicit equations to at least 4th order. The explicit versions
are found by recursively substituting the implicit equation into itself in an
iterative process similar to Picard Iteration. For example, if you have an
implicit equation given by

(66) $\displaystyle y=y_{0}+\Delta tg(y)$

one can recursively apply it twice to obtain the following explicit equation,
which is accurate to the original implicit equation to the 3rd order.

(67) $\displaystyle y=y_{0}+\Delta tg(y_{0}+\Delta tg(y_{0}+\Delta
tg(y_{0})))$

To obtain the 4th order accurate explicit approximation of equations (64) and
(65), they need to be recursively applied at least 3 times. Note that the
explicit equations are 4th order accurate approximations to the original
implicit difference equations and not the original differential equation.

At the end of this process,, we substitute
$(w_{1},w_{2})=(\frac{1}{2}+i\frac{1}{2\sqrt{3}},\frac{1}{2}-i\frac{1}{2\sqrt{3}})$
into the explicit approximation of equation (65) and obtain

$\displaystyle y_{n+1}=$ $\displaystyle y_{0}+\Delta
tf(t_{n},y_{n})+\frac{1}{2}\Delta
t^{2}\bigg{(}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}+\frac{\partial}{\partial
t}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{n},y=y_{n}}$
$\displaystyle+\frac{1}{3!}\Delta
t^{3}\bigg{(}f^{2}{\left(t,y\right)}\frac{\partial^{2}}{\partial
y^{2}}f{\left(t,y\right)}+\left(\frac{\partial}{\partial
y}f{\left(t,y\right)}\right)^{2}f{\left(t,y\right)}$
$\displaystyle+2f{\left(t,y\right)}\frac{\partial^{2}}{\partial y\partial
t}f{\left(t,y\right)}+\frac{\partial}{\partial
t}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}+\frac{\partial^{2}}{\partial
t^{2}}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{n},y=y_{n}}$
$\displaystyle+\frac{1}{4!}dt^{4}\bigg{(}(1+0.096225i)f^{3}{\left(t,y\right)}\frac{\partial^{3}}{\partial
y^{3}}f{\left(t,y\right)}$
$\displaystyle+(4-0.288675i)f^{2}{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}\frac{\partial^{2}}{\partial
y^{2}}f{\left(t,y\right)}+(3+0.288675i)f^{2}{\left(t,y\right)}\frac{\partial^{3}}{\partial
y^{2}\partial t}f{\left(t,y\right)}$
$\displaystyle+3f{\left(t,y\right)}\frac{\partial}{\partial
t}f{\left(t,y\right)}\frac{\partial^{2}}{\partial
y^{2}}f{\left(t,y\right)}+f{\left(t,y\right)}\left(\frac{\partial}{\partial
y}f{\left(t,y\right)}\right)^{3}$
$\displaystyle+(5-0.57735i)f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}\frac{\partial^{2}}{\partial y\partial
t}f{\left(t,y\right)}+(3+0.288675i)f{\left(t,y\right)}\frac{\partial^{3}}{\partial
y\partial t^{2}}f{\left(t,y\right)}$ $\displaystyle+\frac{\partial}{\partial
t}f{\left(t,y\right)}\left(\frac{\partial}{\partial
y}f{\left(t,y\right)}\right)^{2}+3\frac{\partial}{\partial
t}f{\left(t,y\right)}\frac{\partial^{2}}{\partial y\partial
t}f{\left(t,y\right)}+(1-0.096225i)\frac{\partial^{3}}{\partial
t^{3}}f{\left(t,y\right)}$ (68)
$\displaystyle+(1-0.288675i)\frac{\partial^{2}}{\partial
t^{2}}f{\left(t,y\right)}\frac{\partial}{\partial
y}f{\left(t,y\right)}\bigg{)}\and[\Bigg{]}{t=t_{n},y=y_{n}}+O(\Delta t^{5}).$

We compare (8.2) to the first 5 terms of Taylor expansion in equation
(3.2).The error at O($\Delta t^{4}$) is purely imaginary, which demonstrates
we can get 4th order accuracy by extracting the real component of the solution
at real time points (as long as the true solution is real at real time
points).

Now, that we have proved our two-step implicit midpoint method to be 4th order
accurate for nonlinear systems, we test it on the Van Der Pol oscillator,

(69) $\displaystyle\dot{y}=x-\frac{x^{3}}{3}-\frac{\dot{x}}{\mu}$

which can be rewritten as the following system

$\displaystyle\dot{x}$ $\displaystyle=y$ (70) $\displaystyle\dot{y}$
$\displaystyle=\mu(1-x^{2})y-x.$

We set $\mu=10$ and the initial condition as $(x,y)=(2,0)$. We achieve the
expected 4rd order accuracy for the 2-step complex implicit Method (Fig. 14).

Figure 15: The 2-step implicit midpoint method has 4th order accuracy while
the 3-step Backward Euler has 3rd order accuracy for the Van Der Pol
oscillator.

##  9 Discussion

This paper is an exploration of the benefits and trade offs of leaving the
real line to perform numerical wi integration in the complex plane. Our first
goal was to explain the process of developing a complex time-stepper and its
potential benefits in an accessible manner. For both linear and nonlinear
systems, we demonstrate how to derive paths in the complex plane that increase
the order of accuracy of any Runge-Kutta type integrator. The derivation of
these paths for implicit methods involves a more complicated iterative process
and the implicit paths for linear problems show us a natural connection to
Pade approximants. For real valued differential equations, existing Runge-
Kutta methods are less costly since they do not require complex arithmetic.
However, for complex-valued differential equations such as the Schrodinger
equation, the computational cost is similar and the low storage structure of
complex time steppers may be beneficial.

Complex time steppers enable the major advance of breaking the Runge-Kutta
order barrier[15]. Using complex time steps and complex Runge-Kutta
coefficients, we can increase the number of variables available to satisfy the
ever increasing number of order conditions and use fewer function evaluations
to achieve the same order accuracy. For example, a traditional Runge-Kutta
method requires at least 11 function evaluations to obtain 8th order of
accuracy and complex timesteppers could reduce the number of evaluations down
to 8. Whether one would ever need an 8th order method is a different question,
but breaking the order barrier has been a long-standing challenge in the
computational world.

We also showed that complex time stepping can be used to expand stability
regions, allowing larger effective time steps. Expanded stability regions have
been designed using Runge-Kutta methods [21, 25], but complex coefficients in
the optimal stability polynomial may support even larger stability regions in
specific cases. The fact that you can increase the stability region by
changing the integration domain rather than the integration method may also be
appealing from an easy user implementation standpoint. Our preliminary results
about nonlinear stability properties such as SSP, suggest that future
exploration may unearth specific paths for specific methods with desirable
stability properties.

We intend this paper to be introduction to and initiation of the search for
better time integrators that take advantage of paths in the complex plane. We
have mostly focused on creating integrators with increased accuracy and
stability, but there are many other desirable properties. The complex plane
offers an extra dimension on which any integrator can be improved, opening up
a new realm of possibilities when designing time stepping methods.

##  10 Acknowledgements

We would like to thank Emil Constantinescu, Alvin Bayliss and David Chopp for
extensive feedback and valuable suggestions.

## References

  * [1] James Ward Brown and Ruel V Churchill.  Complex variables and applications eighth edition.  McGraw-Hill Book Company, 2009.
  * [2] Norman S Nise.  Control systems engineering.  John Wiley & Sons, 2020.
  * [3] Nicholas J. Higham.  Differentiation with(out) a difference.  SIAM News, 51(5):2, June 2018.
  * [4] George F Corliss.  Integrating odes in the complex plane—pole vaulting.  Mathematics of Computation, 35(152):1181–1189, 1980.
  * [5] Thorsten Orendt, Jürgen Richter-Gebert, and Michael Schmid.  Geometry of numerical complex time integration.  arXiv preprint arXiv:0903.1585, 2009.
  * [6] JE Chambers.  Symplectic integrators with complex time steps.  The Astronomical Journal, 126(2):1119, 2003.
  * [7] François Castella, Philippe Chartier, Stéphane Descombes, and Gilles Vilmart.  Splitting methods with complex times for parabolic equations.  BIT Numerical Mathematics, 49(3):487–508, 2009.
  * [8] Winfried Auzinger, Harald Hofstätter, David Ketcheson, and Othmar Koch.  Practical splitting methods for the adaptive integration of nonlinear evolution equations. part i: Construction of optimized schemes and pairs of schemes.  BIT Numerical Mathematics, 57(1):55–74, 2017.
  * [9] Fernando Casas, Philippe Chartier, Alejandro Escorihuela-Tomàs, and Yong Zhang.  Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations.  Journal of Computational and Applied Mathematics, 381:113006, 2021\.
  * [10] Sergio Blanes, Fernando Casas, and Ander Murua.  Splitting methods with complex coefficients.  SeMA Journal, 50(1):47–60, 2010.
  * [11] Eskil Hansen and Alexander Ostermann.  High order splitting methods for analytic semigroups exist.  BIT Numerical Mathematics, 49(3):527–542, 2009.
  * [12] Fernando Casas and Alejandro Escorihuela-Tomàs.  High order integrators obtained by linear combinations of symmetric-conjugate compositions.  arXiv preprint arXiv:2106.06503, 2021.
  * [13] Tommaso Buvoli and Mayya Tokman.  Constructing new time integrators using interpolating polynomials.  SIAM Journal on Scientific Computing, 41(5):A2911–A2937, 2019.
  * [14] J. Butcher.  Runge-Kutta methods.  Scholarpedia, 2(9):3147, 2007.  revision #91735.
  * [15] John Charles Butcher.  On fifth and sixth order explicit runge-kutta methods: order conditions and order barriers.  Canadian Applied Mathematics Quarterly, 17(3):433–445, 2009.
  * [16] David I Ketcheson, Hendrik Ranocha, Matteo Parsani, Umair bin Waheed, and Yiannis Hadjimichael.  Nodepy: A package for the analysis of numerical ode solvers.  2020\.
  * [17] Yue Guo, Felix Dietrich, Tom Bertalan, Danimir T Doncevic, Manuel Dahmen, Ioannis G Kevrekidis, and Qianxiao Li.  Personalized algorithm generation: A case study in meta-learning ode integrators.  arXiv preprint arXiv:2105.01303, 2021.
  * [18] Mark Huitt Carpenter and Christopher A Kennedy.  Fourth-order 2n-storage runge-kutta schemes.  1994\.
  * [19] Christopher A. Kennedy, Mark H. Carpenter, and R.Michael Lewis.  Low-storage, explicit runge–kutta schemes for the compressible navier–stokes equations.  Applied Numerical Mathematics, 35(3):177–219, 2000.
  * [20] David I Ketcheson.  Highly efficient strong stability-preserving runge–kutta methods with low-storage implementations.  SIAM Journal on Scientific Computing, 30(4):2113–2136, 2008.
  * [21] W Riha.  Optimal stability polynomials.  Computing, 9(1):37–43, 1972.
  * [22] J Douglas Lawson.  An order five runge-kutta process with extended region of stability.  SIAM Journal on Numerical Analysis, 3(4):593–597, 1966.
  * [23] Assyr Abdulle.  On roots and error constants of optimal stability polynomials.  BIT Numerical Mathematics, 40(1):177–182, 2000.
  * [24] Andrei Bogatyrëv.  Effective computation of optimal stability polynomials.  Calcolo, 41(4):247–256, 2004.
  * [25] David Ketcheson and Aron Ahmadia.  Optimal stability polynomials for numerical integration of initial value problems.  Communications in Applied Mathematics and Computational Science, 7(2):247–271, 2013.
  * [26] Matteo Parsani, David I Ketcheson, and W Deconinck.  Optimized explicit runge–kutta schemes for the spectral difference method applied to wave propagation problems.  SIAM Journal on Scientific Computing, 35(2):A957–A986, 2013.
  * [27] Ethan J Kubatko, Benjamin A Yeager, and David I Ketcheson.  Optimal strong-stability-preserving runge–kutta time discretizations for discontinuous galerkin methods.  Journal of Scientific Computing, 60(2):313–344, 2014.
  * [28] Randall J LeVeque.  Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems.  SIAM, 2007.
  * [29] David I Ketcheson, Matteo Parsani, Zachary Grant, Aron Ahmadia, and Hendrik Ranocha.  Rk-opt: A package for the design of numerical ode solvers.  2020\.
  * [30] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor.  Strong stability-preserving high-order time discretization methods.  SIAM review, 43(1):89–112, 2001.