General Disclaimer

One or more of the Following Statements may affect this Document

This document has been reproduced from the best copy furnished by the 

organizational source. It is being released in the interest of making available as 

much information as possible.

This document may contain data, which exceeds the sheet parameters. It was 

furnished in this condition by the organizational source and is the best copy 

available.

This document may contain tone-on-tone or color graphs, charts and/or pictures, 

which have been reproduced in black and white.

This document is paginated as submitted by the original source.

Portions of this document are not fully legible due to the historical nature of some 

of the material. However, it is the best reproduction available from the original 

submission.

Produced by the NASA Center for Aerospace Information (CASI)
,..

F(3) 

^,rr•

t

JNASA-Ch-1753Z±) ASYN;c'TOTI(: AbALYSIS (:F 6d4-153tO

NUStRiCAL riAYE PhLPAGA'ilUN IN FINITE

iA FhhLNCL LvUATICiaS (Md5saC1'usett6 lint. Of

Tech.) 142 i3 nL AJ7/Mt A01 LSLL 1UN Uliclas a

U3/^l JU517

^A'1^1T^

2

A

GAS TURBINE & PLASMA DYNAMICS LABORATORY

MASSACHUEJETTS INSTITUTE OF TECHNOLOGY

CAMBRIDGE, MASSACHUSETTS

0
ASYMPTOTIC ANALYSIS OF

NUMERICAL WAVE PROPAGATION IN

FINITE DIFFERENCE EOUATIONS

by

Michael Giles

W. T. Thompkins; Jr.

GT&PDL Report No. 171 Match 1983

This research, carried out in the Gas Turbine dnd

Plasma Dynamics Laboratory, MIT, was supported by the NASA

Lewis Research Center under Grant No. NAG3-9.
c

2

Abstract

An asymptotic techni q ue is developed for analysing the

propagation and dissipation of wave-like F.)lu *_ions to finite

li`.ference e q uations. It is shown that for each fixed

com p lex frequency there are usually several wave solutions

with different wavenumbers and the slowly varying amplitude

of each satisfies an asymptotic amplitude equation which

includes the effects of smoothly varying coefficients in the

finite difference equation's. The local group velocity

ap p ears in this equation as the velocity of convection of

the amplitude. Asymptotic boundary conditions couplin g the

amplitudes of the different wave solutions are also derived.

A wavepacket theory is developed which preaicts the

motion, and interaction at boundaries, of wavepackets,

wave-like disturbances of finite length. Com p arison with 

numerical experiments demonstrates the success and

limitations of the theory.

Finally an asymptotic global stability analysis is

developed which gives results which agree with other

stability analyses and which can be applied to a wider range

of problems.

C
(

3

Acknow_edgements

This research was supported by NASA Lewis Research

Center under grant No. NAG3-9 with Technical "Monitor Dr. R.

V. Chima, and by a Scholarshi p fcr M. Giles from the Kennedy

Memorial Trust.

Mr.

E'.

t:d
OJ

4

..

Ml.

W4 .

Contents

Abstract 2

Acknowledgements 3

Contents 4

1. Introduccion 7

2. Amplitude Analysis of P.rtial Differential Equations 12

2.1 Fourier Analysis 12

2.2 Asymptotic Amplitucr Equation 16

3. Amplitude Analvsis of :' nite Difference Equations 24

3.1 Fourier Analysis 24

3.2 Asym p totic Amplitu•^4 :a Equation 27

3.3 Examples 33

3.3.1 Trapezoidal scheme 33

3.3.2 Box Scheme 37

3.3.3 Backward Eu1=tr Scheme 39

3.4 Asymptotic Boundar°),

Conditicns 40

3.5 Examines 43

3.5.1 Upstream Boundary 43

3.5.2 Downstream Foundary Space Extrapolation 44

3.3.3 Downstream Boundary 45

Space-time Extrapolation

3.5.4 Downstream Boundary . Box Method 45

4. Ray Theory and Wavepacket-Particle Duality 47

4.1 Pay Theory 47

4.2 Wavepacket-Particle Theory 50

5. Asymptotic Stability and Convergence Analysis I 58

5.1 Theory 58
5

5.2 Example 64

6. Asymptotic Stability and Convergence Analysis II 67

6.1 Theor

6.1.1 Iterative Solution

y

E7

69

6.1.2 Newton-Raphson Solution 71

6.1.3 Stability Domain Method 71

6.2 Exam p le 72

7. Assorted Examples and Further Developments 77

7.1 Instability of Backward Euler with 77

Space-time Extrapolation

7.2 Optimum CFL Number For Trapezoidal Method 80

7.3 Discontinuous CFL Number In Traoezoidal Method 84

8. Assorted L000e Ends 87

8.1 Degeneracy 87

8.2 Eigenmodes and Eigenf.:equencies 90

8.3 Other Asymptotic Ap p roaches 92

9. Wavepacket Test Program 94

9.1 Program Description 90

9.2 Theory. 97

9.3 Numerical Results 103

9.3.1 Traoezoidal Method With Variable 103

CFL Number

9.3.2 Box Method With Variable CFL Nuvbe.r 105

9.3.3 Trapezoidal Method With Space 107

Extrapolation

9.3.4 Trapezoidal Method With Space-time '09

Extrapolation

9.3.5 Trapezoidal Method With Box Boundary 111

Condition

9.3.6 Wavepacket Outflow In Box Method 113

1,^
6

9.3.7 Instability Of Trapezoidal Method With

Space-time Extrapolation

9.3.8 Instability Of Tra p ezoidal Method Due

To Varying Mesh

9.4 Program Listing

10. Conclusions

References

Appendix

A.1 Finite Difference Operator Notation

A.2 Principle Of The Argument

It A.3 Definitions Of Norms And Stability

K.

115

119

121

132

134

135

138

140 

(
7

1. Introduction

Consider the following very sim!Dle problem and

numerical solution. The partial differential equation is

au + c au 0 ( 1 . 1 )

at ax

where c is a positive constant. The domain

considered is O<x<1 . The initial condition is

u(x,0) - exp[-200(x-0.5) 2 ] cos(kx) (1.2)

with k-80. This form of distribution is usually

called a wavepacket. The cos(kx) term defines the

oscillation of a group of waves and the exp[-200(x-0.5)2]

term is an amplitu el - 'envelope'.

The upstream condition is

u(O,t) - 0 ( 1.3)

The solution of this problem is

u(x-ct,0) ct < x < 1

u(x,t) - (1.4)

0 0 < x < ct

F

The numerical solution uses a uniform grid with

computational domain O<j<200 and a trapezoidal scheme.

i- n+1

n- r n+ 1 n n+1 n 1 Tj

j

4` j+1 0+1-1 0 -1

where r = 

cdt ( 1.b)

Ax

In this example r-1. the initial condition is

U^ u(x.,0) (1.7)

J

and the upstream boundary condition is

Et
n _ n

U 200 U199

(1.9)

8

U n

0 (1.8)

0

In addition a numerical boundary condition is

required at the downstream boundary. For this condition

space extrapolation is used.

Figure 1 shows the numerical solution at intervals

E

of 60 time steps with each plot drawn to the same scale.

The first two plots show the initial wavepacket travelling

i downstream in the direction of the phy i,ial characteristic.

Corresponding wavecrests are labelled a-e,

and it can be seen

that the propagation velocity for the wave crests is greater

than for the amplitude envelope. Note for example that the

amplitude maximum lies approximately midway between crests b

and d at n-60 but at n-120 the,maximum is clearly nearer

crest b. At n=180 the numerical disturbance is interacting

with the downstream boundary. The solution appears to be

the sum of two waves, one with the original wavelength , and

one with a very much shorter wavelength. A^- n-240 there is

a reflected wavepacket of wavelength slightly greater than

2. and the plots at n=300,360 show that this wavepacket

travels back up the domain at approximately the same speed

as the origin3L wavepacket. This solution is clearly

numerical and not physical since the analytic, physical

solution moves from left to right across the domain and then

_ out the downstream boundary. The analytic equation does not

have any solutions with waves travelling from right to left.

At n=420 the wavepacket is interacting with the upstream

boundary, and at n-480 there is a reflected wavepacket with

the original wavelength. This completes one cycle. I_f the

solution was continued the wavepacket would travel down to

the downstream boundary and then reflect again into a

C+
9

,NAL Pk'OC i:J

U

URIGINA1- PAGE IS

OF POOR QIjALITY

c

0

U

0

0

n • 60 Wavepocket Disturbance

With Constant Frequency and

Finite Amplitude Envelope.

n = 120 'The Amplitude Envelope

Has Moved at a Slightly Smaller

Velocity Than the Wovecrests.

n = 180 The Wavepocket is

listerocting With tl,e Downstream

Boundary.

n = 240 The Reflected Wovepocket

Has Smaller Amplitude and a

Wovelength Slightly Greater Than 2.

n = 300 The Wovepacket

Travels Upstream.

n = 360 The Wavepocket

Continues To Travel Upstream.

n = 420 The Wovepacket is

Interacting With the Upstream

Boundary.

n = 480 The Reflected Wovepacket

Has the Some Wavelength as

the Original Wavepocket.

0

0 .25 .50 .75 1.00

X

FIGURE 1, NUMERICAL SOLUTION OF CONVECTION EQUATION

i=
l

 I

10

wavepacket with short wavelength, and a decreased amplitude,

travelling upstream.

The qualitative and quantitative prediction of the

behaviour of numerical solutions in problems such as the

above is one of the two objectives of this parer. The

second objective is a global stability analysis

incorporating boundary conditions and smoothly varying

coeff;^ients and predicting both stability and accurate

asymptotic estimates of convergence rates.

To achieve thEse aims a technique is developed to

analyse the approximate time evolution of an amplitude 

modulated wave, i.e. a wave with fixed frequency and a

slowly varying amplituae. Chapter 2 derives the theory for

partial differential equations, while chapter 3 derives the

theory for finite difference equations incorporating

smoothly varying coefficients and boundary conditions. In

the case of dispersi- r e, non-dissipative wave propagation, it

is found that the amplitude is convected at the local group 

velocity, a principle which is well understood in par;:ial

differential equations.

Chapter 4 applies the theory to the motion of

wavepackets which are wave-like disturbances of finite

length and constant frequency such as in the earlier

example. Chapters 5 and 5 der_ve global stability analyses

with different levels of asymptotic approximation. Chapters

7-9 develop further topics and examples including

comparisons between numerical experiments and theoretical

predictions.

Throughout this paper a finite operator notation is

used which greatly simplifies analysis and is a neccessity
11

for general proofs. Since there is no universally accepted

standard notation Appendix t details the notation ujed.

Very little previous p ork appears to have :;een done

along the lines of this paper. The r.oncepc of group velocity

in partial differential equations is well understood and is

explained in many texts 112• The asymptotic approach of

chapter 2 is not common due to the advantages of other

methods but is discussed by Whitham 1• Kentzer 3 has

I

discussed the use of group velocity in analysing finite

difference equations but does not derive a general Equation

for the amplitude or calculate the quantitative effects of

boundary conditions. Vichnevetsky and Bowles s derive the

the group velocity in finite differt-. • e equat^.ons using an

approach which is valid only for constant coefficients.

They also derive amplitude reflection coeffficients at

c

bo^indaries and discuss soma of :he examples given in this

paper.

t
I I 12

2. Amplitude Analysis of Partial Differential Equations

2.1 Fourier Analysis

Consider a homogeneous par-ial differential equation

u(x,t) n 0 < x< 9 , t> 0 (2.1)

where L is a constant linear differential operator.

Iefinel by,

M n

L

-7, C mn t3 x )! (atil (?.2)

m,n

and the coefficients C are constants.

mn

An eigenfunction of the operator L is a function

rt u(x,t) satrsf in ^ g

L u = ku (2.3)

where a is a constant called the eiQenvalue.

An etgenmode is a solution cf the homogeneous

}_ equation (2.7) _.e. -t is an eigenfunction with a,-ger.value

F zero.

3

as exp[i;kx-w1..)] - ik exp[i(kx-ut)] (2.4a)

3

a« exp[i(kx-wt)] - -iu exp(i(kx- wt)]

(2.4b)

3 •mr a I n m n

3x ^atJ exp[i(kx-ut)) - (ik) (-iw) exp[i(kx-,at)] (2.4c)
r

13

OF POOR QUALITY

1

7(a 1mr3 In

=

L C mnlaxj 13t)

exp(i(kx-wt)]

m,n

C (ik)m(-iw)n exp[i(kx-wt) 2.3d)

nn

m,n

Thus exp(i(kx-wt)] is an eigenfunction of

3 3 a M

22

n

