A Conditioning Function for the Convergence of Numerical ODE Solvers and Lyapunov’s Theory of Stability

Divakar Viswanath ^{∗}

22 December 1998

Abstract

For the ordinary differential equation (ODE) ẋ(t) = f(t, x), x(0) = x_{0}, t ≥ 0, x ∈ R^{d}, assume f to be at least continuous in t and locally Lipshitz in x, and if necessary, several times continuously differentiable in t and x. We associate a conditioning function E(t) with each solution x(t) which captures the accumulation of global error in a numerical approximation in the following sense: if x̃(t; h) is an approximation derived from a single step method of time step h and order r then ∥x̃(t; h) − x(t)∥ < K(E(t) + ϵ)h^{r} for 0 ≤ t ≤ T, any ϵ > 0, sufficiently small h, and a constant K > 0.

Using techniques from the stability theory of differential equations, this paper gives conditions on x(t) for E(t) to be upper bounded linearly or by a constant for t ≥ 0. More concretely, these techniques give constant or linear bounds on E(t) when x(t) is a trajectory of a dynamical system which falls into a stable, hyperbolic fixed point; or into a stable, hyperbolic cycle; or into a normally hyperbolic and contracting manifold with quasiperiodic flow on the manifold.

1 Introduction

For the system of ordinary differential equations ẋ(t) = f(t, x), t ≥ 0, x ∈ R^{d}, the initial value problem x(t_{0}) = x_{0}, t_{0} ≥ 0, may not have a solution, and when the solution exists it may not be unique. However, if f(t, x) is continuous in t and not only continuous but also locally Lipshitz in x, there is a unique solution of the initial value problem for t_{0} ≤ t < t_{0} + ϵ, ϵ > 0, which we denote by x(t; t_{0}, x_{0}). We make these assumptions about f(t, x) and only consider solutions x(t; t_{0}, x_{0}) that can be continued till t = ∞. Usually t_{0} = 0, and we denote these solutions by x(t; x_{0}) or simply x(t).

Even at this level of generality, one might ask how accurately x(t; x_{0}) can be approximated by a numerical method. A reasonable first guess is to look at the d × d matrix

E′(t) = \frac{∂x(t; x_{0})}{∂x_{0}},

for 0 ≤ t < ∞; since E′(t) gives the sensitivity of x(t; x_{0}) with respect to x_{0}, one might hope to relate it to the accumulation of global error. However, numerical methods introduce discretization error not just at t = 0 but at every time step of the integration, and therefore, E′(t) proves to be an insufficient concept.

∗Research at MSRI is supported in part by NSF grant DMS-9701755

The following conditioning function E(t) associated with x(t; x_{0}), or more briefly x(t), does capture the accumulation of global error:

E(t) = sup_{v(s)}∥∫_{0}^{t} \frac{∂x(t)}{∂x(s)}v(s)ds∥,  (1.1)

with the supremum taken over continuous functions v : [0, t] → R^{d} with ∥v(s)∥ ≤ 1 for 0 ≤ s ≤ t. Unlike E′(t), E(t) takes into account the sensitivity of x(t) with respect to all x(s), 0 ≤ s ≤ t. In the definitions of E(t) and E′(t) above, we have assumed f(t, x) to be continuously differentiable with respect to t and x. All vector norms in this paper are Euclidean norms and all matrix norms are the corresponding induced norms.

Let x̃(t; x_{0}; h) denote the approximation to x(t; x_{0}) computed by a single step method of order r. The numerical method gives x̃(t; x_{0}; h) at all positive integer multiples of h. For kh < τ < (k +1)h, k = 0, 1, 2, …, we define x̃(τ; x_{0}; h) by following the solution exactly from the initial point t = kh and x = x̃(kh; x_{0}; h) till t = τ. Therefore, x̃(t; x_{0}; h) can be discontinuous only at t = kh, k = 1, 2, … The magnitude of this discontinuity

∥x̃((k + 1)h; x_{0}; h) − x((k + 1)h; kh, x̃(kh; x_{0}; h))∥

is the local discretization error at t = (k + 1)h. We write this discretization error as K_{k+1}h^{r+1}v_{k+1}, where v_{k+1} ∈ R^{d} with ∥v_{k+1}∥ = 1 and K_{k+1} is a non-negative real number. Thus the direction of the local discretization error at the ith step is v_{i} and its magnitude is K_{i}h^{r+1}.

Let us assume that the magnitude of the local discretization error is bounded above by Kh^{r+1} for a constant K > 0. This assumption can hold for example for a Runge-Kutta method of order r if f(t, x) is r + 1 times continuously differentiable with respect to t and x; we discuss this assumption further in Section 2. Then by Theorem 4.5, given ϵ > 0 and T > 0,

∥x(t; x_{0}) − x̃(t; x_{0}; h)∥ < (E(t) + ϵ)Kh^{r}  (1.2)

for 0 ≤ t ≤ T and sufficiently small h. In this bound on the global error, E(t), which is given by (1.1), is independent of details of the numerical method. Further, a bound like (1.2) will not hold for any real-valued function of t which is strictly less than E(t) at some value of t. It is because of these reasons that we call E(t) a conditioning function.

Let us now draw an analogy between E(t) and the absolute condition number of a multivariate function g(x). The formula (see [26])

lim_{δ→0^{+}} sup_{δx<δ} \frac{∥g(x + δx) − g(x)∥}{δ}

is the analogue of (1.1); this absolute condition number measures the sensitivity of g(x) with respect to small changes in x. For a stable numerical method, the error in evaluating g(x) is governed by this absolute condition number in a manner similar to the dependence of global error on E(t) given by (1.2).

Sections 2,3, and 4 define E(t) and develop its properties. We do not begin with (1.1) as the definition of E(t). Section 2 introduces a model of discretization errors which is similar to but more general than the model of Stuart and Humphries [25]. Section 3 defines E(t) using this model. The expression for E(t) in (1.1) is derived in Corollary 6.3 of Section 6. Let us mention the formal similarity of (1.1) to bounds on the global error derived by Iserles and Söderlind [16] using the Alexseev-Gröbner lemma. The results in this half of the paper argue that E(t) is the appropriate vehicle for a study of global errors.

From Section 5 onwards, we undertake the task of relating E(t) to stability properties of x(t; x_{0}). The relation of E(t) to the accumulation of global error as t increases is clear enough from (1.2). If E(t) is bounded by a constant or linearly or by a polynomial of low degree in t, the accurate approximation of the trajectory x(t; x_{0}) can be considered a tractable problem; but if E(t) increases exponentially in t accurate approximation of x(t; x_{0}) is pretty intractable. The standard technique of bounding global errors using the Lipshitz constant is used in convergence proofs of numerical methods [9]. But the bounds on E(t) obtained this way increase exponentially in t and are of hardly any other use.

We upper bound E(t) by making stability assumptions on x(t; x_{0}). Sections 5 and 6 derive linear and constant bounds on E(t) by making stability assumptions on x(t; x_{0}). This connection to stability is somewhat subtle; there are exponentially stable examples with exponentially increasing E(t). However, the work of researchers who followed Lyapunov allows us to clarify and circumvent the difficulties in relating E(t) to stability theory. See Table 1 in Section 6 for a summary.

Let us now make some prefatory remarks about the stability theory of ordinary differential equations. The theory of stability of ordinary differential equations was initiated by A.M. Lyapunov in 1892 [17]. One stream of research which emanates from this remarkable work is about first approximations. Let x(t) = x(t; x_{0}) be a solution of ẋ(t) = f(t, x). If the perturbation y(t) is such that x(t) + y(t) is also a solution of the same equation, then y(t) satisfies

ẏ(t) = f(t, y + x(t)) − f(t, x(t)) = F (t, y),

where clearly F (t, 0) ≡ 0. The solution x(t) is stable in the sense of Lyapunov if for every ϵ > 0 there exists a δ > 0 such that ∥y(0)∥ < δ implies ∥y(t)∥ < ϵ for t ≥ 0. Now (see [22], [5]), x(t) is stable if and only if the zero solution of ẏ(t) = F (t, y) is stable. Thus it is enough to look at zero solutions of systems of the form ẏ(t) = F (t, y), with F (t, 0) ≡ 0.

Lyapunov pointed out an other possible simplification. Since F (t, 0) ≡ 0, if F (t, y) is assumed to be continuously differentiable in y, F (t, y) = A(t)y + o(∥y∥) as y → 0. Here,

A(t) = \frac{∂F (t, y)}{∂y}|_{y=0}
= \frac{∂f(t, x)}{∂x}|_{x=x(t;x_{0})}.  (1.3)

Lyapunov argued that perhaps stability properties of the zero solution of the linear first approximation ẏ(t) = A(t)y might imply stability properties of the zero solution of the nonlinear system ẏ(t) = F (t, y). This program has been carried out by Lyapunov, E. Cotton (1911), O. Perron (1928), I.G. Petrovskii (1934), R. Bellman (1953), and others.

The other stream of research, which also originated with Lyapunov, uses Lyapunov functions V (t, y). If every solution y(t; y_{0}) of ẏ(t) = F (t, y), t ≥ 0, had the property that ∥y(t; y_{0})∥ decreases as t increases, stability of the zero solution would be easy to infer. However, this property does not hold even for stable, linear systems of the form ẏ(t) = Ay, where A is a constant matrix. Lyapunov’s idea was to use a non-negative, real-valued function V (t, y), instead of the norm, such that V (t, y(t)) decreases as t increases. This V (t, y) is required to be related to ∥y∥ in a uniform way for t ≥ 0; precise details depend upon the concept of stability. After Lyapunov, this line was greatly developed by Soviet researchers, including Persidskii, Malkin, Krasovskii, and others, beginning in the 1920s. It was taken up by researchers in the west, including J.L. Massera, T. Yoshizawa, and others, in the 1950s.

Parts of Sansone and Conti [22] and Hale [10] are excellent introductions to the theory of stability of ODEs. There are numerous advanced works; of them we refer to Bellman [1], Malkin [18], and Yoshizawa [27].

Obtaining detailed stability information about solutions of ODEs can be far from trivial. For this reason, applying the theory in Sections 5 and 6 to concrete examples is not a simple matter. In Section 7, we give three applications to dynamical systems. We derive constant or linear upper bounds on E(t) for trajectories falling into stable, hyperbolic fixed points, or into stable, hyperbolic cycles, or into a normally hyperbolic and contracting manifold, with the flow on the manifold being quasiperiodic. These three applications involve the Hartman-Grobman theorem, convergence in phase results for stable cycles, and results of Pugh, Shub, and others about normally hyperbolic flows, respectively. The first of these applications was covered by Stuart and Humphries [25]. The second application is a significant extension of a result due to Cano and Sanz-Serna [4]. The third application appears to be entirely new.

