Griffiths D.F., Higham D.J. Numerical Methods for Ordinary Differential Equations: Initial Value Problems

Подождите немного. Документ загружается.

196 14. Geometric Integration Part I—Invariants

14.2 Linear Invariants

A simple chemical reaction known as a reversible isometry may be written



. (14.1)

Here, a molecule of chemical species X

may spontaneously convert into a

molecule of chemical species X

, and vice versa. The constants k

and k

re-

ﬂect the rates at which these two reactions occur. If these two species are not

involved in any other reactions, then clearly, the sum of the number of X

and

molecules remains constant.

The mass action ODE for this chemical system has the form

(t) = −k

u(t) + k

v(t)

(t) = k

u(t) − k

v(t),

)

(14.2)

where u(t) and v(t) represent the concentrations of X

and X

respectively.

For initial conditions u(0) = A and v(0) = B, this ODE has solution

u(t) =

(A + B)

+ k



A −

(A + B)

+ k



−(k

v(t) =

(A + B)

+ k



B −

(A + B)

+ k



−(k

;

(14.3)

see Exercise 14.1. Adding gives u(t) + v(t) = A + B, showing that the ODE

system automatically preserves the sum of the concentrations of species X

and

Applying Euler’s method to (14.2) gives us

n+1

= u

− hk

+ hk

n+1

= v

+ hk

− hk

Summing these two equations, we ﬁnd that u

n+1

+ v

n+1

= u

+ v

. This shows

that Euler’s method preserves the total concentration, matching the property

of the ODE. This is not a coincidence—by developing a little theory, we will

show that a more general result holds.

Instead of solving the system (14.2) and combining the solutions, we could

simply add the two ODEs, noting that the right-hand sides cancel, to give

(t) + v

(t) = 0; that is,

(u(t) + v(t)) = 0.

We may deduce immediately that u(t) + v(t) remains constant along any solu-

tion of the ODE. Given a general ODE system x

(t) = f(x) with x(t) ∈ R

14.2 Linear Invariants 197

and f : R

→ R

, we say that there is a linear invariant if some linear

combination of the solution components is always preserved:

(t) + c

(t) + ··· + c

(t)) = 0,

where c

, c

, . . . , c

are constants. We may write this more compactly as

x(t) = 0,

where c ∈ R

. Because each x

is given by f

(x), this is equivalent to

(x) + c

(x) + ··· + c

(x) = 0,

which we may write as

f(x) = 0, for any x ∈ R

. (14.4)

The condition (14.4) asks for f to satisfy a speciﬁc geometric property: wher-

ever it is evaluated, the result must always be orthogonal to the ﬁxed vector c.

For Euler’s method, x

n+1

= x

+ hf (x

) and, under the condition (14.4),

we have

n+1

= c

+ hc

f (x

) = c

, (14.5)

so the corresponding linear combination of the components in the numerical

solution is also preserved. More generally, since they update from step to step

by adding multiples of f values, it is easy to see that any consistent linear

multistep method and any Runge–Kutta method will preserve linear invariants

of an ODE; see Exercise 14.2.

Linear invariants arise naturally in several application areas.

– In chemical kinetics systems like (14.2), it is common for some linear combi-

nation of the individual population counts to remain ﬁxed. (This might not

be as simple as maintaining the overall sum; for example, one molecule of

species X

may combine with one molecule of species X

to produce a single

molecule of species X

– In mechanical models, the overall mass of a system may be conserved.

– In stochastic models where an ODE takes the form of a “master equation”

describing the evolution of discrete probabilities, the total probability must

sum to one.

– In epidemic models where individuals move between susceptible, infected and

recovered states, the overall population size may be ﬁxed.

In all cases, we have shown that standard numerical methods automatically

inherit the same property. Now it is time to move on to nonlinear invariants.

198 14. Geometric Integration Part I—Invariants

14.3 Quadratic Invariants

The angular momentum of a free rigid body with centre of mass at the origin

can be described by the ODE system [26, Example 1.7, Chapter IV]

(t) = a

(t) x

(t)

(t) = a

(t) x

(t)

(t) = a

(t) x

(t)











. (14.6)

Here, a

, a

, and a

are constants that depend on the ﬁxed principal moments

of inertia, I

, I

, and I

, according to

− I

, a

− I

, a

− I

Note that, by construction,

+ a

= 0. (14.7)

The quantity x

(t)+x

(t) is an invariant for this system; it remains con-

stant ove r time along any solution. We can check this by direct diﬀerentiation,

using the form of the ODEs (14.6) and the property (14.7):



+ x



= 2 x

+ 2 x

= 2 x

+ 2 x

= 2 x

+ a

)

= 0. (14.8)

So, every solution of (14.6) lives on a sphere; that is,

(0) + x

(0) = R

implies that

(t) + x

(t) = R

for all t > 0. Exercise 14.5 asks you to check that the quantity

(14.9)

is also an invariant for this system. It follows that solutions are further con-

strained to live on ﬁxed ellipsoids. So, overall, any solution lies on a curve

formed by the intersection of a sphere and an ellipsoid. Examples of such curves

are shown in Figure 14.1, which we will return to soon.

14.3 Quadratic Invariants 199

Generally, we will say that an ODE system x

(t) = f (x) has a quadratic

invariant if the function

C x (14.10)

is preserved, where C ∈ R

m×m

is a symmetric matrix. For the rigid body

example (14.6) we have just seen that there are two quadratic invariants, with

C =





1 0 0

0 1 0

0 0 1





and C =





1/I

0 0

0 1/I

0 0 1/I





. (14.11)

By deﬁnition, if (14.10) is an invariant then its time derivative is ze ro, giving

0 =



C x



= x

C x + x

C x

= 2x

C x

= 2x

C f (x),

where we have used C = C

and x

(t) = f (x). This shows that x

C x is a

quadratic invariant if and only if the function f satisﬁes

C f (x) = 0, for any x ∈ R

. (14.12)

The implicit mid-point rule is deﬁned by

n+1

= x

+ hf



+ x

n+1



. (14.13)

This method achieves the commendable feat of preserving quadratic invariants

of the ODE for any choice of stepsize h. To see this, we will simplify the notation

by letting

mid

:= f



+ x

n+1



Then one step of the implicit midpoint rule produces an approximation x

n+1

for which

n+1



+ hf

mid





+ hf

mid



= x

+ 2hx

mid

+ h

mid

From (14.13), we may write

+ x

n+1

−

mid

and so

n+1

= x

+ 2h



+ x

n+1

−

mid



mid

+ h

mid

= x

+ 2h



+ x

n+1



mid

200 14. Geometric Integration Part I—Invariants

But the last term on the right disappears when we invoke the condition (14.12).

Hence,

n+1

= x

conﬁrming that the numerical method shares the quadratic invariant.

By contrast, applying Euler’s method

n+1

= x

+ hf

, (14.14)

we ﬁnd, with f

denoting f (x

), that

n+1



+ hf



C (x

+ hf

)

= x

+ 2hx

+ h

= x

+ h

The extra term h

will be nonzero in general, and hence Euler’s method

does not preserve quadratic invariants of the ODE. In fact, for the two types of

matrix C in (14.11)

we have f

> 0 for any f

6= 0

. So, on the rigid body

problem (14.6) Euler’s method produces an approximation that drifts outwards

onto increasingly large spheres and ellipsoids. This is a more general case of

the behaviour observed in Examples 7.6 and 13.3.

These results are conﬁrmed experimentally in Figure 14.1, which is inspired

by Hairer et al. [26, Figure 1.1, Chapter IV, page 96]. Here, we illustrate the

behaviour of the implicit mid-point rule (left) and Euler’s method (right) on

the rigid body problem (14.6). In each picture, the dots denote the numerical

approximations, plotted in the x

, x

phase space. In this case we have

(0)

+ x

(0)

+ x

(0)

= 1, and the corresponding unit spheres are indicated

in the pictures. The solid curves mark the intersection between the sphere and

various ellipsoids of the form (14.9) with I

= 2, I

= 1 and I

= 2/3. Any

exact solution to the ODE remains on such an intersection, and we see that the

numerical solution provided by the implicit mid-point rule shares this property.

The Euler approximation, however, shown on the right, is s ee n to violate thes e

constraints by spiralling away from the surface of the sphere.

Example 14.1 (Kepler Problem)

The Kepler two-body problem describes the motion of two planets. With ap-

propriate normalization, this ODE system takes the form [26, 61]

(t) = −

+ q

)