3x 3t lax) (atj and L with eigenvalues ik -iw

(ik)(-iw)n, m and

/ C 

mn (ik) m (-iw) n

respectively

i m,n

t

Hence

u(x,t) = exp(i(kx -w*_)] (2.5;

is an exact solution of (2.1) provided

Cmn (ik) m

(-: y ) n 0 (2.6)

_

m,n

This relation between k and w is called the

dis p ersion relation.

Examples of dispersion relations are

+- Surface waves on dee p

water w4 _ ; qkI (2.7a)

s

Acoustic waves J2 c2 k2 ( 2 .7b)

y ;paves prGoaQating along a wavequide w 2 - c 2 (k 2 + k; ) ( 2 .7c ) s,

where q c a-id k o

are constants.

A general solution of (2.1 1 is a superposition of

eigenmodes which in the case of a oarti.al differential

e q uation is ex p ressed as an integral over all the

wavenumbers k of the sum of all the eigeamodes with

wavenumber k.

m N

n

u(:{,t) = J+^ An(k) exa[i(kx-r n t)] dk (2.8)

-Ce n=

t t

1 
14

If the dispersion relation is of order N in W , i.e.

it contains powers of W up to u N

then

W 1 (k) , W 1 (k) , .... , WN(k)

are the N values of W which satisfy the dispersion relation

for a given value of k and the 4 n ('') are the corresponding

constant amplitudes of those eigenmodes.

If A

n 

(k) is non-zero for all k,n then a neccessary and

sufficient condition for u(x,t) to remain bounded and not

increase exponentially is that each eigenmode must remain

bounded. Splitting u into its real and imaginary components

gives ,

W WR + iWI (2.9)

expf-iwt] - exp[-iW R t + WIt] (2.10)

Thus the condition that every eigenmode remain

bounded, and hence a general solution remain bounded, is

,I < 0 for all k,n.

This analysis is lacking in three respects. The

first is that in some situations the initial d_sturbance is

zero exce p t for a finite re g ion and one wants to know the

time evolution of this disturbance, in particular the

p ropagation velocity for the energy. The second failing is

that when the initial-value problem is replaced b y an

initial-value / boundary-value proil lem with boundary

conditions at x = 0,1 there is no eas y way to include the

effect of the boundary conditions in this stability analysis.

The third failing is that exp[i(kx- wt)]

is an eigenfunction

of L only when the coefficients C mn are constant. The

analysis breaks down when the coefficients are non-constant.

ti
15

The resolution of these problems requires the analysis of a

wave of constant frequency with an am p litude which varies

over a characteristic length scale much greater than the

wavelength.

0
16

f.-

2.2 Asymptotic Amplitude Equation

The problem now being considered is

L(x) u(x,t) - 0 < x < - , t > 0 (2.11)

where L(x) is a non-constant linear differential operator

defined bv,

3 L(x) = Cmn(x) (

a x)m (kt) "

(2.12)

m,n

and the coefficients C mn (x) are slowly varying functions of

X.

The theory calculates the approximate evolution of a

wavetrain with waves of a constant frequency u and a slowly

varying amplitude, so u(x,t` is written as

u(x,t) = A(x,t) exp(if(x,t)1 (2.13)

where A(x,t) is the slowly varying amplitude and

Y(x,t) is the phase of the wave which is related to the

frequency u ani wavenumber k by

3Y

at

-u (2.14)

aY = k (2.15)

?x

The frequency u is constant but the wavenumber k

will vary slowly with x because of the slowly varying

coefficients C (x) so the above relations can be intearated

w 

to give ,

x

Y(x,t) J k(^) d^ - ut (2.16)

0

a:

i 
17

,r IGINAL PACE i9

OF POOR QUALITY

To explain the asymptotic approximations which are

made two characteristic length scales L k and L A and one

characteristic time scale T A are defined. L k is the length

scale for variations in k, L p is the length scale for

variations in the amplitude A, and T A is the time scale for

variations in A. Numerical values for L k , L A and T A are

are given by ,

L k min I k / ax l (2.17a)

L A = min A / (2.17b) ax

(

/

T A = min A / t 2. 17c)(

The asymptotic aprr^ximations used in this theory are

L k >> k -1 L A >> k -1 T A » -1

which imply

k

<< k2

2x << Ak

at

<< Au

A Taylor series expansion of A and 'f about a point

( x, ,t, ) gives

A(x,t) = A, + A (x-xo) + - to) ate(t + H.O.T (2.18)

X

(

I

f(x,t) - Y 0 - W(t - t o ) + 1 k(^) dF

X0

X

fo - u(t-to ) + f Ek e + k (^-xe ) + H.O.T j d^

a

xe

o
.J 

iB

OF P' .

'o — u(t — t o ) + k o (x — xo) + 2 ax (x—xo ) 2 + H.O.T. (2. 19)

0

Subscript o denotes terms evaluated at (xo,to).

The H.O.T. , higher order terms, includes terms like

3 2 A AAA a2k

a xe

at e

axe 

A which are O{A(LAk)- 2,A(TAu)-2,A(Lkk) -2]

and are necilected in this asymptotic approximation.

exp (if(x,t)I -

l

= exp (ifo + ik o (x—x o ) — iU(t — to) + 2 aX (X—X, )2 + H.O.T

11 o J

exp ( if 0 + ik o ( x — xo) — iW (t — to) ] ; 1 + t a x o

( x—xo ) 2. I

+ H.O.T. (2.20)

Hence ,

u(x,t) - exp[i?, +iko(x - xo) - iu(t - to)] ( 1 + 2

i 3k

ax 

(x_xo)21 l o

JJ

( Ao + ax 

( x —xo) + at (t—to) 

J 

+ H.O.T.

6 0

exp [i'f 

o + ik o ( x -xo) - iW ( t -to ) 1

• IAo + aX (x-xo) + at (t-to) + 2 A o X (x — xo)2

fil 0 0 0 J

+ H.O.T. (2.21)

To evaluate derivatives of u(x,t) at (xo,to) a

two-variable version of Leibnitz 's rule is used.

( ax) m(at)n[ f(x,t) g(x,t) ]

m n

f 

(m-p)

T ^

f 

q i 1' q 1 ax

(a \ p ^ a_l o f ^ a_\ m-p , a—\ n-q

p

g

L • . (n- ) . ) at) a x

IlatJl

p=0 q=0
19

(NA I

OF POOH ^e^-" 1 ;''

(2.22)

Let

f(x,t) m exp[iY' o + ik 0 (x-x 0 ) - iu(t-t o )] (2.23)

g(x,t) 3 Ao + aA (x - xo) + ato(t-ta) + 2 A,ak (x-xo)2 (2.24)

Then,

(ax)

P

 (at)4 fl = (ik o ) p (-iu) q exp[iY o ] (2.25)

a

g 0 = A 0

( 2.26 )

aq = aA

(2.27)

axa axo

19 =

aA

(2.28) at, ata

zs

iAo ( 2.29)

a x2 0

a X0 

and all other derivatives of g(x,t) evaluated at (x o ,t o )

are zero. Hence ,

(2

a

^

_)m(at)n,1(x,t) = exp[iT0

0

A0 (ik o ) mn + m (iko)m-1(-iu)n + n 

at

(iko)m(-iu)n-1

I

^—x , a

+ i m(2-1) A0 aX( iko) m-2 (-iu) n

+ H.O.T. (2.30) 0

and so ,

O
20

=a.

OF PUCE? JUAL',I tf

L(x) u(x,t) = exp(i?l Cmn ( A(ik) m (-iu) n +

m ,, ,,,,n l

m 

dX (ik) m-1 (-iu) n + a 

3A (ik)m(-iu)n -1 +

.m(2 -1)

 A 

aX 

!ik) m-2 (-iu) n j + H.O . T. (2.31)

To satisfy the homogeneous equation (2.12) the

amplitude A(x,t) must satisfy ,

a ° (k,u,x) A + ai(k,u,x) 3t + az(k,u,x) A + a^(k,u,x) A 3k

ax

= 0 + d.O.T. (2.32)

where

Cmn( x} (ik) m (-iu) n

(2.33a)

m,n

a l (k,u,x) _ T, C mn ( x) n (ik ) m(-iu)n-1

m,n

= 2. a° (2.33b)

a2 (k,u ,x) _ Cmn( x) m (ik)m-1(-iu)n

m,n

-i ak° (2.33c)

and a,(k,u,x) _ Cmn( x) i m(2-1) (ik)m-2(-iu)n

m,n

i a 2 a °

=

2 3k2 (2.33d)

Because of the asymptotic assumptions

A >> u-1 at , k 1 ax , k- A 3x

so ( 2.32 ) can only be satisfied if

fi

t

0
21

a o (k,u,x) 0 (2.34)

This is the dispersion relation between k and u.

k is now a slowly varying function of x due to the slow

variation in the coefficients. Thus the characteristic

length scale L k is related to some characteristic length

scale L„ for variations in the coefficients.

Neglecting the H.O.T. and dividing by a 1 gives the

asymptotic amplitude equation.

at + c ax e A (2.35)

g

where c g - a 2 / a l

(2.36)

and E _ - a, ax / a l

(2.37)

If c is real the left hand side of (2.35) is a

I[ g

Lagrangian-type total time derivative with respect to an

observer moving with velocity c

g

. If the coefficients Cmn

I' ar,: all constant, k is constant, e - 0 and so the amplitude

A is constant along rays moving with velocity c

9 *

In a

wavepacket individual wavecrests move with phase velocity

u/k , usually denoted c P , but the wavepacket, or amplitude

'envelope' , moves with velocity 

cg* 9

For this reason c is

called the group velocity.

Because the group velocity is the propagation

velocity for the amplitude and energy of the wavepacket the

group velocity is often more important than the phase

velocity. One example is the Sommerfeld radiation condition

which states that the waves generated by a fixed source have

a group velocity directed away from the source. In some

f

^_ J
E

i

22

unusual cases the phase velocity of the waves is actuall,r

directed towards the source. A second example is that the

group velocity never exceeds the speed of light which is the

limiting speed of propagation of information while the phase

velocity can exceed the speed of light , as happens in wave

propagation along an electromagnetic waveguide.

To link this derivation of group velocity to other

derivations the dispersion relation (2.34) is differentiated

with x held constant.

da, - a k ° dk + au ° du

0 (2.38

 aa, aa, Hence ( 3w)

ak, x const^ ak au

a 2 / a,

C (2.39)

g

The most common method of showing that (w)x

a

a const

is the group velocity uses the method of stationary phase

which is well explained in the available literature [1,21.

The usual approach is to combine the dispersion relation

the . definition of the group velocity and some physical

principle suc as energy conservation to calculate the

pro p agation of energy. The approach giver, above is not

usually used partly because sometimes the exact partial

iifferential equation is not known and the dispersion

relation has been determined by asymptotic methods (e.g.

water waves) or from empirical data (e.g. seismic waves).

This approach is however suited to analysing finte

difference schemes in which the exact finite difference

equations are known and there _s no general equivalent to
23

the principle of energy conservation.
24

3. Amplitude Analysis of Finite Difference Equations

3.1 Four-er Analysis

Consider a homo g eneous finite difference equation

( L Un 0 (3.1)

J

As explained in the appendix A.1, L can always be

expressed as a sum of step operators,

L = C E E (3.2)

L, mp mx pt

m,p

where the coefficientsare constants, but it is

-mp

often more simply expressed as a polynomial of finite

difference operators written symbolically as,

_

L =_ P(Ex,6x x,E t , 6 t ,Nt) (3.3)

The eigenfunctions and eigenmodes of L are defined

exactl yas ^.n 42.1. The finite operators all have the same

eigenfunctions, exp(i(jm-nQ)J. m and A are related to the

wavenumber k and frequency u of the physical wave being 

modelled by,

r"

9 kdx (3.4)

A uat (3.5)

;^

v

As shown in the appendix A.1

EK exp(i(js - nil)] - exp(im) exp(i(js-nA)] (3.6a)

'. 5x exp(i(js-nA)] - 21 sin($/2) exp(i(jm-na)j (3.6b)

μ x exp(i(j0-na - cos (m/2) exp(i(jm- nA)] (3.6c)

E t

exp(i(j$-nil)] -

exp(-in) exp[i(j0 - n9)j (3.6d)
U exp [ i( j s -nn) J

J

(3.7)

0

25

6 t

mexp - -2i sin( exp[i(j- nn)](3.6e)

V exp [i(jm- nn)] - cos (a/2) exp [i(jo-nn)J (3. 3f)

Emx E Pt exp [ i ( j m -nn ) I - exp [ i (m(P -pil )j exp [ i (j m -nn) ] ( 3 . 1+q )

Thus exp [ i ( j m -nn ) J is an eigenfunction of

E x' 6 x' μ x' E t' 6 t' U and E mx E Pt with eigenvalues

exp(i9), 2i sin( $/2), cos (m/2), exp(-in), -2i sin..(9/2),

cos ( Q/2 ) exp [ i ( m® -p ►1 ) J respectively, and L has ei(renv,.lue

C

mP u exp (i(ms-pn)] or

m,p

P[exp(im),2i sin(O/2),cos(4/2),

exp(-in),-21 sin(n/2),cos(n/2)1

de p ending which expression for L is used.

is

 Lp

C

an exact solution of (3.1) provided

mp exp [ i (m -pit ) m J ^ 0 ( 3.8

This is the dispersion relation between a and n.

Since exp[2nwi] - 1 for all integers n, o + 2n+r is

equivalent to a so only solutions in the ranges

-+r < Re( s ) < a

-R < Re(:2) < +r

need be considered.

If L involves P+1 time leveli and M+1 spatial nodes

the dispersion relation is a polynomial of deg r ee P in

exp (-j.:1 ) and of degree M in exp (im) . Thus for i given m

,:.-iere are P corresponding values of 3, and for a given n
exp[-in+:) - exp( - inA R + nQII (3.11)

26

there are M correspondi n g values of m.

A general solution of (3.'i) for periodic boueary

conditions is a superposition of eigenn,odes.

P

U3 _

)P1Ap

(m) exp(i(jm-nip))

P _ 1

The 0 summation is a summation over all the values

of a which satisfy the periodi-- boundary conditions, and

the p summation is over the P values of A corresponding to

each v..=.lue of s.

If A p ( $) is non - zero for all s,p the, a necessary and

sufficient condition fcr U to remain bounded and not

J

increase exponentially is that each vigenmole must remain

bounded. Splitting Z into its real and imaginary components

gives,

A = it R + ial ( 3. 10 )

.9)

Thus the condition that every eigenmode remain

bounded , and hence a general solution remain bounded , is

.1 < 0 for all s,p.

This analysis is lacking in the same three rrispects

as the analysis of partial differential equations by

e.igenmode exp,insion in the last chapter. The analysis

gives no information about the movement of an initially

localised disturbance, cannot incorporate boundary

conditions or anal y se schemes with non-constant

coefficients.
Z7

3.2 Asymptotic amplitude Eauation

in this section the coefficients C in the

mp

definition of L (3.2) are assumed to be slowly varying

functions of j. The analysis is performed in computational

space with coordinates (j,n) in whicn the grid spacing is

Aj-1, 4n=1. Variations in mesh spacing in physical

coordinates are incorporated directly into the variable

coefficients of the finite difference equations.

The theory calculates the approximate evolution of a

wavetrain with waves of a constant frequency 0 and a slowly

varying amplitude , so U is written as

Un , A(j,n) exp(if(j,n)J (3.12)

where A(j,n) is the slowly varying amplitude and

Y(i,n) is the phase of the wave and is related to the

frequency a and wavenx.ber :0 by

(3.13) ` an

3? 9 (3.14)

which can be integrated to give,

I d^ - nQ (3.15)

0

- As in 12.2 two characteristic length scales, L^ for

variations in 1, and L A

for variations in the amplitude A

and one characteristic time scale T A

for variations in the

t amplitude A , can be defined with numerical values beii.g

given by,

O
28

r PCB'_

Ly - min (1 / a 1

ajl

( 3 . 16a )

L A - min (A / a^ 1 ( 3 . 16b)

T A - min ( A /

a n

( 3 . 16c )

l

The asymptotic approximations are

L 9 >> 1 LA >> 1 ^A > > 1

which imply

am «1

aj

A 3n«A

J J

A Tavlor series expansion of A and 'f about a point

(j,,n,) gives,

SA ?A A(j,+m,n,+o) A,+ m + p

+ H.O.T (3-17)

a j 0 an,

jo+m

Y ( j o

+m , n , yp ) _ ?, - pit + d{

j,

jo+m

= To - P 11 + f (m^ + (`-jo ) + H.O.T J d^

ajo

Jo

a am f, - p11 + -nm, + + H.O.T. {3.18) 2 ajo

Subscript , denotes terms evaluated at (j,,,,)•

The H.O.T. , hiaher order terms, includes terms like

2 2 2 a A a A

a L A which are O(ALA -2 ,AT A -2 ,AL 9 -I I and are

a j 2 3n' a j 2

neglected in this asymptotic approximation.
29

ORIGINAL PAC'- rp

OF POOR QUAL1 i `r'

exp(i f(jo+m ,no+p)]

iM 2J o M Axp [ if, - ipQ + imo o + + H. 0 . T 2 a jo

J

im zam

- exp[ if, - ipm + immo] 1 +

2 ajo

+ H.O.T. (3.19)

Hence

2

U j °o ' aXp ( if o] exp[i(mmo-pil)1 I 1 + 12

+ m

ado

3 3A

Ao + m + p +

ajo ano

H•O.T .

exp(ifo] exp[i(mmo-pil)] A, + m pa ' + p

an

- i2 A.

Jo 0 0

+ H.O.T. (3.20)

Hence,

n c

L 

U j exp[ifJ Cmp(j) exp[i(mm-p1)] .

m,p

^ 2

(A in

A

d- + m 3A + p + 3n +

11 J J

H.O.T. (3.21)

To satisfy the homogeneous equation (3.1) the

amplitude A(j,n) must satisfy,

F ao(9,A,j) A + a l (m ,Q ,7)

aA

+ a 2

(m,f2,j) aA

+ e 3 (a,2,j) A

am

an aj aj

0 + H.O.T. (3.22)

where .

T, C

mp (j) exp [i(mm-pfl)] (3.23a)

M ,, ^,p

a 1 ( m j ) _ T, C ( j ) p exp [ i ( no - PQ

ml,p l

i an la aol (p,j const (3.23b)
O

OF p00K QUALO Y

aj(m,p,j)

T, 

Cm 

P 

(j) m exp(i(m4-nil)]

LLL^^^ 

m 

a ael

-i - (3.23 (-)

m ,j con st

and a^(m,a,j) Cmp(j) 

i2 exp(i(mm-F:3)]

m,p

i ra j a l

- 2 ^3m j J 9, j const ( 3.23d)

In the above derivat:.on of a, ,a l ,a j ,a, the general

shift operator expression for L. (3.2) is used. In

applications it is more convenient to use the finite

operator polynomial expression (3.3). a, is obtained by

replacing each operator with its corresponding eigenvalue

and then a l ,a j ,a, are calculated by differentiating a,.

Because of the asymptotic assumptions

3A 3A am A » 3n aj A 2j

so (3.22) can only be satisfied if

a,(9,.I,]) - 0 + 0 ( L^-1 LA-1 TA-1 t (3.24)

This is the asymptotic form of the dispersion

relation between v and a and will usually be satisfied by

setting a, identically equal to zero. m is now a zl,,wly

varying function of j due to the slow variation in the

coefficients. The characteristic length scale L $ is related

to some characteristic length scale L C for variations in the

ccefficients

Neglecting the H.O.T. and dividing by a, gives the

asymptotic amplitude equation.

ir

O
O

3A 3 

E A

an + rg aj (3.25)

31

ORIGINAL :- -- ,a

OF POOR Q:. A'- r" Y,

where rg = a 2 / a l (3.26)

am and E _ - ( as + a ° ) / a l (3.27)

^ aJ

Differentiating (3.24; with n held constant gives,

i

aa° as °1 da ° _ (a )

dm + ( 3j i

dj m n , j const 0,a const

= 0 + H.O.T. (3.28)

Hence, neglecting the higher order terms:

0

a jj - (aj) $ ,A const

0

/ (am) Q , j const

(3.29)

= i la a ( aj o) m, n const

/ a2

aj E -i a^

a j m , Z const (so )

al a 2 / ( 3 . 30 )

If r is real the left hand side of (3.25) is a

4

Lagrangian-type total time derivative with respect to an

observer moving with velocity r . Thus the amplitude A is

g

beinq convected with velocity r in computational space.

g

Differentiating (3.24) with j held constant gives,

aao da° s(an ) do as°)

+ (am 9,j dm const n,j const

= 0 + H.O.T. (3.31)

Hence, neglecting the higher order terms,

(

s aa°

a

a

9 ) j const (11s ) o n , j const (ail ) m , j const

- a2 / a:

= r (3.32)

9
32

0^ POOR QJALt7Y

Substituting for a and m using (3.5a,b),

a (WA t)

r a (kAx)

At au

-x ak

c At

^x

(3.33)

Thus r is the CFL number corresponding to the group

9

velocity in physical space of the propagating numerical

wave. It is the number of spatial mesh intervals which a

localised disturbance travels in one time step. For the

rest of the paper r

S 

is called the group CFL number.

V

C
33

3.3 Examples

The model problem which is considered is,

at + C( X) u - 0 < x < m (3.34)

ax

Three different methods are anal.ysed.

3.3.1 Trapezoidal Scneme

The trapezoidal scheme is

1 n+1 n ci n+1 n n-1 n

(U. - U.) + (U. +U. ) - (U +U ) = 0 (3.35)

At ] ] 4Axj^ 3 +1 ]+1 j-1 j-1

which can be written using operator notation as

r, 1

It 6 t + 2 ] 62c ut 

J

I Un +2 = 0 (3.36)

 ]

:there the CFL number r is defined as

r. =

 c A t

(3.37)

J ]

and Ax. _ 1 (x, - x.]

)

(3.38 ]

)

] 2 +1 -1

a, is obtained by replacing the operators by their

eigenvalues.

a, _ -2i sin(0/2) + Z 2i sin($) cos(0/2) (:.39)

The dispersion relation is

a, = 0 (3.40)

so tan ( 0/2)= Z s4.n ( m) (3.41)

a 2 , a 2 , a, are obtained by differentiating a,

Opp
34

Jr I GOM

ilI

3a°

a l

 = 1 an

= cos(0 /2) + 2 sin (i) sin(n/2) (3.42)

a2

r cos() cos(n/2) (3.43)

1

a2 a°

a, - 2 i am2

1 3 2

2 3$

_ - 2 sin ($) cos(Z/2) (3.44)

E

Using the dis p ersion relation a l and a, can be

simplified.

a l - cos ( 0/2) + r sin (0) sin(0/2)

= cos(0 / 2) + tan ( n/2) sin(n/2)

cos t (n/2) + sin 2 (n/2) ] / cos(n/2)

= 1 / cos(n/2) (3.45)

a l _ - 2 sin (o) cos(n/2)

- tan(0 / 2) cos(n/2)

-• sin(W2) (3.46)

so r

g

a; / al

r

r 

cos(m) cos 2 (n/2) (3.47)

and E _ -i a, a°l $ ,n const / a t a, a° ! at

zr sin (0) cos(n / 2) ( a^ sin(v) cos(n/2) ]/[r cos($)]
r

OF PC)OR QUAL1Ty

1 3r - —s in = (s) cos' (Q/2) / cos (m ) ( 3 .48 ) 2 33

There are three points of interest

i) sin(n- v) = sin((V)

so for all a there are two corresponding values of

m given by the dispersion relation,

i I satisfying -1r/2 < Re((P l ) < n/2

F and 0 = a - ml

ii ) For real it in the range

0 < I < 2 tan - ' (r/2 )

0 < tan(n/2) < r/2

so 0 < sin(m) < 1

Thus m l and 0, are both real and

0 < $ I < it /2

so r ((P 1 ) > 0

g

and n/2 < 9 2 < it

sor (y= ) < 0 g

riencc for every frequency in the given range there

is one forward travelling wave, travelling in the same

direction as the physical waves being modelled, and one

backward travelli.g wave with wavelength less than 4Ax.

iii) For real A in the range

2 tan -, (r/2) < 0 < ,r

sin(9) > 1

so m l and 0, are complex

Ii C46
0

36

Let m l - 1F/2 + is 

Then sin($,) - cosh(OI)

so m I is real and satisfies

tan(Q/2)

2

 cosh(s

0 2 - r - m 1

- r /2 - is 

These are evanescent waves. If there are

boundaries at j - O,J and the boundary conditions force a

steady oscillation with a frequency in the given range one

wave will decay in amplitude exponentially away from the

boundary at j-0, while the other will decay exponentially

away from the boundary at j-J.

M,r

G 
37

3.3.2 Box Scheme OF pOO(^ Qi_:All ^ d

The box scheme is

1( n+1 n+1 n n l

at J J J J J

c. l

^ ) j ( (U j+1 + Uj +1) - (Un + U^ +1 )

J

 = 0 (3.49)

which may be written in operator notation as

1

Gt ux6t + r j u t a x ) Uj+f = 0 (3.50)

a,= -2i cos ($/ 2) sin ( A/2) + 2ir cos(:2/2) sin($/2) ( 3.51)

The dispersion relation .s

f:an(0/2) - r tan($/2) (3.52)

a

as

, = i ano

= cos($/2) cos(0/2) + r sin ( A/Z,' sin(o/2) (3.53)

a a, a:

a^

= sin($/2) sin(a / 2) + r cos ( n/2) cos(0/2) (3.54)

1 a 2 

a, ao 2 a-^

= 4 cos($/2) sin(A/2) - 4 r cos(n/2) sin(1/2)

_ 6 

ao

0 (3.55)

Thus r a, / a,

4

_

sin($/2) sin(Q/2) + r cos(0/2) cos($/2)

cos(0/2) co3(A/2) + r sin(9/2) sin($/2)
38

tan( */2) tan(n/2) + r

+ r tan(($/2) tan(n/2)

n

_ r _1 + tan 2( */2)

(3.56)

1 + r 2 tang ( s/2 )

and e - 0 (3.57)

There are two points of interest

i) For each real value of n there is one corresponding

real value of * given by the dispersion relation and the

group CFL number r

g 

is real and positive.

ii) When r - 1, r

g 

-r , so waves of all frequenc-es travel

at the same velocity as the physical waves being modelled.

This is because when r-1 the Box scheme reduces to,

UJn+1 - U^ (3.58)

which agrees exactly with the solution of the partial

differential equation,

u(x+ct,t) - u(x,0) (3.59)

^ o
39

PAQ

OF POOR QUALITY

3.3.3 Backward Eule- Scheme

The backward Euler scheme is,

1 (U n+ 1 - U ) +

ci

(Un+ 1 - U rj+ 1

) = 0 (3 - 6 0 ) of ] ] 2 A x ]+1 ]-1

which may be written using operator notation as

0t (7 t + 2

J

r 6 2x ) Un+1 =

J

(3.61)

a,= 1 - exp(in) + it sin(m) (3.62)

f

i The dispersion relation is

exp(inJ - 1 = it sin(m) (3.63)

3 a, al i 30

= exp(in) (3.64)

la,

a,

= r cos O (3.65)

13 2 a,

a,

- _ 2 i

T-1

= 1 3a2

z 3

- Z sin(>t)

2

(3.66)

e

r.

so r = a, / al

4

= r Cos (0 exp( -in) (3.67)

and E

i 

3a,

-1 3

J/ a - / const 1 a 7 a 0 a ,

= lZ sins)) ( i

a

sin(o) J / [ exp(i.1) r cos(m) J

1 3r 3

2 a] sin (a) exp( -in) / cos($) (3.68)
40

,F PGGR QUALITY

3.4 Asymptotic Boundary Conditions

The general solution of

n

L , i1 . 0 ( 3 . 1 )

J ]

is a sum of waves with different constant frequencies A and

slowly varying wavenumber m and amplitude A

M 7

Z

U. Am(j,n) exp(i(

J

V m (^) d{ -nil (3.69)

7

fl m i 0

The outer summation is over different a 'ues rf Q,

and the inner summation is over the M different values of 9

which satisfy the dispersion relation for each A.

For each 11 ,m the amplitude A satisfies its

asymptot : amplitude equation on the interior of the

computationa l_ comain independent of all the other waves.

All the waves of each frequency are however coupled by

boundary conditions.

Suppooe a finite difference boundary condition at

)-J is

B U 

nn

, F (3.70)

where 9 :.s a constant tinite dlfftrence operator which can

be expressed in ooerhtor polynomial form as

5 = P 3 (r x ,5 x ,li x ,F t ,s t ,u t ) (3.71)

and F n ;.s a forcing function which . ::an be expressed as a

sum of inputs of different frequencies.

F 1

f(n) exp(-i7) (3.72)

Q
41

Of PO'R dUALITY

Performing exactly the same asymptotic expansion as

in the derivation of the asymptotic am p litude equation the

boundary condition becomes,

M

ab as ab aA_ i a-b am b 

A M 1 an 

anal - i 3m 33T 2 am t A m ajm

i2 m = 1

J

Exp(i( f O m d^ - nil)] f(52) exp(-ina) + H.O.T. (3.73)

0 A

where b(n,9 m P B [ exp(is ) = m ),21 sin($m/2),cos(am/2),

exp(-ia),-2i sin(A/2),cos(:Z/2) J (3.74)

The coefficients of exp(-inQ) in (3.73) must be

asymptotically equal to zero for each A so,

M ^

ab aA _ ab aA _ i a 2 b am

 ^ b Am + 1 TO anm 1 am aj m 2 amt A m ajm

Lt

J

exp[i co d^ 1 f(Q) +H.O.T. (3.75)

J m

0

This pa p er is primarily concerned with stability

and convergence rates. When analysing perturbations from a f

steady state or constant am p litude oscillation the boundary

condition ftr th perturbation has

f(Q) = 0 (3.76) n

Because the zero order terms will usually dominate

the normal form of the asymptotic boundary conditions is,

1

M

J

b('.I m ) A m (J) exp[i f 9 m d& J 0 (3.77)

m _^ 0
42

The first order terms 3Am

2--Am 3—Om do 3A a are only

] j

important when,

in ) 0 ( T' A -1

LA - ► L * ` )

As explained in §3.1 if the finite difference

operator L in the interior scheme spans M+1 spatial levels

there will be M values of s given by the dis p ersion relation

for a given value of a. If the computational domain is

0<j<J the interior scheme gives finite difference equations

at J-M+1 nodes, so to complete the set of finite difference

equations there must he M finite difference boundary

conditions. Hence for each A the asymptotic amplitude

analysis gives M inde p endent amplitude differential

equations on the interior coupled at the boundaries by M

boundar y conditions.

l
43

ORIGINAL PAGE ;9

OF POOR QUALIV

3.5 Examoles

The same model pro:11em as in 63.3 is considered,

3uau + c 0 (3.34) `3 t ax

0 < x < x 

C(X) > 0

The analytic boundary condition is,

u(O,t) = F(t) (3.78)

For perturbation analysis

u(O,t) - 0 (3.79)

The finite difference scheme using the trapezoid&l

or backward Euler methods on the interior re q uires two

finite difference boundary conditions. For perturbations

the boundary condition at j-0 is,

U0 , 0 (3.80)

The boundary condition at j=J is some form of

extrapolation. Four of the most commonly used are analysed.

3.5.1 Uostream Boundary

U n

0 (3.80) 0

B = 1 (3.81)

so b 1 (3.82)

Hence A 1 (O,n) + A 2 (O,n) = 0 (3.83)

in preparation for the theory developed in chapters

4 and 5 it is useful to define R 0 , the amplitude refraction

coefficient as

C
P

me

44

R O A1/)

(3.84)

A 2 (0 n

so in this example

R O - -1 (3.85)

3.5.2 Downstream 3oun datSpace Extrapolation

The space extrapolat:.cn buur dary condition is

U J UJ-1 (3.86)

3 1 - E-

_ E-fix b

x

(3.87)

so b 1 - exp(-im)

= 2i exp(-im/2) sin(9/2) (3.88)

Hence

2 J

c-1

si:.(mm /2) exp(i jr m^ d& - i 9 m (J)/2 ] Am (J,n) -. 0 (3.89)

n= 1 0

The am^.litude reflection coefficient R is defined J

asp

^ 2(J,n) RJ F.(J,n) (3.90)

so in this example

J

RJ - sln($2/2) eXp

1 f (m l - m t ) d; - 1 ( ,0, ( j )- $2 (j)

z

0

(3.91)

G

a _ O 

,,,, . ... s..^.d. -_ — - - _ - - -
) 

45

UKIGIIV-AL. p^p,^_ITY

OF poop, Q

3.5.3 Downstream Boundary Space-time Extrapolation

The space-time extrapolation boundary condition is

U J UJ-1

(3.92)

B 1 - E E (3.93)

-x -t

b 1 - exp(-is + iA)

= 2i exp(i(Q-m)/2J sin[(*-a)/2] (3.94)

Hence,

2 m -A` J

sin( -2 I exp i I m m d ;'- im^(J)/2 Am(J,n) 0 (3.95) 1.

M- 1 / 0

sin ro ,

 -,71 J

J

R J = -

(

2,1 exp i i f (m 1 -9 2 ) d^ - 2(0 1 (J)-9z (J) )

sing z2.

I` 0

{3.96)

3.5.4 Downstream Boundar., : Box Method

In this example the Box method which was discussed

in 43.3.2 as an interior scheme is now considered as a

downstream boundary condition.

( -1 u x i t + r u t 5 x } UJ = 0 (3.97)

z

s

r

B E_X E_^ t

(

2

u x d t + r u t 6 x ] (3.98)

so b exp(i(a -9)/2] i -2i cos(s /2) sin( ►Z/2) +

2i r cos(A/2) sin((p/2) } (3.99)

Hence,

x^

O
46

ORMINAL F J,,= r^'

OF POOR QUALITY

2

(r sins /2) - tan(f2 /2) cos(m /2)j

m 1

11 m m J

J

Exp[i f A d^ - i (p(J)/2 J A (J, n) 0 (3.100)

m m m

0

R r sin($ 1 /2) - tan(1/2) cos( /2)

= -

J r sin (0 2 /2) - tan( 0/2) cos (m2/2)

J

exp i f (m l -D 2 ) d& - `(m l (J)- 2j ) (3.101;

0

t

t
47

orscwv. r;,, .. t

OF POOR

4. Ray Theory and Wavepacket-Particle Duality

4-1 Ray Theory

In addition to the asymptotic approximations made

in chanter 3 this chapter assumes that for all real

wavenumbers 1, the frequency Q is real for all j and hence

tale group CFL number r

9 

is real.

and

— (

r g ,3¢ j const

3.32)

A Lagrangian - type total time derivative in

computational space is defined by,

- an + 

r (4.1) do a

9 j

so 1-1 _ 2-2+ do an

rg aj

r (4.2)

9

From the asymptotic amplitude equation ( 3.25),

dA = 3A + r aA

do a n g a j

= EA (4.3)

and using ( 3.26) and ( 3.29),

dm _ 3m am

do an + rg a 

a^

a? 1 (a j m ,;1 const / a2

=o

i ( 3

s1

a ^j ) s 4 const (4.4)

iC J
F.P

I L

URIGINAL PAGE Is

OF POOR QUALITY

A general initial value problem for a wave of

fre q uency A and wavenumber s(Q,j) can be solved by

integrating these equations (4.2)-(4.4) with initial

conditions

j(0) - j, (4.5a)

A(0) - A(j,,0) (4.5b)

m(0) _ $(n,j,) (4.5c)

Each value of j generates a ray and all of the

rays together cover the entire domain for n > 0. Figure 2

shows the motion of some typical rays in computational

space.

M!.

n

FIGURE 2, RAYS IN COMPUTATIONAL SPACE

11,(.01, o 
)

49

r is a function of j so at a particular j al'.. the

g

rays have the same slope d . Hence the time separation T of

two rays, illustrated in figure 2 for rays t and 2, remains

constant but the spatial separation L varies as r varies.

g

As explained in chapter 3 if the finite difference

operator spans M+1 spatial nodes then for a particular

value of A there are M values of s which satisfy the

dispersion relation. Define M+ to be the number of

solutions a for which the group CFL number is positive, and

similarly define M_ to be the number of solutions 9 for

which the group CFL number is negative. Let the

computational domain be 0 < j < J as usual. At j=J there

are M+ rays leaving the domain and M_ entering it. The

amplitudes are related through the asymptotic boundary

conditions each of which has the form,

M

Jr

b(Q,^: ) exp( i J m m d^j Am (J,n) - 0 (3.37)

m=^

0

Since the M + amplitudes of the rays leaving the

domain are known and the M_ amplitudes of the rays

entering the domain are unknown there must be M_ boundary

conditions to uniquely determine the amplitudes of the rays

entering domain. Similarly at j=0 there must be M+

^. boundary conditions to uniquely determine the amplitudes of

the rays entering the domain.

IP

4
50

4.2 Wrvepacket-Particle Theory

In terms of ray theory a wavepacket is a ray tube, a

group of rays, along which the amplitude is non-zero. From

the discussion in the last section the time lenqth T A of the

wavepacket remains constant but its spatial length L A will

vary whenever r

g 

has different values at the two ends of

the wavepacket. Provided L A << L 9 all the rays in the ray

tube have approximately the same value for s(A,j) so the

motion of the wavepacket is given by,

11 - rg(^,m(j),j) (4.2) do

and d,D iseal

do al 3j ) m,A const (4.4)

The energy, in physical space, of the wavepacket is

defined as,

x

J

E(n) f ^A(x,tn 2 dx

x

J

r dx

0 d 

Hence the wave energy density in computaticnal

space is defined to be,

o(j,n) - o(j) IA(j,n)l =

(4.7)

where a(j) = dx (4.8)

7

Using (4.3) and the notation that A is the complex

conjugate of A it follows that

)n + 8 (r o) . T( an (aAA) + a (r aAA)

g j 4

L

o 
51

n aA(3n + r Aa ) + aA(an + r aA ) + AA a-(r a)

4 J 9 j j 9

- aA do + aA dA + AA

2

) ( r a

9

M aAEA + aAEA + AA --(r n)

^J 9

I I E + E + 1 a—(r a) P (4.9)

`` a aj 9 J

Hence,

J

dF _ ao 

d

do an J

0

J

f - a (r P) + E + E + 1 a (: a) P dj aj 9 a aj q )

0

J J

a ll

_ - r P +f I E + E + -1 (r

9 a 3 a)

J

P dj (4.10)

9

0 0

If the wavepacket is in the interior of the domain

_ away from the boundaries the energy flux r o at the

9

boundaries j-0,J is zero. Also assuming as before that

L A << L^ then E + E + a —(r i) is approximately constant

J 9

over the length of the wavepacket, so

ll ^r

= + E +

aj(rga)

J J

P dj

do lE

a

0

e+ E+ a ( r a) 1 E ( 4. 1 1)

l a aJ 9 J

Thus equations (4.2), (4.., and (4.11) completely

describe the motion of the -daveparket particle in the

interior of the com p utational domain.

1
52

When the wavepacket reaches the boundary it

interacts with the boundary conditions to produce one or

more reflected wavepacket.s with the same frequency but

different wavenumbers. Th_ only case for which it is easy

to incorporate boundary conditions is when M=2 and

r g (S 1 ,j) > 0

r

g 

(m 2 ,j) < 0

An example of a scheme satisfying this condition is

the trapezoidal method which was applied to the model

convection problem n the introduction.

An additional assumption is that

e r g (* 1 , j) 1, e r g ($ " j) l << J

so that it takes much more than one time step for a

wavepacket to travel from one boundary to the other.

Sup p ose that initially there is one wave p acket with

wavenumber a l as in the introductory example. The

wavepacket travels to the right with position and energy

determined by the equations of motion previously derived

(4.2), (4.4) and (4.11). When the wavepacket :-caches the

boun'_ar y at j-,7 a proportion of the energy E 1 is reflected

into 3 left :.ravelling wavepacket of frequency R, wavenumber

a: and energy E ; . Figure 1 illustrates this interaction.

Equation (4.10) is

J J

dF

do r9G + f

lE

e a a ( z g a)1 A dj (4.10)

0 0 J

The outgoing energy flux is

O
V

at

53

Of POC`:4 QUALITY

r g (s,,J) P,(J,n) r - r g ( s 1 ,J) a(J) IA 1 (J,n) l = (4.12)

The incoming energy flux is

r ( m 2 ,J) o 2 (J,n) - r Q (s.,J) o(J) IA 2 (J,n)l =

(4.13)

g

The amplitude reflection coefficient 

R 

defined by

A2(J,n)

R J A,( J, n)

( 4. 1 4)

is a function of R,s l , o, determined by the asymptotic

boundar y condition.

The energy flux entering the reflected wavepacket Ja factor

r9((D,,J) IRJIi

rg(s1,J)

greater, or less, than the energy flux leaving the incident

wavepacket and so the total energy of the reflected

wave p acket is given by,

r ,J)

L_

(a

9 2 _ I I R J I 2 E 1

( 4. 15 )

r9(01,J)

The reflected wavepacket travels lift according to

.he eq uations of motion for a wave 2 wavepacket until it

reaches --0 where i_ is reflected into a right travelling

wave 1 wave p acket. The reflected energy E 1 is give- by,

E 1 I g 1 _ I I RO 12 E 2( 4. 16 )

r (s2,0)

g

AiiO,n)

where R

A2(0,n) O (4.17)

is determined by the asymptotic boundary condition.
(d) n

1

54

ORIGhIAL . _ Of

OF POOR QUALITY

Figure 3 t41ow s;, ows the particle-like path of the

wavepacket in computational space.

FIGURE 3. WAVEPACKET PATH IN COMPUTATIONAL SPACE 

nt_

In summary the equatio,is for the different parts of

Ar the path are,

ti

(a)

do rg!s1,J) (4.2)

d9, _ 1 3a o`

do ai(dj ) ^,^ cons (4.4)

dE l— 1 3

J 4

IF

C+
55

i

ORIGINAL 4'A-i 

0^ POOR QL AL ;Ty

r (9 2 ,J) _

(b) E2 - ^----- I

! R J I E ; (4.15)

r 9 (0 -"i)

(c) dl = r (m2 .j) do 9

(4.2)

dm 2 i

do a; m,A const

(4.4)

( 3, aj ,) 

dE 2

- 1 a )

E + E + (r a) E 2 4. 11

do z a

)

7 9

I r (a ,0) I

(d) E

1 

_ .^—^ JI I R O I 2 Ez (4.16)

r ( a 2 , OT g

In (a) and (c) the dispersion relation can be used

as a check on the accuracy of the numerical integrat:.on of

the equations or can re p lace the equation for the variation

of 9.

The total number of time steps for a round trip from

0 to J and back again to 0 is

J

N = J^ [ r (9; ,j) ) -1 - [ r ($.,.j) l -1 dj (4.181 

g g

0

The energy growth of a wavepacket travelling from 0

to J is given by

l

1n(E ) =

4

do ;

E dn

r

_ E + E + 1 a

(r a) (4. 19)

a aj g

so d- In( E, ) a- In (E; dj do do

1

i E

-

r E

+ r 1 a a ( r a )

It 7 9

g g
V

56 ORIGiNP

1f OF QOd^

E + E + a

(4.20)

r

aj( In(z g a)

g

Hence, J _ J

ln(E 1 (J)l 1n(E1 (0)) + r

E= E

dj + ln(r a) (4.21) J g

0 g 0

J

r ( m 1 ,J) E + E

so E 1 (J) = g exp f

I dj E 1 (0) (4.22)

rg(0110)

0

r

 g

Similarly the energy growth of the reflected

wavepacket as it travels from J to 0 is

J _

r (9=0) E + E

so E 2 (0) I g exp f - dj E 2 (J) (4.23)

r

9

(s 2 ,J) 0

r g

Combining (4.15),(4.16),(4.22) and (4.23) the round

trip energy amplification factor a is

)

J)

rg a= I

rg(m1,

r (01,0)

I ex

Jp

E- E

d'

r 3

I --

r (^l,J)

I I R

J

12

g C g 1 4

r (92,0)

•

I

J

( E + ^^

d^

I

r(^1,0)

G ( m 2 ,) I

E10 I:

r Q (92,J)

exp

r ^J= r 0 0 a

J

=

R0

R 12 exp Jr

^ t

_ll

E J1 ^E _ll EJ2

dj (4.24)

0

r g r g

Where I Ey E is evaluated at 9. , j

g 1

and IE

E1

is evaluated at 1 21 j

l

r 

g 2

The condition for stability is

a < 1

E CA o 
57

The equivalent average decay rate a is

a _ ln(a) (4.25)

N

When M>2 and there are M+ waves with positive r

g 

and

M_ waves with negative r , one wavepacket with positi ,,r e r

g 9

reachinc j=J produces M + reflected waveoackets with

negative r

g g

, and one wavepacket with negative r reaching

J=0 produces M + reflected wavepackets with positive r

9

Thus the total number of waveoackets increases with time

exponentially. Since each wavepacket has finite lengtr,

this ultimately leads to the problem of determining the

effect of interference between overlap p ing wavepackets. In

general the sum of the energies of two wavepackets is not

equal to the ener7y of the sum of the wavepackets if the

wavepackets are identical the latter is twice the former

while if the second wavepacket has opposite sign to the

first the latter is zero.

C^ o 
58

5. Asymptotic Stabili_v and Convergence Anal}_sj.s I

5.1 Theory

In addition to the as y mptotic approximations in

chapter 3 this chapter assumes that M=Z an y' there is one

boundary condition at both j=O and j=J and if m l is real at

j-0 then a j and m 2 are real cver the whole domain.

Exampled of methods satisfying these conditions are

the trapezoidal method a pp lied to tre model convective

p roblem with variable CPL number r, and the backward Euler

method applied to the model convective problem with constant

CFL number. Methods which do not satisfy these conditions

include Lax-Wendroff and Runge-Kutta type schemes. For

these methods the general stability analysis of chapter 6 is

required.

As explained in 53.1 a standard Fourier series

analysis of the problem with constant coefficients and

p eriodic boundar y conditions shows that the eigenfrequencies

are a(9) wnere 0 is a real wavenumber satisfying the

periodic boundary conditions and 0(0) is the corresponding

frequency g•iien by the dispersion relation.

The common use of a Fourier series analysis to

V

predict the stability of problems with non- p eriodic boundary

conditions implicitly assumes that for real wavenumber s,

:1(9) is a close approximation to an actual eigPnfrequencv.

This chapter follows that assumption, calculates a

F

correction a' to this a($) due to the boundary conditions,

and then determines the validity of the assumption for this

particular class of methods based on the asymptotic errors.

r

J
1-

59

vl,i.GI ` AL PAGE IS

OF POOR QUALITY

For real m 1 ,m, and the corresponding complex Q the

general solution is,

j j

U = A l (j , n ) exp i r $ 1 d^ -nQ + A, (j , n ) exp i J $ 2 d& -n0j

3 0 0

(5.1)

where the amplitudes A l and A, satisfy the asymptotic

am p litude equations

H

an

g m m m

m 

+ (r) 3A = EA m=1,2 (5.2)

and boundary conditions

b 1 ( Q ,m 1 ) A l (O,n) + b 1 ( 11,0 2 ) A 2 (O,n) 0 (5.3) f

^` I J

b. ( a ,a 1 ) exp^i 01 dj A l (J,n)

0

i J

+ b2(a,02) exp^i f 0 2 dj A,(J,n) 0 (5.4)

0

A linear system of finite difference equations with

time-independent coefficients and boundary conditions has

eigenmode solutions of the form,

U^ _ (U g )

j

exp(-inQ ) (5.5)

where A 3

is the comp l ex eigenfrequency and U 5

is

time-independent.

Su pp ose the frequency n in (5.1) is close to Qg.

Define a' by

A B

= A + J1 (5.6)

Thus equating (5.1) and (5.5)

FC 
60

Or!l r, ",I

( U a ) j exp(-inA') OF

 Y_

A,(j,n) exp(i f m l d^ + A 2 (j,n) exp^i J

s = d^ (5.7)

0 l 0

so A l (j,n) _ ?xp(-iaa ') A l (j,0) (5.8)

A 2 (],n) exp(-ini2') A 2 (j,0) (5.9)

and

j 7

( U 3 ) j A,(j,0) exp

l

i r m l dE

1

+ A 2 (j,0) exp i Jr m 2 d^

l 0 0

(5.10)

Substituting (5.8),(5.9) into (5.2) gives

a ^ +

aj [A m (j,0)) = I r 

iS2 '

l Am(j,o) m=1,2 (5.11)

l g )m

which: can be integrated to give

J

r

A m (J,0) Am(0,0) exp J Ir E + iil dj (5.12)

0 l g m J

The boundary conditions (5.3),(5.4) then become in

matrix form

B A 0 (5.13)

(A1(0,0)

where A I (5.14)

(A2(0,0)^

and

B: 1 '31 2

B = (5.15a)

B 2: 8:2

B $ I (5.15b)
61.

ORIGINAL PAGE ',Er

OF POOR QUALM r V

B 12 = b l (A,m 2 ) (5.15c)

J J

B21 = b 2 (S2,m 1 ) exp i m d' ex

E + in' f 5. 1 ] P f ( r dj ( 15d)

0 0

111 g 1

J J

B 22 b2 02) exp(i f 9 2 dj exp f I E r in 

J dj

(5.15e)

00 l g J2

A non-zero solution A of (5.13) exists if, and only

if,

det B = 0 (5.16)

Hence

J

exp ^ia rI^ [ ] 1 - [ r g ((0 2 j) -1 dj

J0

F

bz(a,m2) b,( s2,x t ) r E /E

i

b 2 ('I bl(:2,92) exp l 92 - 9 1 dj exp I(r

1 

-lr 

I dj l 0 ) 0 9l 2 g l

(5.17)

the right hand side of (5.17) can be expressed as

its magnitude

r J

b 2 (Q, m 2 ) b l (n,9 ) p I (E l 

2 

-^E 1 7

f

i b 2 ( `^,s l ) b l (:2,m2) I ex Re Fr J lr I d'

l 0 g g' 1

multiplied by a phase factor exp(if) where

J ,!

arg b2 (2 D l ) b l (il,y 2 ) + f 92 -91 d0 ^ lm 

l^lrgJ: f rgl2dj

l(51.18)

If 0 1 , 9 2 are chosen so that IF/2w is an integer then

(5.17) reduces to
C 

62

n - - N 1n I b2 (n ,m1 ) bl (n,m2 ) I Re f J 1

(=g12-(Egli 

dj l 1 J i

(5.19)

where as defined in chapter 4,

J

N J

( r g (m l ,j) ) 1 - ( r g ( m 2 ,j dj (4.19)

0

The stability criterion is

Im(n a ) - Im(n + p') < 0 (5.20)

Thus the frequency n resulting from a nc • -mal Von

Neumann analysis is corrected by an amount n' due to

boundary conditions and variable coefficients. This

approach, using n as an initial approximation to R 3 the

actual eigenfrequency, is valid provided the asymptotic

errors are small compared to n'.

The asymptotic error is 0( L

C

-2 , T A

-2 ) - 0( L C

-2 , n' 2 )

Now N - O Jr

9

so if r << J then N >> 1

9

and 'hence n' << 1 except near frequencies for which

b 2 : n,m 2 ) b 1 (n ' s )

b t (n, 9 i ) b l ( n ,m!

is zero, or infinite, which usually occurs at n -0. However

these frequencies are heavily damped by the boundary

conditions and so an accurate estimate of their

eigenfrequencies is not essential. This method gives

accurate asymptotic values near the critical frequencies

which are least damped and which therefore determine the

v

k
63

overall spectral radius of the scheme.

C+ 
W7^

64

OR'-GINA! PAGE iS

5 . 2 Examp le OF POOR QUALITY

This example is the Backward Euler method applied to

the model convective problem with constant CFL number r and

space extrapolation at the downstream boundary.

The dispersion relation is

i exp(in) - 1 = it sin(l) (3.63)

so if 91 is real, m 2 = n 1

is also real

b l (A,9) 1 (3.82)

b 2 (0,(D) 2i exp(-im/2) sin(D/2)

(3.88)

b2(A,m2) b l (R,m l ) sin(22/2)

^

,

so =exp [ -1(m 2 -^D 1 )/2) i b2 (+2, 9 1 ) b 1 ( 12 ,m 2

)

sin(m 1 /2 )

cos (0, /2)

= exp[-i( m 2 - 9 1 )/21 sin( ml/2)

t: 'xp[- i($2 -1 1 )/2) cot(1 /2) (5.21)

Since ar = 0 , E=0

Hence Y j -( 9 2 - 41)/2 + J($ 2 - ml)

_ (J - 1/2) (?r - 22 1 ) (5.22)

r,

k 1/2r - n, where n is an integer, implies

E

41 = ?r/2 - 21rn / (2J-1) (5.23)

'"he group CFL nuc•bers are

r g (4 1 ) r cos (4 1 ) exp(

 -

i^j

= r cos((P / ( 1 + -r sin (4 1 ) )

ti
(

65

ORIGINAL PAGv E iS

OF POOR 01JALiTY

r cos(ml )(1 - it sin(m l )) / (1 + r 2 sln 2 (m l )) (5.24)

r g (m 2 ) _ - r cos (a l )(1 - it sin(a i ))/(1 + r 2 sin 2 (o,)) (5.25)

Hence,

t it cos(m ) ( 1 - it sin(m )) rr ll 1 i lnlcot(m1i2)J (5.26)

2J ( 1 + r 2 sin e (o 1 ) 111

exp(-2 Im(A)) _ exp(i:l) I2

( _ 1 + it sin (o, /2 ) 12

i

= 1 + r 2 sin e (D 1 /2 ) ( 5 . 27 )

so Im(A) _ - 2 In ( 1 + r 2 sin 2 ( m l / 2) ) (5.28)

Hence

c

Im( Q + a' ) _ - In [ 1+ r 2 sin 2 ( m l / 2) ) -

r cos( 0, /2 )

In ( cot(a l /2) l (5.29)

2J(1 + r 2 sir. 2 (a l / 2) 

Define the decay rate a to be

o - Im( R + II' ) (5.30)

For small m << 1

C 2 r2a2 + 2J ln(2 /a)

(5.31)

4

Jdm r2a - 2J a -` (5.32)

li

c

Thus do = 0 at $ (2rJ)- (5.33) d 

so min a(a) 4J + 4J 1n18rJ) (5.34)

The spectral radius a is
J

., fY

a max exp(-a)

ma.t ( 1 -a )

4J

+ In(8rJ) J (5.35)

46.2 continues this example provin g rigorously that

every eigenmode is stable, and deriving an asymptotic

expression for the decay rate of the eicenmodes near T,—O,

showing as expected that the decay rate is greater than the

minimum decay rate obtained above and so the above analysis

is valid in calculating the spectral radius.

97
67 C..-

^

QUA'- ► fY

Ur Q

6. Asymptotic Stability and Convergence Analysis 11

6.1 VFeory

This chapter continiies the analysis of chapter 3

without anv additional assumptions or approximations. The

eigenmodes of a linear system of finite difference equations

with time-independent coefficients and boundary conditions

vary exponentially with t±mP Po a general eigenmode czn be

written as

M ] 1

Un exr( -inQ ) Am( j ) exp i Jr m m d& I (5.1 )

J u

t m^ 1 0 J

Note that the amplitudes A are independent_ of n.

m

The time evolution of the ei g enmode is contained solely in

the term exp(-inQ). The complex amplitudes A Pach satisfy

m

their asymptotic amplitude equation,

rAm -

e A m (6.2

g 2

)

j

J

^, r

^_ so A (J) - A (0) exp

J 

i E l dj (6.3)

m m r

0 q m

ldenotes =valuated at sm,j e

l

r,; ) m

4

Now that A (J) is related to A (0) the M asymptotic M m

* boundary condions can be written in tensor summation form as

B km A m - 0 (6.4)

where A r A (0) (6.5)

m m

and

I68

b k (a,s m ) if the k th boundary c o ndition is applied

at j-0

k m

Jf r

1

b k (1,s m ) exp ism dj I (6.6)

0 l * J m J

if the k th boundary condition is applied

at j-J

J

The term b k ('I,* m ) exp f is m dj comes from (2.77)

0

J

and the term exp f "E 

g Jm dj

is the factor relating A (J)

0

lr J m

to A (0) in (6.3).

m

A non-zero solution to ( 6..i) exists if, and only if,

det(B) n 0 (6.7)

Since all the elements of B are implicit functions

of :2 this is the equation which determines the

eigenfrequenc:_es .

If the coefficients are constant all of the A are

m

constant and so there are no asvm p totic errors. The theory

is then exact and is identical to the P-stability analysis

of Beam, warming and Yee [5]. If the coefficients are

variable the asymptotic error is of order 0(LC-`,J-1LC -1)•

The 0(I. 0 -3 ) comes from neglecting second derivatives of Am

and s in the asymptotic amplitude equations. The 0(J-'LC -1)

comes from neglecting the first derivatives of A

m 

and t 

m 

in

the asymptotic boundary conditions.

For all but the , ery simplest problems it will be
n..

L.

69

OF FCj'n ^`^ r^^ITY

4.mpossible to solve (6.7) analytically to obtain the

eiaenfrequencies. Three possible numerical approaches are

outlined below.

6.1.1 Iterative Solution

Suppose a n is a good a p proximation to an

eigenfreauency a and

a

n+1

an aan (6.8)

is to be a better approximation.

The terms in the definition of Bkm which change most

J

rapidly with variations in A Pre exp is

m 

dj since these

0

are oscillatory functions because 9

m 

usually has a real

comooneUt of order 0(1).

as _ 3a)-1

aim (70j m

q M-,

Hence 0m+1 1n + (r g)

m l

 An (6.10)

S+ihscriot m means evaluated at Om,j

S ,iperscript n means evaluated at an

Thus

J J

rJ

1

exp , is dj exp f io n dj expl iAQ (r )-1 dj (6.11 )

0

m

0

m

l0

Ir

9

Define
70

b k ( a ,0m) if the 

kth 

boundary condition

.s a p plied at j-0

r.J J

l

r -1 I km bk (11 , 9m) exp 1 r r 1 n + ism di exp ian f rr m dj

O l g m 0

g

if the k th boundary condition io at j-J (6.12)

Then choose AQ to be the smallest root of

det(B n ) - 0 (6.13)

The method in chapter 5 performs one step in this

iterative procedure. As explained in chapter 5 the

asymptotic error remaining after a correction an is

0( a;l - ^ , L C

-Z ) . For cons- ant coefficients this gives

quadratic convergence to the true eigenfrequency. For

variable coef f icients it will converge to a frequency which

differs From the true eigenfrequency by an asymptotic error

of order 0(r.O - ' ) .

The solution procedure will fail as described in

chapter. 5 near freq , lenc?es for which b k (0,0 m ) is zero since

the fractional variation in b k (A,A m ) is comparable to, or

JJr may exceed, the fractional variation in exp ism dj)

0
71

5.1.'. Newton-Raphson Solution

The Newton-Raphson solution procedure is

^n^1 ^n - det[B(.2n)] _

(6.14)

aZ det [B( ^l n ) ]

Although this seems straightfcrward there will be

problems in practice because det(B) = O has many roots and so

3, det[B(fl n )] has many zeros near which the Newton-Raphson

iterative procedure is badly behaved.

6.1.3 Stabilitv Dom a in Method

If z is defined as

z = exp(iO) (6.15)

then the condition for stability is that det(B) has no zeros

in Izj < 1. In very simple cases this criterion can be

tested analytically (an example is given in sb.2). In more

complicated cases because det[B(z)] is an analytic function

of z the 'Principle of the Argument' method outlined in

appendix A.2 can be used to find whether det[B(z)] has any

zeros in IzI < 1. To find the spectral radius ^ the same

method can be used to find the largest a for which there are

no zeros in 1a1 < X.

i

f 
72

V r`. _. .^ C: t.:a e. 7

6.2 Examole OF POOR QUkUrs

The example is the Tra p ezoidal method applied to

the model convective problem with constant CFL number r and

space extrapolation at the downstream boundary.

The dispersion relation is

tan(Q/2) = Z sin(m) (3.41)

Define K = exp(is) (6.16)

z exp (-iQ ) ( 6 . 1 7 )

1 exp (iQ /2 ) - exp (-ia/ 2 ) Then tan( n/2 exp(in/2) + exp(-iA/2)

1 1 - exp(-i.12)

a —

i 1 + exp (-iQ )

1 1-z

i 1+z

sin( m -I- 23. [ exp(is) - exp(-is) l

= 2i 

(6.19)

Thus the dispersion relation becomes

1-z r

1+Z 4 (6.20)

Note that if K, is one solution then K 2 -K1-, is

zhe other solution. The condition fcr stability is than

IzI < 1 for all eigenmodes.

Let z = R exp(i9) R,a real, R > 0 (6.21)

Then 1-z = 1 - R exp(i9)

1+z 1 + R exp(i9)

_ (1 - R exp(i9); (1 + R exp(-i9)1

[1 + R exp(i9)] [1 + R exp(-i@)]

I(D 
OF PWR QUAL11*Y

U,

73

i) Stability domain for z

1-Z

ii) Stability domain for 1+z

iii) Stability domain for K

FIGURE 4. STABILITY DOMAINS - SHADED REGIONS ARE STABLE

Ec 0
74

ORIG;-INA! FY.r^-

OF FOOD: Q'JA'.; : c

1 - 2iR sin(g) - R2

1 2R cos(9) + R2

(1 - R 2 ) - 2 i R sin(@) (6.22)

[1 + R cos ( @) ] 2 + [R sin(@)]'

Hence R < 1 => Re + Z > 0

R > 1 => Re + Z < 0

Thus if r is )ositive then R I zI < 1 if, and only

i f Re ^+Z > 0 < __> P,[K-K -1 J > 0

Now let K R e.cp (i@ ) ( 6 .23)

K - K -1 - R cos(9) + iR sin(G) - R -1 cos(0) - iR sin(9)

_ (R - R -1 ) cos(9) + i ( R + R -1 ) sin (@) ( 6.24)

so I z I < 1 <=_> Re [ K -K -1 ] > 0

either R > 1 -ir/2 < @ < r/2

or R < 1 n/2 < @ < 31t/2

Figure 4 shows the stability domains for the

different variables.

The upstream Boundary condition has

b l (3.82)

The downstream space extra p olation has

b 2 1 - exp(-iA) (3.88)

Hence B =

^

exp(iJs 1 ) r 1-exp(-i*j)] exp ( iJm1)(1-exp( - i(Dl)l

1 1 l

K J ( 1 - K L) (_K _ 1)1 ( 1 + K) 1

(6.25)
75 _I

Cy

det(B) = 0 implies

K 2J-1, K - 1) - (-1) 1

(K + 1) = 0 (6.26)

r+ 1 _ (_1)J K 2J-1

(6.27)

Consider the twc ins,'. - ..•ility regions I and II marked

in figure A. iii) .

In region I IK + 11 > IK - 1 and SKI > 1

so I 

K+i 

I > I_ 1) J K 2J- 1 1

In region II I K + 11 < 11C - 11 and JKJ < 1

so

I

K

}1 I < I(- 1)JK2J -1I

K-1

Hence there are no possible solutions of (6.32) in

the two regions of instability. Every eigenmode is stable

and so the method is stable. This exam p le is a special

example taken from a more general result, proved by Beam,

'Harming and Yee [5), that the class of A-stable Beam-Warming

multistep methods with qth order space extrapolation at the

downstream boundary is stable.

Continuing this example asymptotic decay rates for

small 91 can be derived.

01 ° 0 > K 1

Let K 1 + d

E q uation (6.31) becomes

(1 5)2J -15

 - (-1) J (2 + 6) 0

Now 1

(1 + J)J = exp(x) J m 

so provided J >> 1 >> 6

(6.28)

(6.29)

(6.30)

I E O
'f

ORIGNAL P!;.._. .

OF POOR QUALITY

(1 + 6) 2J = exv(2J6) (6.31)

Hence exp(2J6) = 2 (-1) J 6 - ' (2 + 6)(1 + 6)

= 2 (-1) J

6

-1 (6.32)

if J is even this has a solution for real 6

2J5 = ln(2/6) (6.33)

so o =

2JJ

In( 2/6 )

1 4J

= 2J 

In 

In ( 2/6 )

ln( J)

In[1n^J)] + terms of order 0( ) (6.34)

2J J

Since K = exp(is,)

= 1 + im 1

(6.35;

m l

-id

i ln(,7)

_ - (5.35)

2J

Linearisinq the dispersion relation about m-u gives

S2 r (P

it In(J) (6-37)

- 2J

The asymptotic decay rate is

r 'n(J)

a (6.38)

2J

If J is odd the smallest o solutions are

6 = r In( J) + Sri (6.39)

2J 2J

and the decay rate is the same_

F ^

A ^^
77 I C

WE I

Assorted Examples and Further Developments

7.1 Ins*_ability of Backward Euler with Soacetime Extrano'_ation

Consider the Backward Euler method applied to the

model convective problem with constant CFL number r on

domain 0 < , < J with space-time extrapolation at the

downstream boundary and J >> r >> 1.

The dispersion relation is

exp(i:1) - 1 - it sin(m) (3.63)

The upstream boundary has

b l = 1 (3.82)

The downstream boundary with space-time

extra p olation has

b, = 1 - exp(i(n-(P)] (3.94)

The ei g enfreauencv e q uation, det(B)=0 reduces to

exp(iJO 1 )(1 - ex j.(i(n- 0 1 )J} = exp(iJ$ 2 )(1 - exp(i(:1-0=)]}

(7.1)

Because m 2 - n-T, this can be written as

exp(21.Js j ) ! - exp(i(n-s, ) ) } - (-1 ) J ( 1 + exp(i(n+m l )J } = 0

(7.2)

This and the dispersion relation form two equations

in the two unknowns Z and 9 1 . Considering only the case in

which J is eve one eigenmode is given by,

Il = n + n n < n (7.3)

so it sin(0, _ -2 -in + 0(ni2) (7.4)

o 
78

m i

1 - — + H.O.T. (7.5)

r r

Substituting into (7.3) with V =0 gives

exp(-4J/r) (1 + exp(2/r)} - 1 + exp(-2/r) - 0 + H.O.T. (7.6)

Expanding the exponentials unde- the approximation

r >> 1 this reduces to,

J4 ln(r) + H.O.T. (7.7)

For a particular even valut of J the value of r

satisfying (7.7) makes this eigenmode asymptotically

neutrally stable. To find whether increasing r makes it

unstable or not (7.2) is differentiated by r.

2iJ dr exp(2iJm 1 ) {i - exp[i(a - ml)]} 

+ exp(2iJm 1 ) {i dz 1 - i d=} exp(i(Q

(i dml + i^} exp(i(:1 + 9 1 )^ - 0 (7.8) it dr

This is evaluated at V-0, J a ln(r) so

ln(r) dml - 1 { dml - dfl } + { dml + do } = 0 +H.O.T. (7.9) dr r dr dr dr dr

Hence d= 1..^.(r) dzl + H.O.T.

2 i ln(r)

_ - + H.O.T. (7.10)

- 2

As r increases from the neutrally stable value,

Im(A) becomes negative so the eigenmode becomes unstable.

T'n us this eigenmode is stable only if

r ln(r) < J

4

In numerical experiments it is found that this is

Q
79

I

the least stable eigenmode and so this condition is both

necessary and sufficient. The condition is asymptotically

equal to an exact stability condition d-rived by Beam,

Warming and Yee ,51.

f

r.

S

(1.1.1
=i

80

Uki6INNL PNGC M

OF POOR QUALITY

7.2 Optimum CFL Number For Trapezoidal Method

Consider the Trapezoidal methcd applied to the

model convective problem with constant CFL number r on

domain 0 < j < J with s p ace extrapolation at the downstream

boundary.

The dis p ersion relation is

tan(a /2) = r sin(m) (3.52) 2

Th o

upstream boundary condition has

1

b 1

(3.82)

The downstream boundary has

b 2 = 2i exp(-im/2) sin($/2) (3.88)

PM

} Following the stability analysis of chapter 5, E=0

since r is constant, and the group CFL numbers are given by,

r = r cos(s,) cos 2 (a/2) (3.47)

g i

rr cos(s 2 ) cos2(3/2)

9 ► 

- r cos(m l ) cos 2 (2/2) (7.11)

L'

J

b 2 (A,s 2 ) b l ("', j ) r

so - ara b (A,a 2 ) b 1 (^,a= ) +

m2 - ^, dj

_. J0

+ J

0 g )2 \

g ^I

(J-1/2) ( $ 2 - m l )

(J-1/2) ( r - 2m 1

) (7.12)

1/21r = n whe y s n is an integer implies

1

U
8i

o; -'ALITY

?f 2rn

2 2J- 1

1n ! b .m, ) b l (n,m i )

v

I

b 2 (Cl,m 2 ) bl

i

N In[ sin 2) / sin((0 2) J

1 In[ cot($,/2) J

(7.13;

J

Re f 1zg12- ( =^)1 dj

(7.14)

J

where N

Jr [rg(s,,j)J-` - [r g (m 2 .j)1 dj (5.20)

0

2J

(7.15)

r cos($,) cos2 (n/2)

The decay rate a was iefinid as

a - Im( 3 + if' ) ( 5 .30 )

a 2J cos(m l ) cos 2 (A/2) ln[cot(ml/2)]

r cos (0 In cat( t,/2)

2J ( 1 + tan 2 (11/2) J

2r cos (9 1 ) In[ cot(t,/2) J (7.16)

J [ 4 r2 sin' (9 1 ) j

Figure 5 shows the decay rate J as a function or r

for various values of t.

When a system of _f inite difference equations is

solved us-, ng a time-independent approach to a steady state

sol , .ition, usually the initial conditions and the final

solution are smooth so the initial error is smooth i.e. the

error is predominantly in the low wavenumber, long

wavelength eigenm..odes. For this problem suppose that the

6c,
N

I

i

In (a^

I ^

= 

1

167

_ 2

16 7T

167T

677

6^

F 0

N

i

m

I

82

161T

-2. 0. 2. a. In i

r rl G.

FIGURE S. LOG-LOG PLOT OF DECAY RATE AGAINST CFL NO.

^f

L qftw o 
2

r - ain($o)

(7.19)

83

PA.

OF r: J9Z QUALA'' .

initial error is in the wavenumber range 0 < $ < $o where $o

is a constant, $o < n/2

3 2r ( sin($) 2r2sin($)cos2($) J ln[cot($/2)]

3$ J 4 + r 2 sin e (^ ) + (4 + r 2 sin e ($) l 2 )

cos($)

< 0 (7.17)

sin($1'4 + r2Sin2($)I

so over the range given above a has a minimum at $,. To

°aximise the overall rate of convergence the CFL number r

can be changed by altering St.

ea - cos(s) ( 4 - r 2 sin 2 ((D) l ln(cot($ /2)j

(7.18)

3r 2J [ 4 - r 

2 s i n 2 ( 0 ) 1

> 0 for 0 < r < 2/sin(o)

for 2/sin(o) < r

so a(m o ) is maximised by choosing

The smoother the initial error the lower the valt•.e

of $o and the higher the optimum CFL number. If the initial 

error is not at all smooth with Zo ap p roaching r/2, the

i

f

optimum CFL number drops Sown towards 2.

i

LE
84 ORIGINAL PAGE M

OF POOR QUALITY

7.3 Discontinuous --FL Number In Trapezoidal Method

Consider the following Trapezoidal method applied to

the model convective problem.

r. r

[ a t + —=-2^ t v x + —2 u AX; U

n+ 2 = 0 (7.20)

Suppose that r is discontinuous

r j < 0

r , j > 0 t

On *..he two sides of the discon—nuity a general

solution of frequency n can be written as,

A,(j,n)exp[i(jV,-nn); + A2 (j,n)exp[i(js 2 -riI)j

j < 0

Ui (7.21)

A,(j,n)exp[i(js,- na)] + A,(3,n)exp[i(jm,-nQ)] j ^> 0

s; and s 2 satisfy the dispersion relation

r

tan(n/2) 2' sin(s) (3.52)

with $2 = 7r — al

and 0 < Re($, ) < x/2

s, and v, satisfy the dis p ersion relation

r

tan(Q/2) 2+ sin(m) (3.52)

with 01, - a - $,

and 0 < Re(s, ) < >r/2
85

ORIGINAL PAGS---' (:i

OF POUR Q'JALMI

The amplitudes A 1 , A 2 are related to A,, A 4 at j-0

by two equations. The first comes from the requirement that

the two expressions for U are equivalent at j=0.

Thus Al + A2 = As + A, ( 7 .22 )

The second requirement is that

d t + 2 r. ut7x + 2 r + 4 6 x ] U0 + ' 0 (7.23)

which implies

( 1

d

+ r μ 0 ] Un+2

- ( d + r a A ] Un+2 (7.24) t - t x 0 t + t x 0

The left hand side involves U at j=0,-1 so the

J

first expression for U n is used. The right hand side

involves U at j = 0,1 so the second expression is used.

c J

Neglecting derivatives of the amplitudes the resulting

equation is,

A 1 { -2i sin(:2/2) + r - cos(n/2)[1-exp(-is,)] }

+ A 2 11 sin(Q/2) + r - cos(A/2)^1-exp(-is 2 )] } 

i = - A,! -2i sin(,'1/2) + r+cos(n/2)[1-exp(-i©,)] } (7.25)

- A,{ -2i sin(A/2) + r + cos(a/2)[1-exp(-is,)] }

Substituting for t_ntQ/2) usin g the dispersion

relations and using

exp(im) - cos(s) + i sin(s) (7.26)

i

this reduces to

r - { A 1 (1-cos(a l )] + A 2 [1-cos(s i )] }

= r+ { A,[1-cos(s 2 )] + A,[1-cos(s.)I } (7.27)

s^
' 86 ORIGINAL

OF POOR QUALI iY

In the general problem the two equations relate the

am p litudes on either side of the discontinuity at j s 0. Now

consider the particular problem in which m l , m 2 , a,, y, are

all real and A, is zero. This is the situation when a

wavepacket with wavenumber m l , and positive group CFL number

travels from j<0 to the interface at j=0 producing a

transmitted wave p acket with wavenumber r 2 and positive group

CFL number and a reflected wavepacket with wavenumber 02 and

negative grou p

CFL number. The reflection coefficient R and

transmisj: on coefficient T are defined by,

r = ;} / A l ( 7.28 )

E T A3 / A, (7.29)

Since A,=O

A l + A 2 = A, ( 7 .30 )

so r- (Ai (1 -cos (m 1 ) ] + A= [ 1 - cos (a 2 ) ] } - r + [ A 1 +A 2 ] ( 1-cos ( a 2 ) ]

(7.31 ) 

r - [1- cos(m l )] - ry(1 -cos(m2)]

Hence ^ _ - ( 7.32 )

r-[1-COS((02)] - r+[1-cos(03)]

aad 1 = 1 + R 

r [1- COS(0 2 )] - r ;1- cos(sj)]

3 r - [1 - COS V 2 )] - r+[1-COS(02)] (7.33)

1r

A'
87

S. Assorted Loose Ends OF POUR 4Jf;tlfY

8.1 Degeneracy

So far in this paper it has been implicitly assumed

that no two ei( • enmodes have the same frequency R and

wavenumber I. This section cons_uers the degenerate case

in which this happens. An example is the Trapezoidal

method. The dispersion relation is

r tan,'Q/2) z sin($) (3.41)

• which has wavenumber solutions m and s 1T These are

1 2 1

ident. :al when

r

which occurs at the degenerate fraglzency A' given by,

tan(A/2) = 2 (8.2)

atl

Note that r ,

q ^s j const

= r cos(a' ) cos y (01/2)

= 0 (8.3)

This is characteristic of degeneracy because in the

} neighbourhool of the double zero

S2 - S2 ^ a( m - m ) a a = constant ( 8 . 4 )

an

so r —

g ax

2a( $

h , 0 at 9 D' (8.5)

For A ¢ 0 1 the general solution for constant r is

0
88

ORiG(yf.'►. PAC r!F

OF POOR QUAL114

Un

J , 

A,exp(i(j0,-na)] + A 2 exp(i(jm 2 -na)] (8.6)

Consider the limit as R approaches Q'

_ l

Hence A l exp [ijm, ] + A 2 [ 2 ]

- exp(ij(0 { A, exp (ij(m,- m')] + A 2 exp (ij(0 2 - m')] }

= exp[ijm'] { A,exp(tj(0,- m')] + A 2 exp(ij($'- $,)] }

2 exp(i 3 o'] { (A,+A:)cos(j(m_-x')] + i(A,- A 2 )sir j(9,—m')] }

= exp [ i j I { A, + A, j }8.9)(

where A, - A, + A 2

(e . 1 0 )

and A, - i(A, - A 2 ) fm, - V') (8.11)

If A, and A, are now considered to be constants in

the limiting process the general form of the degenerate

eigenfunction is,

UnJ , (A3 + A .J ) exp(i(jo'-nil')] (8.12)

Another way of deriving this result is through the

asymptotic am p litude equation. with the amplitude A

tine-independent, r constant, and r

9 

equal to zero the

largest_ term is the second order derivative of A so,

a2 0 (8.13) aj2

which implies A = A, + A, j (8.14)

In all the examples I have analysed I have not yet

is

C 
O

89

found a degenerate eigenfunction which satisfies the

asymptotic boundary conditions and so is a degenerate

eigenmode. The degenerate frequency does however satisfy

the determinant condition

det(B) - 0 (6.7)

because two of the wavenumbers are equal so their columns

are identical and B is singular. This is not a problem in

the stability analysis because the separation of

eigenfrequencies, using the result from chapter 5 , is

approximately,

J -1

21r N 2ir I ^r (m ,j)] - ' - it (^ 1

g

1

g

i .j)] -

dj

0

0 as 2 • 'I'

so there is usually a true eigenfrequency which differs

from the degenerate frequency by less than the overall

asymptotic error.

U
90

8.2 Eigenmodes and Eigenfrequencies

In chapters 5 and 6 it was stated that the

eigenmodes of a linear system of finite difference

equations with time-independent coefficients and boundary

conditions vary exponentially with time. This section

outlines the proof.

Let Un zn V. (8.15)

J J

If the domain is 0 < j < J then there must be J+1

finite difference equations so V. satisfy

7

C V - 0 (8.16)

where V is the J+1 vector of V. elements and C is a (J+1)2

J

matrix whose elements are polynomials of z.

For there to be a non-zero solution requires

det(^:) - 0 (8.17)

This is the equation that determines the eigenvalues

z of the eigenmodes. Apart from the problem of possible

degeneracy the on.y remaining difficulty i_s to show that the

number of eigenmodes equals the number of independent

initial conditions needed to start a numerical solution.

For the Backward Euler method with space extrapolation at

the downstream boundary there are J-1 independent initial

conditions since

U0 = 0 (8.16)

and (JJ , U 0

(8.19)

n:

1

^r
91

L ^•• ^I \J 1' w, w•f T •• '• ^,^

1

-r, z/2 z- 1 r, z/2

-r 2 z,'2 Z-1 r2 z/2

C -r, z/2 z-1 r, z/2

r -rJ-1z/2

z-1 rJ-1z/2

1 

1

(8.20)

so det(C) is a polynomial in z of degree J-1 giving J-1

eigenmodes. Thus any solution can be expressed as a sum of

eigennodes which vary exponentially with time. M.

I 
92

8.3 Other Asymptotic Approaches

One approach which can be used when the variation

in coefficients is extremely small is to let,

n U. A(j,ni exp(i(jm- nil)] (8.21)

3

where n and a are both constant and satisfy

a 0 (S2,m.j) - 0 (8.22}

at some point j.

The asymptotic amplitude equation is the-;

a0A + al 

an 4 a2 34j - 0

(8.23)

where a 0 , a l and a = are defined and calculated as before.

L min A / c A

A j

j

0 (a 2 / a p ) (8.24)

so the fractional error using this method is 0 (a0/a2)2.

If the variations in the coefficients are small this is fine

but: if the variations are (I ( 1 ) the fractional error is O(l)

i.e. the method fails to give accurate asymptotic

approximatiot_s. The asymptotic approach used in this paper

allows total variations in the coefficients of 0(1) and onl.r

requires that the length scale of *_hose variations is much

greater than 1.

Another approach is to set

n

U, A(j,n) exp,i`C(j,n)l

M

with — 9

r ^^

(8.25)

(8.26)
i

93

and

an

A (8.27)

and A, :1, ; are all real and slowly varying. This approacL•

is used extensively in the analysis of water waves and other

partial differential equations with dispersion and very

little dissipation. This approach applied to finite

difference equations would give poor results because

dissipation over one time step is 0(1) so if n is real A

reflects this dissipation and so T A M 0(1). The method

presented in this paper is able to use constant complex A

rather than variab'.e real A as in the above method because

the eigenmodes have constant complex eigenfrequencies

provided the finite difference equations are

time-independent.

F

I
94

S. waveoacket Test P-ogram

9.1 Program Description

The program solves the model ccnvective problem,

8u '

c

au

- 0 (9.1)

3t ex

using a choice of Box or Trapezoidal methods on the domain,

0 < j < 200 (9.2)

and time step range,

0 < n < 400 (9.3)

The CFL ++number

c u t 

r.

= ^x

(9.4)

is specified by the user at j-0,200 and the program

interpolates fcr intermediate values by fitting an

exponential curve through r 0

and

x200'

rr 1

r( j) - r 0 expl20 In( r200/r0)1 (9.5)

Methods 1-3 are different Trapezoidal methods which

are identical if r is constant.. For these methods the

program offers a choirs ^jf space extrapolation, space-time

extrapolation or box condition as the downstream boundary

condition. Method 4 is the Box merhod. The waveoacker.

theory for each of these methods is derived in the next

section.

:he upst-eam boundary condition is,

UO = 0 (9.5)

The initial conditions are given by,

__ Q

'.'
uwIGINAL Fk.GE IS

OF POOR QUALITY 95

n

U 

= A(j,n) exp(i'f(j,n)) (9.7)

where A(j,n) cxp'L- ( j-100)=/200J (9.8)

ar.d Y(j,n) ml (ft,j) (9.9)

k-0

This produces a wavepacket al,out 40 :mesh points

long. The user specifies 7 by specifying the critical CFL

number which is defined as,

r - 2 tar.(a/2) (9.10) crit

so n - 2 tan-` ( r Crl t /2) (9.11)

The significance of the critical CFL number is that

w h E n ,

r - r (9.12) crlt

the group CFL number for the Trapezoidal methods is zero.

Wave propagation can only occur for r > r .

Grit

At each time step the prugram calculates

'experimental' values for the position X(n) and energy =(n)

of the wavepackat. If U was continuous X and F would be

defined as,

x2r00

E(n) - 1 JU(x,t )1 1 dx

/ n

x0

200

fU( j ,n) = dj di

0

(9.13)

^ E o
ni,

96

ORI(;NAL Pi.C►: Ian

OF POOR QUALITY

x2r00

and X(n ) E(n) J IU(x,tn x dx

X

0

200

jU(j,n)I2 x ) dX d (9.14)

E(n) J

0

so since U is discrete: X(n) and E(n) are defined as,

0

E(n) IUjli c j (9.15)

j-0

200

X(n) E(n)

!Un

x a 

(9.16)

]M0

where a j 2(xj+l -xj-1) (9.17)

In the program c -s taken to be _onstant so

variations in r are due to variations in Ax.

c A t Ax - (9.18)

r

so x3+1 xj-1 ° (rj

+,,

i

A physical domain 0<x<1 is used so,

200

ff/ xjrj-1I '7L'i (rj_,) (y.19

1

)

After completing the 400 time steps the program

calculates predicted values for X(n) and E(n) using the

wavepacket equations "erived in the next section with the

experimental values at n-1 as initial r.or-lit.ons.

n %

1
97

t<<' - tat. `.'.G^ t5

OF POOR QUALITY

9.2 Theory

Method 1 is a Trapezoidal method which -s second

order accurate in Ax only when Ax is constant.

[ 6 t

+ 2 μ t 6^ x J III + 0 (9.20)

2 c A t (9.21) r =

xj+l-xj-1

i

( Following the analyses of chapters 3 and 4,

r

i ao = -2i sin(9 / 2) + it cos(Q /2) sin($) (9.22)

a. - cos(.,1/2) + r sin ( A/2) sin(o) (9.23) 2

} a2 = r cos(A/2) cos(s) (9.24)

c_

r

a, _ - 2 cos ( a/2) sin ( V) (9.25)

The dis p ersion relation is

tan(A / 2) a 2 sin($) (9.26)

so a, = 0 (9.27)

and a l - cos(a/2) + sin(n/2) tan(Q/2)

_ ( cos' (n/2) + sin 2 (,1/2) ) / cos(2/2)

sec .1/2 (9.28)

r = a 2 / a l (3.26)

9

r cos(a) cos' ( 1-1/2) 19.29) _.)

do

11

)

'i,s const

a1 (4.v)

e'

o
98 ORiGf^If". Pry is 1;4'

OF POC R QUALITY

_ ar sin

( m) cos 2 ( 12/2) (9.30)

= aj

E = -ia,ll a const 

/ a,"

2 - a o / a i

(3.30)

_ - 1 ar cos

2 ( A/2) sin ( m) tan(m) (9.31)

2 13

Since a is proportional to r-1,

1 as 1 ar

a a^ — r a^

(9.32)

so

a

( r a) r a— ( r /r)

a aj 4 33

g

- r sin.(W) cos 2(A/2) amaj

_ - r sin(a) cos 2 (A/2) i 0 / a2

j

C'

ar 

co g2 ( A/2) sin(o) tan ( m) (9.33)

aj

Hence the ecuations of motion for the wavepacket for

real ;1 and a are 

liz

dr cos(m) cos2(0/2)

(9.34

i = - ar

o a

sin(o) cos2(0/2)

(9.30

do

d 1n(E)

z

j 0

(9.35)

Method 2 is also a Tra p ezoidal method which is

second order accurate in dx only when Ax is constant.

+ r1 ( 

d

+ 7 μ9 + 2-

μ t x7 U.

2 = 0 (9.36)

tX

t _

1.-.
C

99

r. =

cat (9.37)

J+} xJ+1'xJ

a, _ -2i sin(n/2) +Z cos(n/2) I r j+# [exp(im)-1] +

r,-,[1-exp(-im)])

-2i sin(n/2) +2 cos(n/2) Pr ^ + 2 

ar 

)[exp(im)-1]t

ll 7

1

(r - 1 

3r 

)[1-exp(-i(V)] I + H.O.T.

3 2 3 j J

_ -2i sin(0/2) + it cos(n/2) sin(4) +

1 3r cos(n/2)[cos(s)-1! + H.O.T. (9.38) z aI

The H.O.T. are neglected because they are of the

same order of magnitude as other terms already neglected in

the derivation of the asymptotic amplitude equation. To

this same level of as y mptotic accuracy a,, a 2 and a, are

exactly the same as in the analysis of method 1.

The dispersion relation rema'_ns,

tan(n/2) sin($) (9.26)

F so now a, 2 cos(n/2) [cos(m)-1' (9.39)

3i

Afte:- some algebra it follows that the equations of

s,

:notion are,

r cos((P) cos 2 (0/2) (9.40) do

E;

a ' dt

an

1 3z sin(a) cos 2 (n /2)

- z a,

(9.41)

1n(E) =

ar

[1-cos(o)] cos 2 (n/2) (9.42) do 7

method 2 is a Trapezoidal method which is second
100

ORIGMIAL F;-.

OF POOR Q1%_1';%L1 tY

order accurate in Ax for cases with smoothly varying Ax.

(r j +j)

2

(

r 7 - ,) 2

,

( 6 t + r. +r utAx + = +r

ut7

n+

x Jj 0 .4

J + ^! J- J+ j-2

The dis p ersion relation is the same as for methods

1 and 2 so after some algebra,

ao 2.2 cos(A/2) (cos(;)-1; (9.44)

and the equations of motion are,

d1 = r cos(s) cos 2 (A/2) (9.45)

^1 n

dm

3r sin(m) cos 2 (0 /2) (9.46) d 3 - 3 

do In(E) = 2 

ar 

(1- cos(m)1 cos 2 (0 /2) (9.47)

J

When the :raveaackets of methods 1-3 reach the

downstream boudary at j=200 they are reflected into

backward travelling wavepacke*_s. The energy E 2 and the

wavenumber S 2 of the reflected wavepacket are related to 

the energy E 1 and wavenumber $, of the inci3ent wavepacket

1-s the equations,

^V2 - ?r - a 1

(9.48)

ln(E 2 ) = In( E,) + 2 lnj R1 1 (9.49)

where R is the amplitude reflection coefficient.

For s p ace extra p olation (see 43.5.2)

sin( s,/ 2 )

R J I = 9in(m2/2)

tan($,/2) (9.50)

0

R 
101 , r

OF POOR QUAL11"Y

For space-time extrapolation (see 43.5.3)

sin(($,-QU21 IRJ I a

sin(((P2-R)/21

For the box boundary condition (see 43.5.4)

cos($ 1 /2) tan(n/2) - r sin($,/2)

IR J I cos($./2) tan(a/2) - r sin($,/2)

tan' ((P1/2)

after substituting for tan(a/2) using the dispersion

relation and replacing $= byr

(9.51)

(9.52)

The reflection relations at the upstream boundary

are,

$ I - 'r - $ 2 (9.53)

In ( E 1 ) - In( E 2 ) ( 9 .54 )

since R 0 = -1 (see 33.5.1)

Method 4 is the Box methcd discussed in 33.3.2

( u x d t + r

i 

u t d x 1 U^^ 2 = 0 (9.55)

The dis p ersion relation is

tan(A/2) - r tan($/2) (9.56)

R

>E

and the equations of motion for the wsvepacket are,

do

. r ( 1 + tan 2 (m/2) 1 cos 2 (0/2) (9.57)

do

=

-2

;r

tan($/2) cos 2 (Q/2) (9.58) do

r_
102

do In(E) - -2 - tan 2 (0/2) cos 2 (0/2) (9.59)

Wi

v U
103

9.3 Numerical Results

9.3.1 Trapezoidal Method with variable CFI. Number

This example uses,

Method type = 2 ; one of the Trapezoidal methods

Boundary type = 1 ; space extr.-polation

r0 200

= 0.05 r - 0.2 r crit = 0.04

Figure 6 shows R(n) ani ln[i(n)] both predicted and

experimental. This example shows the movement of a

wavepr,-.ket and the change in its energy due to the variation

in the CFL number. The agreement between the predicted and

experimental values is excellent. The energy of the

analytic solution io constant so the wavepacket theory has

successfully predicted almost all of charge in the numerical

energy due to variable Ax.
ORIGI.N., ' :

OF PC t 

9 NUMER I CAL EXIDER I MENT

z WPVEPACKET THEORY

C

X

uo 

00

i

100 200 rl 300 400

104

In(E

Ln

N

i

CND

Ln

n.l

0 100 200 n 300 400

FIGURE 6. EXAMPLE 9.3.1 TRAPEZOIDAL METHOD 2
F-tc; 

105

9.3.2 Box Method with Variable CFL Number

This example uses,

Method type = 4 Box method

r 0 = 0.05

r200

0.2

rcr+t ^ 0.04

Figure 7 shows X(n) and ln[E(n)]. As in 49.3.1 the

agreement between experiment and theory is excellent.

n ;

C

n O
0cn

In(E

C^^

N

Ln

^1J

of 1CppoR No ►s

NUMER I CPL EXPER IME N T

T^^IfV

--WRVEPACKET THECRY

400

uu n 300 400

FIGURE 7, EXAMPLE 9.312 BOX METHOD
107

9.3.3 Trapezoidal Metnod with Space Extrapolation

This example uses,

Method type = 1 ; one of the Trapezoidal methods

Boundary type - 1 ; space extrapolation

r0

1.0 r 200 = 1.0 r 

crit 

W 0.3

Figure 8 shows X(n) and ln(E(n)). This example

illustrates the effect of the downstream boundary reflecting

a wavepacket with reduced energy. Because of the finite

length of the wavepacket the drop in energy is smeared and

X(n) does not quite reach 1.0 .

0
100 200 n 300 400

o^

c0

108 OR;GINAL PAGE IS

OF POOR QUALITY

NUMERICPL EXPERIME NT

WPVEPPCKET THEORY

0

X

0

0

c

n (E

l; .

0 loo 2C0 n 300 400

t^

FIGURE 8. EXAMPLE 9,3,3 TRAPEZOIDAL METHOD 1 WITH SPACE

EXTRAPOLATION

V

T 
iI .

0:

109

9.3.4 Trapezoidal Method with Spacp-time Extrapolation

This example is the same as 19.3.3 except that

Bcundary type - 2 , space-tie extrapolation

Fiqure 9 shows X n) end lr,(Ekni]. The energy of

the reflected wavepacket is less than in 19.3.3. As a

consequence the first orier terms which are neglected in

the asymptotic bounr.-ry conditions are more significant and

so the energy decrease :s more smeared and there is a

larger discre p ant:, between the experimental and predicted

values
110 ORIGINAL PACE 15

OF POOR QUALITY

=—a-^- NUMER I CAL EXPERIMENT

WAVEPACKET TH E09Y

I

0

v

n

Ln

0

c0

roc 20c n 300 4Cc

1:

0

0

in'

c

0

i

0 1c0 200 n 300 400

FIGURE 9. F YIPLE 9.3,4 TRAPEZOIDAL METHOD 1 WITH

SPACE-TIME EXTRAPOLATION

v
111

9.3.5 TraDezDidal Method with Bcx Boundary Condition

This example is the same as 39.3.3 except that,

Boundary type - 3 ; box boundary condition

Figure 10 shows X(n) and in[E(n)]. The energy drop

in this example is three times that in 49.3.3 because,

ln(E 2 ) - in(E1) + 2 ln(I RJ I ) ( 9.49)

tan(m l /2) space extrapolation (9.50)

and IRJI'

tan'(m l /2) space-time extrapolation (9.52)

Thus the box boundary condition increases the

overall convergence rate by factsr 3 with minimal extra

computational effort.
M,-

112

ORIr;INAL PAGE 1,9

OF POOR (QUALI TY

---^-e- NUMER I CAL EXPERIMENT

---INAVEPACKET THE

O

X

Ln I

O

0 I'

I ^ I f

° 0 1C0 200 n 3CC 400

0

c

InIE^

M I

I 

0 100 200 n 300 400

FIGURE 10, EXAMPLE 9,3,5 TRAPEZOIDAL METHOD 1 WITH BOX

BOUNDARY CONDITION

r

f
113

9.3.E waveoacket Outflow in Box Method

This example uses,

:Method *youe = 4 ; Box method

r 30.4

r200 = 0.4

r 0.2

czit

Figure 11 shows X(n) and 1n[E(n)]. Note that when

the wavepacket reaches the downstream boundary the

experimental value for X(n) remains near 1.0 and In(E)

decreases rapidly towards

E!:
114

O-RNCINAL F:ICC :S

OF FOUR QUALITY

$-e--= NUMER I CAL EXFEF 1 M E N T

WAVEFRCKET THEORY N

X

a

O

C ^

1c0 200 n 300 400

f n (E

100 200 n 3c0 400

FIGURE 11. EXAMPLE 9.3,6 BOX METHOD

0

0

y
115

9.3.7 Instability of Trapezoidal Method with Space-time

Extrapolation

This example uses,

Method type - 1 one of the Trapezoidal methods

Boundary.type = 2 space-time extraooiation

r 0 = 3.0 r200 = 3.0 rcrit = 2.1

Figure 12 show; X(n) and ln[E(n);. Both the

theoretical prediction and the numerical result show that

the ener g y increases every time the wavepacket reflects of`

the downstream boundary and so the numerical scheme is

unstable. The CFL stability condition for the case with

constant CFL number is obtained by considering the

amp litude reflection coefficient.

sin l -")

IRJ I ,4 (9.52)

sing 2

The dispersion relation is

tan(n/2) _ 1 sin(m) (9.26)

Now consider the two cases r<2 and r>2 .

a) r < 2

0 < s, < a/2 __> tan (n/2 ) < 1

(ir/2 - U > 0 and (a/2 - 1,) > 0

IT > I ( n/2-0 ) + Or /2-o,) I > I (,r /2-A ) - (,r /2 -0, ) I > 0

rm 2 a rml-^^ 1

__> I sing

2) I

> sin 2 

J
116

OF POOR QiJALIiY

^-^-e- N LI ME R I CRS EXPER I MENT

WRVEPRCKET THEORY

0

X

0

^o 100 200 n 3C0 4CC

' 0 100 20C n 300 40C

FIGURE 12. EXAMPLE 9.3,7 TRAPEZOIDAL METHOD 1 WITH

SPACE-TIME EXTRAPOLATION

Y

i•'I 

n iE

m

cl]

i

I c_

!n

N
117 ORIGINAL PAGE M

OF POOR QUALITY

_>> ^RJ < 1

b) r > 2

For s = n/2 tan(,1/2) > 1 =_> Q > n/2 =_> S2>s

For a = n tan(A/2) = 0 ==> a<$

As 9 2 varies from n/2 to n, 0 - 0 2varies

continuously so at some intermediate value,

'Z = m2

and hence IR J I _

The numerical scheme is thus stable if, and only if

the CFL number is less than 2. This stability condition

has previously been derived by Beam, Warming and Yee (5].

In a numerical experiment an infinite amplitude

reflection coefficient does not occur because the first

order derivatives of the amplitude which are neglected in

the asym p totic boundary conditions become significant. In

fact in all the unstable cases I have tried the acrreem-tnt

between experiment and predict-ion is poor because of the

neglected first order terms in the boundary condition and

the neglected second order terms in the amplitude equation.

The exam p le given is ono of the best. The qualitative

e_`ects of these neglected terms can be understood as

follows;

Since 

an

r 

aA

all the second order terms in the

g 7

2 I I

amplitude equation

an A

ana

a A can be expressed in

I J J

terms of a'A Thus the amplitude equation including

second order terms has the .form,

aA

an + r aA =

2A (9.61)

g 7

a 

]I

L
118 4

where v is a function of R,m and corresponds to an

artificial viscosity. The effect of this artificial

viscosity is to smear the wavepacket increasinq its length

and decreasing its maximum amplitude. This has the largest

effect on X(n) since the longer the wavepacket the further

X(n) must be from the ends of the domain. X(n) still

oscillates approximately in phase with the predictions but

the amplitude of the oscillations decreases steadily. The

effect on tie energy is much smaller.

The downstream boundary condition including firs=

order terms can be written as,

A 2 (j,n) + T 2 arA2(J,n) - R J [A 1 (J,n) T1 an i (J,n)) (9.62)

where T 1 ,T 2 are functions of Z,m. if T 1 and T. are both

small compared to T A then (9.62) is approximately e q ual to

A 2 (J,n +T 2 ) - R J A 1 (j,n +T 1 ) (9.63)

Thus the amplitude is reflected with a delay of

T2-T1 This explains the fact that in almost all the

examples in this chapter the reflected wavepacket lags

rehind the position predicted by wavepacket theory.

`d o
ORIGENAL PAGE 19

119 OF POOR QUALITY

9.3 .8 Instability of Trapezoidal Method Due to Varying M esh

This example uses,

Method type = 3 ; one of the Trapezoidal methods

Boundary type = 1 , space extrapolation

r0 1.0 r m 10.0 r = 0.8

200 grit

Figure 13 shows X(n) and ln[E(n)]. The agreement

between experiment and prediction is good for the energy

but as in 49.3.7 the agreement is poor for X(n) because of

the effect of the second order terms which are neglected in

the asymptotic amplitude equation.. The significance of

this example is that this numerical scheme is stable for

uniform meshes which give constant C F L number r but if the

mesh, and hence r, varies sufficiently as in this example

the scheme becomes unstable. This instability is best

understood by expanding the finite differ..._ equation in

computational space.

ff 2 ( 2

Ir J + 2 1r7-2 n+L

5 t + r. +r t^x + r +r

ut7x J U^ 2 - 0 (9.43)

7 + 2 J'2 J+i J'2

r. 1

so ( a t + -1 u t a 2x + a

μ*5 2 + H.O.T. ] U^ +2 s 0 (9.64)

The *term 3r 

μ t 52 corresponds to a viscous term in

J

computational space. If 

ar 

is positive it corresponds to

J

negative artificial viscosity and so causes the instability

in the above example.

1'i

O
i

c

n lE^

c

^i

c

o 

120 ORIGINAL PAC: IS

OF POOR QUALITY

±---^-^ NUMER I CPL EXPER IMENT

kPVEPACKET THEORY

c

X

^i

0

c

0

^ I II ^ ICI

0 1c0 200 n 300 400

c

q

c icc 200 n 300 400

FIGURE 13. EXAMPLE 9.3.8 TRAPEZOIDAL METHOD 3 WITH SPACE

EXTRAPOLATION

^L 
FM 

121 ; is

OF POOR QUAL I1Y

5.4 Proaram Listing

The program is written in FORTRAN 4 PLUS to be run on a

POP 11-70 with Versatec graphics subroutines.

C •••'• PROGRAM 'NAVE

DIMENSION U(0:202),V(0:202),A(0:200,3),W(0:200),XJ(0:202)

DIMENSION EX(402),EE(402),TX(402),TE(402),T(402)

EQUIVALENCE (TX(1),U(0)),(TX(204),V(0)),(TE(1),A(0,1)),

? (T(1),A(0,3),(T(202),W(0))

EXTERNAL R,J

REAL K,J

COMMON /RCONST/R0,RI,MT

COMMON /X/X(0:200)

C "'••• Input parameters

F

TYPE ♦ ,'INPUT TERMINAL TYPE'

! TYPE • ,'0 VERSATEC'

TYPE *,'3 VT100 WITH GRArHICS'

l' TYPE ','4 CHROMATICS'

TYPE *,'5 VISUAL 500'

ACCEPT *,NT

TYPE

TYPE *,'INPUT =HOD TYPE'

TYPE *,'1-3 DIFFERENT TRAPEZOIDAL METHODS - SEE NOTES'

TYPE • ,'4 BOX METHOD'

t

ACCEPT -,MT

IF (MT.EQ.4) GOTO 1

TYPE ♦ ' '

TYPE *,'INPUT DOWNSTREAM BOUNDARY TYPE'

.^ TYPE *,'? SPACE EXTRAPOLATION'

TYPE *,'2 SPACE-TIME EXTRAPOLATION'

TYPE *,'3 BOX METHOD'

ACCEPT •,MDB

1 TYPE •,' '

TYPE •, 'INPUT CFL 91TMBERS R(0),R(200),RCRIT'

ACCEPT•*,R0,RI,RC

C ***** Jmega definition explained in notes r

OMEGAs2.*ATAN (0.5*RC)

O
,

F,

122

PI-3.14159

JMA.X- 2 0 0

K-400

XJ(0)-0.

DO 2 J1-1,JMAX

OX-1./R(FLOAT(J1)-0.5)

XJ(J1)-XJ(J1-1)+DX

U(J1)-0.

2 V(J1)-0.

DO 3 J1-1,JMAX

XJ(J1)-XJ(J1)/XJ(JMAX)

3 X(J1)-KJ(J1)

C *•"' Initialise U(j,0)+iV(j,0)

PSI-0.

J2-JMAX/2

DO 4 J1-J2-40,J2+40

IF(MT.LE.3) PHI-ASIN(RC/R(FLOAT(J1)))

IF(MT.EQ.4) PHI-2.•A:AN(0.5•RC/R(FLOAT(J1)))

PSI-PSI+PHI

AMP-EXP(-(J1-J2)••2/200.)

U(,^1)-AMP*COS(PSI)

4 V(J1)=A.MP*SIN1PS;)

W(0)-0.

KOUNT-0

9 TYPE *,'NO. OF STEPS TILL NEXT PLuT OF UV

ACCEPT * , `ISTEP

IF (21ST°P.LE.0) 50,0 5

DO 5 KOUNT2-1,NSTEP

KOUNT- )COUNT+ 1

IF ( KO'JNT . GT . M ) GOTO 7

C ••••• Calculate new U+iV

CALL. METHOD(U,A,W,JMAX,R,MT,MDB)

CALL METH0D(V,A,W,JMAX,R,MT,MDB)

C " " • Calculate new X,Loq(E) of wavepacket

S-0.

SJ-0.

XJ(JMAX+1)-XJ(JMAX)

I.-
C

123

DO B ;,1-1,JMAX

P-(U(J1)••2+V(J1)••2)•(XJ(J1+1)-XJ(J1-1))/2.

S-S+P

9 SJ-SJ+P•XJ(J1)

EX (KOUNT) -SJ/S

6 EE(KOUNT)-LOG(S)

5 CALL OUTPTI(XJ ,203,'X ','U ',NT)

COTO 9

C ••••' Calculate predicted X(n),Log(E(n))

C "•'• Initial salu j(1) is passed to prediction subroutine PREI

C ••••• as TX(1)

7 ?'X(1)^J(EX(1))

TEM-EE(1 )

CALL PRED(T,'17X,TE,JMAX,M,RC,R,MT,MDH)

C "•" Plot results

12 TYPE ','PLOT E,X?'

ACCEPT 1000,C

IF(C.EQ.'N ') COTO 1f)

IF(C.EQ.'E ') CALL OUTPT.2(2,EE,TE,402,'N ','LN E',NT)

IF(C.EQ.'X ') CALL 0UTPT2(T,EX,TX,402,'N','X ',NT)

IF(C.EQ.'Y ') COTO 11

COTO 12

11 CG.LL OUTPT3(T,EX,TX,CE,'_'E,402,'N ','X ','LN E',NT)

r

10 CALL PLOT(0.,0.,999)

STOP

1000 FORMAT(A4)

END

C •'••• Function caiculatus R(i)

FC:JC'TION R(J )

COMMON /RCONST/RO,RI,MT

Iq EAL J

R-R0•'•XP(L0G(R1/R0)•J/200.)

P1;TfJ >;:a

E'^1D
124

C ***** °unction calculates X(j)

FUNCTION X(J)

COMMON /X/XJ(0:200)

REAL J

IF(J.GT.200.) J=200.

IF(J.LT.0.) J-0.

J1-INT(J)

X-XJ(J1)+(J- -OAT(J1))*(XJ(J1+1)-XJ(J1))

RETURN

SID

C ***** Function calculates j(X)

FUNCTION J(X)

COMMON /X/XJ(0:200)

REAL J

IF(X.GT.1.) X=1.

IF(X.LT.0.) X=0.

J1-0

1 J1-J1+1

IF(XJ(J1).LT.X) GOTO 1

J-FLOAT(J1)-(XJ(J1)-X)/(XJ(J1)-XJ(J1-1))

RETURN

ND

SUBROUTINE .4ETHOD(U,A,W,JMAX,R,.4T,MDB)

C **'** METHOD sets up the coefficients of the tridiagonal equations

C for the calculation of the new U using .nethod `1T and downstream

C boundary type MOB, if needed. The tridiagonal equations are

solved by TRIO and the riew values of U are returned to T.

DIMaISION U(O:JMAX),A(O:JMAX,3),W(O:JMAX)

IF(MT.EQ.4) GOTO 1

C " "' Set up coefficients for Tra p ezoidal and Backward Euler interior

C **"' schemes

i<

.r
125

00 2 J- 1 , JMA.X-1

GOTO(3,4,5) MT

3 C1=-0.25*R(FLOAT(J))

C2--C1

GOTO 6

4 C1=-0.25*R(FLOAT(J)-0.5)

C2-0._25*R(FLOAT(J)+0.5)

GOTO 6

5 R1-R(FLOAT(J)-0.5)

R2=R(FLOAT(J)t0.5)

C1=-0.5*R1**2/(R1+R2)

C2=0.5*R2**2/(R1+R2)

6 A(J,1)=C1

A(J,2)=1.-C.1-C2

A(J,3)=C2

2 W(J)-U(J)+C1*(U(J)-U(J-1))+C"^*(U(J)-U(J+1))

GOTO(7,8,1) MDB

C ***** Se*_ up coefficients for s pace and space-time extrapolation at

C downstream boundary

7 A(JMAX,1)--1.

A(JM_AX,2) = 1 .

W(JMAX)=0.

GOTO 9

8 A(JMAX,1)=0.

A(,TMAX,2)=1.

W(JMAX)=U(JMA.X-1)

GOTC 9

C ***** Set up coefficients for box method on interior or at downstream

C boundary as appropriate

1 JMIN=JMAX

IF(MT.E2.4)JMIN=1

00 10 J=JMIN,JMAX

C=?(FLOAT(J)-0.5)

A(J,1)=1.-C

A(J,2)=1.+C

A(J,3)=0.

10 W(J)=U(J-1)*A(J,2)+U(J)*A(J,1)

C ***** Se' up coefficients for upstream boundary

9 A(0,2)=1.

A(0,3)-0.

W(0)=0.
14 CALL TRID(U,A,W,JMAX)

RETURN

END

SUBROUTINE TRID(U,A,W,JMAX)

DIMENSION U(O:JMAX),A(O:JMAX,3),W(O:JMAX)

DO 1 J= 1 , JMAX

C C-A(J,1)/A(J-1,2)

Ft

A(J,2)-A(J,2)-A(J-1,3)*C

tf

1 W(J)-W(J)-W(J-1)*C

U(JMAX)-W(JMAX)/A(JMAX,2)

DO 2 J-JMAX-1,0,-1

2 U(J)-(W(J)-U(J+1)*A(J,3))/A(J,2)

RETURN 

END

SUBROUTINE PR%D(T,TX,TE,JMAX,M,RC,R,MT,MDB)

DIMENSION T(M),7X(M),TE(M)

REAL J,K

EXTERNAL R,X

DP.(J) = 100. * (R(J+O.O n 5) -R(J-0.005) )

PI=3.14159

J-TX(1)

T(1)=1.

IF(MT.EQ.4) GOTO 1

C ***** Prediction for trapezoidal schemes

K-ASIN(RC/R(J))

IF(RC.LT.O.) K=PI-K

OM-ATAN(RC/2.)*2.

C1=1./(1.+0.25*RC**2)

4 DO 2 KOUNT-2,.K

T(KOUNT)=FLOAT(KOUNT)

DJ-R(J)*COS(K)*C1

!t

126
o

12 

DK DR(J)*SiN(K)*C1

J=J+0.5*(DJ+R(J+DJ)*COS(K+DK)*C1)

K=K+0.5*(DK-DR(J+DJ)*aIN(Y+[,K)`C1)

TX MOUNT)-X(J)

DE-DR(J+O.S*DJ)*L1*(:.-COS(K-0.5*DK))

IF(MT.EQ.1) TE(KOUNT)-TE(KOUNT-1)

IF(MT.EQ.2) TE(KOUNT)-TE(KCUNT-1)+DE

T_F(MT.EQ.3) TE(KOUNT)=TE(KOUNT-1)+2.*DE

IF(J.GT.O.) GOTO 3

J--J

K=ASIN(RC/R(J))

TX(KOIJNT)-X(J)

3 IF(J.LT.FLOAT(JMAX)) GCTO 2

J-2.*F'-OAT(JMAX)-J

K=PI-ASIN(RC/R(J))

TX(KCUNT)-X(J)

GOTO(6,7,9) MDB

6 TE(KCUNT)=TE(KOUNT)-2.*LOG(TAN(K/2.))

GOTO 2

7 TE(KOUNT)-TE( vuNT)-2.*LOG(ABS(SIN((K-OM)/2.)/SIN((PI-K-OM)/2.)))

F

GOTO 2

8 TE(KOUNT)=TE(KOUNT)-6.*LOG(TAN(K/2.))

2 CONTINUE

RETURN

C ***** Prediction for box scheme

1 T.1= 0.25*RC*RC

T3=0 .5 *RC*T2

K=2.*A:AN(0.5*RC/R(J))

DO 9 KOUNT=2,M

)+ T(KOUNT)-FLOAT(KOUNT)

# IF(J.GE.FLOAI(JMAX)) GOTO 10

I

R1-R(J)

b DJ=(R1+T1/R1)*T2

OK--2.*DR(J)*T3,/R1

R2-R(J+DJ)

J=J+0.5*(DJ+(R2+T1/R2)*T2)

K=K+0.5*(DK+DR(J+DJ)*T3/R2)

TX(KOUNT)-X(J)

TE(KOUNT)-TE(KOUNT-1)-2.*T1*T2*DR(J+0.5*DJ)/(R1*R2)

GOTO 9

r

C

G
0 

10 TX(KOUNT)=1.

TE(KOUNT)-TE(KOUNT-1)

9 CONTINUE

RETURN

END

i SUBROUTINE OUTPTI(X,Y,NPLUS2,C1,C2,NT)

t

DIMENSION X(NPLUS2),Y(NPLUS2),XD(40)

N-NPLUS2-2

NPLUSI-N+1

XL-5.

YL=4.

CAl'_ PLOTS (O , O , NT )

CALL SCALE(X,XL,N,1)

CALL SCALE(Y,YL,N,1)

r- DO 1 I-1 , N/B

1 XD(I)=(X(8.1+1)-X(8'1-7))/X(NPLUS2)

IF(NT.EQ.0) GOTO 2

CALL FACTOR(1 .8 )

CALL PLOT(1.,1.,-3)

CALL AXIS( O.,O.,C1,-4,XL,O.,X(NPLUSI),X(NPLUS2))

CALL AXIS(O.,O.,C2,4,YL,90.,Y(NPLUSI),Y(NPLUS2))

' CALL LINE(X,Y,N,1,1,3)

CALL GRID(0.,0.,1000+N /B,XD,-1,YL,-1)

GOTO 3

! 2 CALL F ACTOR ( 1 .4 )

CALL, PLOT (5 . , 1 . , - 3 )

CALL AXIS(O%,O.,C1,-4,XL,90.,X(NPLUSI),X(NPLUS2))

CALL AXIS(O.,O.,C2,4,YL,180.,Y(NPLUSI),Y("IPLUS2))

Y(NPLUE2)--Y(NPLUS2)

CALI, NEWPEN(2 )

CALL LINE(Y,X,N,1,1,3)

F CALL NEWPEN(1)

CALL GRID(0.,0.,-1,- YL,1000+N/B,XD,-1)

3 CALL PLOT(O.,O.,-999)

RETURN

END

128
129

SUBROUTINE OUTPT2(X,Y1,Y2,NPLUS2,C1,C2,NT)

DIMENSION X(NPLUS2),Y1(NPLUS2),Y2(NPLUS2),Y(4)

N=NPLUS2-2

NPLUSI=N+1

XL=5.

YL=4.

Y( 1 )=1 . X10

Y(2)=-1.E10

DO 1 J= 1 , N

Y(1)=AMIN1(Y(1),Y1(J),Y2(J))

1 Y(2)=AMAX1(Y(2),Y1(J),Y2(J))

CALL PLOTS(O,O,NT)

CALL FACTOR(2.0)

CALL PLOT(1.,0.75,-3)

CALL NEWPEN(2)

CALL SCALE(X,XL,N,1)

CALL SC.ALE(Y,YL,2,1)

CALL A.XIS( 0.,..C1,-4,XL7O.,X(NPLUSI),X(NPLUS2))

CALL AXIS(O.,O.,C2,4,YL,90.,Y(3),Y(4))

Y1(NPLUSI)=Y(3)

Y1(NPLUS2)=Y(4)

Y2(NPLUSI)-Y(3)

Y2(NPLUS2)=Y(4)

CALL NEWPEN(3)

CALL L:NE(X,Y1,N,1,25,1)

CALL LINE(X,Y2,N,1,0,0)

CALL PLOT(O.,O.,-999)

RE ."URN

END

SUBROUTINE OUTPT3(X,Y1,Y2,Z1,Z2,NPLUS2,CX,--Y,CZ,NT)

DIMENSION X(NPLUS2),Y1(NPLUS2),Y2(PIPLUS2),Y(4)

DIMENSION Z1(NPLUS2),Z2(NPLUS2),Z(4)

EQUIVALENCE (Y(1),Z(1))

N-NPLUS2-2

NPLUSI=N+1

XL=4.

YL=2.

ZL=2.
'0

0;2; .: •. - . 3

OF POOR QUALITY

CALL PLOTS(O,O,NT)

CALL FACTOR(1.9)

CALL SCALE(X,XL,N,1)

Y(1)=1 .E10

Y(2)--1.E10

DO 1 J=1,N

Y(1)-A.MIN1(Y(1),Y1(J),Y2(J))

Y(2)-AMAX1(Y(2),Y1(J),Y2(J))

CALL PLOT,3.5,0.5,-3)

CALL SCALE(Y,YL,2,1)

CALL NEWPEN(2)

CALL ACIS(O.,O.,CX,-4,XL,90.,X(NPLUS1),X(NPLUS2))

CALL AXIS(O.,O.,CY,4,YL,180.,Y(3),Y(4))

Y1(NPLUSI)-Y(3)

Y1(NPLUS2)--Y(4)

Y2(NPLUSI)=Y(3)

Y2(NPLUS2)--Y(4)

CALL NEWPEN(3)

CALL LINE(Y1,X,N,1,25,1)

CALL -,INE(Y2,X,N,1,0,0)

CALL 'EWPEN( 1 )

CALL 7-n(0.,0.,2,-1.,4,1.,-21846)

CALL NEW EN(2)

CALL PLOT(-2.3,0.4,3)

CALL PLOT(-2.3,1.0,2)

CALL SYMBOL(-2.23,1.1,0.14,'WAVEPACKET THEORY',90.,17)

CALL PLOT(-2.6,0.4,3)

CALI, SYMBOL(-2.6,0.5,0.08,1,0.,-2)

CALL. SYMBOL(-2.6,0.7,0.06,1,0.,-2)

CALL ::YMBOL(-2.6,0.9,0.08,1,0.,-2)

CALL PLOT(-2.6,1.0,2)

CALL SYMBOL(-2.53,1.1,0.14,'NUMERICAL EXPERIMENT',90.,20)

Z(1)=1.E10

Z(2)=-1.E10

DO 2 J-1 , N

` Z(1)=AMIN1(Z(1),Z1(J),Z2(J))

2 Z(2)-AMAX1(Z(2),Z1(J),Z2(J))

CALL 2I,OT(3.2,0.,-3)

CALL SCALE (Z , ZL , 2 , 1 )

CALL MIS(O.,O., CX , - 4,XL,90.,X(NPLUSI),X(NPLUS2))

CALL AXIS (O.,O.,CZ,4,ZL,180.,Z!3),Z(4))

f Z1(NPLUSI)-Z(3)

Z1(NPLUS2)--Z(4)

Z2(NPLUSI)-Z(3)

Z2(NPLUS2)--Z(4)

CALL NEWPEN(3)

' CALL LINE(Z1,X,N,1,25,1)

ti.

e
131

CALL LINE(22,X,N,1,0,0)

CALL NEWPEN(1)

CALL GRID(0.,0.,2,- 1.,4,1.,-21846)

RETURN

END

1

L

0
I.

132

10. Conclusions

The validity of the asymptotic approach developed in

this paper is demonstrated by the numerical results in

chapter 9. The limitations of the wavepacket theory are due

to the asymptotic approximations involved in treating the

wavepacket as a particle. The stability analyses in

chapters 5 and 6 use fewer a p proximations and so the

asymptotic errors will be substantially smaller. In

particular when the coefficients are constant the analysis

in chapter 6 reduces to the P-stability analysis of Beam,

warming and Yee [51.

The calculation of the asymptotic amplitude equation

and asymptotic boundary conditions for a particular case is

no more difficult than a normal Von Neumann analysis. For

applicable cases the wavepacket theory and the stability

analysis of chapter 5 are straightforward. The general

stability analysis of chapter 6 will usually require

numerical computation. In the more complex cases the main

benefit from this theory will be the insight given by the

asymptotic amplitude equation and boundary conditions. The

amplitude equation gives the group velocities of the

different wavenumbers And the effect of varying

coefficients, which is of great interest since in 2-D

cascade geometries cell lengths can vary by factors of up to

100 in inviscid calculations and 1000 in viscous

calculations. The asymptotic boundary conditions give the

C

amplitude reflection coefficients which provide a practical

criterion for choosing the best numerical boundary

conditions.

There are various possibil4ties for future research

`„e

. 0
133

in this area. Further applications to rela y -vely simple

problems can be done to gain insight into understanding

harder problems and improving boundary conditions

Numerical procedures, such as those suggested in

§6.1.1-6.1.3 , can be developed to solve the equations given

by the stability analysis in chapter 6. Finally t ►,e

asymptotic amplitude equation and boundary conditions can be

extended to 2-D and 3-0.

n 

1C . 
134

References

(11 G.B. Whitham , Linear l+nd Nonlinear Waves , chapter 11

John Wiley 6 Sons (1974)

(21 J. Liqhthill , waves In Fluids , pp 237-260 , Cambridge

Univerity Press (1978)

(3) C.P. Kentzer , 'Group Velocity And Propagation Of

Numerical Errors' , AIAA Paper No. 72-153

(41 R. Vichneve*_sky and J. Bowles , Fourier Analysis of

Numerical Approximations of Hyperbolic Equations , SIAM

Studies In Applied Mathematics (1982)

[5] R.M. Beam , R.F. Warming and H.C. Yee , 'Numerical

Boundary Conditions For Unconditionally Stable Implicit

Difference A p proximations Of Hyperbolic Problems' ,

Paper presented at SIAM 1981 National Meeting, Troy, New

York, June 8-10 1981

(61 M.R. Spiegel , Complex Variables , p127-8 , McGraw-Hill

( 1964 )

6
Appendix

A.1 Finite Difference Onerator Notation

An operator notation for finite difference equations

simplifies the analysi3 of finite difference schemes and is

a necessity for making any general statements and proving

them.

The principal operators are 5 x , central difference,

IA central averaging E

x x x 

, shift operator, 6 , forward

difference and 0 , backward difference. Their definitions

x

are;

5 m 

U. Uj+m/2 Ui

- m/2 (A.1(a))

612U ' (U + U. )

(A.1(b))

7 j +m/2 J- ,n/2

E U U. (A.1(C))

mx ] ]+m 

U. U. - U. (A d) )

mx ] ]+m ]

0 U, U, - U (A.1(e))

mx ] ] ]-m

Usually these definitions will be used with m-1.

The :Hain exceptions are 6 2x which is a node-centered central

difference,

62x U,J 7

' U +1 - U 

J

. -

1

(A.2)

and E

mx 

which can be used to define a general linear

operator,

M

- , U, 1A.3

L

^ j [1 a )

7 jm ]+m

MM 1

1.S

4
^i

136

uR;GINAL PAGE 19

[

M OF POOR QUALITY

so L. ) a

^ t►

E 

m x

(A.4)

L

^

m- 1

When tnere are several independent variables the

subscript on the .f inite operator denotes the direction of

the shift, differencing or averaging. For example if,

U^ = u(x i ,t n ) (A.5)

then 6 x U^ + 0 Uj +1 - Uj (A.6)

and d Un

+f

 . Un +1 - U (A.7)

t 3

The general shift operator expression for a finite

o p erator in 2-D is

L^ _ ) Cmp(]) E mx E pt (A.8)

mop

in applications however this expression can be very

complicated and it is usually sim p ler to express L as a

polynomial in the finite operators. As an example the

operator in 13.3.2 has the polynomial form

L^ - u x d t + r^ 4 t 6 x (A.9)

but in the shift operator form it is,

_ Imo - 1 -r _-z

L 3 2 E x/2 E t/2 + 2 E -x/2 E t/2 2 Ex/2E-t/2

1+r E_ Ex/2--t/2

(A. 10 )

?art of the advanta g e in using o p erator notation

when analysing finite difference schemes arses 'because all

of the finite operators have the same a_genfunction which in

2-0 is exp ( i(34-nA)I.

f^

O
137 o^,CTIAL PAGE 1$

OF POOR QUALI'Y

6 exp(i(j0-nA)) exp(i((j+}) o - nQ)I - exp(i((j- 1 )m- nil) )

21 sin(0 /2) exp(i(j0-nQ)j (A.11(a))

a

x

exp(i(j0- na)) = 2 (exp(i((J+J)*-nQ)j e + exp(i((j-J)0-n:l)J^

= cos(4/2) exp(i(j0 - nQ )) (A.11(b))

E

x 

exp(i(j9-nn)) - exp(i((j +1)0 - nQ)j

exp(i9) exp(i(j0 -nil)]

(A.11(c))

x 

exp(i(j0- nQ)j - exp(i((j+1;0-nQ)j - exp(i(j0-nn)j

( exp(is) - 1 ) exp(i(j9-nn)) (A.11(d))

7 % exp(i(j0 - nQ)j = exp(i(j0-nQ)) - exp(i((j- 1)0-nQ)j

exp(-i0) } exp(i(j0-nQ)j

(A- 11(e))

and similarly,

6 t exp[ik, - nil)] _ -21 sin( /2) exp(i(j0-nil',) (A.11(f))

u t exp[i(j0-nQ)j cos(n/2) exp[i(je-nit)) (A. 11(q))

E . 

= exp(-in) exp(i(j0-nQ)j (A.11(h))

a 

exp(i(j0•-nA)j = { exp(-i:) - t } exp(i(ja-nA)) (A.it(i))

o f exp(i (j0-n0)j _ { 1 - exp(i2) } exp[ i (j0-nil)] (A.11(j))

v

(
138

A.2 ?ri.nciDl.e of the Argument

Let f(z) be an analytic complex function with

simple poles in a region of the complex z-plane and let C

be a closed curve in tl^e region. Then the number of zeros

of F minus the number of poles of f lying inside C is equal

to,

1 1 f , ( z)

271 J

Z

^(Z)

C

= 1

2 _ i [ ln,f) 1C

1

= 2 1T [ arg(f) i C (A.12)

[ ) C denotes the change as Z goes round C

anticlockwise.

arg(f) is defined by,

f(z) = R exp(i9) R,3 real R>0 (A.13)

arg(f ) = A (A. 14)

with the restriction that a must vary continuously as z goes

round C.

The jroc.t is given in many standard text$ on complex

analysis., e.g. [61. Tnis provides a very simple test when

d

considering stability problems in which it is sufficient to

know whether there ara anv zeros in a critical region

without knowing their exact position. This is the basis of

z the Nyquist criterion in control theory stability analysis.

' The test can also be performed numerically relative]!

easily. The step size Iz in going round C is decreased, if

1: -

1

, (^

-o-m^

-

-.1-

- - --

-AN," --^ ^ - - (- ^
139

1

necessarv, until 142rg(f)I < E. since ( arg(f) ] cis

an integer there is no rounding error. The only possible

error is if the magnitude o1? Aarg(f; over one step lies in

the range 2nw -E < Garg(f) < 2nw +E

for some integer n

other than zero. Decreasing E reduces the chance of an

error at the expense of increased computation. c-i/6 should

be adequate in most case.

rs' 
140

A.3 Definitions of Norms and Stabilit

The nn^.ms used i% this p aper are generalised L.

norms. For a continuous func t ion u(x,t) defined on 0<x<X

the norm H u(t)II is defined by,

X

Ilu(t)II 2 f Ju(x,t)1 2 a(x) dx (A.15)

0

where a ( x) is a p ositive non-zero function.

For a discrete function Un defined on C<j<J the J

norm 11U * 11 which is a function of n is defined by,

J

IjUn,12 jUn12 a j(A.16)

0

where a is a positive non-zero function.

j

The s t ability used in this paper is asymptotic

Liapounov stability which is defined as,

Given 6>0 there exists e>0 such that

IIu(0)iI < E __> a) (Iu ( t)II < d

anal b) 1 1 u(t) I I + 0 as t-m

Condition a) is Liapouncv stability which limits

how large an initially small disturbance can become.

Condition b) is asymptotic stability which s p ecifies that a

i.

si:fficiently small initial disturbance must ultimately tend

to zero.

For linear systems of e q uations an equivalent

definition is

a) There exists M>0 such that H U (t)II < M IIu(0)iI
1 

141

and b) J j u (t) j i • 0 as tThe corresponding de°ir,irions for a discrete

function are,

Given 5>0 there exists E>O such that

Ilu 0 ^^ < E -_> IIUn,I < 6

and I j U n j j • 0 as n•m

and for linear discre*.e systema,

a) There exists M>') such that IjU n I, < M IIUO11

ar.d b) I I Un I - 0 a • • a--

I

(