DIFFERENCE METHODS FOR ORDINARY D I F F E R E N T I A L

EQUATIONS

WITH APPLICATIONS TO P A R A B O L I C

EQUATIONS

by E u s e b i u s Jacobus

Doedel

B . S c . , U n i v e r s i t y of B r i t i s h C o l u m b i a , M . S c , U n i v e r s i t y of B r i t i s h C o l u m b i a ,

1971 1973

A T h e s i s S u b m i t t e d i n P a r t i a l F u l f i l m e n t of the R e q u i r e m e n t s for the Degree of D o c t o r of P h i l o s o p h y i n the

Institute of

Applied Mathematics

We accept this t h e s i s as c o n f o r m i n g to the required

standard

T H E UNIVERSITY OF BRITISH C O L U M B I A January,

1976

In

presenting

an

advanced degree

the I

Library

further

for

this

shall

agree

thesis

in

at

University

the

make

that

it

freely

permission

s c h o l a r l y p u r p o s e s may

by

his

of

this

representatives. thesis

for

partial

financial

is

for

The

of

University

t

of

gain

^l/^/jC

British

2075 W e s b r o o k P l a c e V a n c o u v e r , Canada V6T 1W5

of

Columbia,

British

by

shall

Columbia

f

for

extensive the

understood

written permission.

Department

of

available

be g r a n t e d

It

fulfilment

requirements

reference copying of

Head o f

that

not

the

I

agree

and

be a l l o w e d

that

study.

this

thesis

my D e p a r t m e n t

copying or

for

or

publication

without

my

ABSTRACT

The of

first

c h a p t e r o f the

finite difference

nth

order

thesis

approximations

ordinary differential

collocation procedure

to

a method

of u n d e t e r m i n e d

of

these finite

minants states

certain normalization factor compact

difference

The

new

difference

conditions Kreiss based also

are

i s applied

expressed

given.

treated i n detail.

to i n v e s t i g a t e

The

basic

the d e t e r theorem

equations with only be

improved

by

of known,

as

Approximations stability

number

one

suitable well

as

to b o u n d a r y

t h e o r y of H.

the s t a b i l i t y of f i n i t e d i f f e r e n c e A

a

T h i s i s the c a s e f o r

Several examples

are

A

as

p r o v i d e d o n l y that

for difference

upon these a p p r o x i m a t i o n s .

O. schemes

of n u m e r i c a l e x a m p l e s

are

given.

In the f i r s t linear

the s e c o n d chapter

first

'included and is

be

o r d e r of c o n s i s t e n c y m a y

approximations also

can

a r e consistent,

c h o i c e of the c o l l o c a t i o n p o i n t s .

based

which i s equivalent

s m a l l dimension.

and

in linear

It i s s h o w n that the c o e f f i c i e n t s

d o e s not v a n i s h .

equations

c o l l o c a t i o n point.

with polynomials,

approximations

that t h e s e a p p r o x i m a t i o n s

construction

This construction is

coefficients.

of m a t r i c e s of r e l a t i v e l y

w i t h the

to b o u n d a r y v a l u e p r o b l e m s

equations.

upon a l o c a l

difference

is concerned

can

order

c h a p t e r i t i s s h o w n how be

extended

to i n i t i a l v a l u e p r o b l e m s

ordinary differential

the w e l l - k n o w n

the c o n s t r u c t i o n m e t h o d

stability

equations.

Specific

of

for systems examples

theory f o r these d i f f e r e n c e

of

are

equations

summarized.

It i s t h e n s h o w n how linear first

parabolic discretizing

partial

these d i f f e r e n c e methods

differential

i n space

by

a

e q u a t i o n s i n one

suitable method

from

may space

be

applied

variable

the f i r s t

to

after

chapter.

-111-

The

stability

investigated

of

using

e f f e c t of v a r i o u s The

such difference an

approximations

Numerical

for parabolic

eigenvalue-eigenvector

r e l a t i o n of t h i s a n a l y s i s

indicated.

schemes

to the

to the

examples

are

analysis.

boundary

stability also

equations i s In p a r t i c u l a r ,

conditions

theory

included.

of

J.

is M.

the

considered. V a r a h is

-ivTABLE

OF

CONTENTS

ABSTRACT TABLE

i

OF CONTENTS

iv

ACKNOWLEDGMENTS

v

Chapter I.

II.

BOUNDARY VALUE

PROBLEMS

1

1.

Introduction

1

2.

The C o n s t r u c t i o n of F i n i t e Difference A p p r o x i m a t i o n s

6

3.

Effect of the C h o i c e of C o l l o c a t i o n Points

19

4.

E x a m p l e s of F i n i t e Difference A p p r o x i m a t i o n s

21

5.

F i n i t e Difference A p p r o x i m a t i o n s to B o u n d a r y Conditions

29

6.

E x a m p l e s of A p p r o x i m a t i o n s to B o u n d a r y Conditions

34

7.

The S t a b i l i t y of F i n i t e Difference Schemes

37

8.

Numerical Examples

46

INITIAL

VALUE

PROBLEMS

1.

Introduction

2.

Initial Value Problems

54 54 in O r d i n a r y D i f f e r e n t i a l

Equations

55

3.

E x a m p l e s of F i n i t e Difference A p p r o x i m a t i o n s

58

4.

The S t a b i l i t y of F i n i t e Difference A p p r o x i m a t i o n s to Initial Value Problems

65

5.

Initial Boundary Value Problems

71

6.

The S t a b i l i t y of F i n i t e Difference A p p r o x i m a t i o n s to I n i t i a l B o u n d a r y Value P r o b l e m s

76

Numerical Examples

91

7.

BIBLIOGRAPHY

117

-v-

ACKNOWLEDGMENTS

The

author

and

guiding this

her

patience.

w i s h e s to thank

research. Finally,

the

J . M.

Varah for stimulating

H e i s a l s o t h a n k f u l to h i s w i f e

he i s i n d e b t e d

e x c e l l e n t t y p i n g of the m a n u s c r i p t proofreading

Professor

manuscript.

Adrienne for

to M r s . V i v i a n Davies

a n d to D r . F a u s t o

f o r her

Milinazzofor

-1-

Chapter

BOUNDARY

1.

VALUE

PROBLEMS

INTRODUCTION

In to

