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

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

66 5. Linear Multistep Methods—II

1. Adams–Bashforth (1883): These have ﬁrst characteristic polynomials

ρ(r) = r

− r

k−1

which have a simple root r = 1 and a root of multiplicity (k − 1) at r = 0.

Methods in this family are also explicit, so

n+k

− x

n+k−1

= h



k−1

n+k−1

+ ··· + β



and the k coeﬃcients β

, . . . , β

k−1

are chosen so that C

= C

= ··· =

k−1

= 0 and give, therefore, a method of order p = k. This is the

highest possible order of a zero-stable explicit k-step method (part 3 of

Theorem 5.7) and is attained despite imposing severe restrictions on its

structure.

We have already encountered the k = 1 member—Euler’s method—and

the AB(2) method (4.10) for k = 2. With k = 3 the AB(3) method is

n+3

− x

n+2



23f

n+2

− 48f

+ 5f



, (5.9)

which has order p = 3 and error constant C

2. Adams–Moulton (1926): These are implicit versions of the Adams–Bashforth

family, having the form

n+k

− x

n+k−1

= h



n+k

+ ··· + β



There is one more coeﬃcient than the corresponding Adams–Bashforth

method so we can attain order k + 1. For k = 1 we have the trape-

zoidal rule (4.7). The two-step Adams–Moulton method (AM(2)—see Ex-

ercise 4.12) has order 3 and not the maximum possible order (4) for an

implicit two-step method implied by Theorem 5.7. For k = 3, the AM(3)

method

n+3

− x

n+2



n+3

+ 19f

n+2

− 5f

n+1

+ f



has order 4 and error constant C

= −19/720.

3. Nystr¨om method (1925):s These are explicit methods with k ≥ 2 having

ﬁrst characteristic polynomial ρ(r) = r

− r

k−2

, so take the general form

n+k

− x

n+k−2

= h



k−1

n+k−1

+ ··· + β



They have the same number of free coeﬃ cients as the AB family and are

chosen to achieve order p = k, the s ame order as the corresponding AB

method. The two-step version is the sec ond-order mid-point rule (4.21).

5.3 Analysis of Errors: From Local to Global 67

4. Milne–Simpson (1926): These are the implicit analogues of Nystr¨om meth-

ods:

n+k

− x

n+k−2

= h



n+k

+ ··· + β



The simplest member of the family is Simpson’s rule (see Exercise 4.13),

which has k = 2 and order p = 4. The high degree of symmetry present

in the structure of the method allows it to attain the maximum possible

order for a zero-stable two-step method (unlike AM(2) above).

5. Backward diﬀerentiation formulas (BDFs—1952). These form a generaliza-

tion of the backward Euler method (see Table 4.2). The simplest possible

form is chosen for the second characteristic polynomial consistent with the

method being implicit: σ(r) = β

. Thus,

n+k

+ α

k−1

n+k−1

+ ··· + α

= hβ

n+k

and the (k+1) free coeﬃcients are chosen to achieve order k—not the higher

order (k+2) that can be achieved with the most general implicit zero-stable

k-step LMMs. Despite this, they form an important family because they

have compensating strengths, as we shall see in Chapter 6.

It may be shown (see Iserles [39, Lemma 2.3], for instance) that the ﬁrst

characteristic polynomial, ρ(r), of these methods is given by

ρ(r) =

j=1

k−j

(r − 1)

, c =

j=1

. (5.10)

Zero-stability is easily checked for each k—the BDF(2) method of Exer-

cise 4.10 is stable, as are other members up to k = 6; thereafter they are

all unstable.

5.3 Analysis of Errors: From Local to Global

The convergence of Euler’s method was analysed in Theorem 2.4 for the special

case when it was applied to the IVP

(t) =λx(t) + g(t), 0 < t ≤ t

x(0) = 1.

(5.11)

In this section we shall indicate the additional features that have to be ac-

counted for when analysing the convergence of k-step LMMs with k > 1. To

