On theoretical upper limits for valid timesteps of implicit ODE methods

1.   Introduction

Many mathematical models take the form of a system of ordinary differential equations (ODEs) for a vector of unknowns ${\mathit{\boldsymbol{q}}}(t) \in \mathbb{R}^m$ subject to boundary data:

Standard transformations reduce systems that are higher order, non-autonomous, or subject to interior-point data to this first-order autonomous form with boundary data at the cost of increased system dimension ^[9]. When $\tau_i \equiv \tau_0$, $i = 1, 2, \ldots, m$, we have an initial-value problem (IVP); otherwise, we have a two-point boundary-value problem (BVP).

In practice, numerical methods for the solution for such ODEs involve successive approximations at successive timesteps and are either implicit or explicit. Implicit methods typically involve the iterative solution, at each time step, of systems of nonlinear algebraic equations, generally a theoretically infinite process with a potentially non-unique or non-existent result. Explicit methods, in contrast, can be implemented directly, generally a theoretically finite process with a unique result. The iterative solution process of an implicit method can incur a significant run-time cost, but the use of such methods may result in greater overall efficiency or fidelity. For example, the increase in the timestep afforded by an implicit method when solving stiff ODEs typically offsets the increased cost per step. Also, when integrating a Hamiltonian system, an implicit method may be used to arrange that the simulation itself preserves energy or is symplectic ^[3,23].

The existence and uniqueness theory for IVPs is much more decisive than for BVPs. IVPs have unique solutions under mild assumptions that are typically satisfied in practice, whereas BVPs may have from zero to uncountably many solutions. However, when an implicit method is involved in approximating the solution of an IVP, the possibility emerges of divergence or convergence to one of multiple solutions. Convergence to spurious solutions is well recognized (^[5,11] and references therein), particularly in the numerical solution of BVPs (^[9,16,17,18] and references therein). Less attention is typically given to the context of solving IVPs (^[14,15,19] and references therein), where there is theoretically a unique solution and where the presence of a "good" initial guess is taken for granted. In this article, we consider the context of an implicit IVP method that has multiple solutions at a given timestep and how to choose from among them, as opposed to a qualitatively incorrect numerical solution of an IVP or BVP.

Standard methods for error estimation and control via timestep selection tend to adjust timesteps such that, in practice, any ambiguity arising from multiple solutions is avoided, but they are not usually specifically designed to do so. But there are two specific scenarios in which constant timesteps are often used in practice: simulations of symplectic systems using a symplectic method ^[3] and simulations of large physical systems. Efforts toward adaptive symplectic methods have been made, but they tend to be specialized and require a significant amount of user judgment ^{[20,21,22,33]}. Software packages for the simulation of large physical systems, especially on distributed architectures ^{[29,30,31,32]}, often use constant timesteps, because of the large relative expense of estimating the error and changing the timestep. However, at constant timestep, a simulation may unexpectedly encounter a time-localized region of complex dynamics, and convergence can easily be construed as nominal when in fact it is not. A large number of independent simulations, e.g., explorations of a parameter space, may not be feasible or efficient at small constant timestep. It may not be easy to automate the detection of anomalous behaviour in the absence of convergence failure, and a manual inspection may not be feasible.

Numerical methods for solution of ODEs can also be classified according to how many past steps are stored and used in computing the next step. In general, the next step can be a function of $k$ past steps, leading to multi-step methods. If only the current state is stored, then the method is one-step. Although a multistep method may be regarded as a one-step method on a Cartesian product of state spaces ^[28], the theory of multistep methods is complicated by the possibility of uncontrolled growth of the error in the past states ^[25]. This is not the focus of this article; we consider only one-step methods.

To formalize, the time-advanced state $\boldsymbol q^{n+1} \approx {\mathit{\boldsymbol{q}}}(t_{n+1})$ after one step of a one-step implicit method is obtained from the given current state $\boldsymbol q^n \approx {\mathit{\boldsymbol{q}}}(t_n)$ by solving a generally nonlinear equation of the form

