NUMERICAL INTEGRATORS THAT CONTRACT VOLUME

ROBERT I. MCLACHLAN AND G.R.W. QUISPEL

Abstract. We study numerical integrators that contract phase space volume even when the ODE does so at an arbitrarily small rate. This is done by a splitting into two-dimensional contractive systems. We prove a sufficient condition for Runge-Kutta methods to have the appropriate contraction property for these two-dimensional systems; the midpoint rule is an example.

1. Introduction. What is a dissipative system? In physics, the term usually refers to possession of a scalar function (such as energy) which decreases in time, and one speaks of, e.g., the dissipative pendulum, ẍ = − sin x−εẋ, for which \frac{d}{dt}(\frac{1}{2}ẋ^{2}−cos x) = −εẋ^{2} ≤ 0. (See [1, 6] for some general formulations of such systems.) In dynamical systems, it usually refers to a decrease of phase space volume in time, as in the “dissipative Hénon map” (x, y) ↦ (y, 1 + bx − ay^{2}), with Jacobian determinant −b— phase space area decreases if |b| < 1. Another example is the famous Lorenz system, which contracts volume at a constant rate. In the numerical analysis of ODEs, it has been used to describe systems that decrease some norm of the solution, either in the sense that \frac{d}{dt}∥x∥^{2} < a − b∥x∥^{2} for some a, b > 0, or \frac{d}{dt}∥x∥^{2} < 0 for all ∥x∥ > R > 0 [11].

In the field of geometric integration, much work has been done in maintaining the preservation of a conserved quantity (first integral) [3, 6, 7], the decrease of a dissipated quantity (a Lyapunov function) [6, 7], or the preservation of phase space volume [2, 5]. Here we look at the missing case, and study how to maintain the property of contracting phase space volume.

Consider the ODE

ẋ = f(x), x ∈ ℝ^{n} (1)

with solution x(t) and Jacobian (first variation) A(t) = ∂x(t)/∂x(0) which evolves according to

Ȧ = F A, A(0) = I,

where F (x) = df(x) is the derivative of the vector field f. We have

\frac{d}{dt} det A = det A tr (A^{−1}Ȧ)
= det A tr F

so that phase space volume contracts, is preserved, or expands when tr F < 0, tr F = 0, or tr F > 0 for all x, respectively. tr F is the divergence or trace of the vector field f. Strongly contractive systems are those for which there is a b such that tr F < b < 0. In this case any consistent numerical integrator will be contractive for small enough time step h. Therefore we concentrate on weak contraction (defined below), which is a closed property and more difficult to preserve. It turns out that requiring contractivity for all h > 0 and all contractive f is prohibitively difficult, which leads to the following definition. We consider one-step methods x_{n+1} = g(x_{n}) with Jacobian A = dg(x_{0}).

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

Definition 1. The ODE (1) is (weakly) contractive if tr F ≤ 0 for all x. An integrator is (weakly) contractive if for any matrix norm ∥ · ∥ and all L > 0 there is a time step h^{∗} > 0 such that | det A| ≤ 1 for all 0 < h < h^{∗}, for all x, and for all f such that ∥F ∥ < L and tr F ≤ 0.

That is, there might be stiffness problems (for large L, h^{∗} might be small), but the time step needed to preserve contractivity should not tend to zero as tr F → 0. Note that a contractive integrator as defined here is not necessarily volume-preserving when the ODE is, nor is the relative amount of contraction necessarily correct as tr F → 0. These would be true if we added the requirement ln(det A)/h tr F → 1 uniformly as tr F → 0 uniformly, for all fixed h < h^{∗}. The midpoint rule (see Proposition 3, below) satisfies this, for example.

Since there are no known linearly covariant volume-preserving schemes in more than two dimensions [2], we expect that the same is true here, and we immediately consider systems in two dimensions.

2. Dissipative schemes in two dimensions.

Example 2. Euler’s method is not contractive in two dimensions. We have x_{n+1} = x_{n} + hf(x_{n}) so A = I + hF. In two dimensions,

det A = det ( \begin{matrix}1 + hF_{11} & hF_{12} \\ hF_{21} & 1 + hF_{22} \end{matrix} ) = 1 + h tr F + h^{2} det F.

