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

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

56 4. Linear Multistep Methods—I

h =0.1 h =0.01 h =0.001

=0.544 x

=0.938 x

=1.070

=0.199 x

=0.567 x

=0.535

=1.735 x

=2.384 x

=3.205

= −6.677 x

= −6.810 x

= −10.159

=37.706 x

= 39.382 x

= 56.697

=−197.958 x

=−193.017 x

=−277.788

Table 4.3 The numerical solutions for Example 4.11

Since f

= −x

and f

n+1

= −x

n+1

the LMM leads to the diﬀerence equa-

tion

n+2

= −4(1 + h)x

n+1

+ (5 − 2h)x

,n=0, 1, 2,...,

which is used to compute the values given in Table 4.3. In each case we see

that the numerical solution oscillates wildly within just a few steps, making

the method completely worthless.

Behaviour of the type seen in the preceding example had been e xperienced

frequently ever since LMMs were ﬁrst introduced in the 19th century, but

it took until 1956 for the Swedish numerical analyst Germund Dahlquist to

discover the cause—this example was use d in his original work and reapp ears

in his book with A. Bj¨orck [16]. He proved that for a method to be convergent

(and therefore useful) it not only had to be c onsistent (have order ≥ 1), but it

also had to satisfy a property known as zero-stability. This is addressed in the

next chapter.

We conclude this chapter by extending two-step LMMs so as to incorporate

further terms from the history of the numerical solution. This generalization

will allow us the opportunity of summarizing the main ideas to date.

4.3 k-Step Methods

The one- and two-step methods that have been discussed thus far generalize

quite naturally to k-step methods. The most general method takes the form

n+k

+α

k−1

n+k−1

+···+α

= h

�

n+k

+β

k−1

n+k−1

+···+β

�

(4.22)

and is of implicit type unless β

= 0, when it becomes explicit. It uses k past

values of the pair x

n+j

(j =0:k − 1), as well as f

n+k

in the implicit case,

in order to calculate x

n+k

. The coeﬃcients have been normalized so that the

coeﬃcient of x

n+k

is α

= 1. This method has ﬁrst and second characteristic

4.3 k-Step Methods 57

polynomials (see Deﬁnition 4.6)

ρ(r) = r

+ α

k−1

+ ··· + α

σ(r) = β

+ β

k−1

+ ··· + β

(4.23)

and an associated linear diﬀerence operator deﬁned by (see Deﬁnition 4.2)

z(t) ≡

j=0

z(t + jh) − hβ

(t + jh). (4.24)

Taylor expanding the right-hand side about the point h = 0 and collecting

terms in powers of h leads to (see (4.17))

z(t) = C

z(t) + C

(t) + ···+ C

(p)

(t) + C

p+1

(p+1)