where $h$ is the given timestep. For example, the backward Euler method has ${\mathit{\boldsymbol{F}}}({\mathit{\boldsymbol{y}}}; {\mathit{\boldsymbol{x}}}, h) = {\mathit{\boldsymbol{y}}} - {\mathit{\boldsymbol{x}}} - h\, {\mathit{\boldsymbol{f}}}(\mathit{\boldsymbol{y}})$.

We assume, for all $\boldsymbol x$, that

where $\boldsymbol 1$ is the identity matrix and subscripts of $\boldsymbol F$ denote partial derivatives. Then, by the implicit function theorem, for any fixed $\boldsymbol x$, there is a unique smooth solution $\boldsymbol y(h)$ defined for sufficiently small $h$, and $\boldsymbol y = \boldsymbol x+h\boldsymbol f(\boldsymbol x)+O(h^2)$, as is required for consistency.

Two solutions $\boldsymbol y_1(h)$ and $\boldsymbol y_2(h)$, defined on open intervals containing $h = 0$ and satisfying the condition that $\boldsymbol F_{\boldsymbol y}\bigl(\boldsymbol y_i(h); \boldsymbol x, h\bigr)$ is nonsingular, are equal on the intersection of their domains (the set on which $\boldsymbol y_1(h) = \boldsymbol y_2(h)$ is nonempty, closed, and open by the implicit function theorem, and the intersection of intervals is connected). Therefore, there is a maximal such solution, which we call the principal solution branch, and there is generally a critical timestep $h_c$ after which the principal branch ceases to exist. The condition that the solution set of $F(\boldsymbol y; \boldsymbol x, h) = \boldsymbol 0$ is smooth in the space $\{\boldsymbol y, h\}$ is weaker than the condition that it defines $\boldsymbol y$ as a function of $h$. Solutions may be continuously connected after a fold bifurcation, for example, where the solution manifold turns backwards from the direction of increasing timestep ^[26,27]. Continuing through such a bifurcation leads to multiple solutions at smaller timesteps than the critical one at which the bifurcation occurs. These solutions co-exist with the solutions from the principal branch, and they may persist even as the timestep approaches zero.

A simple example demonstrating the existence of a critical timestep $h_c$ at which a fold bifurcation occurs is the application of the backward Euler method to the scalar IVP

With the backward Euler method, the update equation (1.1) can be written as

having solutions

Evidently, there are two solutions for $h < h_c = (4x)^{-1}$, a single solution at $h = h_c$, and no solutions for $h > h_c$; see Figure 1. We note the existence of two solutions for all $0 < h < h_c$. By the Newton–Kantarovich Theorem ^[24], convergence to either solution is possible for an appropriately chosen initial guess.

The principal solution branch contains the initial condition at zero timestep, which is exactly correct, and it contains all solutions near zero timestep that continuously emanate from the initial condition. If a more complicated bifurcation occurs, say a pitchfork, as opposed to a fold (e.g., bifurcations of the solutions of $y^3-(h-h_c)y = 0$ and $y^2-(h-h_c) = 0$, respectively), then there are multiple solutions beyond the critical timestep, none of which can be continuously and monotonically continued back to the initial condition without passage through the bifurcation itself. There may be two successive fold bifurcations, resulting in an "S"-shaped solution manifold, with timesteps larger than $h_c$, for which there is a unique solution. But in that situation, the two additional solutions at timesteps $h\lesssim h_c$ would be rejected in favour of that on the principal branch, and for $h\gtrsim h_c$, the unique solutions are near an inconsistent one and so would also be rejected.

Therefore, if there are multiple solutions at a given timestep smaller than the critical one, then the principal branch should be chosen, and the timestep cannot be validly increased to be larger than the critical one, irrespective of whether or not there are unique solutions there. Simulations should generally not be carried on from solutions that are not on the principal branch, regardless of whether the timestep used is larger or smaller the critical timestep, because they tend to result in non-negligible perturbations to the solution. So as to have a concise term, we make the following definitions.