Let us mention the similarity of our analysis to the asymptotic analysis of global error of Henrici [12] [13] and several following researchers, including Gragg [8]. Theorem 6.2, in particular, is implicit in Henrici’s work. Beginning with Dahlquist [6], it has been known that the use of one-sided Lipshitz constants can sometimes give meaningful bounds on the global error. We show in Section 8 that the use of one-sided Lipshitz constants fits naturally into our framework. Section 8 also considers variable time stepping, multistep methods, and other issues.

To summarize briefly, this paper consists of three parts. The first part, Sections 2, 3, and 4, derives a conditioning function E(t) and associates it with the global error in numerically approximating the solution of an ordinary differential equation. The second part, Sections 5 and 6, relates E(t) to the stability theory of ordinary differential equations. The third part, Section 7, applies this theory to dynamical systems.

2 A Model for Discretization Errors

The model of discretization error which we now present is close to the way discretization errors are made by single step methods with constant step sizes. Stuart and Humphries [25] model discretization errors of single step methods in a similar manner.

Let α(h) be a continuous, strictly increasing function of h for h ≥ 0. Assume also that α(0) = 0. Then an approximation x̃_{α}(t; x_{0}; h) to x(t; x_{0}) is defined as follows:

x̃_{α}(0; x_{0}; h) = x_{0}
x̃_{α}(nh; x_{0}; h) = x(nh; (n − 1)h, x̃_{α}((n − 1)h; x_{0}; h)) + hα(h)v_{n} n ≥ 1,  (2.1)

where v_{n} ∈ R^{d} can be any vector with ∥v_{n}∥ ≤ 1. In words, the approximate solution at t = nh, n ≥ 1, is obtained by exactly propagating the point x̃_{α}((n − 1)h; x_{0}; h) at t = (n − 1)h under ẋ(t) = f(t, x) till t = nh, and then adding the discretization error or discontinuity hα(h)v_{n}, where ∥v_{n}∥ ≤ 1. For (n − 1)h ≤ t < nh,

x̃_{α}(t; x_{0}; h) = x(t; (n − 1)h, x̃_{α}((n − 1)h; x_{0}; h)).  (2.2)

Since v_{n} can be any vector with ∥v_{n}∥ ≤ 1, this actually defines a whole family of approximate solutions which we denote by X̃_{α}(x_{0}; h). The exact solution x(t) is the only member of this family which is continuous. Figure 1 gives an example of an approximate solution.

We now address how single step numerical methods are related to approximations from the family X̃(x_{0}; h). In our description of single step methods, the discretization error at every step was taken to be K_{i}h^{r+1}v_{i}. We now assume that the K_{i} are bounded by a constant K which does not depend upon h or i. This can be proven in some circumstances; see [25]. Besides, if there is no such K, the numerical method will not in practice behave as if it were of order r. Thus Figure 1: The thick lines show an exact solution of the equation ẋ(t) = sin t + (cos t)x and an approximation to it with h = 0.5 and α(h) = 2h for.

when the order of accuracy of the numerical method for approximating x(t; x_{0}) is r, we can take α(h) = Kh^{r}, and there will always be an approximation in the family X̃_{α}(x_{0}; h) which is the same as the trajectory of the numerical method.

The notations we have established so far will be adhered to in the rest of the paper. For the system, ẋ(t) = f(t, x), f(t, 0) is not necessarily zero. The solution with x(t_{0}) = x_{0} is denoted by x(t; t_{0}, x_{0}), and by x(t; x_{0}) when t_{0} = 0. The approximations to x(t; x_{0}) obtained as in (2.1) and (2.2) are denoted by x̃_{α}(t; x_{0}; h). The family of approximations is X̃(x_{0}; h). The same notations apply for ẏ(t) = F (t, y) with the xs changed to ys. But here F (t, 0) ≡ 0, and ỹ_{α}(t; h), t ≥ 0, is an approximation to the zero solution, and Ỹ_{α}(h) is a collection of those approximations. When we speak of a linear first approximation, it is always obtained as in (1.3). Sometimes we omit the word linear. This same equation is sometimes called the equation of first variation or simply the linearization. Let us note that when we speak of stability, it is always stability in the sense of Lyapunov and his followers; in particular, it is not numerical stability in the sense of Dahlquist.

3 Definition of E_{α}(t)

Let us recall that α(h) is assumed to be a strictly increasing, continuous function of h for h ≥ 0 which is zero at zero; for example, α(h) can be Kh^{r} for a positive integer r. Assume the functions f(t, x) and F (t, x) to be defined for 0 ≤ t < ∞ and x ∈ R^{d}, to be continuous in t, and locally Lipshitz in x. Besides, F (t, 0) ≡ 0. Assume the solution x(t; x_{0}) of the initial value problem ẋ(t) = f(t, x), x(0) = x_{0}, to be continuable till t = ∞. Only some of these assumptions are restated in theorems that follow.

The global error e_{α}(t; x_{0}; h) is defined as follows:

e_{α}(t; x_{0}; h) = sup_{x̃_{α}∈X̃_{α}}∥x̃_{α}(t; x_{0}; h) − x(t; x_{0})∥.

We later use this to define E_{α}(t).

Instead of saying that x̃_{α}(t; x_{0}; h) exists for 0 ≤ t ≤ T, we say that x̃_{α}(t; x_{0}; h) can be continued till t = T. We say that every approximation x̃_{α}(t; x_{0}; h) can be continued till T if every x̃_{α}(t; x_{0}; h) exists for 0 ≤ t ≤ T with any allowed choice of discontinuities at t = kh. Lemma 3.1 introduces h_{0}(T, r) and L(T, r).

Lemma 3.1. Assume, as usual, that f(t, x) is continuous in t, t ≥ 0, and locally Lipshitz in x, x ∈ R^{d}. Let x(t; x_{0}), t ≥ 0, be the unique solution of the initial value problem ẋ(t) = f(t, x), x(0) = x_{0}. Then,

(i) There exists a constant L(T, r) > 0 such that

∥f(t, x_{1}) − f(t, x_{2})∥ ≤ L(t, r)∥x_{1} − x_{2}∥

for 0 ≤ t ≤ T and ∥x_{i} − x(t; x_{0})∥ ≤ r, for any T > 0, r > 0, and i = 1, 2.