(t) + O(h

p+2

in which C

= ρ(1), C

= ρ

(1) − σ(1), and the coeﬃcients C

, C

, . . . are

linear combinations of β

, β

, . . . , β

, α

, . . . , α

(α

= 1). The method has

order p if

= C

= ··· = C

= 0

and the ﬁrst non-zero coeﬃcient C

p+1

is the error constant.

The general implicit (explicit) k-step LMM has 2k + 1 (2k) arbitrary coeﬃ-

cients; so, provided the linear relationships are linearly independent, we would

expect to be able to achieve order 2k with implicit methods and (2k − 1) with

explicit methods. We do not oﬀer proofs of these assertions since it will turn

out in the next chapter that, because of the type of instability observed in Ex-

ample 4.11, convergent methods cannot, in general, achieve such high orders.

EXERCISES

4.1.

Distinguish the implicit methods from the explicit methods in Ta-

ble 4.2.

4.2.

Consider the problem of using the trapezoidal rule to solve the

IVP x

(t) = −x

(t), x(0) = 1. Show that it is necessary to solve

a quadratic equation in order to determine x

n+1

from x

and that

an appropriate root can be identiﬁed by making use of the property

that x

n+1

→ x

as h → 0.

Hence ﬁnd an approximate solution at t = 0.2 using h = 0.1.

It is common to deﬁne the error constant to be C

p+1

/σ(1)—we shall call this the

scaled error constant. See the footnote on page 69 for an explanation.

58 4. Linear Multistep Methods—I

4.3.

Consider the IVP x

(t) = 1 + x

(t), x(0) = 0. When the backward

Euler method is applied to this problem the numerical solution x

at time t

= h is deﬁned implicitly as the solution of a quadratic

equation. Explain carefully why it is appropriate to choose the root

given by

1 +

√

1 − 4h

4.4.

The backward Euler method is to b e applied to the IVP x

(t) =

x(t), with x(0) = 1. Use an argument similar to that in Exer-

cise 4.2 to explain carefully why this leads to

n+1



h +

+ h



for a convergent process.

4.5.

Write down the linear diﬀerence operator L

associated with the

backward Euler method and, working from ﬁrst principles, ﬁnd the

precise form of the leading term in the LTE, i.e. ﬁnd its order and

error constant.

4.6.

Investigate the consistency of the LMM

n+2

− x

h(3f

n+1

− f

4.7.

What values should the parameters a and b have so that the following

LMMs are consistent:

(a) x

n+2

− ax

n+1

− 2x

= hbf

(b) x

n+2

+ x

n+1

+ ax

= h(f

n+2

+ bf

4.8.

If β

were to be kept free in Example 4.10, say β

= θ, show that

this would lead to the one-parameter family of order 1 methods:

n+2

− x

= h(θf

n+1

+ (2 − θ)f

What is the error constant of this method? What value of θ produces

the highest order?

4.9.

Show that the LMM x

n+1

= x

+ 2hf

is not consistent and that,

when used to solve the IVP x

(t) = 1 for t ∈ (0, 1], x(0) = 0, it

leads to x

= 2nh. Hence , by writing down the GE x(t

) −x

when

= 1, show that the method is not convergent.

Exercises 59

4.10.

Determine the order and error constant of the LMM

n+2

− 4x

n+1

+ x

= 2hf

n+2

, (4.25)

which is known as the “two-step backward diﬀerentiation formula”:

BDF(2).

4.11.

Determine the two-step LMM of the form

n+2

+ α

n+1

− ax

= hβ

n+2

that has maximum order when a is retained as a free parameter.

4.12.

Show that an LMM of the form

n+2

− x

n+1

= h(β

n+2

+ β

n+1

+ β

)

has maximum order when β

= 5/12, β

= 2/3, and β

= −1/12.

What is the error constant of this method (known as the two-s tep

Adams–Moulton method (AM(2))?

4.13.

???

Write down the linear diﬀerence operator L

associated with

Simpson’s rule

n+2

− x

h(f

n+2

+ 4f

n+1

+ f

Use the Taylor expansions of z(t+2h) and z(t) about the point t+h:

z(t + 2h) = z(t + h) + hz

(t + h) +

(t + h) + . . . ,

z(t) = z(t + h) − hz

(t + h) +

(t + h) − . . . ,

with similar expansions for z

(t + 2h) and z

(t), to show that it has

order p = 4 (the maximum possible for a two-step method) with

error constant C

= −1/90.

What is the advantage of using these non-standard expansions?

4.14.

???

Show that the general coeﬃcient C

(m > 0) in the expansion

(4.17) is given by

j=0



−

(m − 1)!

m−1



, (4.26)

where α

= 1 is the normalizing condition.

Use this expression to verify the error constants found in Exam-

ples 4.4 and 4.10.

60 4. Linear Multistep Methods—I

4.15.

Show that the conditions (4.15) for consistency are equivalent to

the equations L

z(t) = 0 with the choices z(t) = 1 and z(t) = t.

Hence prove that L

z(t) = 0 whenever z(t) is a linear function of t.

4.16.

If x

(t) = g(t), t ∈ [a, b], with x(a) = 0, show that x(b) =

g(t) dt.

Use the trapezoidal rule with h = 1/2 to compute an approximate

value for the integral

dt and compare with the exact value.

4.17.

???

Complete the details of the proof of the second part of The-

orem 4.8. Explain, in particular, why it is necessary to invoke

l’Hˆopital’s rule.

4.18.

Determine the quadratic polynomial

p(t) that satisﬁes the condi-

tions

p(t

n+1

) = x

n+1

, p

n+1

) = f

n+1

, p

) = f

Show that x

n+2

= p(t

n+2

) leads to the AB(2) method (4.10).

4.19.

Find a quadratic polynomial p(t) that satisﬁes the conditions

p(t

) = x

, p

) = f

, p

n+1

) = f

n+1

Show that x

n+1

= p(t

n+1

) leads to the trapezoidal rule (4.7).

The construction of p(t) can be simpliﬁed by seeking coeﬃcients A, B, and C

such that p(t) = A + B(t − t

n+1

) + C(t − t

n+1

)

Linear Multistep Methods—II:

Convergence and Zero-Stability

5.1 Convergence and Zero-Stability

Means of determining the coeﬃcients in LMMs were described in Chapter 4

and criteria now need to be established to identify those methods that are

practically useful. In this section we describe some of the behaviour that should

be expected of methods (in general) and, in subsequent sections, indicate how

this behaviour can be designed into LMMs.

A basic requirement of any method is that its solutions should converge

to those of the corresponding IVP as h → 0. To formalize this concept, our

original speciﬁcation in Deﬁnition 2.3 must be tailored to accommodate the

additional starting values needed for k-step LMMs.

In order to solve the IVP

(t) = f (t, x(t)), t > t

x(t

) = η



(5.1)

over some time interval t ∈ [t

, t

], we choose a step-size h, a k-step LMM

n+k

+α

k−1

n+k−1

+···+α

= h



n+k

+β

k−1

n+k−1

+···+β



, (5.2)

and starting values

= η

, x

= η

, . . . , x

k−1

= η

k−1

. (5.3)

Springer Undergraduate Mathematics Series, DOI 10.1007/978-0-85729-148-6_5,

D.F. Griffiths, D.J. Higham, Numerical Methods for Ordinary Differential Equations,

62 5. Linear Multistep Methods—II

The values x

, x

k+1

, ..., x

are then calculated from (5.2) with n = 0 : N − k,

where Nh = t

− t

Our aim is to develop methods with as wide a range of applicability as

possible; that is, we expect the methods we develop to b e capable of solving

all IVPsof the form (5.1) that have unique solutions. We have no interest in

methods that can solve only particular IVPsor particular types of IVP.

In addition to the issues related to convergence that we previously discussed

in Section 2.4 (for Euler’s method) and Section 3.4 (for TS(p) methods), the

main point to be borne in mind when dealing with k-step LMMs concerns the

additional starting values (5.3). These will generally contain some leve l of error

which must tend to zero as h → 0, so that

lim

h→0

= η, j = 0 : k − 1. (5.4)

In practice, the additional starting values x

, . . . , x

k−1

would be calculated

using an appropriate numerical method. For example, condition (5.4) would be

satisﬁed if we used k − 1 steps of Euler’s method.

Deﬁnition 5.1 (Convergence)

The LMM (5.2) with starting values satisfying (5.4) is said to be convergent if,

for all IVPs (5.1) that possess a unique solution x(t) for t ∈ [t

, t

lim

h→0

nh=t

∗

−t

= x(t

∗

) (5.5)

holds for all t

∗

∈ [t

, t

The next objective is to establish conditions on the coeﬃcients of the general

LMM that will ensure convergence. Theorem 4.8 shows that consistency is a

necessary prerequisite, but Example 4.11 strongly suggests that, on its own, it

is not suﬃcient. Recall that consistency implies that ρ(1) = 0 and ρ

(1) = σ(1):

j=0

= 0,

j=0

jα

j=0

, (5.6)

with our normalizing condition α

= 1.

Example 5.2

Explain the the non-convergence of the two-step LMM

n+2

+ 4x

n+1

− 5x

= h(4f

n+1

+ 2f

)

Although it is natural to choose η

= η, the given starting value for the ODE,

this is not necessary for convergence.

5.1 Convergence and Zero-Stability 63

that was observed in Example 4.11 by applying the method to the (trivial)

IVP x

(t) = 0 with initial condition x(0) = 1. Use starting values x

= 1

and x

= 1 + h (recall that the deﬁnition of convergence requires only that

→ x(0) as h → 0).

The method becomes, in this case,

n+2

+ 4x

n+1

− 5x

= 0. (5.7)

This is a two-step constant-coeﬃcient diﬀerence equation whose auxiliary equa-

tion is

+ 4r − 5 = (r − 1)(r + 5).

The general solution of the diﬀerence equation (see Appendix D) is

= A + B(−5)

where A and B are arbitrary constants. The initial condition x

= 1 implies

that A + B = 1 and x

= 1 + h leads to

1 + h = A + B(−5).

These solve to give B = −h/6 and A = 1 − h/6, and so

= 1 +

h[1 − (−5)

It is the presence of the term (−5)

that causes the disaster: supp ose , for

instance, that t = 1, so nh = 1, then

h|(−5)

| =

→ ∞ as h → 0.

When h = 0.1, for example, n = 10 and

≈ 10

, so the divergence of x

is,

potentially, very rapid.

The auxiliary equation of (5.7) is the ﬁrst characteristic polynomial, ρ(r),

of the LMM. It has the property ρ(1) = 0, i.e. r = 1 is a root of ρ(r) for all

consistent methods. This means that the ﬁrst characteristic polynomial of any

consistent two-step LMM will factorize as

ρ(r) = (r − 1)(r − a)

for some value of a (in the previous example we had a = −5, which led to the

trouble). A method with this characteristic polynomial applied to x

(t) = 0

would give a general solution

= A + Ba

which suggests that we s hould restrict ourselves to methods for which |a| ≤ 1

so that |a|

does not go to inﬁnity with n. This turns out to be not quite

suﬃcient, as the next example shows.

64 5. Linear Multistep Methods—II

Example 5.3

Investigate the convergence of the three-step LMM

n+3

+ x

n+2

− x

n+1

− x

= 4hf

when applied to the model problem x

(t) = 0, x(0) = 1 with starting values

= 1, x

= 1 − h, and x

= 1 − 2h.

Establishing consistency of the method is left to Exercise 5.7.

The homogeneous diﬀerence equation x

n+3

+ x

n+2

− x

n+1

− x

= 0 has

auxiliary equation

ρ(r) = (r − 1)(r + 1)

and, therefore, its general solution is

= A + (B + Cn)(−1)

(see Appendix D). With the given starting values, the solution can be shown

to be

= 1 − h + (−1)

(h − t

), (5.8)

while the e xact solution of the IVP is, of course, x(t) = 1. Thus, for example,

at t = t

∗

= 1, x(1) = 1 while, |x

− 1 + h| = 1 − h for all values of h with

∗

= nh = 1. The GE, therefore, has the property |x(1) −x

| → 1 and does not

tend to zero as h → 0. Hence, the method is consistent but not convergent.

This leads to the following deﬁnition.

Deﬁnition 5.4 (Root Condition)

A polynomial is said to satisfy the root condition if all its roots lie within or

on the unit circle, with those on the boundary being simple. In other words,

all roots satisfy |r| ≤ 1 and any that satisfy |r| = 1 are s imple.

A polynomial satisﬁes the strict root condition if all its roots lie inside the

unit circle; that is, |r| < 1.

Thus, ρ(r) = r

− 1 satisﬁes the root condition, while ρ(r) = (r − 1)

does not.

Deﬁnition 5.5 (Zero-Stability)

An LMM is said to be zero-stable if its ﬁrst characteristic polynomial ρ(r)

satisﬁes the root condition.

We say that λ is a simple root of ρ(r) if λ − r is a factor of ρ(r), but (λ − r)

not.

5.2 Classic Families of LMMs 65

All consistent one-step LMMs have ﬁrst characteristic polynomial ρ(r) = r −1

and so satisfy the root condition automatically, which is why this notion was

not needed in the study of Euler’s method or the Taylor series methods. We

are now in a position to state Dahlquist’s celebrated theorem.

Theorem 5.6 (Dahlquist (1956))

An LMM is convergent if, and only if, it is both consistent and zero-stable.

The main purpose of this theorem is to ﬁlter out the many LMMs that fail its

conditions; we are left to focus our attention on those that pass—the convergent

methods. An explanation of why both conditions together are suﬃcient is given

in Section 5.3 in the context of a sc alar linear IVP.

Zero-stability places a signiﬁcant restriction on the attainable order of

LMMs, as the next theorem attests—recall that the order of an implicit (ex-

plicit) k-step LMM could be as high as 2k (2k − 1).

Theorem 5.7 (First Dahlquist Barrier (1959))

The order p of a stable k-step LMM satisﬁes

1. p ≤ k + 2 if k is even;

2. p ≤ k + 1 if k is odd;

3. p ≤ k if β

≤ 0 (in particular for all explicit methods).

Proofs of Theorems 5.6 and 5.7 were originally given in the landmark papers

of Dahlquist [14, 15]. They may also be found in the book of Hairer et al. [28].

Had we been armed with this barrier theorem earlier, the method described

in Example 4.11 could have been dismissed immediately, since it violates part

3 of the theorem (it is explicit with k = 2 and p = 3).

5.2 Classic Families of LMMs

The ﬁrst four families of classical metho ds we describe below are constructed by

choosing methods that have a particularly simple ﬁrst characteristic polynomial

that will automatically be zero-stable. The remaining coeﬃcients are chosen so

that the resulting method will have maximum order.