get across the key concepts, it is suﬃcient to look at the behaviour of two-step

methods.

68 5. Linear Multistep Methods—II

As discussed in Section 4.3, LMMs are constructed from Taylor expansions

in such a manner that their associated linear diﬀerence operators have the

property

z(t) = C

p+1

(p+1)

(t) + ··· ,

where z is any (p + 1)-times continuously diﬀerentiable function.

The LTE, denoted by T

n+2

, is deﬁned to be the value of the left-hand side

of this expression when we replace z by x, the exact solution of our IVP. So,

for our two-step method we deﬁne, at t = t

n+2

= L

x(t

). (5.12)

When x(t) is a (p + 1)-times continuously diﬀerentiable function we have

n+2

= C

p+1

(p+1)

(t) + ··· ,

so that T

n+2

= O(h

p+1

). Let y

= x(t

) (for all n) denote the exact solution of

the IVP (5.11) at a typical grid p oint t

. The general two-step LMM applied

to the IVP (5.11) leads to

n+2

+ α

n+1

+ α

hλ



n+2

+ β

n+1

+ β



+ h(β

g(t

n+2

) + β

g(t

n+1

) + β

g(t

)). (5.13)

With Deﬁnition 4.2 for L

and equation (5.12) we ﬁnd that the exact solution

satisﬁes the same equation with the addition of T

n+2

on the right:

n+2

+ α

n+1

+ α

hλ



n+2

+ β

n+1

+ β



+ h(β

g(t

n+2

) + β

g(t

n+1

) + β

g(t

)) + T

n+2

. (5.14)

Subtracting (5.13) from (5.14), the GE e

= x(t

) − x

≡ y

− x

is found to

satisfy the diﬀerence equation