So det A ≤ 1 for all h if det F ≤ 0, and det A ≤ 1 for

h ≤ \frac{− tr F}{det F}

if det F > 0, so small contractivity can require a small time step to be captured.

Note that since det A = det(I + hF) = ∏(1 + hλ_{i}), where λ_{i} are the eigenvalues of F, Euler’s method is contractive in n dimensions on systems with bounded negative eigenvalues. We look at this further in Section 3.

Proposition 3. The midpoint rule, x_{n+1} = x_{n} + hf(x̄), x̄ = (x_{n} + x_{n+1})/2, is contractive in two dimensions.

Proof. We have

A = (I − \frac{1}{2}hF (x̄))^{−1}(I + \frac{1}{2}hF (x̄)),

so

det A = \frac{1 + he + h^{2}d}{1 − he + h^{2}d}

NUMERICAL INTEGRATORS THAT CONTRACT VOLUME 3

where e = \frac{1}{2} tr F (x̄), d = \frac{1}{4} det F (x̄). Thus (det A)^{2} ≤ 1 if

(1 + he + h^{2}d)^{2} ≤ (1 − he + h^{2}d)^{2}

or

e (1 + h^{2}d) ≤ 0

Since e ≤ 0, this is true for all h if e = 0 (the well-known result that the midpoint rule is area-preserving, or symplectic), for all h if d ≥ 0, or for h < 1/\sqrt{−d} if d < 0.

Proposition 3 can be generalized as follows. Proposition 4. The symplectic Runge-Kutta methods with b_{i} > 0 for all i are contractive in two dimensions. Proof. For the terminology, see [9]. Our proof closely follows their proof of symplecticity. An s-stage Runge-Kutta method is defined by

X_{i} = x_{n} + h ∑_{j=1}^{s} a_{ij}f(X_{j}), (2)

x_{n+1} = x_{n} + h ∑_{j=1}^{s} b_{j}f(X_{j}), (3)

and is symplectic if b_{i}b_{j} −b_{i}a_{ij} −b_{j}a_{ji} = 0 for all i and j. Note that in two dimensions, A^{T} JA = J det A, where J = ( \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ), so we evaluate the left hand side. Let D_{i} = df(X_{i}) = F (X_{i})dX_{i} =: F_{i}A_{i}. Differentiating (3) gives

\begin{aligned}A^{T} JA & = (I + h ∑_{i} b_{i}D_{i})^{T} J(I + h ∑_{j} b_{j}D_{j}) \\ & = J + h ∑_{i} b_{i} (JD_{i} + D_{i}^{T} J) + h^{2} ∑_{i,j} b_{i}b_{j}D_{i}^{T} JD_{j}. \end{aligned}

Differentiating (2) gives

A_{i} = I + h ∑_{j} a_{ij}D_{j} (4)

or

JD^{i} = AT i JD^{i} − h ∑_{j} a^{ij}DT j JD^{i} .

Inserting,

A^{T} JA = J + h ∑_{i} b_{i} (A_{i}^{T} JD_{i} + D_{i}^{T} JA_{i}) + h^{2} ∑_{i,j} (b_{i}b_{j} − b_{i}a_{ij} − b_{j}a_{ji}) D_{i}^{T} JD_{j}.

The last term is zero because of the assumption on the coefficients b_{i}, a_{ij}. Now D_{i} = F_{i}A_{i}, so

A^{T} JA = J + h ∑_{i} b_{i}A_{i}^{T} (JF_{i} + F_{i}^{T} J)A_{i}

= J + h ∑_{i} b_{i}A_{i}^{T} JA_{i} tr F_{i}

= J(1 + h ∑_{i} b_{i} det A_{i} tr F_{i})

so

det A = 1 + h ∑_{i} b_{i} det A_{i} tr F_{i}.

From (4), det A_{i} is bounded and equal to 1 + O(h). Using b_{i} > 0 and tr F_{i} ≤ 0 gives the result.

The assumption b_{i} > 0 is necessary. Suppose there are s = 2 stages with b_{1} > 0 and b_{2} < 0. Then the vector (b_{i}) lies in the fourth quadrant, and all we know of the vector (tr F_{i}) is that it lies in the third quadrant. In regions where the trace varies relatively quickly, the angle between these two vectors can be less than \frac{π}{2}, leading to (b_{i}) · (tr F_{i}) > 0 and det A > 1.