Definition 1.1. Given an initial condition $\boldsymbol x$, a solution $\boldsymbol y$ with timestep $h$ is called consistent if $(\boldsymbol y, h)$ is on the principal solution branch.

Definition 1.2. The smallest timestep $h$ such that there is a bifurcation of the principal branch is called the critical timestep $h_c$. Timesteps $h > h_c$ are called invalid. By definition, any solutions obtained with invalid timesteps cannot be on the principal branch.

Another basic example relevant to Lagrangian systems that shows the existence of a fold bifurcation for a critical timestep $h_c$ involves the Lagrangian

which models the dynamics of a nonlinear spring that has a stiffness proportional to its linear displacement from equilibrium. Applying the first-order Variational Taylor method (described in Section 3 below) yields an update equation of the form

to be solved for $y$ and $\dot{y}$. The solution to this update equation is

which again shows the potential for multiple solutions. Taking, for example, the initial state $x = 1$, $\dot{x} = 0$, it is straightforward to show that (1.2) has two real solutions for $h < h_c = \sqrt{2+2\sqrt{2}}$, a single solution at $h = h_c$, and no real solutions for $h > h_c$.

2.   Numerical continuation of the implicit solution

We consider one-parameter families of solutions to (1.1) in the hyperplane defined by the vector $({\boldsymbol y}, h)$ for fixed ${\boldsymbol x}$. For practical computation, we parameterize these families by arclength, which monotonically increases throughout the computation.

Suppose we have two points on the solution curve $({\boldsymbol y}_0, h_0)$ and $({\boldsymbol y}_1, h_1)$ with known tangent vector ${\mathit{\boldsymbol{\mathcal T}}}_0$ at the first point. Our goal is to find a next point on the solution curve in the same arclength direction following $({\boldsymbol y}_1, h_1)$. This pseudo-arclength continuation is accomplished in a predictor-corrector fashion. The prediction is performed via a linear approximation to the curve at the point $({\boldsymbol y}_1, h_1)$ and the correction by computing the minimum-norm solution to the system with rank-1 deficiency.

The tangent vector $\mathit{\boldsymbol{\mathcal T}}_1$ at $({\boldsymbol y}_1, h_1)$ satisfies

and determines $\mathit{\boldsymbol{\mathcal T}}_1$ up to a scalar factor. To preserve the direction of the orientation of the branch, we require

Using a natural decomposition $\mathit{\boldsymbol{\mathcal T}}_i = (\mathit{\boldsymbol{\mathcal T}}_i^{(\boldsymbol y)}, {\mathcal T}_i^{(h)})$, $i = 0, 1$, we write the above as a single system

the solution to which uniquely determines $\mathit{\boldsymbol{\mathcal T}}_1$.

Using a step of length $\Delta s$, we create a predictor

to approximate the next point on the curve. The name pseudo-arclength comes from the fact that $\Delta s$ measures arclength along the tangent line.

From the predictor, a Newton-like method is applied to obtain the next solution point on the curve $(\boldsymbol y_2, h_2)$. Our specific choice of algorithm is the Gauss–Newton method. This particular approach for path-following was first used in ^[10]. It is identified as being the Gauss–Newton method in ^[8]. A practical description is given in ^[2] and summarized as follows.

We seek a solution to

such that $\|(\Delta\boldsymbol y, \Delta h)\|$ is minimal. Because $\boldsymbol F$ is nonlinear, we set up the Newton iteration as $(\boldsymbol y^0, h^0) = (\hat{\boldsymbol y}, \hat{h})$, and $\boldsymbol y^{k+1} = \boldsymbol y^k + \Delta\boldsymbol y$, $h^{k+1} = h^k+\Delta h$, where ($\Delta\boldsymbol y$, $\Delta h$) is the minimum-norm solution to