(ii) There exists a constant h_{0}(T, r) such that for 0 < h < h_{0}(T, r) every approximation x̃_{α}(t; x_{0}; h) can be continued till t = T, and satisfies ∥x̃_{α}(t; x_{0}; h) − x(t; x_{0}∥ < r, for 0 ≤ t ≤ T.

Further, for 0 < h < h_{0}(T, r), e_{α}(t; x_{0}; h) ≤ te^{L(T,r)t}α(h).

Proof. Similar to proof of Theorem 3.4.6 in Stuart and Humphries [25]. Similar estimates using the Lipshitz constant are found in [9] and other places.

For the zero solution of ẏ(t) = F (t, y), y(0) = 0, where F (t, 0) ≡ 0, e_{α}(t; h) is defined as follows:

e_{α}(t; h) = sup_{ỹ_{α}∈Ỹ_{α}}∥ỹ_{α}(t; h)∥,  (3.1)

where Ỹ_{α}(h) is the family of approximations to the zero solution with time step h.

Proposition 3.2. Assume h < h_{0}(T, r) for r > 0, T > 0, and that 0 ≤ t ≤ T. The e_{α}(t; x_{0}; h) of the solution x(t; x_{0}) of ẋ(t) = f(t, x), x(0) = x_{0}, and the e_{α}(t; h) of the zero solution of ẏ(t) = F (t, y), y(0) = 0, where F (t, y) = f(t, y + x(t; x_{0})) − f(t, x(t; x_{0})), are the same.

Proof. It is enough to show that members x̃_{α}(t; x_{0}; h) of X̃_{α}(x_{0}; h) and members ỹ_{α}(t; h) of Ỹ_{α}(h) can be matched so that x̃_{α}(t; x_{0}; h) = x(t; x_{0}) + ỹ_{α}(t; h).

It is well known (see [5]) that

x(t + τ; t; x(t; x_{0}) + δ) = x(t + τ; x_{0}) + y(t + τ; t, δ)

for τ ≥ 0. Hence, if x̃_{α}(kh; x_{0}; h) = x(kh; x_{0}) + ỹ_{α}(kh; h) then x̃_{α}(t; x_{0}; h) = x(t; x_{0}) + ỹ_{α}(t; h) for kh < t < (k + 1)h. Further, the discontinuities of x̃_{α}(t; x_{0}; h) and ỹ_{α}(t; h) at t = kh, where k is a positive integer, can be exactly matched. Therefore, the approximate solutions can be matched as desired.

The E_{α}(t) which corresponds to the zero solution of ẏ(t) = F (t, y), y(0) = 0, t ≥ 0, is defined as follows:

E_{α}(t) = lim sup_{h→0} \frac{e_{α}(t; h)}{α(h)}.  (3.2)

By Lemma 3.1, E_{α}(t) ≤ te^{Lt}, where L = L(T, r) for some T > t and r > 0. For a nonzero solution x(t; x_{0}) of ẋ(t) = f(t, x), x(0) = x_{0}, t ≥ 0, E_{α}(t) is defined to be the same as the E_{α}(t) for the zero solution of ẏ(t) = f(t, y + x(t; x_{0})) − f(t, x(t; x_{0})), y(0) = 0, t ≥ 0.

As in the stability theory of ODEs, we can and do confine ourselves to an analysis of the zero solution of ẏ(t) = F (t, y) without any loss of generality. From here on, assume L(T, r) to be such that

∥F (t, y_{1}) − F (t, y_{2})∥ ≤ L(t, r)∥y_{1} − y_{2}∥

for 0 ≤ t ≤ T, ∥y_{1}∥ ≤ r, ∥y_{2}∥ ≤ r, where T > 0 and r > 0. Also, assume h_{0}(T, r) to be such that for 0 < h < h_{0}(T, r) every approximation ỹ_{α}(t; h) can be continued till t = T and satisfies ∥ỹ_{α}(t; h)∥ < r.

4 Properties of E_{α}(t)

This section is devoted to properties of E_{α}(t). We show that E_{α}(t) is continuous (Theorem 4.3) and independent of the choice of α(h) (Theorem 4.4) so that it can be written as E(t). Theorem 4.5 and its corollary relate E(t) to the accumulation of global error when x(t; x_{0}) is approximated.

In the proofs in this section, we repeatedly use Theorem 10.1 of [9]. That theorem allows us to bound the divergence of two solutions of a differential equation using a Lipshitz constant.

A sequence of inequalities are often combined in a way we now describe. Let e_{0} = 0, and e_{i} ≤ fe_{i−1} + r_{i−1} for 1 ≤ i ≤ n. Then, we have,

e_{n} ≤ (f^{n−1}r_{0} + f^{n−2}r_{1} + ⋯ + r_{n−1}),

assuming f ≥ 0. In fact, most of the time r_{i} will all be equal. In that situation, e_{n} ≤ r_{0}(1 + f + ⋯ + f^{n−1}). And when f ≥ 1, we can get e_{n} ≤ f^{n−1}(r_{0} + ⋯ + r_{n−1}), or e_{n} ≤ nf^{n−1}r_{0} if the r_{i} are all equal.

Lemma 4.1. Let h < h_{0}(t + s, r), r > 0. Then

e_{α}(t + s; h) ≤ e^{Ls}e_{α}(t; h) + (s + h)e^{Ls}α(h),

where L = L(t + s, r) and s ≥ 0.

Proof. Let ỹ_{α}(t; h) be an approximation to the zero solution of ẏ(t) = F (t, y), y(0) = 0.

Let t + h_{l} be the first multiple of h after t, t + s − h_{r} the last multiple before t + s, and let t + s − h_{r} = t + h_{l} + nh. Using Theorem 10.1 of [9] and taking into account the discontinuities at multiples of h, we have

∥ỹ_{α}(t + h_{l}; h)∥ ≤ e^{Lh_{l}}∥ỹ_{α}(t; h)∥ + hα(h),
∥ỹ_{α}(t + h_{l} + kh; h)∥ ≤ e^{Lh}∥ỹ_{α}(t + h_{l} + (k − 1)h; h)∥ + hα(h), 1 ≤ k ≤ n,
∥ỹ_{α}(t + s; h)∥ ≤ e^{Lh_{r}}∥ỹ_{α}(t + nh; h)∥.

Combining these inequalities, we get

∥ỹ_{α}(t + s; h)∥ ≤ e^{Ls}∥ỹ_{α}(t; h)∥ + (s + h)e^{Ls}α(h).

The proof is easy to complete using (3.1).

Lemma 4.2. With the same assumptions about h and L as in the previous lemma,

e_{α}(t + s; h) ≥ e^{−Ls}e_{α}(t; h),

where s ≥ 0.

Proof. Let ỹ_{α}(τ; h), 0 ≤ τ ≤ t, be an approximate solution. Continue it from τ = t to τ = t + s exactly as a solution of ẏ(t) = F (t, y). Then,

∥ỹ_{α}(t + s; h)∥ ≥ e^{−Ls}∥ỹ_{α}(t; h)∥.

The proof is now easy to complete using (3.1).

Theorem 4.3. The E_{α}(t) of the zero solution of ẏ(t) = F (t, y), y(0) = 0, is continuous for t ≥ 0.

Proof. From Lemmas 4.1 and 4.2, we get

e_{α}(t − δ; h)e^{−Lδ} ≤ e_{α}(t; h) ≤ e^{Lδ}e_{α}(t − δ; h) + e^{Lδ}(δ + h)α(h),

when 0 < t − δ < t, and when t < t + δ,

(e_{α}(t + δ; h) − e^{Lδ}(δ + h)α(h))e^{−Lδ} ≤ e_{α}(t; h) ≤ e^{Lδ}e_{α}(t + δ; h).

Divide these two inequalities by α(h) and use (3.2) to deduce continuity of E_{α}(t).

Theorem 4.4. The E_{α}(t) of the zero solution of ẏ(t) = F (t, y), y(0) = 0, is the same for any α(h), with α(0) = 0 and α(h) continuous and strictly increasing for h ≥ 0.

Proof. Let β(h) be another function like α(h), which is continuously increasing for h ≥ 0 and zero at zero. We will show that E_{β}(t) = E_{α}(t).

Take L = L(t, r) for some r > 0. Take h_{α} and h_{β} small enough that the approximate solutions in the families Ỹ_{α}(h) and Ỹ_{β}(h) are not only continuable till t but stay within a radius r of zero. For ϵ small enough, there exist unique h_{α} and h_{β} such that α(h_{α}) = β(h_{β}) = ϵ. Assume this relationship between h_{α}, h_{β}, and ϵ. Then h_{α} → 0 and h_{β} → 0 as ϵ → 0.

It is enough to show that

lim_{ϵ→0} \frac{|e_{α}(t; h_{α}) − e_{β}(t; h_{β})|}{ϵ} = 0.

In fact, it is enough to show that for every approximate solution ỹ_{α}(t; h_{α}) there is an approximate solution ỹ_{β}(t; h_{β}) such that

∥ỹ_{α}(t; h_{α}) − ỹ_{β}(t; h_{β})∥ ≤ ϵf(t)g(ϵ),  (4.1)

where f(t) is finite and g(ϵ) → 0 as ϵ → 0. It will also have to be shown that a given ỹ_{β}(t; h_{β}) can be approximated by some ỹ_{α}(t; h_{α}) in the same manner; but the proof of this is gotten by transposing α and β.

Given a ỹ_{α}(t; h_{α}), we construct a ỹ_{β}(t; h_{β}) as follows. ỹ_{β}(τ; h_{β}) is determined by the differential equation ẏ(t) = F (t, y), except when τ = ph_{β}, p being a positive integer. At points τ = ph_{β}, introduce a discontinuity of magnitude smaller than or equal to h_{β}β(h_{β}) in such a way as to get as close to ỹ_{α}(ph_{β}; h_{α}) as possible.

Let h_{β} = (m + f)h_{α}, m ≥ 0 being an integer, and 0 ≤ f < 1. Assume that the number of integer multiples of h_{α} in (ph_{β}, (p + 1)h_{β}] is k_{p}. Let e_{p} = ∥ỹ_{α}(ph_{β}; h_{α}) − ỹ_{β}(ph_{β}; h_{β})∥.

Then e_{0} = 0 and

e_{p+1} ≤ e^{h_{β}L}e_{p} + k_{p}h_{α}α(h_{α})e^{h_{β}L} − h_{β}β(h_{β})
= e^{h_{β}L}e_{p} + k_{p}h_{α}ϵe^{h_{β}L} − h_{β}ϵ.

The first two terms can be derived from Theorem 10.1 of [9]. The third term is because of the discontinuity introduced into ỹ_{β}(τ; h_{β}) at τ = ph_{β}. If this bound on e_{p+1} is negative, we can take e_{p+1} = 0. If n = ⌊\frac{t}{h_{β}} ⌋, and there are k_{n} multiples of h_{α} in (nh_{β}, t], then

∥ỹ_{α}(t; h_{α}) − ỹ_{β}(t; h_{β})∥ ≤ e^{Lh_{r}}e_{n} + k_{n}h_{α}α(h_{α})e^{Lh_{r}},

where h_{r} = t − nh_{β}.

Let p be the highest among {0, 1, …, n} such that e_{p} = 0. Combining the inequalities for e_{i}, p ≤ i ≤ n, we get

e_{n} ≤ e^{L(n−p)h_{β}}(h_{α}ϵe^{h_{β}L}(k_{p} + ⋯ + k_{n−1}) − (n − p)h_{β}ϵ).

But k_{p} + ⋯ + k_{n−1} being the number of multiples of h_{α} in (ph_{β}, nh_{β}] is bounded above by (n − p)h_{β}/h_{α} + 1. Let t′ = t − h_{r}. Then,

e_{n} ≤ e^{Lt′}((n − p)h_{β}e^{h_{β}L}ϵ + h_{α}ϵe^{h_{β}L} − (n − p)h_{β}ϵ)
≤ max(t′e^{Lt′}, e^{Lt′})((e^{h_{β}L} − 1) + h_{α}e^{h_{β}L})ϵ.

Finally,

∥ỹ_{α}(t; h_{α}) − ỹ_{β}(t; h_{β})∥ ≤ max(e^{Lt}, te^{Lt})((e^{h_{β}L} − 1) + h_{α}e^{h_{β}L} + k_{n}h_{α})ϵ.

Since k_{n} is the number of multiples of h_{α} in (t − h_{r}, t], k_{n}h_{α} < h_{β} + h_{α}. Thus the bound above satisfies all conditions required of (4.1), and the proof is complete.

From here on, we drop the subscript α from E_{α}(t), ỹ_{α}(t; h), Ỹ_{α}(h), x̃_{α}(t; x_{0}, h), X̃_{α}(x_{0}; h), and e_{α}(t; h).

Theorem 4.5. Let E(t) be associated with the zero solution of ẏ(t) = F (t, y), y(0) = 0. Given T > 0 and ϵ > 0, there exists h_{0} > 0 such that h < h_{0} implies

sup_{ỹ∈Ỹ}∥ỹ(t; h)∥ = e(t; h) ≤ (E(t) + ϵ)α(h)

for 0 ≤ t ≤ T.

Further, if θ < E(t) there exists an ϵ > 0 such that e(t; h) > (θ + ϵ)α(h) for arbitrarily small h.

Proof. Let ϵ(t; h) = e(t; h)/α(h) − E(t).

Take h < h_{0}, where h_{0} to begin with is smaller than h_{0}(T, r) for some r > 0, and L = L(T, r). By definition of E(t) (3.2),

lim sup_{h→0} ϵ(t, h) = 0  (4.2)

for 0 ≤ t ≤ T.

To prove the first part, given ϵ > 0 we find an h_{0} > 0 such that h < h_{0} implies ϵ(t; h) < ϵ for 0 ≤ t ≤ T. In fact, it is enough to show that there exists an h_{t} > 0, which depends on t, and an open neighbourhood of t in [0, T] such that for τ in that neighbourhood and h < h_{t}, ϵ(τ; h) < ϵ. If the neighbourhoods of t_{1}, …, t_{n} give a finite cover of the compact interval [0, T], we can take h_{0} to be the minimum of h_{t_{1}}, …, h_{t_{n}}.

Clearly,

ϵ(t + δt; h) = ϵ(t; h) + \frac{e(t + δt; h) − e(t; h)}{α(h)} − (E(t + δt) − E(t))

≤ ϵ(t; h) + |\frac{e(t + δt; h) − e(t; h)}{α(h)}| + |E(t + δt) − E(t)|.

By (4.2), ϵ(t; h) < ϵ_{1} for h small enough. By continuity of E(t), |E(t + δt) − E(t)| < ϵ_{2} for |δt| small enough.

Using Lemmas 4.1 and 4.2, and the bound tE^{Lt}α(h) on e(t; h) from Lemma 3.1, we get

\frac{|e(t + δt; h) − e(t; h)|}{α(h)} ≤ max(te^{Lt}(1 − e^{−Lδt}), te^{Lt}(e^{Lδt} − 1) + (δt + h)e^{Lδt}).

for δt > 0, and for δt < 0,

\frac{|e(t + δt; h) − e(t; h)|}{α(h)} ≤ max(te^{Lt}(e^{−Lδt} − 1), te^{Lt}(1 − e^{Lδt}) + (δt + h)).

Therefore, we can choose |δt| and h small enough to ensure |e(t + δ; h) − e(t; h)|/α(h) < ϵ_{3}. We need only take ϵ_{1} + ϵ_{2} + ϵ_{3} < ϵ to find the desired h_{t} and the neighbourhood of t.

The second half of the theorem is immediate from the definition of E(t).

Corollary 4.6. Let E(t) be associated with the solution x(t; x_{0}) of ẋ(t) = f(t, x), x(0) = x_{0}, where f(t, x) is continuous in t and locally Lipshitz in x in an open neighbourhood of the solution (t, x(t; x_{0})), t ≥ 0. Given T > 0 and ϵ > 0, there exists an h_{0} > 0 such that

∥x̃(t; x_{0}; h) − x(t; x_{0})∥ < (E(t) + ϵ)α(h)

for 0 ≤ t ≤ T, h < h_{0}, and any approximation x̃(t; x_{0}; h) in the family X̃(x_{0}; h).

So far, we have defined the notion of an approximate solution (see Figure 1), and an E(t) for every solution of an ODE. Theorem 4.5 and its Corollary 4.6 relate E(t) to the accumulation of global error. In the upper bound (E(t) + ϵ)α(h) for the global error, the details of the numerical method have been pushed into α(h); the conditioning function E(t) is something intrinsic to the solution x(t; x_{0}). Figure 2 shows E(t) for some simple examples.

As shown in Figure 2a, the bound on global error given by E(t) may be pessimistic in practice. Our model of approximations allows arbitrary discontinuities whose norms are bounded by hα(h), while for any given numerical method the discretization errors are fixed. Our analysis does not say how small an h is good enough for the bound on the global error to be governed by E(t) as in Theorem 4.5. Will an h in a practical computation be small enough? We answer this in Henrici’s words [13]: usually if a stepsize is small enough to yield an accurate solution, it is also small enough that an asymptotic formula gives a correct indication of the size of the error. Numerical experiments in [3] and [4] lend ample support to this statement.

5 E(t) for Linear Systems

Our investigation of the relationship between E(t) and stability properties of the exact solution begins with the linear system ẏ(t) = A(t)y, y(0) = 0. The relationship is not as simple as one might wish. There are both asymptotically stable examples with exponentially increasing E(t) and

(a) (b)

Figure 2: (a) The thick line is the E(t) for ẏ(t) = −(2(t + 1) + (t + 1)^{−1})y, y(0) = 1. The two thin lines below it are global errors of forward and backward Euler divided by 8h; here h = .01. (b) E(t) for ẏ(t) = λy. The exact E(t) are obtained using Theorem 5.1.

unstable examples with linearly bounded E(t). However, Theorems 5.5 and 5.6 give conditions for E(t) to be bounded by a constant or to be linearly bounded. Linear systems are an important class of problems by themselves. What is more, they can be used for understanding the E(t) of nonlinear systems. We will see this in Theorem 6.2 and in Section 6.2.

In ẏ(t) = A(t)y, the d × d matrix A(t) is assumed to be continuous for t ≥ 0. For such linear systems, it is easy to show that y(t) is a linear function of y(0) for t ≥ 0 [22]. Thus y(t) = Y (t)y(0) for Y (t) ∈ R^{d,d}, where Y (t) is called the principal fundamental matrix of the linear system. The matrix Y (t) itself is continuous in t and always nonsingular. Moreover, y(t) = Y (t)Y^{−1}(s)y(s) for t ≥ s. Although it does not appear to be directly related, let us mention the beautiful derivation of the Magnus series for Y (t) by Iserles and Nørsett [15].

For scalar linear systems ẏ(t) = a(t)y, a(t) ∈ R, we have the following theorem.

Theorem 5.1. The E(t) for the zero solution of ẏ(t) = a(t)y, y(0) = 0, is given by

E(t) = e^{g(t)} ∫_{0}^{t} e^{−g(s)}ds,

where

g(t) = ∫_{0}^{t} a(τ)dτ.

Proof. The fundamental matrix, which is scalar in this situation, is given by Y (t) = e^{g(t)}. Since Y (t) is always positive, the optimal choice of v(s) in Theorem 5.2 is v(s) ≡ 1.

Let us now consider the concepts of stability put forward by Lyapunov [17]. The definitions will be stated for nonzero solutions x(t; x_{0}).

Definition 1. The solution x(t; x_{0}) of ẋ(t) = f(t, x), x(0) = x_{0}, is stable if given any ϵ > 0 there exists a δ > 0 such that ∥x′_{0} − x_{0}∥ < δ implies ∥x(t; x′_{0}) − x(t; x_{0})∥ < ϵ for t ≥ 0. In fact, stability implies that given ϵ > 0, there exists a δ(τ) > 0 for every τ ≥ 0 such that ∥x(τ; x′_{0})−x(τ; x_{0})∥ < δ(τ) implies that ∥x(t; x′_{0}) − x(t; x_{0})∥ < ϵ for t ≥ τ. Let us emphasize that δ(τ) can depend on τ.

Definition 2. The solution x(t; x_{0}) is asymptotically stable if given ϵ > 0 there exists a δ(τ) > 0 for every τ ≥ 0 such that ∥x′_{τ} − x(τ; x_{0})∥ < δ(τ) implies not only that ∥x(t; τ, x′_{τ}) − x(t; x_{0})∥ < ϵ for t ≥ τ but also that ∥x(t; τ, x′_{τ}) − x(t; x_{0})∥ → 0 as t → ∞.

Implicit in the definitions is an assumption about the existence of solutions which begin near the solution x(t; x_{0}). Obviously, asymptotic stability implies stability. For the scalar, linear problem ẏ(t) = a(t)y, y(0) = y_{0}, a necessary and sufficient condition for asymptotic stability is g(t) → −∞ as t → ∞, where g(t) = ∫_{0}^{t} a(s)ds. However, the following examples show that both these concepts of stability are insufficient for bounding E(t).

Example 5.1. We will show using the scalar, linear problem ẏ(t) = a(t)y that even when a solution is asymptotically stable the E(t) associated with it can increase at an arbitrarily high rate. Given a rate r(t), consider a continuously differentiable function g(t), t ≥ 0, g(0) = 0, such that

1. g(t) ≤ −t for all t ≥ 0,

2. e^{g(k)} ∫_{0}^{k} e^{−g(s)}ds > r(k) for k = 1, 2, 3, …

For the linear system, take a(t) = g′(t). The first condition ensures asymptotic stability of the zero solution, and the second condition implies E(k) > r(k) for positive integers k. Such a g(t) is easy to construct. Take g(k) = −k for k = 0, 1, 2, . . .. For k − 1 < t < k, k ≥ 1, define g(t) so that g(t) ≤ −t and

∫_{k−1}^{k} e^{−g(s)}ds ≥ r(k)e^{k}.

This can be carried out for any continuous r(t), for example r(t) = e^{t}.

Example 5.2. On the other hand, there are unstable solutions with linearly bounded E(t). Consider the scalar, linear system ẏ(t) = \frac{α}{t+1}y, t ≥ 0. For this ODE, y(t) = (1 + t)^{α}y(0) implying instability of the zero solution for α > 0. Yet, for 0 < α < 1 E(t), which is (1−α)^{−1}(1+t)^{α}((1+t)^{1−α}−1), is linearly bounded. For α = 1, E(t) is (1 + t) log(1 + t).

Example 5.3. In Example 5.1, |a(t)| will be unbounded. Does asymptotic stability of ẏ(t) = a(t)y, t ≥ 0, imply a linear bound for E(t) if |a(t)| is bounded? The answer is no; E(t) can still increase exponentially in t. We sketch the construction of a g(t) to show this. First take g_{1}(t) = −t and g_{2}(t) = −2t. Take g(t) = g_{2}(t) for 0 ≤ t ≤ t_{1}, and let g(t) increase monotonically till g(t_{2}) = g_{1}(t_{2}) for t_{2} ≥ t_{1}, and then let g(t) decrease monotonically till g(t_{3}) = g_{2}(t_{3}) for t_{3} ≥ t_{2}. Repeat the same construction from t_{3} onwards with t_{4}, t_{5}, and t_{6} in place of t_{1}, t_{2}, t_{3}, and so on. The construction may be arranged so that

1. if g(τ) = g_{1}(τ) then g(t) = g_{2}(t) for f_{1}τ ≤ t ≤ f_{2}τ for any fixed 0 < f_{1} < f_{2} < 1,

2. |a(t)| = |g′(t)| is bounded.

It is easy to check that E(τ) ≥ e^{−τ/2}(e^{2f_{2}τ} − e^{2f_{1}τ}), for τ such that g(τ) = g_{1}(τ). Therefore, for f_{2} > 1/2, E(t) increases exponentially. Let us note that the linear system in this example has a negative Lyapunov exponent of −1.

Notions of stability needed for bounding E(t) either by a constant or linearly are introduced after Theorem 5.2. Theorem 5.2 is comparable to Theorem 3.1 of Cano and Sanz-Serna [4].

Theorem 5.2. The E(t) of the zero solution of ẏ(t) = A(t)y, y(0) = 0, is given by

E(t) = sup_{v(s)}∥∫_{0}^{t} Y (t)Y^{−1}(s)v(s)ds∥,

where the supremum is over all continuous functions v(s) : [0, t] → R^{d} with ∥v(s)∥ ≤ 1 for 0 ≤ s ≤ t. As before, Y (t) is the principal fundamental matrix of ẏ(t) = A(t)y and A(t) is continuous.

Proof. As usual, let the discontinuity of ỹ(t; h) at t = kh be hα(h)v_{k}. Let n = ⌊t/h⌋ and h_{r} = t−nh. Then

ỹ(kh; h) = Y (kh)Y^{−1}((k − 1)h; h)ỹ((k − 1)h; h) + hα(h)v_{k},

for 1 ≤ k ≤ n, and ỹ(t; h) = Y (t)Y^{−1}(nh)ỹ(nh; h). Combine these equalities to get,

ỹ(t; h) = ∑_{k=1}^{n} Y (t)Y^{−1}(kh)hα(h)v_{k},  (5.1)

This expression is also used in the proof of Theorem 6.2.

To prove,

sup_{v(s)}∥∫_{0}^{t} Y (t)Y^{−1}(s)v(s)ds∥ ≤ E(t),  (5.2)

it is enough to find a ỹ(t; h) for every v(s) such that

∫_{0}^{t} Y (t)Y^{−1}(s)v(s)ds = \frac{ỹ(t; h)}{α(h)} + η(h),  (5.3)

where η(h) → 0 as h → 0. Take the discontinuity of ỹ(t; h) at t = kh to be hα(h)v(kh). Use (5.1) and it is apparent that ỹ(t; h)/α(h) approximates the Riemann integral ∫_{0}^{t} Y (t)Y^{−1}(s)ds. The standard convergence theorem for Riemann integrals gives η(h) → 0 as h → 0.

To show (5.2) in the reverse direction, it is enough to find a v(s) for a given ỹ(t; h) so that (5.3) holds. First, consider a discontinuous ṽ(s) defined by ṽ(s) = v_{k} for kh ≤ s < (k + 1)h. Clearly, ∥ṽ(s)∥ ≤ 1. Then,

∫_{0}^{t} Y (t)Y^{−1}(s)ṽ(s)ds = \frac{ỹ(t; h)}{α(h)} + δ(h)t,  (5.4)

where δ(h) = sup_{0≤s≤t}∥Y (t)Y^{−1}(s + h) − Y (t)Y^{−1}(s)∥. Since Y (s) is continuous, δ(h) → 0 as h → 0. Lusin’s theorem [21] is a basic result for approximating measurable functions by continuous functions. It guarantees a continuous v(s) such that ∥v(s)∥ ≤ ∥ṽ(s)∥ and μ(v(s) ≠ ṽ(s)) < h, where μ is the Lebesgue measure. If M = sup_{0≤s≤t}∥Y (t)Y^{−1}(s)∥, then

∫_{0}^{t} Y (t)Y^{−1}(s)v(s)ds = ∫_{0}^{t} Y (t)Y^{−1}(s)ṽ(s)ds + 2mh, |m| ≤ 2M.  (5.5)

Together, (5.4) and (5.5) complete the proof.

Corollary 5.3.

E(t) ≤ ∫_{0}^{t}∥Y (t)Y^{−1}(s)∥ds.

Corollary 5.4. For 0 ≤ δ ≤ t ,

E(t) ≥ ∥∫_{0}^{δ} Y (t)Y^{−1}(s)ds∥.

The definitions of uniform stability and uniform asymptotic stability that follow seem to have been introduced by Malkin [18]. Theorems which deduce the stability of a nonlinear system from its linear first approximation usually (always?) assume the linear first approximation to be uniformly stable or uniformly asymptotically stable [22]. The uniformity assumptions are not explicitly stated sometimes, for example in [1] and [5]. In these cases, the A(t) in ẏ(t) = A(t)y is either constant or periodic, which means that stability implies uniform stability and asymptotic stability implies uniform asymptotic stability. Uniformity assumptions are natural in the theory of Lyapunov functions as well [27]. We will find them useful for bounding E(t).

Definition 3. The solution x(t; x_{0}) of ẋ(t) = f(t, x) is uniformly stable if for every ϵ > 0 there exists a δ > 0 such that ∥x(τ; x_{0}) − x′_{τ} ∥ < δ for τ ≥ 0 implies ∥x(t; x_{0}) − x(t; τ, x′_{τ})∥ < ϵ for t ≥ τ.

Definition 4. The solution x(t; x_{0}) is uniformly asymptotically stable if it is uniformly stable and the choice of δ in the previous definition can be made in such a way that ∥x(t; x_{0})−x(t; τ, x′_{τ})∥ → 0 as τ → ∞ in a uniform way; i.e. given ϵ′ > 0 there exists T_{ϵ′} such that ∥x(t; x_{0}) − x(t; τ, x′_{τ})∥ < ϵ′ for all t > τ + T_{ϵ′} and x′_{τ} satisfying ∥x(τ; x_{0}) − x′_{τ} ∥ < δ.

Theorem 5.5. If the zero solution of ẏ(t) = A(t)y, y(0) = 0, is uniformly stable, its E(t) is linearly bounded; i.e., E(t) ≤ Kt for some K > 0 and 0 ≤ t < ∞.

Proof. Uniform stability of the zero solution is equivalent to boundedness of ∥Y (t)Y^{−1}(s)∥ for t ≥ s ≥ 0 [22] [27]. If ∥Y (t)Y^{−1}(s)∥ ≤ K for t ≥ s ≥ 0, Corollary 5.3 implies E(t) ≤ Kt.

Theorem 5.6. If the zero solution of ẏ(t) = A(t)y, y(0) = 0, is uniformly asymptotically stable, its E(t) is bounded by a constant; i.e. E(t) < K for some K > 0 and 0 ≤ t < ∞.

Proof. Uniform asymptotic stability of the zero solution is equivalent to ∥Y (t)Y^{−1}(s)∥ < Me^{−ν(t−s)} for ν > 0, M > 0, and t ≥ s ≥ 0 [22] [27]. Again, we can use Corollary 5.3 to complete the proof.

Theorem 5.5 implies that E(t) for the solution of ẋ(t) = Ax, x(0) = x_{0}, is linearly bounded if all the eigenvalues of A have negative or zero real parts, and the ones with zero real part are simple. If all the eigenvalues of A have strictly negative real parts then, in fact, E(t) is bounded by a constant by Theorem 5.6. The necessary stability properties of the zero solution of ẋ(t) = Ax are verified in numerous places including [22].

6 E(t) for Nonlinear Systems

This section gives two approaches to the analysis of E(t) of nonlinear systems. Theorem 6.2 proves that the E(t) of the solution of a nonlinear system and the E(t) of the zero solution of its first approximation are the same. So one approach is to look at the linearized problem. The other approach is to directly make stability assumptions about the solution of the nonlinear system (Theorems 6.5 and 6.6). Both these approaches are illustrated in Section 7.

Lemma 6.1. Assume that F (t, 0) ≡ 0, that F (t, y) has continuous first order partial derivatives with respect to y, and that F (t, y) is continuous with respect to t, t ≥ 0. Then,

F (t, y) = A(t)y + g(t, y),

where A(t) = \frac{∂F (t,y)}{∂y} |_{y=0} and g(t, y) = o(∥y∥) as y → 0 uniformly over compact intervals of t.

Proof. See [5].

As noted in the introduction, the following theorem is implicit in the work of Henrici [12] [13]. It might be possible to generalize it to partial differential equations.

Theorem 6.2. As in Lemma 6.1, let F (t, 0) ≡ 0, let F (t, y) have continuous first order partial derivatives in y and be continuous in t. Then the zero solution of

ẏ(t) = F (t, y), y(0) = 0,

and the zero solution of

ẏ(t) = A(t)y, y(0) = 0,

where A(t) = \frac{∂F (t,y)}{∂y} |_{y=0} have the same E(t).

Proof. Let F (t, y) = A(t)y + g(t, y) as in Lemma 6.1. If needed, take h < h_{0}(t, r) for some r > 0. To reduce clutter in the proof, denote ỹ(kh; h) by ỹ_{k} and Y (kh) by Y_{k}; Y (t) is the principal fundamental matrix of ẏ(t) = A(t)y. Let n = ⌊t/h⌋ and h_{r} = t − nh.

Routine application of the variation of constants formula [9] [22] gives

ỹ_{k} = Y_{k}Y_{k−1}^{−1}ỹ_{k−1} + ∫_{(k−1)h}^{kh} Y_{k}Y^{−1}(s)g(s, ỹ(s; h)) ds + hα(h)v_{k},

for k = 1, 2, …, n, and

ỹ(t; h) = Y (t)Y_{n}^{−1}ỹ_{n} + ∫_{nh}^{t} Y (t)Y^{−1}(s)g(s, ỹ(s; h)) ds.

Here, ỹ_{0} = 0 and Y_{k}Y^{−1}(s) = I + O(h) for (k − 1)h ≤ s ≤ kh. Further, ∥ỹ(s; h)∥ ≤ te^{Lt}α(h) for 0 ≤ s ≤ t by Lemma 3.1, and ∥g(s, y)∥ = o(y) as y → 0 uniformly over s ∈ [0, t]. Therefore, we can write Y_{k}Y^{−1}(s)g(s, ỹ(s; h)) = Cα(h)η(h)u(s) for 0 ≤ s ≤ t, a constant C which can depend on t but not on h, and ∥u(s)∥ ≤ 1, where u(s) ∈ R^{d}; here η(h) → 0 as h → 0. The expressions for ỹ_{k} and ỹ(t; h) become,

ỹ_{k} = Y_{k}Y_{k−1}^{−1}ỹ_{k−1} + Cα(h)η(h) ∫_{(k−1)h}^{kh} u(s)ds + hα(h)v_{k},
ỹ(t; h) = Y (t)Y_{n}^{−1}ỹ_{n} + Cα(h)η(h) ∫_{nh}^{t} u(s)ds.

Combining these equalities,

\begin{aligned}\frac{ỹ(t; h)}{α(h)} = h ∑_{k=1}^{n} Y (t)Y_{k}^{−1}v_{k} & + Chη(h) ∑_{k=1}^{n} Y (t)Y_{k}^{−1} ∫_{(k−1)h}^{kh} u(s)ds \\ & + Cη(h) ∫_{nh}^{t} u(s)ds. \end{aligned}

Both the 2nd and 3rd terms vanish as h → 0 because of η(h). Thus the nonlinear term g(t, y) has no effect on E(t). The proof can be formally completed using (5.1).

Corollary 6.3. Let f(t, x) be continuous in t and have continuous first order partial derivatives with respect to x. The E(t) of the solution x(t; x_{0}) of ẋ(t) = f(t, x), x(0) = x_{0} , and the E(t) of the zero solution of ẏ(t) = A(t)y, y(0) = 0, where A(t) = \frac{∂f(t,x)}{∂x} |_{x=x(t;x_{0})}, are the same. In fact,

E(t) = sup_{v(s)}∥∫_{0}^{t} \frac{∂x(t)}{∂x(s)}v(s)ds∥,

where the supremum is taken over continuous functions v(s) satisfying ∥v(s)∥ ≤ 1.

Proof. For a proof of the above formula for E(t), note that if Y (t) is the principal fundamental matrix of ẏ(t) = A(t)y, then Y (t)Y^{−1}(s) = \frac{∂x(t)}{∂x(s)}; for a proof see any basic book on ODEs. Plug this into Theorem 5.2 to get the formula.

For comparison with Theorem 6.2, we state Theorem 33 of Chapter IX of Sansone and Conti [22]. This result is typical of Lyapunov’s theory of first approximation.

Theorem 6.4. Assume the zero solution of ẏ(t) = A(t)y, t ≥ 0, is uniformly asymptotically stable. If g(t, y) is continuous in t and y, and g(t, y) = o(∥y∥) uniformly over t ≥ 0 as y → 0, then the zero solution of ẏ(t) = A(t)y + g(t, y) is also uniformly asymptotically stable.

The assumption about g(t, y) being o(∥y∥) uniformly over the semi-infinite interval t ≥ 0 in Theorem 6.4 is quite stringent. A nonlinearity of the form g(t, y) = t[y_{1}^{2}, …, y_{d}^{2}]^{T} does not satisfy that assumption. But the theorem does not hold if the assumption is weakened to what is known about g(t, y) from Lemma 6.1; for a counterexample, see Bellman [1, p. 87]. A stringent assumption about g(t, y) is not required in Theorem 6.2 because how small an h we take to get convergence of e(t; h)/α(h) to E(t) can depend upon t. In definitions of stability, in contrast, we need one small perturbation which stays “close” till t = ∞.

The rest of this section is about bounding E(t) by making stability assumptions about the solution of a nonlinear system. For linear systems, this program was carried out in Section 4. We introduce a technique for error analysis of approximate solutions which uses Lyapunov functions. But first an example to point out some difficulties.

Example 6.1. Let us first consider the zero solution of ẏ(t) = y − e^{t}y^{3}, y(0) = 0, t ≥ 0. Its E(t), by Theorems 6.2 and 5.1, is e^{t} − 1. But we show that the zero solution is actually uniformly asymptotically stable. It is even exponentially asymptotically stable.

Figure 3 is the portrait of trajectories of ẏ(τ) = y − e^{t}y^{3}. The portrait for y ≤ 0 is a reflection about y = 0. Thus we can restrict ourselves to trajectories which are always in the upper half plane. Every point (t, y) with y ≥ e^{−t/2} is on a trajectory that is pointed downwards.

We now verify uniform stability of the zero solution using the last observation. Given ϵ > 0, choose t_{ϵ} so that e^{−t_{ϵ}/2} < ϵ. By continuity properties, we can choose a δ > 0 so that if τ ≤ t_{ϵ} and Figure 3: The portrait of trajectories of ẏ(t) = y − e^{t}y^{3}. All solutions tend towards the curve ẏ(t) = e^{−t/2}.

y_{τ} ≤ δ, the trajectory through (τ, y_{τ}) stays below ϵ till t_{ϵ}. The maximum possible height (along y) of a trajectory beginning at (τ, y_{τ}), τ ≥ t_{ϵ} and y_{τ} < δ, is bounded by the larger of e^{−τ/2} and |y_{τ} |. Since e^{−τ/2} < ϵ and δ < ϵ, uniform stability is verified.

The verification of uniform asymptotic stability will be sketchy. We base it on the following facts:

1. The solution of ẏ(t) = y − e^{t}y^{3}, y(0) = y_{0} tends to zero as t → ∞ for 0 ≤ y_{0} ≤ 1,

2. Further, y(t; y_{0}) ≤ y(t; 1) for t ≥ 0, 0 ≤ y_{0} ≤ 1,

3. And, 0 ≤ y(t + τ; τ; y_{0}) ≤ y(t; 0; y_{0}) for any y_{0} ≥ 0, τ ≥ 0, t ≥ 0.

The proof of item 1 involves a bit of elementary work which we do at the end. Item 2 is trivial. For item 3, think of y(t; 0, y_{0}) as the solution of ẏ(t) = y − e^{t}y^{3}, y(0) = y_{0}, and of y(t + τ; τ, y_{0}) = z(t) as the solution of ż(t) = z − e^{t+τ} z^{3}, z(0) = y_{0}, and use a differential inequality [11, p. 27].

Now let T_{ϵ} be such that ∥y(t; y_{0})∥ < ϵ for t ≥ T_{ϵ} and y_{0} = 1. Then ∥y(t + τ; τ, y_{0})∥ < ϵ for any τ ≥ 0, t ≥ T_{ϵ} and |y_{0}| ≤ 1. Thus, the zero solution is uniformly asymptotically stable.

To prove item 1, it is enough to verify that the solution with the initial condition y(0) = 1 satisfies y(t) < 2e^{−t/2} for t ≥ 0. This is obvious from the portrait of trajectories: every trajectory that cuts the curve y = 2e^{−t/2} is in the downward direction; therefore any trajectory that starts below that curve has to stay below it for t ≥ 0. In fact, y(t) < 2e^{−t/2} for the trajectory with y(0) = 1 implies that the zero solution is exponentially asymptotically stable, if we adopt Yoshizawa’s definition of exponential asymptotic stability [27].

Here we have an example which is uniformly asymptotically stable, yet has an E(t) which increases exponentially with t. This is possible because the approximations ỹ(t; h) can introduce errors at all points kh while the definitions of stability allow just one perturbation.

The next two theorems are nonlinear analogues of Theorem 5.5 and 5.6. The proofs this time rely on the theory of Lyapunov functions. Let us note that Lyapunov functions V (t, y) are always assumed to be continuous in t and locally Lipshitz in x. V′_{F} (t, y) is the rate of increase of V (t, y) along a solution of ẏ(t) = F (t, y) which goes through (t, y). More precisely, if y(t + τ; t, y) is such a solution, then

V′_{F} (t, y) = lim sup_{δ→0^{+}} \frac{V (t + δ, y(t + δ; t, y)) − V (t, y)}{δ}.

If V′_{F} (t, y) ≤ aV (t, y) then V (t + δ, y(t + δ; y, t)) ≤ e^{aδ}V (t, y), and if V′_{F} (t, y) ≤ 0 then V (t + δ, y(t + δ; y, t)) ≤ V (t, y); these two facts can be inferred from differential inequalities [9] [11].

We say that F (t, y) is uniformly Lipshitz in a neighbourhood of zero if ∥F (t, y_{1}) − F (t, y_{2})∥ < L∥y_{1} − y_{2}∥, for a constant L > 0 and any y_{i} with ∥y_{1}∥ < r, ∥y_{2}∥ < r where r > 0; L is the same constant for any t ≥ 0.

Theorem 6.5. Let F (t, y) be uniformly Lipshitz in y in a neighbourhood of 0. Assume that the zero solution of ẏ(t) = F (t, y), y(0) = 0, is exponentially stable in the sense that

∥y(t; t_{0}, y_{0})∥ ≤ Ke^{−c(t−t_{0})}∥y_{0}∥

for ∥y_{0}∥ < r, t_{0} ≥ 0, t ≥ t_{0}, c > 0, and K > 0. Then E(t) of the zero solution is bounded above by a constant.

Proof. Stability assumptions in the theorem imply existence of Lyapunov function with following properties (Yoshizawa [27, p. 97], Hale [10]):

1. ∥y∥ ≤ V (t, y) ≤ C∥y∥, where C > 0 is a constant,

2. |V (t, y_{1}) − V (t, y_{2})| ≤ L∥y_{1} − y_{2}∥,

3. V′_{F} (t, y) ≤ −qcV (t, x) for some 0 < q < 1.

The domain of V (t, y) is t ≥ 0 and ∥y∥ ≤ r for r > 0.

Take h < h_{0}(t, r). Because of item 3, V (t, ỹ(t; h)) decreases at least by a factor e^{−qch} along the approximate solution when t increases from kh to (k + 1)h; on that interval of t the approximate solution follows the exact solution till the discontinuity at t = (k +1)h. The discontinuity can cause an increase in V (t, y) of at most L times its magnitude by item 2. Therefore,

V (kh; ỹ(kh; h)) ≤ e^{−qch}V ((k − 1)h, ỹ((k − 1)h; h)) + Lhα(h),

for k = 1, 2, …, n and

V (t, ỹ(t; h)) ≤ e^{−qch_{r}}V (nh, ỹ(nh; h)).

Combining these inequalities, we get

V (t, ỹ(t; h)) ≤ e^{−qch_{r}}Lα(h)(\frac{1 − e^{−nqch}}{1 − e^{−qch}} )
≤ Kα(h).

That K above can be a constant independent of h and n can be deduced from basic calculus using h < h_{0}. Now, by item 1, ∥ỹ(t; h)∥ ≤ Cα(h) implying a constant upper bound for E(t).

The difficulty in proving Theorem 6.5 directly using the norm ∥·∥ is that when K > 1 the discretization error might actually be amplified by a factor greater than 1 over any given time step. Since we have to make the worst possible assumption at every time step, the final bound on E(t) obtained this way will actually be exponential in t when K > 1. The proof of Theorem 6.5 uses a carefully constructed Lyapunov function to get around this difficulty.

Theorem 6.6. Assume as in the previous theorem that F (t, y) is uniformly Lipshitz with respect to t in a neighbourhood of y = 0. If the zero solution of ẏ(t) = F (t, y), y(0) = 0, t ≥ 0, is uniformly asymptotically stable, then E(t) ≤ Kt for some constant K.

Proof. Stability assumptions in this theorem imply the existence of a Lyapunov function with the following properties: (Hale [10, Theorem 4.2, Chapter X], Yoshizawa [27])

1. ∥y∥ ≤ V (t, y),

2. V (t, 0) ≡ 0,

3. V′_{F} (t, y) ≤ 0, where V′_{F} (t, y) is defined as in the previous proof,

4. |V (t, y_{1}) − V (t, y_{2})| ≤ K∥y_{1} − y_{2}∥ for some constant K > 0.

The domain of definition of V (t, y) is the same as in the previous proof.

The proof is similar to that of Theorem 6.5, but this time

V (kh, ỹ(kh; h)) ≤ V ((k − 1)h, ỹ((k − 1)h; h)) + Khα(h)

for k = 0, 1, …, n − 1 and

V (t, ỹ(t; h)) ≤ V (nh, ỹ(nh; h)).

Combining these inequalities, we have V (t, ỹ(t; h)) ≤ Ktα(h). As before, ∥ỹ(t; h)∥ ≤ Ktα(h), which this time implies that E(t) ≤ Kt.

Corollary 6.7. Let x(t; x_{0}) be the nonzero solution of the initial value problem ẋ(t) = f(t, x), x(0) = x_{0}. Assume f(t, x) is uniformly Lipshitz in x in a neighbourhood of x(t; x_{0}). Then uniform asymptotic stability of x(t; x_{0}) implies that its E(t) is linearly bounded.

Let us call attention to the necessity of making a uniform Lipshitz assumption about F (t, x) or f(t, x) in Theorem 6.6 or Corollary 6.7. Example 6.1, which is uniformly asymptotically stable, does not satisfy the uniform Lipshitz assumption and has an E(t) which increases exponentially. We do not know if Theorem 6.6 is still true if the assumption of uniform asymptotic stability is weakened to just uniform stability. If such a theorem were true, its wider applicability might be of use. Table 1 summarizes all of Sections 5 and 6 except Theorem 6.2.

7 Three Applications to Dynamical Systems

We give three examples to illustrate the applicability of our methods for bounding the accumulation of global error.

7.1 Hyperbolic Sinks of C^{1} Dynamical Systems

Let p be a fixed point of a C^{1} dynamical system ẋ(t) = f(x); i.e. let f(p) = 0. Then p is a hyperbolic sink, if all the eigenvalues of \frac{∂f}{∂x} |_{x=p} have strictly negative real parts. The following theorem can be derived using Chapter 6 of [25]. We give the theorem here because our method of proof is different.

Theorem 7.1. Let x(t; x_{0}) be a trajectory of the dynamical system ẋ(t) = f(x), f ∈ C^{1}(R^{d}), which falls into a hyperbolic sink p as t → ∞. Then its E(t) is bounded above by a constant.

Stability E(t) can increase exponentially even if ∥A(t)∥ is bounded

Asymptotic stability E(t) can increase exponentially even if ∥A(t)∥ is bounded; Example 5.3

Linear Problems Uniform stability E(t) must be linearly bounded;

Theorem 5.5

Uniform asymptotic stability E(t) must be bounded by a con-

stant; Theorem 5.6

Uniform stability E(t) can increase exponentially Uniform stability with uniform Lipshitz assumption Not known if E(t) must be linearly bounded

Uniform asymptotic stability E(t) can increase exponentially;

Example 6.1

Nonlinear problems Uniform asymptotic stability with uniform Lipshitz assumption E(t) must be linearly bounded; Theorem 6.6

Exponential stability as in Theorem 6.5 Not known if E(t) must be linearly bounded

Exponential stability as in Theorem 6.5 with uniform Lipshitz assumption E(t) must be bounded by a constant; Theorem 6.5

Table 1: Summary of part of Sections 5 and 6. The second column is the stability assumption about the solution; the last column says what is known about the conditioning function E(t) corresponding to such a solution.

Proof. Without loss of generality, take p = 0. By [20, p. 150], there is a neighbourhood U_{0} of 0 such that x_{0} ∈ U_{0} implies

∥x(t; x_{0})∥ < ce^{−at}

for constants a > 0, c > 0, and for t ≥ 0.

Using continuity properties of solutions of differential equations (or openness of the basin of attraction), we can assume an open neighbourhood U of {x(t; x_{0})|t ≥ 0} in R^{d} and a compact set K containing U such that K and consequently U are both contained in the basin of attraction of p = 0. Since all trajectories beginning in K enter U_{0} in a finite amount of time, which depends only on K, we can assume

∥x(t; x_{0})∥ < Ce^{−at}

for constants a > 0, C > 0, for t ≥ 0, and for any x_{0} ∈ K. We can also take K and U to be invariant under the flow.

For any trajectory x(t; x_{0}), t ≥ 0, with x_{0} ∈ U, there exists r > 0 such that if ∥x(t; x_{0})−x_{1}∥ < r, t ≥ 0 then x_{1} ∈ K. For such an x_{0}, if ∥δ∥ < r,

∥x(t; x_{0}) − x(t; τ, x(τ; x_{0}) + δ)∥ < 2Ce^{−at}

for t ≥ τ ≥ 0. If x_{0} ∈ /U but is in the basin of attraction of p, the trajectory x(t; x_{0}) enters U in a finite amount of time. So we can get the same kind of bound as above by adjusting C and a if necessary.

To apply Theorem 5.5 and deduce boundednes of E(t), we need only verify uniformity of the Lipshitz condition on f(x) for x ∈ K. This is trivial since K is compact and f(x) is C^{1}.

7.2 Hyperbolic, Attracting Cycles of C^{1+ϵ} Dynamical Systems

Let x(t), t ≥ 0, be a periodic orbit of the C^{1} dynamical system ẋ(t) = f(x) in R^{d}. Let T > 0 be its period so that x(t + T) = x(t). Denote the set of points on this orbit by γ.

The characteristic multipliers of the cycle γ can be defined in two ways. One is to pick a point p ∈ γ, take a cross-section Σ at p, define a Poincaré map for Σ, and then define the characteristic multipliers as the (d−1) eigenvalues of the linearization of the Poincaré map at p. The other way is to consider the linear first approximation ẏ(t) = A(t)y on the cycle γ. Obviously, A(t + T) = A(t) for t ≥ 0. The Floquet numbers of this linear system can also be used to define characteristic multipliers. For a lucid account of these matters, see Robinson [20].

The cycle γ is hyperbolic and attracting if all its characteristic multipliers are strictly less than 1 in magnitude.

Theorem 7.2. Let x(t; x_{0}), t ≥ 0, be an orbit of a C^{1+ϵ} dynamical system ẋ(t) = f(x) in R^{d} which falls into a hyperbolic, attracting cycle γ as t → ∞. Then its E(t) is linearly bounded from above.

Let us first prove the following lemma. If any one solution of a linear system is uniformly stable, so is every other solution. So we might speak of the linear system itself as being uniformly stable.

Lemma 7.3. Assume x_{0} ∈ γ so that x(t; x_{0}) is a periodic orbit. Let its linear first approximation be ẏ(t) = A(t)y, t ≥ 0. If γ is hyperbolic and attracting, ẏ(t) = A(t)y is uniformly stable, and the E(t) associated with x(t; x_{0}) is linearly bounded.

Proof. Uniform stability of ẏ(t) = A(t)y is an easy consequence of the characteristic multipliers of γ being strictly less than 1. See Chapter IX of [22]. The linear bound on E(t) is implied by Theorem 5.5 and Corollary 6.7.

Lemma 7.3 is contained in a different form in the work of Cano and Sanz-Serna [4]. But let us note that Theorem 7.2 goes beyond Lemma 7.3 in a significant way. In practice, it is highly unlikely that x_{0} itself is on the cycle γ. But it is often easy to find x_{0} so that x(t; x_{0}) falls into a cycle γ.

The following lemma is known as the Dini-Hukuhara-Caligo theorem. It is Corollary 1 of Chapter IX of [22]. Its proof, which we omit, is short and simple, and illustrative of an important technique in stability theory.

Lemma 7.4. Assume the linear system ẏ(t) = A(t)y, t ≥ 0, is uniformly stable. Assume also that B(t), t ≥ 0, is continuous with ∫_{0}^{∞}∥B(t)∥dt < ∞. Then the linear system ẏ(t) = (A(t) + B(t))y is also uniformly stable.

Proof of Theorem 7.2. By Hartman [11, p. 254], there exists a point x′_{0} ∈ γ such that

∥x(t; x_{0}) − x(t; x′_{0})∥ < ce^{−at},  (7.1)

for constants a > 0, c > 0, and for t ≥ 0. This is called convergence in phase [20].

Let ẏ(t) = A(t)y, where A(t) = \frac{∂f}{∂x} |_{x=x(t;x′_{0})}, be the first approximation along x(t; x_{0}). By Lemma 7.3, this linear system is uniformly stable.

Let ẏ(t) = (A(t) + B(t))y, where A(t) + B(t) = \frac{∂f}{∂x} |_{x=x(t;x_{0})} be the first approximation along x(t; x_{0}). The estimate (7.1) for convergence in phase implies

∥B(t)∥ < c_{1}e^{−a_{1}t},

for constant a_{1} > 0 and c_{1} > 0. This is because both x(t; x_{0}) and x(t; x′_{0}) stay within a compact region of R^{d}, and f(x) is C^{1+ϵ} (this where we need the ϵ in C^{1+ϵ}). By Lemma 7.4, ẏ(t) = (A(t) + B(t))y is also uniformly stable.

Since A(t)+B(t) gives the linearization of x(t; x_{0}), the proof is easy to complete using Theorem 5.5 and Corollary 6.7.

7.3 Normally Contracting Manifolds with Quasiperiodic Flows

Let us introduce the notation ϕ_{t} for the flow induced on R^{d} by ẋ(t) = f(x). With this notation x(t; x_{0}) = ϕ_{t}x_{0}. Let V be a compact C^{1} manifold which is invariant under this flow. We consider the situation when the flow on V is differentiably conjugate to quasiperiodic flow on a torus, and V is normally hyperbolic and contracting, or briefly, normally contracting. We now explain the two italicized concepts in this paragraph.

A torus T^{n} is n copies of the circle S^{1}. If the angle on the ith circle is parameterized by θ_{i}, a quasiperiodic flow on T^{n} is of the form θ_{i}(t) = (θ_{i}(0) + α_{i}t) mod 2π. In fact, the flow is periodic if the α_{i} are all mutually commensurable. Denote this flow by ψ_{t}.

When we say that the flow ϕ_{t} is differentiably conjugate to quasiperiodic flow on a torus, we mean that there exists a C^{1} homeomorphism h : V → T^{n} such that h(ϕ_{t}x) = ψ_{t}(hx) for x ∈ V and t ≥ 0.

To define normal contractivity [14] [20], associate a direct sum decomposition T_{x} ⊕ N_{x} of R^{d} with every x in V. In this splitting T_{x} is the tangent space of V at x, and N_{x}, the normal space, varies continuously with x. If the N_{x} can be chosen so that

∥Π_{N_{y}} \frac{∂y}{∂x}|N_{x}∥ < ce^{−μt},

where y = ϕ_{t}x, the matrix inside the norm is the restriction of the derivative \frac{∂y}{∂x} to act from N_{x} to N_{y}, c > 0, and μ > 0, then V is normally contracting. Usually, the definition of normal contractivity comes with an other assumption which says contraction in the normal direction dominates any contraction on the manifold V. But since we have assumed that the flow on V is differentiably conjugate to quasiperiodic flow on a torus, this other assumption can be dropped.

Obviously, the tangent spaces T_{x} are invariant under the derivative map \frac{∂ϕ_{t}x}{∂x}. It is actually possible to choose N_{x} so that they too are invariant under the derivative map [14] [20]. We take this to be the case. So the derivative map \frac{∂ϕ_{t}x}{∂x} maps T_{x} to T_{ϕ_{t}x} and N_{x} to N_{ϕ_{t}x}.

Theorem 7.5. Let x(t; x_{0}) be a trajectory of the C^{1+ϵ} flow ẋ(t) = f(x) which falls into a normally contracting and invariant manifold V. Assume that the flow on V is differentiably conjugate to quasiperiodic flow on a torus. Then the E(t) of x(t; x_{0}) is linearly bounded.

The above theorem generalizes Theorem 7.2. Its proof is exactly analogous. We begin with a lemma about trajectories that begin on V.

Lemma 7.6. Let x_{0} ∈ V so that x(t; x_{0}) stays on V for t ≥ 0. Make the same assumptions about V as in Theorem 7.5. Then the linearization ẏ(t) = A(t)y along x(t; x_{0}) is uniformly stable, and the E(t) of x(t; x_{0}) is linearly bounded.

Proof. The principal fundamental matrix of the linearization in the lemma is given by Y (t) = \frac{∂ϕ_{t}x_{0}}{∂x_{0}} , which is the derivative map. Therefore, it is enough if we show that ∥\frac{∂ϕ_{t}x}{∂x}∥ is bounded by a constant for any x ∈ V and t ≥ 0.

We already know that the maps induced by the derivative map between tangent spaces and between normal spaces are bounded in norm because of differentiable conjugacy to flow on a torus and normal contractivity, respectively. Since the tangent spaces and normal spaces are both invariant under the derivative map, it is enough if we show that the angle between T_{x} and N_{x} (in the sense of the CS decomposition) is bounded away from 0. That this angle is bounded away from 0 is implied by the compactness of V.

The proof of Theorem 7.5 can be completed exactly like the proof of Theorem 7.2 using the following result about convergence in phase.

Theorem 7.7. As in Theorem 7.2, let x(t; x_{0}) be a trajectory of the C^{1} flow ẋ(t) = f(x) which falls into a normally contracting and invariant manifold V, and let the flow on V be differentiably conjugate to quasiperiodic flow on a torus. Then there exists x′_{0} ∈ V such that

∥x(t; x_{0}) − x(t; x′_{0})∥ < ce^{−at},

for positive constants c and a, and t ≥ 0.

Proof. This theorem can be deduced from Theorem 4.1 and the remark following its proof in Hirsch, Pugh and Shub [14]. See in particular part (a) of that theorem about stable manifolds and part (g) about conjugacy to linearized flows.

8 Three Remarks

This last section is a collection of parenthetical remarks about matters which are related to E(t) and which we have not investigated in detail.

(i) Multistep methods and variable time stepping. The model for discretization errors in Section 2 is adapted to single step methods with constant step sizes. For linear multistep methods with constant step sizes, we believe the accumulation of global error can be worse but not better (after excluding some trivial cases) than indicated by E(t).

For variable time stepping, it is sometimes but not always true that the ratio of the largest to the smallest time step is bounded by a constant [23] [24]. In this case at least, the effect of variable time stepping is to improve global errors by a constant factor over the indication given by E(t).

(ii) One sided Lipshitz conditions. Let us briefly summarize the approach to global error analysis using one sided Lipshitz conditions for the linear system ẏ(t) = A(t)y. Since ∥y(t)∥^{2} = y^{T} (t)y(t), we have

\frac{d∥y(t)∥^{2}}{dt} = y^{T} (t)(A^{T} (t) + A(t))y(t).

If λ(t) is the maximum eigenvalue of A^{T} (t) + A(t), then

\frac{d∥y(t)∥^{2}}{dt} ≤ λ(t)∥y(t)∥^{2}.

Thus upper bounds for ∥y(t)∥, and hence for ∥Y (t)∥ where Y (t) is the principal fundamental matrix of the linear system, can be written down in terms of λ(t). These can be plugged into Theorem 5.2 to get bounds on E(t) and hence the accumulation of global error. For a detailed account, see [9].

In our view, one sided Lipshitz conditions are basically a way to get a handle on stability by looking at the evolution of the norm of y(t). This is a far less general approach to stability than the two methods of Lyapunov we have used (this deficiency of norms is something Lyapunov must have realized). It is also of far lesser applicability; we do not see a way to derive any of the results in Section 7 by using one sided Lipshitz conditions.

(iii) Numerical symplectic integrators and Hamiltonian problems.

Special methods which preserve properties like symplecticity or energy conservation or volume conservation in phase space show milder increase of global errors with t than general methods in numerically integrating Hamiltonian problems. Calvo and Sanz-Serna [3] showed using careful numerical experiments that the accumulation of global error with t was linear for a symplectic method but quadratic for a Runge-Kutta method when the methods were applied to Kepler’s problem. A full explanation came after [2] and [7] in the paper [4] by Cano and Sanz-Serna. There are numerical experiments in Quispell and Dyt [19] which are unexplained.

The conditioning function E(t) is not suitable for studying the accumulation of global error in these special methods. The class of local discretization errors allowed by the model in Section 2 is too general. The discretization errors made by these special methods are restricted in some way; for example, the discretization errors of an energy preserving method cannot take the approximation out of the constant energy manifold. But is it possible to define conditioning functions which are adapted to these special methods? We think it might be possible. The concepts of stability needed to bound such conditioning fuctions will also, no doubt, be different.

9 Acknowledgements

It is a pleasure to thank Arieh Iserles and Mike Shub for encouragement and enlightening discussions. I have benefited from talking to several visitors to M.S.R.I., Berkeley, in Fall 98 about this work. I thank Timo Eirola and Andrew Stuart for suggestions which directly influenced the presentation of this paper.

References

[1] R. Bellman, Stability Theory of Differential Equations, McGraw-Hill, New York, 1953.

[2] M.P. Calvo and E. Hairer, Accurate long-term integration of dynamical systems, Appl. Numer. Math. 18 (1995), 95-105.

[3] M.P. Calvo and J.M. Sanz-Serna, The development of variable-step symplectic integrators, with application to the two-body problem, SIAM J. Sci. Comput. 14 (1993), 936-952.

[4] B. Cano and J.M. Sanz-Serna, Error growth in the numerical integration of periodic orbits, with application to Hamiltonian and reversible systems, SIAM J. Numer. Anal. 34 (1997), 1391-1417.

[5] E.A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, McGraw-Hill, New York, 1955.

[6] G. Dahlquist, Stability and error bounds in the numerical integration of ordinary differential equations, Trans. of the Royal Inst. of Techn., Stockholm, Sweden, Number 130.

[7] D.J. Estep and A.M. Stuart, The rate of error growth in Hamiltonian-conserving integrators, Z. Agnew. Math. Phys. 46 (1995), 407-418.

[8] W.B. Gragg, Repeated extrapolation to the limit in the numerical solution of ordinary differential equations, Thesis, Univ. of California, 1964; see also SIAM J. Numer. Anal. 2 (1965), 384-403.

[9] E. Hairer, S.P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I, Springer-Verlag, New York, 1980.

[10] J. Hale, Ordinary Differential Equations, John Wiley, New York, 1969.

[11] P. Hartman, Ordinary Differential Equations, John Wiley, New York, 1973.

[12] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley, New York, 1962.

[13] P. Henrici, Error Propogation for Difference Methods, John Wiley, New York, 1963.

[14] M.W. Hirsh, C.C. Pugh and M. Shub, Invariant Manifolds, Lecture notes in Mathematics 583, Springer-Verlag, New York, 1977.

[15] A. Iserles and S.P. Nørsett, On the solution of linear differential equations in Lie groups, Phil. Trans. of Royal Soc. A, to appear.

[16] A. Iserles and G. Söderlind, Global bounds on numerical error for ordinary differential equations, J. Complexity 9 (1993), 97-112.

[17] A. M. Lyapunov, Problème générale de la stabilité du mouvement, Comm. Soc. Math. Kharkov 2 (1892), 3 (1893); Ann. Fac. Sci. Toulouse 9 (1907), 204-474; Ann. of Math. Studies 17, Princeton University Press, Princeton, 1949.

[18] J. G. Malkin, Theory of Stability of Motion, Atomic Energy Commision Tech. Rep. 3352, 1956.

[19] G.R.W. Quispel and C.P. Dyt, Volume-preserving integrators have linear error growth, Physics Letters A 242 (1998), 25-30.

[20] C. Robinson, Dynamical Systems. Stability, Symbolic Dynamics, and Chaos, CRC press, Boca Roton, 1995.

[21] W. Rudin, Real and Complex Analysis, 3rd Edition, McGraw-Hill, New York, 1987.

[22] G. Sansone and R. Conti, Non-linear Differential Equations, Macmillan, New York, 1964.

[23] D. Stoffer and K. Nipp, Invariant curves for variable step size integrators, BIT 31 (1991), 169-180.

[24] A.M. Stuart, Probabilistic and deterministic convergence proofs for software for initial value problems, Numerical Algorithms 14 (1997), 227-260.

[25] A.M. Stuart and A.R. Humphries, Dynamical Systems and Numerical Analysis, Cambridge University Press, Cambridge, UK, 1996.

[26] L.N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, Philadelphia, 1997.

[27] T. Yoshizawa, Stability Theory by Liapunov’s Second Method, Mathematical Society of Japan, 1966.

Mathematical Sciences Research Institute (MSRI) 1000 Centennial Drive Berkeley, CA 94720. divakar@cs.cornell.edu