These methods actually preserve area when tr F = 0. In fact, this is not necessary for contractivity in two dimensions, because we can allow a small amount of “numerical contractivity” even as tr F → 0; away from tr F = 0 the inherent contractivity of the ODE contributes. It turns out that only methods of order 2, 3, 6, 7, …, can achieve this.

Lemma 5. Let R(z) be the linear stability polynomial of a consistent Runge-Kutta method. In two dimensions, the method is contractive on linear ODEs if there is a u^{∗} > 0 such that

R(u)R(−u) ≤ 1, R(iu)R(−iu) ≤ 1

for all 0 ≤ u < u^{∗}.

Proof. In n dimensions, a Runge-Kutta method on linear problems ẋ = F x has derivative A = R(hF). Therefore det A = ∏_{i} R(hλ_{i}) = 1 + h tr F + O(h^{2}), so the method is contractive if tr F < 0. If tr F = 0 we have to examine det A in more detail. In n = 2 dimensions, there are only two such cases: the eigenvalues can be (u, −u) or (iu, −iu). This gives the result.

We note that the result also applies to nonlinear problems with 1-stage methods, since then F (x) is evaluated at only a single point.

Proposition 6. Let the method have order p, so that R(z) = e^{z} + az^{p+1} + bz^{p+2} + O(z^{p+3}). If 4|(p + 1) and a < 0, or if 4|(p + 2) and b < a, then the method is contractive on linear problems in two dimensions.

Proof. We expand

\begin{aligned}R(z)R(−z) − 1 & = e^{−z}(az^{p+1} + bz^{p+2}) + e^{z}(a(−z)^{p+1} + b(−z)^{p+2}) + … \\ & = az^{p+1}(1 − (−1)^{p}) + (b − a)z^{p+2}(1 + (−1)^{p}) + … \end{aligned}

The leading term must be negative for z = u and for z = iu, so it must be a fourth power. If p is even, the leading term is z^{p+2} so 4|(p + 2) and we need b < a; if p is odd, the leading term is z^{p+1} so 4|(p + 1) and we need a < 0.

An example is any 3-stage, 3rd order Runge-Kutta, which has R(z) = 1+z + \frac{1}{2}z^{2} + \frac{1}{6}z^{3}.

This result can be extended to more dimensions. For example, a longer calculation shows that Proposition 6 holds with p = 3 in three dimensions. We are not sure how it extends to nonlinear systems. It seems that if the eigenvalues of F are varying rapidly, contractivity could be lost.

3. More than two dimensions. For systems in more than two dimensions, we generalize the volume-preserving method of Feng and Wang ([2]; see also [5]). We write the ODE as a sum of two-dimensional contractive systems (i.e., ones for which ẋ_{i} = 0 except for two indices i), apply a contractive method to each term, and compose the resulting maps with positive time steps. Since contractivity is a semigroup property, we can build a contractive integrator of order 1 or 2 in this way [5]. This relies on the following proposition.

Proposition 7. Any C^{r+1} contractive ODE is the sum of two-dimensional C^{r} contractive ODEs.

Proof. Consider ẋ = f(x), F = df. We shall write f in the form f_{i} = ∑_{j} ∂_{j}L_{ij} (where ∂_{j} = ∂/∂_{x_{j}}.)

Let s_{ij}(x) be n^{2} functions with

s_{ij}(x) + s_{ji}(x) ≥ 0 and ∑_{i,j=1}^{n} s_{ij}(x) = 1

for all x. Let

S_{ij} = ∫∫ s_{ij}(x) tr F (x) dx_{i} dx_{j}

where any values of the indefinite integrals can be taken. Let

f̃_{i} = f_{i} − ∑_{j} ∂_{j}S_{ij} = f_{i} − ∑_{j} ∫ s_{ij} tr F dx_{i},

so that

tr df̃ = (1 − ∑_{i,j} s_{ij}) tr F = 0.

Thus, f̃ is traceless and can be written as

f̃_{i} = ∑_{j} ∂_{j}A_{ij}