The minimum-norm solution is obtained by solving the matrix system

and constructing

with

Quadratic convergence of the above algorithm is proven in ^[8], and in our experiments, we iterate solutions to a tolerance $\|\boldsymbol F\| < F_{\text{tol}} = 10^{-9}$, where $\| \cdot \| = \| \cdot \|_2$.

An initial step length $\Delta s_0$ must be chosen to initialize the continuation. For subsequent steps along the continuation curve, we choose the next step length according to

where $\epsilon$ is a user-defined tolerance for the absolute error in $(\hat{\boldsymbol y}, \hat h)$; we set $\epsilon = 100\, F_{\text{tol}}$.

3.   The double pendulum

We applied pseudo-arclength continuation as described in Section 2 to the double pendulum system of two bobs connected by a frictionless pin joint, the first moving on a circle with fixed center, the second moving on a circle centered at the first, and both moving in the presence of gravity; see Figure 2. The system is Lagrangian (and so symplectic), with Lagrangian function

in terms of the angular coordinates $(\alpha, \beta)$ shown in Figure 2a, and $g = 9.81$. The system evolution is in accord with Euler–Lagrange differential equations of the form

For initial conditions, we use the near-vertical arrangement as in Dharmaraja et al. ^[1],

We compute a reference solution trajectory, $\boldsymbol q_{\text{ref}}$; details are given at the end of this section. In the resulting motion, the system quickly goes through a fast loop, as seen in the lower right of Figure 2b. Events such as this, where the dynamics change quickly relative to the initial or overall dynamics, seem to readily induce the bifurcations in (1.1) in which we are interested. Accordingly, we chose a point on the trajectory that is just before the loop for our $\boldsymbol x$ value, specifically the point corresponding to time $t = 0.9$.

The fold points we seek are a property of the implicit method used for the integration as defined by $\boldsymbol F(\boldsymbol y; \boldsymbol x, h)$. We now investigate the properties of these fold points for various implicit integration methods.

We first consider a Variational Taylor (VT) method ^[6,7], which uses the first-order Taylor expansion of trajectories $\boldsymbol x \rightarrow \boldsymbol x +t\dot{\boldsymbol x}$ to obtain the discrete Lagrangian,

The resulting one-step method is obtained by finding the critical points of the discrete action $L_h(\boldsymbol y)+L_h(\boldsymbol x)$ subject to the constraints

where $\boldsymbol a_1$ and $\boldsymbol a_2$ are constants, the values of which are not required for the method implementation. The first two of these constraints are discrete analogues of the fixed-endpoint constraint for the continuous variational principle. The last, which connects the discrete trajectories, is a discrete analogue of the continuous constraint that the derivative of configuration is velocity. The resulting constrained optimization problem is equivalent to discrete Euler–Lagrange equations, and, in the case at hand, the associated Lagrange multipliers may be explicitly eliminated to obtain equations of the form (1.1). Such methods extend to any order, but they become much more complicated, and it suffices for our purpose here to consider only the first-order method ^[12,13].

In addition to this, we investigate the fold points of a variety of implicit Runge–Kutta methods of the form

Specifically, using the usual Butcher tableau notation ^[4], we use the methods below (in row order): backward Euler (BE); trapezoidal rule (TR); the trapezoidal rule backward differentiation formula 2 split-step method, with optimized split-step size $\gamma = 2-\sqrt{2}$ as in ^[1] (TRB); third- and fifth-order Radau IIA methods (R3, R5); and fourth- and sixth-order Gauss–Legendre methods (GL4, GL6). In the experiments described below, the nonlinear systems associated with any implicit method were solved using a classical Newton iteration.

A reference solution for this problem was computed using timestep $h_{\text{ref}} = 2\times 10^{-5}$ with the VT method. This reference solution was compared to that obtained from a VT integration with timestep size $h_{\text{ref}}/2$. The state of the system at $t = 2.0$ is given in Table 1, where we see agreement of 8 significant digits between these solutions (Table 1).