(1 − hλβ

n+2

+ (α

− hλβ

n+1

+ (α

− hλβ

= T

n+2

, (5.15)

with the starting values e

= 0 and e

= x(t

)−η

, which may be assumed to be

small (the deﬁnition of convergence stipulates that the error in starting values

must tend to zero as h → 0). This is the equation that governs the way “local”

errors T

n+2

accumulate into the “global” error {e

}. Usage of the terms “local”

and “global” reﬂects the fact that the LTE T

n+2

can be calculated locally (for

each time t

) via Equation (5.12), whereas the GE results from the overall

accumulation of all LTEs—a global process.

In order to simplify the next stage of the analysis, let us assume that the

LTE terms T

n+2

are constant: T

n+2

= T , for all n. As discussed in Appendix D,

the general solution of (5.15) is then comprised of two components:

5.4 Interpreting the Truncation Error 69

1. A particular solution, e

= P , which is constant. Substituting e

n+2

n+1

= e

= P into (5.15) and using the consistency condition 1+α

+α

0, we ﬁnd that

P =

hλσ(1)

;

so, if T = O(h

p+1

), it follows that P = O(h

2. The general solution of the homogeneous equation. If the auxiliary equation

(1 − hλβ

+ (α

− hλβ

)r + (α

− hλβ

) = 0

has distinct roots, r

6= r

then this contribution is

+ Br

where A and B are arbitrary constants.

Thus, the general solution for the GE is, in this case,

= Ar

+ Br

+ P. (5.16)

The constants A and B are determined from the starting values. As h → 0,

the roots r

and r

tend to the roots of the ﬁrst characteristic polynomial ρ(r).

If the root condition in Deﬁnition 5.4 is violated, then their contribution will

completely swamp the GE and lead to divergence. The LTE contributes to the

GE through the term P = O(h

), so consistency (p > 0) ensures that this

tends to zero with h.

Consistency of an LMM ensures that the local errors are small, while zero-

stability ensures that they propagate so as to give small GEs.

5.4 Interpreting the Truncation Error

We extend the discussion in the previous subsec tion so as to give an alternative—

and perhaps more intuitive—interpretation of the LTE (5.12) of an LMM. To

simplify the presentation we will again restrict ourselves to the case of a lin-

ear ODE solved with a two-step method, but the conclusions are valid more

generally.

The presence of the factor σ(1) in the denominator induces many writers to deﬁne

the error constant of an LMM to be C

p+1

/σ(1) and not C

p+1

as we have done. We

shall refer to C

p+1

/σ(1) as the scaled error constant; it is seen on page 89 to be a key

feature of Dahlquist’s Second Barrier Theorem.

70 5. Linear Multistep Methods—II

Supp ose we were to compute x

n+2

from (5.13) when the past values were

exact—that is, x

n+1

= y

n+1

and x

= y

(recall that y

n+j

= x(t

n+j

))—and

we denote the re sult by ex:

n+2

+ α

n+1

+ α

hλ



n+2

+ β

n+1

+ β



+ h(β

g(t

n+2

) + β

g(t

n+1

) + β

g(t

)).

Subtracting this from (5.14) gives

(1 − hλβ

)(y

n+2

− ex

n+2

) = T

n+2

1. In the explicit case (β

= 0), T

n+2

= y

n+2

− ex

n+2

: the LTE is the error

committed in one step on the assumption that the back values are exact.

2. In the implicit case (β

6= 0), the binomial expansion:

(1 − hλβ

)

−

= 1 + hλβ

+ O(h

) = 1 + O(h)

and the fact that T

n+2

= O(h

p+1

) for a method of order p shows that

(1 + O(h))T

n+2

= y

n+2

− ex

n+2

and hence

n+2

= y

n+2

− ex

n+2

+ O(h

p+2

So, the leading term in the LTEis the error committed in one step on the

assumption that the back values are exact.

The assumption that back values are exact is known as the localizing as-

sumption (Lambert [44, p. 56]).

EXERCISES

5.1.

Investigate the zero-stability of the two-step LMMs

(a) x

n+2

− 4x

n+1

+ 3x

= −2hf

(b) 3x

n+2

− 4x

n+1

+ x

= ahf

Are there values of a in part (b) for which the method is convergent?

5.2.

Show that the order of the LMM

n+2

+ (b − 1)x

n+1

− bx

h [(b + 3)f

n+2

+ (3b + 1)f

]

is 2 if b 6= −

1 and 3 if b = −1.

Exercises 71

Show that the method is not zero-s table when b = −1, and illustrate

the divergence by applying it to the IVP x

(t) = 0, x(0) = 0 with

starting values x

= 0 and x

= h.

5.3.

Show that the LMM

n+2

+ 2ax

n+1

− (2a + 1)x

= h((a + 2)f

n+1

+ af

)

has order 2 in general and express the error constant in terms of a.

Deduce that there is one choice of the parameter a for which the

method has order 3 but that this method is not zero-stable. How is

this method related to Example 4.11?

What value of a leads to a convergent method with smallest possible

error constant?

5.4.

Consider the one-parameter family of LMMs

n+2

+ 4ax

n+1

− (1 + 4a)x

= h



(1 + a)f

n+1

+ (1 + 3a)f



for solving the ODE x

(t) = f (t, x(t)).

(a) Determine the error constant for this family of methods and

identify the method of highest order.

(b) For what values of a are members of this family convergent?

Is the method of highest order convergent? Explain your answer

carefully.

5.5.

Prove that any consistent one-step LMM is always zero-stable.

5.6.

Supp ose that the ﬁrst characteristic polynomial ρ(r) of a convergent

LMM is a quadratic polynomial. Prove that the most general form

of ρ is

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

for a parameter a ∈ [−1, 1).

5.7.

For Example 5.3, verify that

(a) the method is consistent;

(b) the numerical solution at time t

is given by the expression (5.8);

This example illustrates the conﬂicting requirements of trying to achieve maxi-

mum order of consistency while maintaining zero-stability. When b = −1 the charac-

teristic polynomials ρ and σ have a common factor—such methods are called reducible

LMMs.

72 5. Linear Multistep Methods—II

5.8.

Prove that the method

n+3

+ x

n+2

− x

n+1

− x

= h(f

n+3

+ f

n+2

+ f

n+1

+ f

)

is consistent but not convergent.

5.9.

Investigate the convergence, or otherwise, of the method

n+3

− x

n+2

+ x

n+1

− x

h(f

n+3

+ f

n+2

+ f

n+1

+ f

5.10.

Prove that σ(1) 6= 0 for a convergent k-s

tep LMM.

5.11.

Consider the two-parameter family of LMMs

n+2

− x

= h (β

n+1

+ β

Show that there is a one-parameter family of convergent methods of

this type.

Determine the error constant for each of the convergent methods and

identify the method of highest order. What is its error constant?

5.12.

???

Apply the LMM x

n+2

− 2x

n+1

+ x

= h(f

n+1

− f

) to the IVP

(t) = −x(t) with x(0) = 1 and the starting values x

= x

= 1

(note that these satisfy our criterion (5.4)). Show that the exact

solution of the LMM in this case is x

= A + B(1 − h)

. Find the

values of A and B and discuss the convergence of the metho d by

directly comparing the expressions for x

and x(t

Extend your solution to the initial values x

= 1, x

= 1 + ah.

Relate your conclusions to the ze ro-stability of the method.

5.13.

Consider the family of LMMs

n+2

+ (θ − 2)x

n+1

+ (1 − θ)x

h ((6 + θ)f

n+2

+ 3(θ − 2)f

) ,

in which θ is a parameter.

Determine the order and error constant for the method and show

that both are independent of θ.

For what range of θ values is the method convergent? Explain your

answer.

5.14.

Determine the coeﬃcients in the three-step BDF(3) method (see

page 67) and calculate its error constant. Verify that it has order

p = 3 and that the coeﬃcients give rise to the ﬁrst characteristic

polynomial (5.10) with k = 3. Examine the zero-stability of the

resulting method.

Exercises 73

5.15.

The ﬁrst characteristic polynomial of the general BDF(k) is given

by (5.10). Use consistency of the method to derive an expression for

the second characteristic polynomial σ(r).

Verify that, when k = 2, these characteristic polynomials give rise

to the BDF(2) method given in Exercise 4.10.

5.16.

???

Use the starting values e

= 0 and e

= x(h) − η

to determine

the constants A and B in (5.16) and hence show that

= (x(h) − η

)

− r

+ P



1 −

1 − r

− r

1 − r

− r



If the ﬁrst characteristic polynomial has roots r = 1 and r = a,

then it may be shown that r

= 1 + O(h) and r

= a + O(h), where

−1 ≤ a < 1 and, consequently, r

− r

= 1−a+O(h) is bounded away

from zero as h → 0. It then follows that e

= O(h

) if P = O(h

)

and the error in the additional starting value is of the same order:

= x(h) + O(h

Linear Multistep Methods—III:

Absolute Stability

6.1 Absolute Stability—Motivation

The study of convergence for LMMs involves the limit

→ x(t

∗

) when n → ∞ and h → 0 in such a way that t

= t

+ nh = t

∗

where t

∗

is ﬁxed value of the “time” variable t.

Thus, convergent methods generate numerical solutions that are arbitrarily

close to the exact solution of the IVP provided that h is taken to be suﬃ-

ciently small. Since non-convergent methods are of little practical use we shall

henceforth assume that all LMMs used are convergent—they are consistent and

zero-stable.

Being in possession of a c onvergent method may not be much comfort in

practice if one ascertains that h needs to b e particularly small to obtain results

of even modest accuracy. If h < 10

−6

, for instance, then more than 1 million

steps would be needed in order to integrate over each unit of time; this may

be entirely impractical. In this chapter we begin the inve stigation into the

behaviour of solutions of LMMs when h is not arbitrarily small. We give a

numerical example before pursing this question.

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

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