I

this

chapter

the l i n e a r

(1.1)

the c o n s t r u c t i o n

differential

Ly(x) = y

( n )

(x)

of f i n i t e d i f f e r e n c e

approximations

equation

+

n

j

1

a (x)y k

( k )

( x ) = f(x), 0 < x
J

> -

1

L e t h = max

= 1 } with J

h..

mesh-

It w i l l a l w a y s

be

j

h f o r some

c o n s t a n t K.

F o r a n y f u n c t i o n w(x)

w(x.). J

difference approximations

the p o i n t x = x^ a r e a s s u m e d

to the d i f f e r e n t i a l

e q u a t i o n (1.1) a t

to h a v e the f o r m

s.

J

(1. 2)

L u.^ h

I i=-r.

0 ,

Q

J

where

f. i s s o m e

approximation

to f..

J precisely

Thus

J - n + 1 such

equations.

J approximation i s said

The

m e s h S,

J T h e w i d t h co.. of e a c h

to be c o m p a c t .

f o r s o m e i n t e g e r co max b

local

J

truncation

error

there a r e

Im-

a p p r o x i m a t i o n i s d e f i n e d a s co. s r . + s. + 1.

co. < co j — max

for each

J

finite

difference

If co. = n + 1 t h e n the

J The assumption

that d o e s

i s made

that

n o t d e p e n d o n S, . h r

i s g i v e n b y T . = L ^ y ^ - f..

If t h e r e i s a

constant C

n of h such that |T.| 0, J s

independent

high as possible,

then the difference approximation

n

as s

i s said to be consistent

the order of c o n s i s t e n c y i s equal to n s.

with the differential equation and To

and

specify y(x) completely it i s n e c e s s a r y

conditions, which a r e a s s u m e d to take the

to i m p o s e boundary

form

n (0) k

(1.1a)

B (0)y(0)=

£

k

b ^ .(0)y (0) = b ( 0 ) ,

l

.(l)u

J + i

=b (l).

l < k < n - n

k

0

.

i=-r (l) k

H e r e b (0) k

The

and b ( l ) k

a r e approximations to b (0) k

and b ( l )

respectively.

k

definitions of truncation e r r o r and c o n s i s t e n c y are

similar

the previous definitions and the equations are compact if s (0) = ( 0 ) n

k

r (l)

= n (l).

k

k

order (1.2)

A finite difference

s_ if the approximate and (1.2a,b) and

k

to and

scheme is said to be (accurate) of

solution Uy

|y^-u. | < K ^ h

S

for

(0 < j

< J),

is uniquely defined by

some constant

that does not

depend on h. N u m e r o u s texts and papers a r e partly or entirely c o n c e r n e d with finite difference approximations to (1.1) books by Collatz (I960) and K e l l e r

and ( l . l a , b ) .

In p a r t i c u l a r ,

the

(1968) need to be mentioned h e r e . The

general theory of difference approximations to systems of boundary value p r o b l e m s is contained i n papers by G r i g o r i e f f (1970) and K r e i s s (1972). An

extensive survey of finite difference methods for boundary value

p r o b l e m s can be found i n a paper by K e l l e r

(1975).

V a r i o u s finite difference approximations of the f o r m (1.2) the l i t e r a t u r e .

(1.3)

To the second o r d e r differential

y"(x)

+ a^xjy'fx) + a°(x)y(x)

= f(x)

equation

,

the s i m p l e s t difference approximation is perhaps given by

(1.4)

appear

in

-4H e r e D. u. = (u h J J+1 uniform, i . e . , h. = h.

-u. ,) / 2h and the m e s h i s The

order of consistency i s equal to two,

as i s

the order of a c c u r a c y of a difference scheme based upon (1.4), p r o v i d e d that the approximations

to the boundary conditions are consistent.

F r e q u e n t l y i t i s d e s i r a b l e to use equations for which the o r d e r of consistency i s m o r e than two.

One

way

of obtaining such equations i s by

allowing the width of the approximation to be greater than n+1. example,

(1.5)

The

a fourth o r d e r finite difference approximation to (1.3) i s

L u. h

For

EE

width of this approximation i s equal to five.

If a difference a p p r o x i -

mation i s not compact then an entire difference scheme cannot be upon that equation alone. used.

here.

N e a r the boundaries other equations need to be

In (1.5) for example,

p e r i o d i c i t y i s assumed.

based

L^u. i s not defined i f j = 1 or j = J-1 unless

Hence other approximations need to be

E x a m p l e s of these are given by Shoosmith (1975).

used

In this chapter

such extra boundary conditions w i l l be d i s c u s s e d i n Sections 7 and

8.

It is not n e c e s s a r y to i n c r e a s e the width of the approximation to obtain higher o r d e r equations, for one can also r e q u i r e m o r e " l o c a l i n f o r m a t i o n " to be available.

In this case each difference equation

involves the values of the coefficient functions and the inhomogeneous t e r m at m o r e than one point.

F o r example,

to

(1.6)

y"(x) + a°(x)y.(x) = f(x)

a fourth order approximation

-5i s given by

(1.7)

D*u.

+ K ( f ) (a°u.) = K ( f ) f . , h

h

with

K (c)w. h

The

c o n s t r u c t i o n and

s

i|£. w.^

+ cw.

+ f w. 1

£

+ 1

a n a l y s i s of such higher o r d e r

mations i s also included i n this

.

compact a p p r o x i -

chapter.

U s u a l l y the statement that a c e r t a i n difference approximation given order Systematic

of consistency can be justified by

unless one

methods such as G a l e r k i n , Ritz and

e.g.,

V a r g a (1966),

Shampine (1972), Strang and

F i x (1973) and

Swartz (1973),

Schultz (1973)).

compare v a r i o u s compact equations,

with

piecewise