3/2

, q

(t) = p

(t) = −

+ q

)

3/2

, q

(t) = p

(14.15)

And, more generally, for any symmetric positive deﬁnite matrix C.

14.3 Quadratic Invariants 201

Fig. 14.1 Dots show numerical approximations for the implicit mid-point rule

(left) and Euler’s method (right) on the rigid bo dy problem (14.6). Intersections

between the unit sphere and ellipsoids of the form (14.9) are marked as solid

lines. (This ﬁgure is based on Figure 1.1 in Chapter IV of the book by Hairer

et al. [26])

where q

(t) and q

(t) give the position of one planet moving in a plane relative

to the other at time t.

This system has a quadratic invariant of the form

− q

, (14.16)

which corresponds to angular momentum. This can be checked by direct dif-

ferentiation:

− q

) = q

+ q

− q

= p

−

+ q

)

3/2

− p

+ q

)

3/2

= 0.

However, the system also has a nonlinear invariant given by the Hamiltonian

function

H(p

, p

, q

) =



+ p



−

+ q

, (14.17)

(see Exercise 14.11), so called because the original ODEs (14.15) may be written

in the form

(t) = −

∂

∂q

H(p

, p

, q

), q

(t) =

∂

∂p

H(p

, p

, q

(t) = −

∂

∂q

H(p

, p

, q

), q

(t) =

∂

∂p

H(p

, p

, q

(14.18)

Hamiltonian ODEs are considered in the next chapter

We will consider problems where H : R

→ R, but the same ideas extend to the

more general case of H : R

→ R.

202 14. Geometric Integration Part I—Invariants

14.4 Modiﬁed Equations and Invariants

Supp ose now that the ODE system x

(t) = f (x), with x(t) ∈ R

and f :

→ R

, has a general nonlinear invariant deﬁned by some smooth function

F : R

→ R, so that F(x(t)) remains constant along any solution. By the

chain rule, this means that, for any x(t),

0 =

F(x(t)) =

i=1

∂F(x(t))

∂x

(f(x(t)))

We may write this as

(∇F(x))

f(x) = 0, for any x ∈ R

, (14.19)

where ∇F denotes the vector of partial derivatives of F. This is a direct gen-

eralization of the linear invariant case (14.4) and the quadratic invariant case

(14.12).

Supp ose we have a numerical method of order p that is able to preserve the

same invariant. So, from any initial value x(0), if we take N steps to re ach time

= Nh then the global error satisﬁes

− x(t

) = O(h

), (14.20)

and, because the invariant is preserved, we have

F(x

) = F(x(0)). (14.21)

Now, following the ideas in Chapter 13, we know that it is generally possible

to construct a modiﬁed IVP of the form

(t) = f (y) + h

g(y), y(0) = x(0), (14.22)

such that the numerical method approximates the solution y(t) even more

accurately than it approximates the original s olution x(t). In general, by adding

the extra term on the right-hand side in (14.22), we are able to increase the

order by one, so that

− y(t

) = O(h

p+1

). (14.23)

We will show now that this modiﬁed equation automatically inherits the in-

variance property of the ODE and numerical method. This can be done through

a contradiction argument. Suppose that the modiﬁed equation (14.22) does not

satisfy the invariance property (14.19). Then there must be at least one point

∈ R

where (∇F(x

))

g(x

) 6= 0. By switching from F to −F if

neces-

sary, we may assume that this nonzero value is positive and then, by continuity,

there must be a region B containing x

for which

(∇F(x))

g(x) > C, for all x ∈ B,

14.5 Discussion 203

where C > 0 is a constant. If we choos e x

as our initial condition and let [0, t

]

be a time interval over which y(t) remains inside the region B, then,

F(y(t)) = h

(∇F(y(t)))

g(x) > Ch

, for t ∈ [0, t

After integrating both sides with respect to t from 0 to t

, we ﬁnd

|F(y(t

)) − F(y(0))| > Ct

. (14.24)

On the other hand, since the method preserves the invariant, we have F(x

) =

F(y(0)), so

F(y(t

)) − F(y(0)) = F(y(t

)) − F(x

Because F is smooth and we are working in a compact set B, the diﬀerence on

the right may be bounded by a multiple of the diﬀerence in the arguments

So, using (14.23), we have

F(y(t

)) − F(y(0)) = O(h

p+1

However, this contradicts the bound (14.24). So we conclude that the modiﬁed

equation must preserve the invariant.

14.5 Discussion

It is natural to ask which standard numerical methods can preserve quadratic

invariants. In the case of Runge–Kutta methods, (9.5)) and (9.6), there is a

very neat result. Cooper [10] proved that the condition

+ b

= b

, for all i, j = 1, . . . , s, (14.25)

characterizes successful methods. Of course, it is also possible to construct ad

hoc methods that deal with speciﬁc classes of ODEs with quadratic, or more

general, invariants.

The result that a modiﬁed equation inherits invariants that are shared by

the ODE and the numerical method can be established very generally for one-

step methods, and the same principle applies for other properties, including

reversibility, volume preservation, ﬁxed point preservation and, as we study in

the next chapter, symplecticness. The proof by contradiction argument that we

used was taken from Gonzalez et al. [22], where those other properties are also

considered.

More precisely, we may assume that a Lipschitz condition |F(x)−F(y)| ≤ Lkx−

yk holds for x, y ∈ B.

204 14. Geometric Integration Part I—Invariants

The modiﬁed equation concept plays a central role in the modern study

of geometric integration. By optimizing over the number of extra terms in the

modiﬁed equation, it is possible to show that the diﬀerence between a numerical

solution and its closest modiﬁed equation can be bounded by C exp(−D/h)

over a time interval 0 ≤ t ≤ E/h, where C, D and E are constants. This is

an exponentially small bound that applies over arbitrarily large time; a very

rare phenomenon in numerical ODEs!

Establishing the existence of a modiﬁed

equation with the same structure as the original ODE is often a key step in

proving further positive results, such as mild growth of the global error as a

function of time.

Finally, we should add that not all prop e rties of an ODE and numerical

method are automatically inherited by a modiﬁed equation; see Exercise 14.12

for an illustration.

EXERCISES

14.1.

Conﬁrm that (14.3) satisﬁes the ODE system (14.2). Also check

that this solution gives u(t) + v(t) ≡ A + B. What happens to this

solution as t → ∞? Explain the result intuitively in terms of the

reaction rate constants k

and k

14.2.

By generalizing the analysis for the Euler case in (14.5), show

that any consistent linear multistep method and any Runge–Kutta

method will preserve linear invariants of an ODE. You may assume

that the starting values for a k-step LMM satisfy c

= K, j = 0 :

k − 1, for some constant K when c is any vector such that equation

(14.4) holds.

14.3.

Show that the second-order Taylor series method, TS(2), from

Chapter 3, applied to (14.2) takes the form



n+1







+ h



1 −

h(k

+ k

)





−k





Conﬁrm that this method also preserves the linear invariant.

14.4.

For the system (14.2) show that u

(t) = Au(t), where A may be

written as the outer product

A =



−1





, −k



We should note, however, that the bound does not involve the solution of the

original ODE. Over a long time interval the modiﬁed equation might not remain close

to the underlying problem.

Exercises 205

Deduce that A



−k

− k



j−1

A (j = 1, 2, 3, . . . ) and hence gen-

eralise the previous exercise to the T S(p) method for p > 2.

14.5.

By diﬀerentiating directly, as in (14.8), show that the quadratic

expression (14.9) is an invariant for the ODE system (14.6).

14.6.

For (14.6) conﬁrm that (14.12) holds for the two choices of matrix

C in (14.11).

14.7.

Show that the implicit mid-point rule (14.13) may be regarded as

an implicit Runge–Kutta method with Butcher tableau (as deﬁned

in Section 9.2) given by

14.8.

???

By following the analysis for the implicit mid-point rule, show that

the implicit Runge–Kutta method with Butcher tableau

−

√

−

√

called the order 4 Gauss method, als o preserves quadratic invariants

of an ODE system for any step size h.

14.9.

The backward Euler method produces x

n+1

= x

+ hf (x

n+1

instead of (14.14). In this case show that

= x

n+1

+ h

when (14.12) holds for the ODE. Hence, explain how the picture

on the right in Figure 14.1 would change if Euler’s method were

replaced by backward Euler.

14.10.

Show that the invariant (14.16) may be written in the form

















where C ∈ R

4×4

is symmetric, so that it ﬁts into the format required

in (14.10).