Atkinson K. An Introduction to Numerical Analysis

~-I

I

STIFF DIFFERENTIAL EQUATIONS AND

THE

METHOD

OF

LINES

411

Table

6.19

Coefficients of

BDF

method (6.9.6)

p

{3

ao a!

a2

aJ

a4

1

2

4 1

2

-

3

3 3

6

18

9

2

3

11

12

48

36

16

3

4

25

60

300 300

200

75

12

5

137

60

360

450

400

225

72

6

147

147 147

xn+I•···•

xn-p+l

{see

(3.1.5))

.. Use

Combining with (6.9.3) and solving for Y(xn+

1

),

p-1

Y(xn+l)

='=

L

ajY(xn-)

+

h/3/(xn+l•

Y{xn+

1

))

j=O

The p-step

BDF

method

is

given by

p-1

Yn+l

= L

ajYn-j

+

h/3/(xn+l•

Yn+l)

j=O

as

10

147

(6.9.4)

{6.9.5)

{6.9.6)

The coefficients for the cases of p = 1,

...

, 6 are given in Table 6.19. The case

p

= 1 is simply the backward Euler method, which was discussed following

(6.8.54) in the last section. The truncation error for (6.9.6) can be obtained from

the error formula for numerical differentiation, given in (5.7.5):

(6.9.7)

for some

xn-p+l

,:5;

~n

,:5;

xn+l·

The regions of absolute stability for the formulas

of

Table 6.19 are given in

Gear

(1971, pp.

215-216).

To create these regions,

we

must find all values

hA

for

which

llj(h:\)

I<

1

j=0,1,

...

,p

(6.9.8)

I

l

I

•

I

.....

----

--

-

_____

I

412 NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS

where the characteristic roots rj(h'A) are the solutions of

p-1

rP

= L

ajrp-l-j

+

hA.f3rP

j-0

{6.9.9)

It

can

be

shown that for p = 1 and p =

2,

the BDFs are A-stable, and that for

3

~

p

~

6, the region of absolute stability becomes smaller

as

p increases,

although containing the entire negative real axis

in each case.

For

p

;;:::

7, the

regions

of

absolute stability are not acceptable for the solution

of

stiff problems.

For

more discussion of these stability regions, see Gear (1971, chap. 11) and

Lambert (1973, chap.

8).

There are still problems with the BDF methods and with other methods that

are chosen solely on the basis of their region

of

absolute stability. First, with the

model equation

y'

=

A.y,

if Real('A)

is

of large magnitude and negative, then the

solution Y(x) goes to the zero very rapidly, and

as

!Real('A)i increases,

the convergence to zero of

Y(

x)

becomes more rapid.

We

would like the same

behavior to hold for the numerical solution of the model equation {

y,

}.

But with

the A-stable trapezoidal rule, the solution [from (6.5.24)]

is

[

h'A

]"

1+-

Y, = h

2

A

Yo

1--

2

n;::.:O

If

!Real(A)f is large, then the fraction inside the brackets is near to

-1,

andy,

decreases to 0 quite slowly. Referring to the type

of

argument used with the Euler

method in (6.8.45)-(6.8.50), the effect

of

perturbations will not decrease rapidly

to zero for larger values of

A.

Thus the trapezoidal method may not be a

completely satisfactory choice for stiff problems.

In

comparison the A-stable

backward Euler method has the desired behavior.

From

(6.8.55), the solution of

the model problem

is

[

1

..

] "

Yn

=

1

_

hA

Yo

n;::.:O

As

IAI

increases, the sequence

{y,}

goes to zero more rapidly. Thus the

backward Euler solution better reflects the behavior of the true solution of the

model equation.

A second problem with the case of stability regions is that

it

is

based

on

using

constant A

and

linear problems. The linearization (6.9.1) is often valid, but not

always.

For

example, consider the second-order linear problem

y"

+ ay' +

(1

+ b · cos(27rx})y =

g(x)

x;::.:O

(6.9.10)

in which

onl!

coefficient is not constant. Convert it

to

the equivalent system

Y2

=

-(1

+ b · cos{27rx))y

1

-

ay

2

+

g(x)

(6.9.11)

STIFF DIFFERENTIAL EQUATIONS AND THE METHOD

OF

LINES 413

We

will assume a >

0,

1 b 1 < 1.

The

eigenvalues

of

the homogeneous

equation

[g(x)

=

0]

are

-a±

Ja

2

-

4[1 + b

·cos

(2'1Tx)]

;\.

=

---'----------

2

(6.9.12)

These

are

negative real numbers

or

are complex numbers with negative real parts.

On

the

basis

of

the stability theory for the

constant

coefficient

(or

constant

A)

case, we would assume

that

the

effect

of

all perturbations

in

the

initial

data

would

die

away

as x

--+

oo.

But

in fact,

the

homogeneous

part

of

(6.9._10)

will

have

unbounded

solutions.

Thus

there will

be

perturbations

of

the

initial values

that

will

lead

to

unbounded

perturbed

solutions

in

(6.9.10). This calls

into

question

the

validity

of

the

use

of

the

model

equation

y'

= ;\.y +

g(x).

Its use

suggests

methods

that

we

may

want

to

study

further,

but

by

itself, this

approach

is

not

sufficient to encompass

the

vast variety

of

linear

and

nonlinear

problems.

The

example

(6.9.10) is taken from Aiken (1985, p. 269).

Solving

the

finite difference method

We

illustrate

the

problem

by

considering

the

backward

Euler method:

Yn+1

=

Yn

+ hj(xn+1•

Yn+1)

If

the

ordinary

iteration formula

(j+l)

- +

hif(

(j)

)

Yn+I

-

Yn

xn+I• Yn+I

is used,

then

n:::::O

_

U+l),;,

h

aj(xn+l•

Yn+I) ( _

Ul]

Yn+l

ay

Yn+l Yn+l

For

convergence,-we would need

to

have

(6.9.13)

(6.9.14)

(6.9.15)

But

with

stiff equations, this would again force h

to

be

very small, which we are

trying

to

avoid.

Thus

another

rootfinding

method

must

be

used

to

solve for Yn+l

in

(6.9.13).

The

most

popular

methods for solving (6.9.13) are

based

on

Newton's

method.

For

a single differential equation,

Newton's

method

for

finding Yn+l is

{6.9.16)

for

j::::: 0. A

crude

initial guess is

y~~

1

=

Yn•

although generally this

can

be

improved

upon.

With

systems

of

differential equations

Newton's

method

be-

comes

very expensive.

To

decrease

the

expense,

the

matrix

some

z,;,

Yn

(6.9.17)

;

i

I

.

..

---

.... 1

414 NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS

is used for all j and for a number of successive values of

n.

Thus Newton's

method [see

Section

2.11]

for solving the system version of (6.9.13) is approxi-

mated by

{6.9.18)

YHil>

=

y~{\

+

a<n

for j

~

0.

This amounts to solving a number

of

linear systems with the same

coefficient matrix. This can be done much more cheaply than when the matrix is

being modified (see the material in Section 8.1).

The

matrix in (6.9.17) will have

to be

updated

periodically,

but

the savings will still be very significant when

compared to an exact Newton's method.

For

a further discussion of this topic,

see Aiken

(1985, p.

7)

and

Gupta

et al. (1985, pp. 22-25). For a survey of

computer codes for solving stiff differential equations, see Aiken

(1985, chap. 4).

The method

of

lines Consider the following parabolic partial differential equa-

tion problem:

U{O,

t) = d

0

{t)

U(x,O)

=

f(x)

O<x<1

O~x~1

t>O

t

~

0

{6.9.19)

( 6.9.20)

{6.9.21)

The

notation

U,

and

Uxx

refers to partial derivatives with respect to t and

x,

respectively.

The

unknown function U(x,

t)

depends on the time t and a spatial

variable

x.

The

conditions (6.9.20) are called boundary conditions, and (6.9.21)

is

called

an

initial condition. The solution U can

be

interpreted as the temperature

of

an insulated rod of length

1,

with U(x,

t)

the temperature

at

position x

and

time

t;

thus (6.9.19) is often called the heat equation. The functions G, d

0

,

d

1

,

and

f are assumed given and smooth.

For

a development of the theory of

(6.9.19)-(6.9.21), see Widder (1975) or any standard introduction to partial

differential equations. We give the

method

of

lines for solving for

U,

a numerical

method that has become much more popular in the past ten to fifteen years.

It

will also lead to the solution of a stiff system

of

ordinary differential equations.

Let

m > 0

be

an integer, and define 8 =

1/m,

j=0,1,

...

,m

We discretize Eq. (6.9.19) by approximating

the

spatial derivative. Recall the

formulas

(5.7.17) and (5.7.18) for approximating second derivatives. Using this,

8i

a

4

U(

~

1

•

t)

12

ax

4

for j = 1, 2,

...

, m -

1.

Substituting into (6.9.19),

( )

_·

U(xJ+l•

t)-

2U(x

1

,

t) +

U(x

1

_

1

,

t) ( )

U,

x

1

,

t -

82

+ G x

1

,

t

8

2

a

4

u(~

1

.t)

12

ax

4

1~j~m-1

{6.9.22)

STIFF DIFFERENTIAL EQUATIONS AND THE METHOD OF LINES 415

Equation (6.9.I9)

is

to

be approximated at each interior node point

xj.

The

unknown

~j

E

fxj-I•

xj+d·

Drop the

final

term in (6.9.22), the truncation error in the numerical differenti-

ation. Forcing equality in the resulting approximate equation,

we

obtain

(6.9.23)

for j

=I,

2,

...

,

m-

I.

The functions

u/t)

are intended to be approximations

of

U(xp

t ), I

~

j

~

m -

I.

This

is

the method

of

lines approximation to (6.9.I9),

and it

is

a system of m - I ordinary differential equations. Note that u

0

(t)

and

um(t), which are needed in (6.9.23) for j = 1 and j =

m-

1, are given using

(6.9.20):

(6.9.24)

The initial condition for (6.9.23)

is

given by (6.9.21):

I

~j

~

m-

1

(6.9.25)

The name method

of

lines comes from solving for

U(x,

t)

along the lines

(xj.

t),

t

~

0,

1

~

j

~

m - I, in the

(x,

t)-plane.

Under suitable smoothness assumptions

on

the functions d

0

,

d

1

,

G,

and

f,

it

can be shown that

M.ax

ju(xj, t) -

u/t)

I~

Cr~

2

Os;Js;m

Os; cs; T

(6.9.26)

Thus to complete the solution process,

we

need only solve the system (6.9.23).

It

will

be convenient

to

write (6.9.23) in matrix form. Introduce

-2

I 0 0 0

1

-2

I 0

I

0 I

-2

1

A=

82

0

(6.9.27)

I

-2

1

0 0 0

1

-2

The matrix A

is

of order m -

I.

In the definitions of u and

g,

the superscript T

indicates matrix transpose, so that u and g are column vectors of length m -

I.

!

I

j

416 NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS

U:sing these matrices, Eqs. (6.9.23)-(6.9.25) can

be

rewritten as

u'(

t)

=

Au(

t) +

g(

t)

u(O)

= u

0

(6.9.28)

If

Euler's method is applied, we have the numerical method

(6.9.29)

with

tn

= nh

and

vn

= u(tn). This is a well-known numerical method for the heat

equation, called the simple explicit method.

We

analyze the stability

of

(6.9.29)

and

some other methods for solving (6.9.28).

Equation

(6.9.28) is in the form of the model equation, (6.9.1),

and

therefore

we

need

the

eigenvalues

of

A to examine the stiffness

of

the system.

It

can be

shown

that

these eigenvalues are all real and are given

by

1

~j

~

m

-1

(6.9.30)

We

leave

the

proof

of

this to Problem 6 in Chapter 7. Directly examining this

formula,

(6.9.31)

A =

-stn

2

=-

-4

. (

(m-

1)'11) .

-4

m-1

82

2m

82

with

the

approximations valid for larger m.

It

can

be

seen

that

(6.9.28) is a stiff

system

if

8 is small.

Applying (6.9.31) and (6.8.49) to the analysis

of

stability

in

(6.9.29), we must

have

j=1,

...

,m-1

Using (6.9.30), this leads to the equivalent statement

4h (

j'IT

)

0 <

2

sin

2

- < 2

8

2m

1~j~m-1

This will

be

satisfied if

4hj8

2

~

2,

or

{6.9.32)

If

8 is

at

all small, say 8 = .01, then the time step h must

be

quite small to have

stability,

;

... ·l

········-···.!

STIFF DIFFERENTIAL EQUATIONS AND THE METHOD

OF

LINES

417

In

contrast

to

the restriction (6.9.32) with Euler's method, the

backward

Euler

method

has

no

such restriction since

it

is A-stable.

The

method

becomes

(6.9.33)

To

solve this linear problem for

Vn+l•

(6.9.34)

This

is a tridiagonal system

of

linear

equations

(see Section 8.3).

It

can

be

solved

very rapidly,

with

approximately 5m

arithmetic

operations

per

time step, exclud-

ing

the

cost

of

computing

the right side

in

(6.9.34).

The

cost

of

solving the

Euler

method

(6.9.29) is

almost

as great,

and

thus

the

solution

of

(6.9.34) is

not

especially time-consuming.

Example Solve the

partial

differential

equation

problem

(6.9.19)-(6.9.21) with

the functions

G,

d

0

,

d

1

,

and

f determined

from

the

known

solution

U = e

-.lr

sin

('IT

X)

O~x~1

t~O

(6.9.35)

Results for Euler's

method

(6.9.29)

are

given

in

Table

6.20,

and

results for

the

backward

Euler

method

(6.9.33)

are

given

in

Table

6.21.

For

Euler's method,

we

take m = 4, 8, 16,

and

to

maintain

stability, we

take

h = S

2

/2,

from

(6.9.32).

Note

this leads

to

the

respective time steps

of

h

~

.031,

.0078, .0020.

From

(6.9.26)

and

the

error

formula

for Euler's

method,

we

would

expect

the

error

to

be

proportional

to

8

2

,

since

h = 8

2

/2.

This

implies

that

the

Table6.20

The method

of

lines:

Euler's

method

Error

Error Error

t

m=4

Ratio

m=

8

Ratio

m =

16

1.0

3.89E-

2

4.09

9.52E-

3

4.02

2.37E-

3

2.0

3.19E-

2 4.09

7.79E-

3

4.02

1.94E-

3

3.0

2.61E-

2

4.09

6.38E-

3

4.01

1.59E-

3

4.0

2.14E-

2 4.10

5.22E-

3 4.02

1.30E-

3

5.0

1.75E-

2 4.09

4.28E-

3 4.04

1.06E-

3

Table

6.21

The

method

of

lines: backward

Euler's

method

Error

t

m=4

m = 8

m =

16

1.0

4.45E-

2

l.lOE-

2

2.86E-

3

2.0

3.65E-

2

9.01E-

3

2.34E-

3

3.0

2.99E-

2

7.37E-

3

1.92E-

3

4.0

2.45E-

2

6.04E-

3

1.57E-

3

5.0

2.00E-

2

4.94E-

3

1.29E-

3

·;·

'

)

'

I

.J

418

NUMERlcAL

METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS

error should decrease by a factor of 4 when m

is

doubled, and the results in

·Table 6.20 agree. In the table, the column Error denotes the maximum error at

the node points

(xj,

t), 0

5.j

5.

n, for the given value

oft.

For

the solution of (6.9.28) by the backward Euler method, there need no

longer

be

any connection between the space step 8 and the time step

h.

By

observing the error formula (6.9.26) for the method of lines and the truncation

error formula

(6.9.7) (use p =

1)

for the backward Euler method,

we

see that the

error in solving the problem

(6.9.19)-(6.9.21) will be proportional to h + 8

2

•

For

the unknown function U of (6.9.34), there

is

a slow variation with t. Thus for the

truncation error associated with the time integration, we should be able to use a

relatively large time step

h

as

compared to the space step

8,

in order to have the

two sources of error be relatively equal in size. In Table

6.21,

we

use h =

.1

and

m

==:'

4,

8, 16. Note that this time step

is

much larger than that used in Table 6.20

for Euler's method, and thus the backward Euler method

is

much more efficient

for this particular example.

For

more

discu~sion

of the method of lines, see Aiken (1985, pp. 124-148).

For

some method-of-lines codes to solve systems

of

nonlinear parabolic partial

differential equations, in one and

two

space variables, see Sincovec and Madsen

(1975)

and

Melgaard and Sincovec (1981).

6.10 Single-Step

and

Runge-Kutta Methods

Single-step methods for solving

y'

=

f(x,

y)

require only a knowledge of the

numerical solution

Yn

in order to compute the next value Yn+l· This has obvious

advantages over the p-step multistep methods that use several past values

{Jn,

...

,

Yn-p+d,

since then the additional initial values

{y

1

,.

•.

,

Yp-d

have to

be computed by some other numerical method.

The best known one-step methods are the

Runge-Kutta

methods. They are

fairly simple to program, and their truncation error can be controlled in a more

straightforward manner than for the multistep methods. For the fixed-order

multistep methods that were used more commonly in the past, the

Runge-

Kutta

methods were the usual tool

for· calculating the needed initial values for the

multistep method. The major disadvantage of the

Runge-Kutta

methods

is

that

they use many more evaluations of the derivative

f(x,

y)

to attain the same

accuracy, as compared with the multistep methods. Later

we

will mention some

results on comparisons of variable-order Adams codes and fixed-order

Runge-Kutta

codes.

The most simple one-step method

is

based on using the Taylor series. Assume

Y(

x)

is r + 1 times continuously differentiable, where

Y(

x)

is

the solution of the

initial value problem

y'

=

/(x,

y)

(6.10.1)

Expand Y(x

1

)

about x

0

using Taylor's theorem:

h'

hr+l

Y(x

1

)

=

Y(x

0

)

+

hY'(x

0

)

+ · · ·

+-

y<r>(x

0

)

+ ( )

y<r+~>a)

(6.10.2)

r! . r + 1 !

.....

~

.....

J

SINGLE-STEP AND

RUNGE-KUTIA

METHODS 419

for some x

0

~

Eo~

Xp

B:)'

dropping the remainder term,

we

have an approxima-

tion for

Y(x

1

),

provided

we

can calculate

Y"(x

0

);

••••

y<r>(x

0

).

Differentiate

Y'(x)

=

f(x,

Y(x)) to obtain

Y"(x)

=

jx(x,

Y(x))

+ fv(x,

Y(x))Y'(x),

Y" =

fx

+ /y/

and proceed similarly to obtain the higher order derivatives of

Y(

x

).

Example Consider the problem

y(O)

= 1

with the solution

Y(x)

=

1/(1

+

x).

Then

Y"

=

-2YY'

= 2Y

3

,

and (6.10.2)

with r = 2 yields

We drop the remainder to obtain an approximation of

Y(x

1

).

This can then be

used in the same manner to obtain an approximation for

Y(x

2

),

and so on. The

numerical method

is

n~O

(6.10.3)

Table 6.22 contains the errors in this numerical solution

at

a selected set of node

points. The grid sizes used are

h = .125 and h = .0625, and the ratio of the

resulting errors·

is

also given. Note that when h is halved, the ratio

is

almost

4.

This can be justified theoretically since the rate of convergence can be shown to

be

O(h

2

),

with a proof similar to that given in Theorem

6.3

or Theorem 6.9,

given later in this section.

The Taylor series method can

give

excellent results. But it

is

bothersome to use

because of the need to analytically differentiate

f(x,

y).

The derivatives can be

very difficult to calculate and. very time-consuming to evaluate, especially for

systems of equations. These differentiations can be carried out using a symbolic

Table 6.22

Example of the Taylor series method (6.10.3)

h =

.0625

h =

.125

X

Y~r(x)

Y(x)-

Y~r(x)

Y(x)-

Y~r(x)

Ratio

2.0

.333649

-3.2E-

4

-1.4E-

3

4.4

4.0

.200135

-1.4E-

4

-5.9E-

4

4.3

6.0

.142931

-7.4E-5

-3.2E-

4

4.3

8.0

.111157

-4.6E-

5

-2.0E-

4

4.3

10.0

.090941

-3.1E-

5

-1.4E-

4

4.3

420 NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS

manipulation language on a computer, and then it is easy to produce a Taylor

series method. However, the derivatives are still likely to be quite time-consuming

to

evaluate,

and

it

appears that methods based

on

evaluating

just

f(x,

y)

will

remain more efficient.

To

imitate the Taylor series method, while evaluating only

f(x,

y),

we

tum

to the

Runge-Kutta

formulas.

Runge-Kutta

methods The

Runge-Kutta

methods are closely related to the

Taylor

series expansion

of

Y(x)

in (6.10.2),

but

no differentiations

of

f are

necessary

in

the use of the methods.

For

notational convenience, we abbreviate

Runge-Kutta

to RK. All

RK

methods will be written in the form

Yn+l

=

Yn

+ hF(xn,

Yn>

h;

/)

n;;::;O

(6.10.4)

We

begin with examples of the function

F,

and

will later discuss hypotheses for

it. But

at

this point, it should

be

intuitive that we want

F(x,

Y(x),

h;!)

~

Y'(x)

=

f(x,.Y(x))

(6.10.5)

for all small values

of

h. Define the truncation error for (6.10.4) by

n;;:: 0 (6.10.6)

and

define Tn(Y) implicitly by

Rearranging (6.10.6),

we

obtain

n;;:: 0 (6.10.7)

In

Theorem 6.9, this will be compared with (6.10.4) to prove convergence

of

{

Yn}

toY.

Example 1. Consider the trapezoidal method, solved with

one

iteration using

Euler's method as the predictor:

n

;;::

0 ( 6.10.8)

In

the notation of (6.10.4),

F(x,

y,

h;

f)=

Hf(x,

y)

+

J(x

+

h,

Y +

hf(x,

y))]

As

can

be

seen in Figure 6.9, F is an

average

slope

of

Y(x)

on

[x, x +

h].

2.

The following method is also based

on

obtaining an average slope for the

solution

on

[xn,

xn+d:

(6.10.9)

Atkinson K. An Introduction to Numerical Analysis

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