( F o r a d i s c u s s i o n of these,

C i a r l e t Schultz and V a r g a (1967),

deBoor and

appear i n f r e -

c o n s i d e r s finite element

collocation procedures

as finite difference methods.

but no

a

T a y l o r expanding.

methods of d e r i v i n g finite difference equations

quently i n the l i t e r a t u r e however,

polynomials

simply

has

see

R u s s e l l and

as w e l l as the books of B i r k h o f f and

Gulati: (1974)

attempt i s made to m e t h o d i c a l l y

construct such compact difference approximations. a p r o c e d u r e for d e r i v i n g difference approximations

Swartz (1974) d i s c u s s e s that r e s e m b l e s

the

p r o c e d u r e given i n this thesis, but no complete a n a l y s i s i s given.

The

m a i n purpose of the present

c l a s s of finite difference approximations their o r d e r

of consistency.

chapter

i s then to d e r i v e a l a r g e

to the p r o b l e m (1.1) and

T h i s i s done i n the next two

to give

sections.

Section 4 a number of examples i s given, including known as w e l l as

In new

-6difference equations. approximations

I n S e c t i o n 5 the c o n s t r u c t i o n of f i n i t e

to the b o u n d a r y c o n d i t i o n s i s c o n s i d e r e d .

t h e s e a p p e a r i n S e c t i o n 6.

g i v e n i n this chapter.

THE To

C O N S T R U C T I O N OF

u p o n the

i n S e c t i o n 8 the

are given.

FINITE D I F F E R E N C E

APPROXIMATIONS

c o n s t r u c t a f i n i t e d i f f e r e n c e a p p r o x i m a t i o n o f the f o r m

the d i f f e r e n t i a l e q u a t i o n (1.1) one i n d e x j l e t z , < z_ < • • • < z 1 2 m J

satisfying

Finally,

| z . . - z , | >—-

h,

be

may r

i # k\

proceed

as f o l l o w s .

For

(1.2) to given

p o i n t s i n the i n t e r v a l Tx. ,x.. 1, j - r j+s L

(To

J

s i m p l i f y the n o t a t i o n the

subscript

K

1

j d e n o t i n g d e p e n d e n c e o n the p a r t i c u l a r d i f f e r e n c e a p p r o x i m a t i o n w i l l omitted here,

so that,

e.g.,

exceeding

P_j d e n o t e the s p a c e d=r

+ s + m - l .

be

z. . b e c o m e s z.). J> 1

Let

of

(1972) i s

of d i f f e r e n c e s c h e m e s b a s e d

r e s u l t s of s o m e n u m e r i c a l e x p e r i m e n t s

2.

Examples

I n S e c t i o n 7 the t h e o r y of K r e i s s

a p p l i e d to i n v e s t i g a t e t h e s t a b i l i t y difference approximations

difference

1

of a l l p o l y n o m i a l s w i t h d e g r e e

Assume

t h a t f o r s o m e i n t e g e r i,

not (-r < l< s ) ,

SL

a p o l y n o m i a l p (x) e P ^ (2. 1) P (x )_ =

is determined u j+i

uniquely by

i +

r < i < s , (2.2)

Lp (z.)

e q u a t i o n s (2.2) w i l l b e

the e q u a t i o n (2.3)

equations

t

f(z.) 1 < i < m

(The

the r + s + m

referred

to as

. "collocation"

equations.)

Then

relates

the q u a n t i t i e s j ^ »

(-r i=l e

j + i

L

z

of p(x) i t f o l l o w s

°-

=

that

^

which

truncation

LQP(X)=1,

+

r

is a contradiction.

So

i f E + 0 then n + 1 or

_ d. / d. = 1 + 0(h), the s a m e l ' l

Since

h is small

Pj j_ = 0 a n d

d^

i=-

_ i=l d. a r e n o n z e r o . l

more d^,

m

ip

d

is true

f o r the

enough.

error

T. i s d e f i n e d

as

J m r

where

j -

V j

s

~ I i ( i> i=l e

conditions

= I Vj+i i = -r

z

y(x) i s the s o l u t i o n

boundary

- I i i=l e

£

(

z

i

•

)

to the d i f f e r e n t i a l e q u a t i o n

(1. l a , b).

d T.

f

m

y

Taylor

t-i i=-r

for

some

the

quantity

i

xd+1

y

d.(x.,.-x.) l j+i j

between

square brackets

follows such

w i t h the u n d e r l y i n g that t h e r e

equals

zero

By

f o r e a c h k.

d. a n d e. a n d the b a s i s 1

functions

construction F r o m the w (x), 1

1

v

assumption

a r e constants

,

v

s. b e t w e e n x. a n d x., . a n d t- b e t w e e n x. a n d z. . i J J+i i J i

d e f i n i t i o n s o f the c o e f f i c i e n t s together

1

v _ i _ r / ^+1 - > e. — , , , ... L(z.-x.) Li I (d+1)! I y i=l

i

[S }

— . ,, ,, , — (d+1)!

C^,

that ^

< tu < h,

and

(1 n.

Assume

m-1 r+s ~| |~ ( x - t , ) 1 T ( x - x .

k=i

k=o

K

J

"

that E + 0

, ) and a s s u m e r + K

s*** -f s T in. t h a t LpW

( z ^ = 0,(1 < i < m ) . T h e n the o r d e r

difference approximation If m

(2.4) i s a t l e a s t e q u a l t o r + s.+ m

- n + 1.

- 1 a n d the d i f f e r e n c e a p p r o x i m a t i o n i s c o m p a c t , i . e . , r + s = n,

~ r+s+m n then w (x) = ~| ["(x-x. k=0 J

, ). _

r

I n t h i s c a s e t h e r e i s only one

+

collocation point f o r which a higher m

of c o n s i s t e n c y of the f i n i t e

order

of c o n s i s t e n c y i s o b t a i n e d .

= 1 and r + s > n then t h e r e a r e r + s + 1 - n p o s s i b l e c h o i c e s

point.

( A s s u m i n g E ± 0.)

-parameter

F o r the c a s e m

f a m i l y o f p o i n t s z^, (1 < i < m ) ,

of c o n s i s t e n c y i s o b t a i n e d . i n t h e d e f i n i t i o n of w these 4.

of t h i s

> 1, r + s = n t h e r e i s a f o r w h i c h the i m p r o v e d

(m-1) order

The p a r a m e t e r s i n q u e s t i o n a r e the points t ^ (x) i n t h e o r e m (3. 1).

Particular

e x a m p l e s of

s p e c i a l c o l l o c a t i o n p o i n t s w i l l be i n c l u d e d i n the n e x t s e c t i o n .

EXAMPLES

E x a m p l e 4. 1. two

If

OF

FINITE DIFFERENCE

L e t n = 1, r = 0 a n d s = 1;

point finite difference approximation

(4.1)

Take m

APPROXIMATIONS

i . e . , we

w a n t to c o n s t r u c t a

to the e q u a t i o n

L y ( x ) = y'(x) + a°(x)y(x) = f(x)

= 1 and z

i

= x^ ± +

= _-(x.. + X j

+

1

).

T h e n the b a s i s f u n c t i o n s w ( x ) a s 1

d e f i n e d b y (2.13) a n d (2.14) a r e w°(x) = - ( x - x . .) / h., a n d w ( x ) = ( x - x . ) / h . , J+1 J J' J T h e d i f f e r e n c e a p p r o x i m a t i o n h a s the f o r m X

•22where

the c o e f f i c i e n t s

d ,

a n d e^ c a n b e

Q

obtained

from

(2. 10) a n d

(2. 11), v i z .

A

0

T

d

Q

= Lw

1

= Lw

,

"I

V

( )

=

Z l

,

1

+

2

0 a. ± , +

and

1/ \ 1 , 1 0 ( ) = h. — + 2 a, , i 1+1

T

A

d

ZjL

J

The

difference

approximation

(4.2)

{u

c a n then be w r i t t e n as

j r j j j+i V j+i u

)/h

+ia

(

u

•

) =

+

which is recognized Keller least

(1969). equal

as the c e n t e r e d

According

to one.

That

follows

from

w (x) =

(x-Xj)(x-x.. ).

2

to t h e o r e m the o r d e r

(3.1),

because

equation,

or Box

(2.1) the o r d e r

scheme

of

of c o n s i s t e n c y i s a t

of c o n s i s t e n c y a c t u a l l y z^ i s the c r i t i c a l

equals

two

p o i n t of

+ 1

An

equally simple

d e r i v e d when the

theorem

Euler

form

taking m

difference

equation,

= 2, w i t h z, = x. a n d z_ = x. , 1 j 2 j+1

d_.u. + d . u. , . = e . f. + e f . ... 0 j 1 J+1 1 J 2 j+1

1 E

2

This

equation

' , :

0

With w

the t r a p e z o i d a l m e t h o d , i s

2 (x) = ( x - x _ . ) ( x - X j ^ ) / h. the c o e f f i c i e n t s a r e

Lw

0

(z^)

Lw

2 (z^) d

Lw°(z ) 2

Lw (z ) 2

2

l

=

E

Lw

1 (z^)

Lw

2 (z^)

Lw

1 (z )

Lw

2 (z )

2

2

has

-23-

e

_

1 2 = __" L

w

L^y(x) = y

(_.) z

i

a

(x).

n

d

e

2

-1 ~ET "

=

L i W

The equation

2 ^ i)' z

w

i

t

n

= Q

E

L

^ l^

w

2

z

~ 0 L

2 ^ 2^

W

Z

A N C I

thus obtained i s

(u., .-u.) / h. + |(a°u.+a°, ,u. .) = _-(f'.+f. ,). J+ 1 J J J J J+1 J+1 J j+l' 2 V

B y theorem (2. 1) the o r d e r

of consistency of this difference equation i s

equal to two.

Next, consider the general two point difference approximation m

= 2.

That such a difference approximation

f r o m theorem (2.4). i s equal to two.

i s always consistent follows

T h e o r e m ( 2 . 1 ) i m p l i e s that the o r d e r

From

with

of consistency

theorem (3. 1) i t follows that the o r d e r

becomes ~ 3

three i f the collocation points coincide with the c r i t i c a l points of w (x-t, )(x-x.)(x-x. , ,). 1 J J+1

(x) =

L e t t, = x. + r i h . , z. = x. + £,h. and z_=x. + £_h.. 1 J J 1 J 1 J 2 j 2 j

T h e n the special values

of ^

i

and £

2

that give a third o r d e r

the roots of p (£) = 0, where p(£) = £(£-l)(£-n).

formula are

These roots a r e

rz—

f

t

_

5

-

V+ 1 ± v ••

••

f) -T7 + 1

3

If n = •§• then the roots a r e £ = \ ± ^

. (0))

l0, 0 —

k„ < j < r - l ,

y

X

0 -

J

n

-

i=l

0 H

of unknowns.

to m a t c h the n u m b e r o f

equation and given b y

i=0

(7. 5)

i n order

that t h e s e e x t r a

s.

(7.4)

i . e . , i f r + s > n then r + s - n e x t r a

j

m

1 j,i J i i=-r. d

U

+

J

i{z

e

V

55

= I ),i j,i i=l

]

J"B+I "" > T

T

T n

k ' " ' J-n+k ' T

T

Q

T

T

k
«

1 < k
0

i=0 H e r e N = r + s - n a n d t h e c o e f f i c i e n t s a^ a r e t h e s a m e characteristic polynomial

c(t) i n ( 7 . 6 ) .

a s those of the

I f (7.10) s u b j e c t

to (7.11) a d m i t s

a n o n t r i v i a l s o l u t i o n then i t i s e a s i l y s e e n that the p o l y n o m i a l s ((r+s-n) / 2

=k

dependent.

H e n c e we

Theorem

(7. 5).

Q

< j< r-l),

and

< j < m a x N.-N), a r e l i n e a r l y

h a v e s h o w n the f o l l o w i n g

L e t the homogeneous p r o b l e m c o r r e s p o n d i n g

( l . l a , b ) only have the t r i v i a l be consistent.

p^(t), (o

c^(t),

Assume

to (1.1) a n d

s o l u t i o n a n d l e t the d i f f e r e n c e s c h e m e (7.6)

t h a t c(t) i s s y m m e t r i c .

Also

suppose that the

d e g r e e o f c(t) i s e v e n a n d that c(t) i s s t r i c t l y d i a g o n a l l y d o m i n a n t positive coefficients. boundary conditions

L e t the c h a r a c t e r i s t i c p o l y n o m i a l s be r e l a t e d a s i n E q u a t i o n ( 7 . 1 4 ) .

c..(t) a n d P j ( t ) d e f i n e d

s t a n t K i n d e p e n d e n t of h s u c h that

Example

7.6.

of the e x t r a

I f the p o l y n o m i a l s

b y (7. 8) a n d (7. 15) r e s p e c t i v e l y a r e l i n e a r l y

pendent then the d i f f e r e n c e s c h e m e i s s t a b l e .

with

Hence there

inde-

exists a con-

| | e ^ | | < K | | T ^ | |.

2 L e t c(t) = a ^ t + ( l - 2 a Q ) t + a ^ w i t h

d e g r e e of t h e c h a r a c t e r i s t i c p o l y n o m i a l

- co < a ^
(t) = " ( t - t ) / A t , (oV) = ( t - t _ ! ) / 0

0 < i < 2 — —

X

n

n

A

t

a

n

d

n

^ (t) = (t-t _ )(t-t )/A 2

n

1

n

i t i s f o u n d that

A-t^W

" " A

fc

B

^

•

2 t

,

•60-

L

From

AT

2,1

A

Li

(t -i) n

A(t ) n v

(3.1) i t f o l l o w s t h a t E , A ( t

,) = E _ A ( t ).' T a k e E , =1. T h e n E. i s 2 n 2 1 e q u a t i o n E , A ( t ,) = A ( t ). C l e a r l y , i t i s m o r e 1 n-1 n n-1

1

the

AT

2,2

s o l u t i o n of t h e m a t r i x

v

1

d i f f i c u l t n o w to g i v e e x p l i c i t r e p r e s e n t a t i o n s was

the c a s e f o r s c a l a r equations

a l w a y s obtain the c o e f f i c i e n t s

of the c o e f f i c i e n t s E^ than

i n C h a p t e r I.

Of course,

one c a n

numerically.

I n t h e s p e c i a l c a s e w h e r e A ( t ^) a n d A ( t ) c o m m u t e , n

n

i ti s possible

to s t a t e e x p l i c i t l y w h a t the c o e f f i c i e n t s o f t h e a b o v e f i n i t e d i f f e r e n c e equation

are.

In particular,

2 I i t i s f o u n d that E ^

*

h

D

T h u s the d i f f e r e n c e equation , n n-1. ) A (u -u At

which i s recognized consistency

Example

i f A i s a constant -i

=

AT

t h e n w i t h E_> =

fB(t ,) a n d D = - f - A - | B ( t ). * n-1 0 At ^ n n

A

m a y be w r i t t e n as

^ ( B ^ u ^ B ^ ^ u

1

1

-

1

)

= |(f(t ) n

m^j)

+

as the u s u a l t r a p e z o i d a l m e t h o d .

of t h i s e q u a t i o n i s a l s o e q u a l

3.2.

matrix,

W i t h p = 2, |_L = 1 a n d

T h e o r d e r of

to t w o .

= t

^ + € A t the d i f f e r e n c e approxi-

m a t i o n to (2.1) b e c o m e s

|(2M)u -2(e-Du - + l ( 2 e - 3 ) u n

A

^

At

n

1

n

2

B(^ ) 1

n-1

ie(e-i)u-e(e-2)u- + n

x

+ l(£-l)(£-2)

( T h e m e s h i s u n i f o r m i n t h i s e x a m p l e . ) T h e o r d e r of c o n s i s t e n c y o f t h i s approximation order

i s equal

to t w o .

of c o n s i s t e n c y i n c r e a s e s

From

T h e o r e m (2.1) i t f o l l o w s t h a t t h e

to t h r e e i f o n e t a k e s

£ ^- \^3 =

±

~ 1-58.

•61If £ = 2 t h e n the a b o v e f o r m u l a

v

reduces

to

3 n _ n-1 , i n-2 — u - 2u + 2

B(t

U

n' A t

)u

= f(t ) , n

n

n

w h i c h i s one o f G e a r ' s b a c k w a r d d i f f e r e n t i a t i o n f o r m u l a s . p.

(See G e a r

(1971),

217.)

Example

3. 3.

Consider

the c o n s t a n t c o e f f i c i e n t

Aw

(3. 2)

(t) - B w ( t ) = f(t) ,

w(0)

(3. 2a)

where A

= w

i s nonsingular and where A

example when A

= I.

problem

0 < t < T

,

,

and B

commute.

L e t p. = 2 a n d p = 1.

This occurs for

The construction procedure

of

S e c t i o n 2 a g a i n l e a d s to t h e m a t r i x e q u a t i o n (3.1) g i v e n i n E x a m p l e (3.1), now

f o r g e n e r a l c h o i c e of t h e p o i n t s

Let E B

commute,

E

= -L co (£ ) A 2

1

0

i t i s found

= L co (q)A Z

2

2

0

+ ^{t,^

a n d t,^.

B.

Then,

u s i n g the f a c t t h a t A

that

- co (?, )B,. 2

1

D _ j _ = E {L co°(c: )A - u°(t )B}+ E {L co°(t: )A-co (C )B} , 0

1

and

D

Q

-

0

i

E ^ L Q W ^ C ^ A

>1

- W ^ ^ B }

2

+

0

E.,{ L.

2

Qbi

1

[t )A-( \t , )B}, >2

w h e r e LQto(t) = co'(t).

In p a r t i c u l a r

i f one u s e s the G a u s s i a n

2

points

li

2

and

-62-

= (t

n - 1

+ t ) / 2 - AtN/3/6 n

and

C = 2

(^.i+tj/

2

+ At^3/6

t h e n the d i f f e r e n c e e q u a t i o n b e c o m e s

{^A

(3. 3)

2

-±AB

- §B }u "

= {|A

+

2

n

^I-B}

1

+ {^A -lAB 2

f( ) 4- { ± A Cl

|

+

- ^ - B

2 B

}

} u

n

f( ) &2

T h e o r d e r of c o n s i s t e n c y o f t h i s e q u a t i o n i s e q u a l to f o u r . f(t) = 0 t h e n the f i n i t e d i f f e r e n c e a p p r o x i m a t i o n

(3. 3) c a n a l s o be

f r o m a particular Pade approximation

to t h e s o l u t i o n w(t) =

(3.2),

o r i g i n a l l y due to P a d e

(3.2a).

These approximations,

been extensively studied i n connection with their (See

B

If obtained w^

of

(1892) h a v e

stability properties.

e. g. , V a r g a (1961).)

Example

3.4.

The preceding

finite difference approximation

example suggests

that the g e n e r a l

p-step

to (3. 2) w i t h p. c o l l o c a t i o n p o i n t s h a s the

form

i.=-p

where

v

=

1=1

0

v=0

a n d (3^ a r e c o n s t a n t s that d e p e n d o n the c h o i c e of the p o i n t s t,^,

(1 < SL < p.). JL

T o v e r i f y t h i s a n d to s h o w h o w found, (3. 5)

c o n s i d e r the s c a l a r

the c o e f f i c i e n t s a

equation

L w ( t ) = w'(t) - X w(t) = f(t) ,

v

JL

a n d |3

may

be

-63The

finite difference approximation

to (3.5) i s g i v e n b y

n+4

(3.6)

l=-p

1=1

where

(3.7)

e.

i-il

P+ // 1

T

|JL+4+l

X

(£ )

Lw

and

6^ = J

(3.8)

L

H + i w

(C )e k

,

k

-P< k< 0

k=l The n o r m a l i z i n g factor

S i s g i v e n b y (2.4).

Considering

(3.5),

(3.7) a n d

(3. 8) i t i s o b s e r v e d that the c o e f f i c i e n t s h a v e the f o r m u.-l

[X

v=0

for c e r t a i n constants In order i n equation

a

v=0

4

and 6

4

to s h o w that t h e s e c o n s t a n t s

(3.4),

note that b y d e f i n i t i o n

a r e t h e s a m e a s the c o n s t a n t s

•64-

Y

e La |

k )

) = 0 ,

(^

p+l < k < p

+ HL-1

4=1 i.e., .-1 = o ,

4=1 for

v=0

a l l X and f o r p + l < k < p

+ [ i - l

So q^(X.)= 0 a n d t h e r e f o r e ° x ( C ) = 0 f o r a l l m a t r i c e s

C.

k

In

p a r t i c u l a r w i t h C = A B one obtains X

H

u-1 0 .

L Q C O ^ ^ J I - C O ^ C ^ A ^ B

4=1

v=0

V

U p o n m u l t i p l i c a t i o n b y A ^ a n d u s i n g the a s s u m p t i o n that A a n d B one

commute

obtains p.

u-1 L co (^)A-oo (^)B k

k

0

Q

4=1

,

p + l < k < p +[i-l .

v=0

Thus i f one defines .-1 v=0

t h e n t h e l a s t u.-2 r e l a t i o n s i n t h e m a t r i x Further,

the f i r s t

equation

(2.3) a r e s a t i s f i e d .

p + l r e l a t i o n s o f (2.3) a r e a u t o m a t i c a l l y s a t i s f i e d b y

letting P-

D -

S

I

E

L cop ^(c: )A(c: )+

v

0

V

v

p+i w

(c )B(c: ) V

v

T h e r e f o r e i t has

b e e n shown,

that i n c a s e

p r o b l e m (3.2), w i t h c o m m u t i n g A coefficient matrices

4.

THE

and

S T A B I L I T Y OF

INITIAL V A L U E

and

B,

of the c o n s t a n t c o e f f i c i e n t one

c a n e x p l i c i t l y g i v e the

of the d i f f e r e n c e a p p r o x i m a t i o n ( 2 . 2 ) .

FINITE D I F F E R E N C E APPROXIMATIONS

TO

PROBLEMS

In this section some w e l l - k n o w n concepts of finite d i f f e r e n c e a p p r o x i m a t i o n s

i n the s t a b i l i t y

analysis

to s y s t e m s o f f i r s t o r d e r o r d i n a r y

d i f f e r e n t i a l equations w i l l be

summarized.

d i f f e r e n c e a p p r o x i m a t i o n one

o n l y has

To

d e f i n e s t a b i l i t y of a f i n i t e

to c o n s i d e r the f o r m

this

difference

i

e q u a t i o n t a k e s w h e n a p p l i e d to the e q u a t i o n w (t) = 0 . scalar function.

T h u s f o r the t y p e of a p p r o x i m a t i o n s

H e r e w(t) i s a discussed i n Section

2 the r e l e v a n t d i f f e r e n c e e q u a t i o n i s 0 Y,

(4.1)

i =

&i u

n

+

i

= 0 ,

p < n < N ,

-p

where

5

i

(-1Y =

v^v v

p+1

v •• V^X

w i t h L Q C o ( t ) = to'(t).

T h e m e s h i s t a k e n to be u n i f o r m i n t h i s ( At)

n

= A t = 1/N,

(1 < n < N).

The

11

satisfies

so t h a t

a p p r o x i m a t i o n i s s a i d to b e

p r o v i d e d that f o r a l l g i v e n i n i t i a l data u , (4.1)

section,

(0 < n
1.

for which

are

o n l y i f £ > 1.

2

are

6^ = ( 4 - 4 £ ) / 2 A t

e q u a t i o n (4.2)

| r j | < 1 i f and

s t a b l e i f and

of G e a r ,

the c o e f f i c i e n t s

b

n

= 0 ,

r o o t s of the c h a r a c t e r i s t i c

let p = 2 and

6

n

S Q = ( 2 £ - 3 ) / 2 At,

g i v e n by

o b t a i n e d i f % - 1 + -|- N/3~

Now In

n-2

= (2£-3)/(2£-l).

2

the a p p r o x i m a t i o n the

u

the

- 8

2

1

-3(e +e ) + 6 e e 1

2

x

2

+ 2.

and

6

The

characteristic

equation

°"Q( " )

*i § For If £

the a p p r o x i m a t i o n

2

-9(

6

= 0 has

2

.

roots n ^ = 1 ,

e + e) + 6 x

strongly

>

2

and

2

+ 2

stable i t i s n e c e s s a r y

p r o v i d e d that ^ 5

and

+14

2

x

to be

that

-

Y

-3(e + e ) + 6 ^ e

2

= 2 t h e n t h i s i s the c a s e

then i t i s n e c e s s a r y

b

= -

Q

if




that

2 "J*

^

|i? |< !• 2

^2

=

14 needs — < ^

^ < 2.

-68Another

w e l l - k n o w n c o n c e p t i s that of s t a b i l i t y

difference approximation. eigenvalues

with large

approximations stable,

accuracy. do

A t to be

this

eigenvalues

from

requirement

then a l a r g e

this

may

shortcoming.

under

w h e r e w(t) i s a

consideration

scalar

function.

For

not be

necessary

Aw'(t) = Bw(t),

where w

p r o v i d e d that A

difference approximations

^B

investigating

analysis

a complete

to the s c a l a r

of d i f f e r e n c e

approximations the

A

and

that

difference

f o r the

are

constant

set of e i g e n v e c t o r s .

e q u a t i o n w'(t)

of

(t) - \w(t) = 0,

valid B

be

e f f e c t o f the

j_ ^w(t)= w

remains

to

for reasons

to a p p l y the f i n i t e

i s a v e c t o r and

has

those

to the p r o b l e m

The

class

has

s m a l l f o r the a p p r o x i m a t i o n

u p o n the s t a b i l i t y i t i s c u s t o m a r y

approximations

matrices,

very

finite

(2.1) the m a t r i x A

T h u s i t i s of i n t e r e s t to c h a r a c t e r i z e

not s u f f e r

system

negative r e a l part,

requires

even though

If i n the s y s t e m

r e g i o n of a

= X.w(t) h a v e

Finite the

form

0

Y

where

the c o e f f i c i e n t s

6

6^

n+l

u ' * n

L

are

given

= 0 ,

P < n < N

,

by

(-l) "1

L

) A.

F o r stability i t i s again n e c e s s a r y

that the r o o t s r j . o f the

U.

characteristic

0

e q u a t i o n q(rj) =

At Y

8.r]

P

= 0 satisfy

the r o o t c o n d i t i o n

i=-P From

(4.4) i t f o l l o w s

that q(r?) h a s

the

form

(4.3),

(4.3a).

-69-

'

(4.5)

q(n)

=

3 =o L

m

2

1

I

c

A

'

S " J

m

I

h

A

c

h

A

h h

=

c

0

-

Hence each component cl^ j

equation

v

h

A

is the diagonal m a t r i x of eigenvalues.

(6.5)

II

h

p.

h h L v

v

=0

0

satisfies

Substitution into (6.4),

p.

4=-p

of

c^.

gives

AT

where

2

3 one c a n w r i t e u ^ =

I % 4=-p v=0

h

Then

| } 2

"',

P < n< N = T / A t . +1 ^ c5* P ) and let h,j + 1

Define the p X p m a t r i x A ( z ) by

-82-

• p-l P (z) p

"Pp-2 (z)

( z )

( Z )

P p

p

P (z) p

1

0

0

0

1

0

A(z) =

1

v , ^ a m z' (z) = 7. v=0 v= one step difference equation

0 1, a n d i n a d d i t i o n ,

n e a r z = 1 m u s t be s a t i s f i e d .

s h o w n b y V a r a h h o w one c a n c h a r a c t e r i z e e i g e n v a l u e s E s s e n t i a l l y the same c h a r a c t e r i z a t i o n i s obtained where the eigenvalue values and

problem L ^ v ^ = ^ ^ n

z w i t h R e ( z ) > 0.

(6.4).

v n

n

z with

It i s a l s o | z | > 1.

when analyzing the case

°^ t h i s s e c t i o n h a s

T h i s w i l l be c o n s i d e r e d

again

eigensystem analysis.

p u r p o s e of the next e x a m p l e i s t o i n d i c a t e w h y i t i s d e s i r a b l e

at l e a s t the e n t i r e negative

E x a m p l e 6.3. each eigenvalue

(6.6)

eigen-

i n E x a m p l e (6.3)

to r e s t r i c t t o d i s c r e t i z a t i o n s i n t i m e w h i c h have a s t a b i l i t y r e g i o n includes

some

It i s n o t c l e a r h o w t h e c o n d i t i o n o n t h e s p e c t r u m o f Q r e l a t e s

to the c u r r e n t

The

Further

°f t h e d i f f e r e n c e s c h e m e Q g j ^ =

h a v e no e i g e n v a l u e s

c o n d i t i o n on t h e s p e c t r u m of Q

sufficient

by V a r a h , a r e f i r s t that

( T h u s t h e w o r k o f W i d l u n d c a n be u s e d t o c h e c k t h i s c o n d i t i o n . ) the

c a n be

T h e s e effects have a l s o been investigated by

V a r a h f o r the s c a l a r d i f f u s i o n equation, systems,

A n advantage i s

real axis.

It f o l l o w s f r o m t h e g e n e r a l of the eigenvalue

Ly

= y"

theory

of K r e i s s

problem

= Xy

that

,

y(0) = y'(l) = 0

,

(1972), t h a t

-86is approximated by an eigenvalue of the discrete eigenvalue problem ( 6

-

- W h

Vh

7 )

.

as h —- 0, provided that L^u-^ = K^f^ defines a consistent and stable finite difference approximation to L,y=f. values of (6.6) are given by \

(See Section (I. 7) J

= - L"

1 '

5

fc

k

- *» ^

Since the eigenf o l l o w s

are eigenvalues X.^ of (6.7) such that \^ -*• -co as h — 0.

t h

a t there

Hence the

restriction to discretizations i n time with an unbounded region T 0

Theorem (6.2) is desirable.

If T

a

a

as i n

were bounded then a restriction on the

size of the time step would have to be imposed.

F o r example, consider

the case where L, u. = -— (u. ,-2u.+u. ,,) and K. u. = u. , (1 < i < J-1). J 2 j-1 j j+r h j J - 1

The

J

h

discrete boundary conditions are

UQ

=

0 and

(Uj-Uj_^)/h

= 0.

Then the

eigenvalues of (6.7) are easily found to be given by

Thus i f the region T

a

includes only a part (-(3,0) of the negative r e a l axis Q

then At must satisfy At < -*r h ' — 4 The

2 in order to guarantee that At L . e T . h, x a to

stability of particular finite difference approximations w i l l be

investigated i n the next two examples as w e l l as i n numerical examples i n the next section.

Since one can always apply a strongly

A(a)~stable

method to the system (6.3), it is necessary only to study the effect of various spatial difference approximations on the eigenvalues of the eigenvalue problem (6.7).

Specifically, it i s important to check that this dis-

crete eigenvalue problem does not allow eigenvalues X.^ that tend to 00 in the right half plane and that the condition number -("v"^) is bounded. K

-87-

Example 6 . 4 .

-fAlthough stability and

approximation

imply that the eigenvalues of ( 6 . 6 ) are approximated

consistency of the spatial difference

those of the discrete eigenvalue p r o b l e m ( 6 . 7 ) ,

by

it need not be true that

a l l eigenvalues of ( 6 . 7 ) converge to eigenvalues of

(6.6).

C o n s i d e r for example the case where L ^ and

are as i n the

previous example, but where the boundary condition at x = 1 i s a p p r o x i mated by

=

(2-c) h

i s consistent with y (1) = 0 ,

T h i s approximation

0

p r o v i d e d that c + 2 .

The

4

order of consistency i s equal to one.

If c =-j

(See also E x a m p l e ( 1 . 8.4').)

is equal to two.

F o r this d i s c r e t i z a t i o n the

i n ( 6 . 3 ) i s the identity and the operator L ^ can be w r i t t e n as

operator = Pjjl^,

where " -2

P

then the order of consistency

1

and

h = 2-c

1

-2 1

i

~

Let v ^ = JL ~

v^.

T h e n the eigenvalue p r o b l e m T ^ v ^ = X ^ K ^ v ^ becomes

JL ~

~

P £ -J-hPj^h

=

^h h' v

w

^

c

^

with r e a l eigenvalues.

admits The

a complete

matrices

set of o r t h o n o r m a l

and

containing the eigenvectors

~ v ^ ^ and v ^ ^ r e s p e c t i v e l y by column are then r e l a t e d by JL

^ )2I U h"2 2

K ( V

h

)

"

V

h

2

V

h

2 -

" 'h"2

11

"h

"2

eigenvectors

II 7
1 it follows that

there is a root z ^ of (6.9) s u c h that z ^ — c - i as h -*• 0.

The

-89c or responding

eigenvalue

is obtained

from

(6.8)

and

asymptotically

given

by

h z

h

Thus

X^

—>

equivalent

+

h —» 0.

to the t e c h n i q u e

cretization of the

oo as

(see

h (c-l)

2

Varah

2

h

This

eigenvalue

of V a r a h w h e n a p p l i e d to the

(1970), p.

33, a n d

discrete boundary conditions

of the a m p l i f i c a t i o n m a t r i x , eigenvalues

of K ^ L ^ .

computation is essentially

(1971).)

There

i s r e l a t e d d i r e c t l y to the

whereas

If the

Varah

current

here

these

effects are

dis-

the

effect

eigenvalues

r e l a t e d to

d i s c r e t i z a t i o n i s c o m p l e t e d by

applying

the a

i

multistep

m e t h o d to the

system

K^u^

= J-^^ >

n

m e t h o d i s s t r o n g l y A(a)-stable,

but

in its stability

resulting scheme

c < 2.

The

rule.

(See

simplest

c h e c k e d by

example

half - plane.

applying

f o r w h i c h the

e n o u g h h,

co h a v e t h i s Example

6.5. (6.7)

2.

this

That this

i f this

multistep

p o s i t i v e r e a l a x i s i s not

scheme

m e t h o d i s the stability

is unstable

g ^ = v^,

trapezoidal

region

v^

There

d i f f e r e n c e s c h e m e w i l l be a l l m e t h o d s that a r e

only i f

i s the

for c > 2 can

where

positive eigenvalue

In f a c t ,

contained

i s stable i f and

m e t h o d the

i n i t i a l data

complete

even if c >

d

of s u c h a m u l t i s t e p For

i t to the

a s s o c i a t e d w i t h the

in time

version

t h e n the

E x a m p l e (1.4. 1).)

entire negative

vector

region,

the

a

i s the

be

eigen-

are discretizations stable f o r

small

s t r i c t l y s t a b l e at

property. Again consider

the

eigenvalue

problem

with

-4-

L, v. = n J h

2

(v. , - 2v. 3

3

+

v..,) 3

(6.6)

and

its discrete

-90-

and K, v . = -r— v . - + -pr v . + -pr v . , , h j 12 j-1 12 j 12 j+1

Thus this eigenvalue p r o b l e m a r i s e s i n the s t a b i l i t y a n a l y s i s when u s i n g the fourth o r d e r s p a t i a l difference a p p r o x i m a t i o n d i s c u s s e d i n E x a m p l e (1.4.4).

L e t the a p p r o x i m a t i o n to the boundary c o n d i t i o n y ' ( l ) = 0 have the

form

2

(6.10) 1

=

_

i

b

=

u J

+

i

0

J

r

It i s now difficult to check ^ ( V ^ ) ,

so only the eigenvalues w i l l be c o n s i d e r e d .

P r o c e e d i n g as i n the p r e v i o u s e x a m p l e , i t i s found that the

eigenvalues

a r e g i v e n by

x. =

(6.11)

for

2

ll

1 1

1 2

2

( Z

-

1 ) 2

(z* +

10Z

+

1)

each z w h i c h i s a root of 0

2

(6.12)

i=Note that i f z i s a root of

" z"

then so i s \ .

(6.12)

j( ) = . Z z

1

Now

equation

(6.12)

=

) =

0.

P j

n o m i a l of the d i s c r e t e boundary c o n d i t i o n

c

( J + 1 )

_

r

b. z

(6.10)

The c h a r a c t e r i s t i c p o l y i s g i v e n by

J

J

has a root z s a t i s f y i n g | z | > i + e , J l a r g e , i f and

only i f the c h a r a c t e r i s t i c equation C j ( z ) = 0 has a root i f z i s a root of

(6.12)

with | z |

= 1 ,

|z[ > 1 .

Further,

then an easy computation shows that

the c o r r e s p o n d i n g eigenvalue i s r e a l a n d n e g a t i v e .

Thus i t i s sufficient

t o c h e c k t h a t n o n e of t h e r o o t s o f t h e c h a r a c t e r i s t i c

e q u a t i o n of the d i s -

c r e t e b o u n d a r y c o n d i t i o n e x c e e d s one i n a b s o l u t e v a l u e .

As

i n the p r e v i o u s e x a m p l e the s e c o n d o r d e r

h

(

2

J

U

"

2

u

J-l

2 J-2

+

u

)

=

approximation

°

s a t i s f i e s t h e e i g e n v a l u e c o n d i t i o n b e c a u s e t h e r o o t s of i t s c h a r a c t e r i s t i c equation a r e

= 1 and z^

the t h i r d o r d e r

approximation

6S