where the matrix A is antisymmetric and as smooth as f [2, 5]. Therefore

f_{i} = ∑_{j} ∂_{j}(A_{ij} + S_{ij}) (5)

or L = A + S.

For an explicit splitting, we take the diagonal elements s_{ii}(x) = 0. Then ẋ = f is the sum of the following n(n − 1)/2 two-dimensional ODEs:

\begin{aligned}ẋ_{i} & = ∂_{j}L_{ij} \\ ẋ_{j} & = ∂_{i}L_{ji} \\ ẋ_{k} & = 0 for k ≠ i, j \end{aligned}

for each pair (i, j) of indices from 1 to n. Each is contractive because each A piece is traceless and each S piece has trace (s_{ij} + s_{ji}) tr F ≤ 0.

One degree of smoothness is lost in this splitting, because each S piece depends on tr F.

An interesting solution is obtained by taking s_{ii} = 1/n, s_{ij} = 0 for i ≠ j, and

nL_{ij} = ∫ f_{i} dx_{j} − ∫ f_{j} dx_{i} + δ_{ij} ∫∫ tr F dx_{i} dx_{i}.

However, a more practical decomposition is to take the same S but A_{ij} = 0 for |i − j| > 1; this gives the minimum of n − 1 two-dimensional ODEs.

Although the above proof is constructive, it may be possible to find a more convenient splitting by ad hoc methods, in some cases leading to an explicit contractive integrator.

Firstly, if f is the sum of integrable contractive vector fields, then their flows can be composed to give a contractive integrator for f. For example, the Lorenz system,

ẋ = ( \begin{matrix}−σ & σ & 0 \\ ρ & −1 & 0 \\ 0 & 0 & −β \end{matrix} ) x + ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & −x_{1} \\ 0 & x_{1} & 0 \end{matrix} ) ∇ (\frac{x_{2}^{2} + x_{3}^{2}}{2} ) ,

is the sum of a linear, contractive part and a Poisson, non-contractive part, each of which may be solved exactly, giving an integrator with exactly correct contractivity.

Secondly, it may be possible to use a simpler method, such as Euler, on some of the pieces. Here are some criteria which allow this.

Proposition 8. Euler’s method is contractive in n dimensions if there is a b such that tr(F^{2}) > b > 0. This condition is equivalent to ∥S∥^{2} > ∥A∥^{2} + b, where F = A + S, A = −A^{T}, S = S^{T}, and ∥ · ∥ is the Frobenius (sum of squares) norm. This condition is satisfied if all the eigenvalues of F are bounded away from the sectors \frac{π}{4} < |θ| < \frac{3π}{4}; in particular, if they are all real and bounded away from zero.

Proof. Let λ_{i} be the eigenvalues of F. For Euler’s method we have

ln det A = ln det(I + hF)

= ln ∏_{i}(1 + hλ_{i})

= ∑ ln(1 + hλ_{i})

= h ∑_{i} λ_{i} − \frac{1}{2}h^{2} ∑_{i} λ_{i}^{2} + O(h^{3})

= h tr F − \frac{1}{2}h^{2} tr(F^{2}) + O(h^{3}).

If there is a b such that tr(F^{2}) > b > 0, this is less than 0 for all small enough h, i.e., the method is contractive. Splitting F into its symmetric and antisymmetric parts,

tr(F^{2}) = ∑_{i,j} F_{ij}F_{ji} = ∑_{i,j}(S_{ij} + A_{ij})(S_{ij} − A_{ij}) = ∥S∥^{2} − ∥A∥^{2},

giving the second part of the proposition. Now tr(F^{2}) = ∑_{i} λ_{i}^{2}, and if each λ_{i} is outside the specified sectors, then each real eigenvalue or complex conjugate pair of eigenvalues gives a positive contribution to this sum, giving the last part of the proposition.

Note that the eigenvalues of elliptic or nearly elliptic fixed points lie near the imaginary axis—right in the middle of the “bad” sector. Perhaps this was only to be expected.