To demonstrate the potential effects of convergence to non-principal-branch solutions from an implicit one-step method, we show inconsistent solutions obtained with the VT method and timestep size $h = 0.1225$ in Figure 3. The nonlinear solver used in computing the solution at the next timestep converges in all steps displayed, yet the inconsistent solutions appear markedly different from the reference solution of Figure 2. There is even a difference in trajectories with inconsistent solutions when all that is changed is the initialization of the nonlinear solver, in this case, from the solution at the previous timestep to a linear extrapolation based on the solution at the previous two timesteps; this implies the existence of at least two non-principal solution branches at some point along the time integration.

All of the methods tested exhibit fold points before $h = 0.33$, as seen in Figure 4. Four of the methods (BE, TR}, TRB, and R3) have two fold points, resulting in a continuous connection between the $h = 0$ solution and the $h = 0.35$ solution. Three of the methods (R5, GL4, and GL6) have a single fold point, beyond which we were unable to find other solutions. The solution obtained with the VT method exhibits two separate solution branches. The main branch originating from the $h = 0$ solution folds back, with norm approaching infinity as $h$ returns to zero. A second branch extends beyond the fold point, allowing for solutions for larger $h$ that are not continuously and monotonically connected to the $h = 0$ solution.

The relative error in the trajectories, computed as

is visualized for all solution curves in Figure 5. We notice that the computed solutions near the fold points for each of the methods have relatively large (typically $\mathcal{O}(1)$) errors in this instance and that the methods are out of their regions of asymptotic convergence. Passing through the fold points in arclength, the error generally increases as the timestep decreases.

4.   Discussion

Bifurcations, as a function of the timestep, are to be expected in the solutions of the update equations defined by implicit numerical methods for nonlinear IVPs. As demonstrated, such bifurcations can occur in purely mathematical examples, such as the backward Euler method applied to the differential equation $\dot q = q^2$, as well as in more physically interesting models, such as the double pendulum.

For simulations where the accuracy of single trajectories is important and variable step sizes are used, standard step-size control may reduce the timestep to help prevent convergence to inconsistent solutions. Such convergence, however, may persist even at small timesteps, depending on how the iterative solver for a given method is initialized. Accordingly, convergence of an implicit method solver alone is unreliable for assessing whether or not a solution is consistent or whether a simulation should be carried on. In fact, it may be possible to identify a critical timestep size $h_c$, corresponding to the first point at which a bifurcation of the principal solution branch occurs, such that timesteps $h > h_c$ can be considered to be invalid and solutions obtained are by definition inconsistent.

In other instances, such as the long-time integration of symplectic systems, where timesteps are often constant or more generally where consideration of backward error is a driving motivation, specific trajectories may not be as important as the general behaviour of the system. Even though the $\mathcal{O}(1)$ relative errors we observe in trajectories at fold points in the double pendulum example imply that the solutions may be inconsistent for timesteps smaller than the critical timestep, it is helpful to be aware of the presence of fundamental method- and state-dependent limitations on the size of the valid timesteps. Simulations with timestep sizes that are larger than the critical timestep cannot be expected to generate reasonable and robust results in practice.

We used pseudo-arclength continuation to follow the solution branches of the numerical solution to the double pendulum problem using several common implicit numerical methods. Pseudo-arclength continuation can be used in this fashion as a post hoc solution validation technique. That is, given a simulation, we can use the pseudo-arclength continuation procedure described to verify that the solutions at each timestep lie on the principal solution branch. In principle, this can be done in parallel with the computations of the time-advanced state. It would be interesting to develop such a tool to monitor timesteps of implicit methods and to flag if an inconsistent solution is generated.

Acknowledgements

This work was funded by the National Science and Engineering Research Council of Canada under grant number 228090-2013. The authors thank the anonymous referees for their insightful comments.

Conflict of interest

All authors declare no conflicts of interest in this paper.