Experts will recognize the last part of Proposition 8 as the appearance of an order star of a Runge-Kutta method [4] (the set {z : |R(z)| < |e^{z}| where R(z) is the method’s linear stability polynomial). For linear problems, or nonlinear problems with 1-stage methods, a method is more contractive than the flow of the ODE if hλ_{i} lies in the order star of the method for each eigenvalue λ_{i}. However, this seems rather restrictive so we do not explore further.

Proposition 9. There are explicit contractive integrators.

Proof. Let f be any contractive vector field with ∥F ∥ < L. Because eigenvalues vary continuously and can only become imaginary when two eigenvalues meet, and because symmetric matrices have real eigenvalues, there is a symmetric, traceless matrix M with distinct eigenvalues such that the derivative of f_{1} := f − Mx has real eigenvalues. Let f_{2} := Mx and split f = f_{1} + f_{2}. f_{1} is contractive and admits an explicit contractive integrator (e.g. Euler’s method, see Proposition 8); f_{2} is traceless and can be solved explicitly. Composing these maps gives the result.

We close with some open questions we hope to report on in the future.

1. Are there explicit contractive integrators of any order? (Proposition 9 constructs a first order method.) There are if one only demands linear contractivity. The order cannot be increased by composition, because the adjoint of Euler’s method—backward Euler—is not contractive for f_{1}.

2. The present method reduces to the volume-preserving method of Feng and Wang [2] when the vector field f is traceless. There is another approach to volumepreserving integration due to Quispel [8] and to Shang [10], which does not rely on a splitting at all; moreover, it has a generalization to systems preserving non-Euclidean measures, which we have not even considered here. Can this approach be carried over to the contractive case?

3. The splitting used in the proof of Proposition 5 writes f = a+b where tr(da) = 0 and b ≡ 0 when tr(df) ≡ 0. Are there splittings with the property that b(x) = 0 when tr(df(x)) = 0? If so, they could be used for systems in which tr(df(x)) changes sign on a compact hypersurface; the interior would then be invariant and one could construct an integrator which preserved it and was contractive there. This was done for the case of dissipation of scalar functions in [7].

Acknowledgements. The authors thank John Butcher and all the organizers of ANODE for financial support and for providing the atmosphere in which this paper was conceived, and the MSRI where it was concluded. We also thank the Marsden Fund of the Royal Society of New Zealand for their financial support.

References

[1] A. Bloch, P.S. Krishnaprasad, J.E. Marsden, and T.S. Ratiu, The Euler-Poincaré equations and double bracket dissipation, Comm. Math. Phys. 175 (1996), 1–42.

[2] K. Feng and D.-L. Wang, Dynamical systems and geometric construction of algorithms, in Contemporary Mathematics 163, eds. Z. Shi and C. Yang, 1994.

[3] O. Gonzalez, Time integration and discrete Hamiltonian systems, J. Nonlinear Sci. 6(5), 449– 467 (1996).

[4] A. Iserles and S.P. Nørsett, Order Stars, Chapman & Hall, London, 1991.

[5] R.I. McLachlan and G.R.W. Quispel, Generating functions for dynamical systems with symmetries, integrals, and differential invariants, Physica D 112 (1998) 298–309.

[6] R.I. McLachlan, G.R.W. Quispel, and N. Robidoux, Geometric integration using discrete gradients, Phil. Trans. Roy. Soc. A, to appear.

[7] R.I. McLachlan, G.R.W. Quispel, and N. Robidoux, A unified approach to Hamiltonian systems, Poisson systems, gradient systems, and systems with Lyapunov functions and/or first integrals, Phys. Rev. Lett., to appear.

[8] G.R.W. Quispel, Volume-preserving integrators, Phys. Lett. 206A (1995) 26-30.

[9] J.-M. Sanz-Serna and M.P. Calvo (1994), Numerical Hamiltonian Problems, Chapman & Hall.

[10] Shang Zai-jiu, J. Comp. Math. 12 (1994), 265.

[11] A.M. Stuart and A.R. Humphries, Dynamical systems and numerical analysis, CUP, 1996.

Robert I. McLachlan, Mathematics, Massey University, Private Bag 11–222, Palmerston North, New Zealand

 E-mail address: R.McLachlan@massey.ac.nz

G.R.W. Quispel, Mathematics Department, LaTrobe University, Bundoora, Melbourne 3083, Australia

 E-mail address: R.Quispel@latrobe.edu.au