Stability Theorems Concerning High Order Explicit Algorithms for the Linear Advection Equation

College of Aeronautics Report No. 8517

October 1985

UJCHTV/

Kluy

_Ve

"••LfT

Stability Theorems Concerning High Order

Explicit Algorithms for the Linear Advection Equation

by M. Lobo

.10

College of Aeronautics

Cranfield Institute of Technology

Cranfield, Bedford MK43 OAL, England

Cranfield

College of Aeronautics Report No. 8517

October 1985

Stability Theorems Concerning High Order

Explicit Algorithms for the Linear Advection Equation

by M. Lobo

College of Aeronautics

Cranfield Institute of Technology

Cranfield, Bedford MK43 OAL, England

ISBN 0 947767 14 2

£7.50

"The views expressed herein are those of the authors alone and do not necessarily represent those of the Institute."

linear advection equation involves the calculation of the roots of a

polynomial of degree N-1. A theorem is proved to show that if the

algorithm is of the nth order (n<N), then the associated 'stability polynomial'

must have a multiple root of order [—0—] so that, in effect, the analysis need

only be carried out for a polynomial of degree N - [-y-]. This technique for

reducing the complexities involved in stability analysis is illustrated for

families of 4- and 5- point algorithms. Moreover in the special case of

'ideal algorithms' (for which n = N - 1 ) , a complete stability analysis is carried

through, a theorem being proved that such algorithms ar^jlways stable in their

1. Introduction

In the theory of numerical solutions to hyperbolic partial differential

equations, the linear advection equation,

u + au = 0, u(x,0) = f(x), -°°<x<oo (1.1)

is a useful model equation on which to test algorithms for despite its

trivial analytic solution,

u(x,t) - f ( x - a t ) , (1.2)

the attempt to generate numerical solutions leads to many of the difficulties

encountered when trying to solve equations representing physically realistic

situations. One of the pioneer researchers into numerical algorithms for

solving Eq. (1.1) was S.K. Godunov who deduced certain basic results including

the theorem which now bears his name (Godunov (1959)). More recently, a

fairly detailed study was also carried out by Roe (1981).

Any attempt to assess the value of a particular algorithm should take

into account the range of Courant numbers for which the algorithm could

yield stable solutions. During the Second World War, J. von Neumann developed

a technique for testing stability which has now become the standard tool of

the stability analyst; if an initial line of discontinuities is expressed in

terms of a finite Fourier series, the stability of the algorithm would depend

on whether the mangitude of the errors - that is, the absolute value of the

Fourier series - could be held in check. More explicity, if

U. . , = y c u. . (1.3)

(where i & j enumerate the grid points along x & t respectively and m varies

over a prescribed limit) denotes a numerical algorithm for Eq. (1.1), its

stability would depend on whether the associated 'amplification factor'

g(G)

- I ^/^^)^^

(1.4)

m

has an absolute value not exceeding unity for all values of the frequency 0.

If the algorithm (1.3) has N- point support (by this it is implied that

all the coefficients c are identically zero except for a given range of N

consecutive values of the subscript; some of the c may be zero inside the

stated range), the square of the absolute value of the Fourier series (1.4)

can be converted into a polynomial of degree (N-1) in OosQ; this is proved

in Theorem 1 of this paper. Thus, analysis of the stability of an N- point explicit algorithm for Eq. (1.1) involves the calculation of the roots of a polynomial of degree (N-1). Since the coefficients of this polynomial are, in general, rather complicated functions of the Courant number, the stability analysis could turn out distractingly complex even for algorithms with as little as 4- point support. Itisnot surprising, therefore, that researchers in this area have tended to concentrate more on deriving general results of a theoretical nature such as the establishment of a relationship between stability and accuracy. Among recent papers which deal with this problem, mention must be made of Engquist & Osher (1981) and of Iserles & Williamson (1984). In this paper, by contrast, the emphasis is solely on techniques for analysing stability. The most important conclusions derived herein form the statements of two

theorems which have been designated Theorem A and Theorem B. The statement of Theorem A is that the 'stability polynomial' associated with an algorithm of order n must have a multiple root of order [-0-] (where [ ] denotes the

greatest integer function). As a consequence of this theorem, it only remains to investigate the behaviour of a polynomial of degree N - 1 - [ - y - ] , so that the stability analysis is considerably simplified. The statement and proof of Theorem B form a complete study of the special case of ideal algorithms - which are defined as algorithms having the maximum possible order obtainable from the supporting points (in other words, an ideal algorithm is an algorithm for which n = N - 1 ) .

The main body of this paper is divided into three sections. Section 2 at first sets out to prove several relatively easy theorems; all these results

'add up' to yield the proof of Theorem A. The 'use' of the conclusions arrived at in Theorem A is illustrated in Section 3 where the stability of various families of 4- and 5- point algorithms are analysed. The special case of ideal algorithms forms the subject of Section 4 'culminating' in the proof of

2. Preliminary Theoretical Results and the Proof of Theorem A

It is appropriate to begin this section by proving the statement made

in the introduction - that stability analysis for an algorithm with N- point

support involves the investigation of the behaviour of a polynomial of degree

(N-1).

Theorem 1

If the algorithm (1.3) has N- point support, the absolute value of the

square of the amplification factor, as defined by (1.4), can be expressed as a

polynomial of degree (N-1) in Cos 0.

Proof

If r and s denote the minimum and maximum values of m for which c is

m

non-zero (s-r=N-l), the amplification factor g may be expressed in the

form

. s

g(e) = I c {Cos(m0) + /^Sin(m0)} (2.1)

m=r

so that

lg(0)l' = { I c Cos(m0)P + { I c Sin(m0)}2 (2.2)

m=r m=r

which, with the aid of some of the fundamental trignometric identities,

reduces to

lg(0)l' =

I I

c c Cos|(p _ q)gj_ (2.3)

p=r q=r ^ ^

The terms in Eq. (2.3) may be rearranged to assume the form

N-1

|g(0)r = I d.Cos(j0) (2.4)

j=0 J

where

and 3_j

'j ^ \ l V m + j <^°^ 1<^'<N-1) (2.6)

Recalling that

Cos (j0) = T^(Cos0) (2.7)

where T. is the Chebyshev polynomial of degree j , it follows, by substitution

Ü

of (2.7) in (2.4),that|g|^ can be expressed as a polynomial of degree (N-1) in Cos 0.

Remark

Although |g|^ can be expressed as a polynomial in Cos 0, it can equally conveniently be expressed as a polynomial (of the same degree) in

a = Cos 0 - 1; in fact this representation turns out to be more amenable to analysis. Consequently, if we define the polynomial,

N-1

6(a) = |g{Cos"'(a+l)}M = I d.T.(a+l) . (2.8)

j-0 J J

stability of the numerical algorithm would require that G(a) have a value not exceeding one for all values of a in the closed interval [-2,0].

Eq. (2.8) may be expressed more compactly in the form

G(a) = I 6 a'' (2.9)

k=0 ^

where

and t. , is the coefficient of a in the polynomial expansion of Cos j0 in a.

J > K

k

t

From the theory of Chebyshev polynomials, it is fairly easy to express t. . as an explicit function of j & k but, for the sake of completeness, it will be indicated, in the next theorem, how this expression may be derived

from elementary principles. Theorem 2

Each coefficient t. , (j>0, k>0) is a polynomial of degree k in j^ or,

J ,K more specifically, tj^O = ^' (2.11) and, if k > 1 , j^(j^-l)(j^-2)...(i^-(k-l)^) 1 . . . . .

rr~. 3 . 5 . . . (2k-T) • Ti ^^ ^-^

L o if k>j.

(For ease of notation, t. . may sometimes be simply denoted as t . . ) .

Proof By d e f i n i t i o n

J k

Cos j0 = Z t . , a .

k=0 ^'^

(2.13)

S u b s t i t u t i o n of Eq. (2.13) i n t o the well-known t r i g n o m e t r i c i d e n t i t y

Cos j0 - 2 Cos 0 Cos (j-l)0 - Cos (j-2)G (2.14) yields j ,. J-1 .. J-2 (2.15) I t , a = 2(a+l) I t. , .a - I t. „ ^ a J-2

I

k=0 ^'" k=0 ^ "'" k=0

and, after equating like powers of a, the following recurrence relations emerge

^j,0 ^ ^^j-1,0 " ^j-2,0 if j >2 (2.16)

^j,k = 2(tj.i^,.l . t._^^,) - t._,^, if j > 2 , k > l , (2.17)

Eq. (2.16) subject to the initial conditions t«Q= 1 and t,Q = 1 yields the solution (2.11). Eq. (2.11) now forms one of the boundary conditions for the recurrence relation (2.17); the other two boundary conditions are

'0,0 1, t

0,k

0 i f k > 0 (2.18)

h . o = i' ^1,1 = 1' ti,k = o if k > l

(2.19)

The solution of the recurrence relation (2.17) subject to the boundary conditions (2.11), (2.18), (2.19) is precisely (2.12).

Remark

Eq. (2.12) may be expressed more compactly in the form .(j + k - 1 ) ! 2 \ .. .

^ j , k = < (2.20)

'- 0 if k > j

is impor'

polynomial of degree k in j^ whose constant term is zero if k ^ L

but for the purposes of Theorem A, it is important to be aware that t. . is a

At this point it is appropriate to insert a fundamental relationship

between the coefficients c of the algorithm (1.3) which was worked out by

Roe (1981). As a consequence of this relationship, another will be derived

(Theorem 3) which, in conjunction with Theorem 2,will quickly yield the proof

of Theorem A.

Theorem (Roe (1981))

If the algorithm

has an order of accuracy equal to n, then the relationship

I m\

= (-v)^ (2.22)

m=r

holds for all values of k between 0 and n (here v = aAt/Ax is the Courant

number associated with the grid proportions). Conversely if (2.22) holds

for all k between 0 and n, the algorithm will be of the nth order.

Using concepts from linear algebra. Roe gave a neat proof of this

theorem; he showed in fact that if the function u, on which the algorithm

(2.21) is applied, is a polynomial of degree n, then (2.21) would yield

the exact solution of the linear advection equation (1.1) if and only if

conditions (2.22) hold. He further demonstrated that for more general

functions (assumed to be differentiable at least (n+1) times), conditions

(2.22) implied that the truncation error would have a leading term of order

(Ax)" (Theorems 1 and 2 of Roe (1981)).

Theorem 3

If (2.21) is an algorithm of order n, then

I

c c (p-q)'< = 0 (2.23)

(P,q) P "^

for all even positive values of k such that k<n. (The summation

I

is

understood to extend over all distinct pairs (p,q) in the range ^P'"^'

i^<P<s, r<q<s; if the algorithm has N- point support, there are totally NCa such

pairs. Note that unless k is even

I

is not properly defined.)

Proof

I

CpC (p-q)^

-Ï I I

c c (p-q)^

(p,q)

V

^ p=r q=r

^

^

- Ï I I

c c I(-l)l(kC )p^-^l

p=r q=r P ^ 1=0 '

-

Ï I

{(-l)^kC )[

I c/-^]\_ I cj-]}

1=0 ' p=r P q=r ^

= 1 I {(-l)\kCJ(-v)''"^-v)^} (from (2.22)

1=0

k

^

1

=

H - v r I

(-l)'(kC,) = 0. (2.24)

1=0 '

From Theorem 3, it follows, a fortiori, that

s s

I

I

V (p-q)'' = 0 (2.25)

p=r q=r ^ ^

for all positive (even or odd) values of k such that k<n. In fact if k is odd

Eq. (2.25) is trivially satisfied and the reasoning laid forth in the proof

of Theorem 3 is superfluous. However, for the purposes of Theorem A, it is

important to note that the identity holds when k is even.

Remark

Theorem 3 provides an alternative proof of Godunov's Theorem which

states that there are no monotonie algorithms with higher than first order

accuracy. (An algorithm is called monotonie if whenever the data lu- .} is

a monotonie increasing or decreasing function of i, the solution iu.. ,} is

monotonie in the same direction). In fact Godunov (1959) showed that an

algorithm is monotonie if and only if none of the e are negative (the proof

has also been given by Roe (1981)); if an algorithm has second or higher order

accuracy, there will be at least oneeven value of k for which Eq. (2.23) must

hold so that (neglecting the trivial ease of an algorithm in which all but one

of the c„ vanish - in which case it would reduce to an identity) at least one

m •''

of the e must be negative and the algorithm cannot be monotonie.

With the aid of the results deduced so far, we are now in a position to

prove Theorem A.

Theorem A

If G(a) is the polynomial representation of the amplification factor

associated with an algorithm of order n, then G(a) - 1 has a root of order

[•^ip—] at a = 0. (Here [ ] denotes the greatest integer function).

Proof

In other words, it must be shown that 6n = 1 and that 61, - 0 for

l<k<[^]. Now

N-1

6Q = j

d.t.Q

(from (2.10))

N-1

=

I

d. (from (2.11))

j=0

^

s

=

I c ' + 2 I

c c (from (2.5) & (2.6))

m=r

(p,q)

^ ^

s

I

m=r

ll

cj^ = 1, (2.26;

the last equality following from Eq. (2.22) with k=0. For l<k<[^],

N-1

K

= I d.t . ; (2.27)

this follows from Eq. (2.10) and the fact that t«. = 0 for k>0. For

j>l, k>l. Theorem 2 states that t., can be expressed as a polynomial of

degree k in j^ with constant term zero. Writing this as

^jk =

^l^\/^

(2.28)

and subsituting in (2.27)

N-1 k

j

k N-1

6. =

I

d..(

I a..p)

^

j=l J 1=1

^^

1=1 •" j = l J

..21)

k

= 2

I

a^ (

I

c c (p-q)^^) (from (2.6)) . (2.29;

1=1 (p,q) P ^

For each term, the index 1 is a positive integer such that 21<n; consequently, Theorem 3 states that the expression within the brackets vanishes for each value of 1. This shows that ^i, = 0 for l<k<[Y] and completes the proof of Theorem A.

It need hardly be emphasized that the conclusions of Theorem A provide an approach for drastically reducing the complexities involved in stability

analysis. A measure of this 'reduction in complexity' is best illustrated if we consider the case of second order algorithms having 3- point support. For

such algorithms. Theorem A states that the second order polynomial G(a) -1 has a double root at a = 0 so that analysing the behaviour of G(a) effectively implies analysing the behaviour of a constant! This constant is precisely

62 =_ tasda = 2d2 = 4c e . What has in fact been proved is

Theorem 4a

A second order algorithm (of the type (2.21)) with 3- point support is stable if and only if the product of its two extreme coefficients does not exceed zero.

It is possible to go a step further. The coefficients of the

second order algorithm, spanning the grid-points (i+r, i+r+1, i+r+2) on the jth time level, may be determined by solving the system of equations (2.22) with m = r, r+1, r+2 and k = 0,1,2. In fact one obtains

c = Hr+l+v)(r+2+v)

^r+1 = -(r+v)(r+2+v) (2.30) ^r+2 = H>"+v)(r+l+v)

and the product c c ^ will not exceed zero if and only if

-(r+2)<v<-r. (2.31) In other words, we have proved

Theorem 4b

A second order algorithm (of the type (2.21)) with 3- point support is stable for all Courant numbers in the CFL range.

It ought to be mentioned at this point that Theorem 4b is not really a 'new' result because all second order algorithms with 3- point support are essentially equivalent to each other (one can be obtained from another by an appropriate translation of the Courant number ) and therefore the stability

properties of such algorithms may be studied by considering their prototype

- the 'Lax-Wendroff' (for which r = -1, s = 1) and it is well known that

the 'Lax-Wendroff' is stable for all Courant numbers in its CFL range [-1,1].

But the manner in which Theorem 4b was proved does serve to illustrate the

importance of the conclusions derived in Theorem A.

In the next section, the stability criteria for more general

families of 4- and 5- point algorithms will be analysed with the aid of

the theory presented in this section.

3. Application of the Theory for Analysing the Stability of Three Specific Families of Algorithms

Roe (1981) proved that algorithms of order n with N- point support form a family describable by (N-n-1) parameters, each parameter being a function of the Courant number v (this can in fact be directly deduced from Eq. (2.22) because (n+1) of these equations would have to be satisfied by N coefficients). While studying the properties of such a family, it is useful to establish a relationship between the stability criteria and the parameter(s); for practical applications, one can then select a parameter yielding an algorithm with

desirable stability properties. In this section, the stability of the one-parameter family of second order algorithms with support over the four points

(i-2, i-1, i, i+1), the one-parameter family of third order algorithms with support over the five points (i-2, i-1, i, i+1, i+2) and the two-parameter family of symmetric second order algorithms with support over the five points

(i-2, i-1, i, i+1, i+2) will be analysed with this view in mind. It ought to be mentioned that the first of these three problems has already been

investigated by P.L. Roe, J. Pike and P.K. Sweby (see, for example. Pike (1984)); here the analysis is carried out mainly with the object of illustrating how

the theory of the preceding section may be applied to advantage.

A convenient parameter to choose, for each of the one-parameter families, is X = 1-CQ. The coefficients c , for the second order 4- point family, can now be determined by solving the system of equations,

Cj, + C , + C , = A -2c_2 - c_^ + c^ = -V (3.1) 4c_2 + c_, + c, = v2 yielding

c_2 = i(v^ - A)

c_i = 1^(1 -^) + ^ (3.2)

^1 =4v(l4v)+^A

and, of course C Q = 1 - X. (3.2a)

For an algorithm with 4- point support, the function G(a) assumes the form

3

I

k=0

but, from Theorem A, <SQ = 1 and 5, = 0, so that

G(a) = I 6|^a^ (3.3) G(a) = 6 a' + 6 a^ + 1 3 2 = ( f d.t.Jct^ + ( I^d.t.Ja^ + 1 (from (2.10)) j=o ^ J-^ • j=0 -J -^^ = 4d3a' + (12d3 + 2d2)a^ + 1. (3.4)

Now if we define the linear function

L(a) = [G(a) - l]/a^

= 4d3a + (12d3 + 2d2), (3.5)

the condition for stability becomes L(a)<Ofora c[-2,0]. Now if d3>0, L(a) attains a maximum at a = 0; if d3<0, the maximum occurs at a = -2. Two conditions therefore arise

d3>0 and 6d3 + d2<0 (3.6a) d3<0 and 2d3 + d2<0 . (3.6b) Recalling that d3 = 2c_2C, ; d2 = 2(C_2CQ + c_,c,), (3.7) s u b s t i t u t i o n of Eqs. (3.2) in (3.6) y i e l d {|(3 - V) - A } { X - V M > 0

4 ( 2 v - l ) A + v2(3v^ - 8 v + l)<0

and

(3.8a)

{^(3 - V) - A H A - V ^ } < 0

{l{3 + v2)-AHA-iv^}k 0 .

(3.8b)

The shaded region in Fig. (la) depicts the portion of the (A,v) plane where conditions (3.8a) are satisfied and the shaded region in Fig. (lb) depicts the portion where conditions (3.8b) are satisfied. The two results are combined in Fi g, (Ic).

Important members of this family of algorithms are the 'Lax-Wendroff' (which spans the 3 points i-1, i, i+1) and the 'Warming-Beam' (which spans the 3 points i-2, i-1, i). The respective values of the parameter A may be obtained by setting c_2 = 0 or c, = 0 yielding A = v^ (Lax-Wendroff) and

A = ^-(3-v)(Warming-Beam). Substitution into conditions (3.8) quickly reveals that these algorithms are stable for Courant numbers in the intervals [-1,1] and [0,2] respectively. It is noteworthy that though the curves A = v^ and

A = ^(3-v) form boundaries to the shaded areas in Figs.(la) and (lb), they are both in the interior of the stability regime depcited in Fig. (Ic).

A member of the family having a special status is its unique third order representative. Its coefficients c (m = -2, -1, 0, 1) may be obtained by solving Eqs. (2.22) for k = 0, 1, 2, 3. It turns out that A = .^v(l +2v-v^) yielding the stability regime [0,1] for the Courant number.

We now turn our attention to the third order family of algorithms supported by the 5- point network (i-2, i-1, i, i+1, i+2). With the

parameters A defined as before by A= 1 - C Q , the coefficients c (m = -2,-1,1,2) form the solutions to the system

'-_2 + c_, + c, + C2 = A (3.9; -2c_2 - c_, + c, + 2C2 = -V 4c_2 + c_, + c, + 4C2 = ^^^ -8c_2 - c_, + c, + 8C2 = -^^ that is c_2 = ^ [ - ^ ^ + 2v2 + v^ - 2A] c 1 = i [4v - v^ - v' + 4A] (3.10) Cj = -g- [-4v - v^ + v^ + 4A] ^2 " 17 1^^ "^ 2v^ - v^ - 2A] .

For this family, G(a) is a polynomial of degree 4 with S Q = 1 and 6, = 0; it may be quickly verified that

62 = 40d^ + 12d3 + 2d2

63 = 32d4 + 4d3 (3.11)

so that if we define the quadratic function Q(a) = [G(a)-l]/a2

= Sd^a^ + (32d^ + 4d3)a + (40d^ + 12d3 + 2d2), (3.12)

stability would require that Q(a)<0 for ae[-2,0]. It is convenient to first derive a necessary condition for this by evaluating Q(a) at a = 0 and a = -2; in fact Q(0) = 2[20d4 + 6d3 + d2] (3.13a) Q(-2)= 2[4d4 + 2d3 + d2]. (3.13b) Recalling that d^ = 2c_2C2 d3 = 2(c_2C, + c_,C2) (3.14) d2 = 2 ( C _ 2 C Q + c_^c^ + CQC2),

the conditions Q(0)<0 and Q(-2)<0 y i e l d , r e s p e c t i v e l y , on s u b s t i t u t i o n of (3.10)

^ v - ( 5 - v 2 ) < A ( 3 . 1 5 a )

lv^<A< l-(v2 + 3) . ( 3 . 1 5 b )

The subset of the (A,v) plane where inequalities (3.15) hold is depicted in Fig. (2a). Of the three disconnected regions S_, S Q , S , only S Q lies wholly within the CFL range. Moreover, itmay be checked, albeit somewhat tediously,

that

Q(-l) = 2 [Sd^ + 4d3 + d2]

is always positive in the regions S_ and S so that one need only concentrate on S Q .

In fact it is true that a parameter A restricted to S Q suffices to ensure stability. If the maximum value of Q(a) in [-2,0] does not occur at either of its end points, then Q"(a)<0, that is c_2C2 < 0 or

v 2 - 1-^(1 - v2) < A < v2 + Lyi(l - v ^ ) . (3.17) The shaded area in Fig. (2b) depicts the subset of S Q where inequality (3.17)

doesnothold thereby providing one sufficiency condition. Another may be obtained by determining the subset of S Q for which Q'(0) >0,an inequality which reduces to

A < i(2 + 5v- -v-*), (3.18)

depicted by the shaded area in Fig. (2c). The two shaded areas overlap and together make up the region S Q (Fig. (2d)).

The unique fourth order representative of the family is of particular interest. By solving Eqs. (2.22) for k = 0,1,2,3,4, one obtains its

coefficients c (m = -2,-1,0,1,2) and A turns out to be ^(5v^ - v"*) which though stable forv £[-1,1] - the maximum range possible for a member of the family - nevertheless borders on the instability zone. This is a consequence of the fact that Q(0) becomes identically zero if the algorithm is fourth order - the statement of Theorem A for a fourth order algorithm.

At this point, a brief comparison of the stability regimes depicted in Figs. (Ic) and (2d) might be of interest. It is noteworthy that though the analysis of the behaviour of Q(a) is more complex than the analysis of the behaviour of L(a), the stability regime which finally emerges in the former case is relatively simple and bounded by only two curves. This appears to be a consequence of the symmetric nature of the family spanned by the 5- point network. By contrast, the stability regime for the asymmetric 4- point family involves three curves on the boundary and one of these appears as two disjointed segments. It is rather curious that the curve A= •j(v^ + 3) forms part of the boundary of the stability regime for each of these families.

Another striking distinction is that while the regime in Fig. (Ic) can admit parameters yielding algorithms that would be stable for all Courant number in the CFL range, this is not true of the regime in Fig. (2d). This appears to be a consequence of the restrictions generated by the additional complexities of the higher order family and if one seeks algorithms which are stable in the full range [-2,2], it would be necessary to consider the two-parameter family

of second order algorithms in the 5- point network (i-2, i-1, i, i+1, i+2) to which we now turn our attention.

The stability analysis which follows will not be as comprehensive as in the preceding two cases; rather the aim will be to determine a subclass of the two-parameter family of second order algorithms with support over the 5- point network (i-2, i-1, i, i+1, i+2) all of whose members are stable in the full CFL range [-2,2]. Other desirable properties are symmetry aboutv= 0 and exactness atv= 0. If A- and AQ denote the even and odd parts of the coefficient C2 (as a function of v ) , then

C2 = A^ + AQ (3.19) and symmetry would require that

c_^- A^ - AQ (3.20)

Ap and AQ may now be considered the two parameters of the second order family and the remaining coefficients c_,, C Q , C, may be determined by solving the system of equations

c_, + CQ + c, = 1 - 2A -c_^ + c^ = -V - 4AQ (3.21) c_, + c, = v2 - 8A^ yielding c 1 - 2 - 4A^ +_ v(l+ v) 4, 2A ^^AQ C Q = 1 - V ^ + 6A^ * Ci = - ^^^^^- 4A^ - 2AQ . (3.22)

Note that Eqs. (3.19), (3.20), (3.22) define an algorithm which is 'Lax-Wendroff' plus a fourth order dissipation (A ) and a third order dispersion ( A Q ) .

For stability, it is necessary and sufficient that the quadratic Q(a)

defined by (3.12) should not exceed zero fora G [ - 2 , 0 ] ; it is convenient

to first derive the necessary conditions, Q(-2) <0, Q(-l) <0, Q(0) <0, which yield respectively on substitution of (3.14), (3.19), (3.20) and (3.22)

- i v ^ d - v ^ ) + 2(1 -2v2)A^ + 16Ap2 < O ( 3 . 2 3 a )

- | v 2 ( l - v ^ ) +2A^ + 4VAQ < O . (3.23c)

From Eq. (3.23b), i t f o l l o w s , a f o r t i o r i , t h a t

[ ^ E " 1 ( ^ - ^ ' ^ ] '

-^

A

• (^-^'^^

Eqs (3.23a) and (3.24) may be expressed more simply as

^(v^ -1) < A^ <iv2 (3.25a)

and

- i(2-v^) < A^ < iv^ (3.25b)

respectively.

For practical purposes, it is desirable to keep the parameters simple and it would be useful to know if an assumption of the form

A^ = a^ + bgV^ (3.26a)

AQ = a^ + b^v^ (3.26b)

could yield a wholly stable algorithm for some particular choice(s) of a„, b„, a„, b„. The condition for exactness at v = 0 would require that

e e 0 0

a = O and (assuming a = 0) conditions (3.25) applied at v = 2 together 1

imply that b = 3 oi^

A^ = ^v^ (3.27) Substitution of (3.27) into the left-hand side of (3.23a) (which is simply

Q(-2))yields Q(-2) = 0. It follows that a necessary condition for stability

is Q'(-2) < 0. In conjunction with condition (3.23c) (that Q(0) i 0 ) , this

necessary condition would also become sufficient. The condition Q'(-2) < 0 assumes the form

On s u b s i t u t i o n of ( 3 . 2 7 ) , the t h i r d and f o u r t h terms cancel each other and i n e q u a l i t y (3.28) may be r e w r i t t e n as

- ^ < A Q < 0 i f v > 0 (3.29a)

0 < A p < - ^ i f v < 0 . (3.29b)

Moreover c o n d i t i o n ( 3 . 2 3 c ) , a f t e r s u b s t i t u t i o n of ( 3 . 2 7 ) , reduces simply to

AQ < - ^ i f v > 0 (3.30a)

AQ > - ^ i f v < 0 . (3.30b)

Combining (3.27), (3.29) and (3.30), we obtain the necessary and sufficient

conditions for the algorithm to be stable in [-2,2], that is

r - ^ < A Q < - ^ i f v > o

A = -gv^ andj (3.31)

l - X ^ < A Q < - ^ i f v < 0 .

Note that in the derivation of (3.31), though it has been assumed that

(3.26a) holds, it is not essential that A Q must be of the form (3.26b).

The two extremes

^n = - j

and A Q = -y^ provide interesting special

cases. If A Q = - J, c 1 = c, = 0 and the algorithm is simply a 'Lax-Wendroff'.

If A Q = - yg, substitution in (3.19), (3.20) and (3.22) yield

^-2 8^ ^ 1 6 ^

^-1 = h - h"

<^o

^1

•=2

=

l - > ^

= - iv . |v=

•

5"^

-

iV'

(3.32)

an algorithm which, apart from being stable in the full CFL range, also has the

desirable properties of simplicity, symmetry and exactness at v = -2,0,2.

To conclude this section, it may be noted that the algorithm defined

by (3.32) is a special case of a class of algorithms studied by Hall (1985)

Hall investigated M-step algorithms, each step having 3- point support;

to reduce computational demands the dissipative term A^u. . was evaluated

only on the first step. For M = 2, this class of algorithms assumes the

form

^-j^r ^ j - vVi,j^^i^^^'^,j

^,j+i = ^,j -"Vi,j+r V^'^^j'

(3.33)

where a., 6.. and

$2

^""^ parameters which may be varied to provide algorithms

satisfying desirable properties. If expressed in the form of (1.3), the

class of algorithms (3.33) yield the coefficients,

a 6

c_2 = > ' + -T^'

^-1 " I ^ ^2^' • ^1^'

a-,

c. = 1 - ( 4 + 26jv2

^ ^ ^ (3.34)

c^ =- J +

B^v'

+

B^v'

a

B,

C2 4V 2 ^

It may be easily checked that any algorithm of this class is always first

order accurate and the condition for second order accuracy is

a^ + B2 =

J.

(3.35)

Moreover, the class of algorithms is symmetric with

a^ B,

^E " T"^^ ^"^ ^0 " ' T"^'-

(

^

-

^

^

^

The preceding analysis now enables one to select the free parameters

a, and 6, to ensure stability over the CFL range [-2,2]. In fact we must

have

The condition for second order accuracy (3.35) then yields B2 = 0 and

the two step algorithm (3.33) may be rewritten as

^j^i =^,j-l^Vi.j^r^^^"ij

(3.38)

which, if expressed in the form (3.34), is precisely the same as (3.32).

In practice. Hall tested the algorithms for nonlinear equations in

two dimensions so that the parameters could only be selected by trial and

error. For these cases, he obtained stability for Courant numbers with

modulus up to 1.9.

4. The Special Case of Ideal Algorithms

An algorithm will be called ideal if its order of accuracy is the maximum that its supporting points can offer. Examples of ideal algorithms are the 'Lax-Wendroff', the 'Warming-Beam', the unique third order

representative of the family of algorithms whose stability is depicted in Fig. (Ic) and the unique fourth order member of the family whose stability is depicted in Fig. (2d).

An ideal algorithm of order n will be called normal if its CFL range is symmetrically situated about the interval [0,1] if n is odd and about the interval [-1,1] if n is even. This interval will be called the Central range of the normal algorithm Among the examples of ideal algorithms cited in the previous paragraph, the 'Lax-Wendroff' and the third and fourth order algorithms are all normal but the 'Warming-Beam' is not.

In general, the Central range of an ideal algorithm of order n will be defined as the interval of length 1 if n is odd and 2 if n is even -about which the CFL range of the algorithm is symmetrically situated. More

specifically, if the algorithm has support over the (n+1) points (r, r+1 ... r+n), its Central range extends from -[r + (•'^)] and - [r + (^^)] if n is odd and from -[r + (5- + 1)] and -[r + (7 - 1)] if n is even.

Since all ideal algorithms of a fixed order are essentially equivalent (see the remark following Theorem 4 b ) , the properties of an ideal algorithm can be studied by considering the normal algorithm of the same order. In the analysis which follows, it will sometimes be found convenient to assume that the algorithm is normal, but all results thereby deduced will of course hold for all ideal algorithms of that order (with an appropriate translation of the Courant number).

The primary aim of this seciton is to investigate the stability of ideal algorithms. As in Section 2, certain preliminary results will first be

deduced - all of which lead up to the statement and proof of Theorem B. The chain of reasoning required to complete the proof of Theorem B is somewhat tricky and involves excursions into fields of mathematics which are somewhat far removed from the essence of this paper. For this reason - and in order to be able to present the significant results first - part of the analysis will be relegated to an appendix.

Theorem 5

If an ideal algorithm of order n is expressed in the form s

u. . 1 = V c u.

1,j+l L m i+m,j

(where, of course, s = r+n), then the coefficients c are given by

m

^''{\>) n C m-r r ^ ' | m-r i

v+m

(4.1) (4.2)

where v denotes the Courant number and

'lJ^(v) = (v + r)(v + r + 1) ... (v + s - l)(v + s) . (4.3)

(Here C denotes the combinataric function and should not be confused with the coefficients c).

Proof

Roe (1981) showed that ideal algorithms must yield exact solutions for all integral values of the Courant number within the CFL range (in fact if V is an integer, say v = -p, conditions (2.22) are satisfied by Cp = 1 and all other c = 0 so that the algorithm is exact). In other words u. . 1 must be obtained from the (n+1) values of u. . by Lagrange

interpolation for which each c must be defined by (4.2). Theorem 6

In the case of an ideal algorithm, the coefficients 6.(k>0), as defined by (2.10), may be expressed in the form

6, =21

_I

_nC„

s-m

„ ^ m-r v+m , ^ „ .m=r j=r-m ( I (-1)'^*T.. .. nC, j,k m-r+j (4.4) where

j,k

(J+k-1)

U-k)!

2^

T2kr

if j is positive (0 < k <; j)

-J.k

if j is negative (0<k<-j) (4.5)

for all other combinations

(Note that if j is positive, T . . = t. ./j and that although t. . is

J J "^ J >l^ J > K

always an integer, T . . need not necessarily be one).

Proof By s u b s t i t u t i n g (4.2) in ( 2 . 6 ) , we get f o r j ^ l 2

'i-'

^^A-^V

n!

-1 m=r

s-j 2(m-r)+j n C ^ nC^,^^^

^ ^ ' v+m v+m+j

,(-1)"

r , S

^ ; ( v )

_ 2

I nC,

nC

1

-J m=r

m-r m-r+j 'v+m

v+m+j^ (4.6)

S u b s t i t u t i n g (4.6) i n (2.10) (and r e c a l l i n g t h a t i f k > 0 , the lower l i m i t of summation i n (2.10) may be replaced b y l ) , w e g e t , a f t e r rearranging terms

>^(v) Y

• s

I nC

/ j I

m=r

m-r v+m

~ [l ( - 1 ) - ^

j=r-m

"'^m-r^J *

I (-l)'V"Vr.jl

S-m

I

j=l

(4.7)

which, in view of the definition of T .- >., is precisely (4.4).

J ,K

Remark

From Theorem A, we know of course that 6. must vanish if 1 <k < [•2], but it would be interesting to verify this directly from Eq. (4.4). In fact T.. .

J , K

is a polynomial of degree (2k-l) in j, so that, taking into account the combinatoric identity

I (-l)^t^nC. = 0 if 0 < 1 <n, t=0 ^

(4.8)

it follows that 5|^ must vanish if 2k-l <n. Since the formulation (4.4) for 6. holds only when k > l , the result follows.

At this point, it must be remarked that the formulation (4.4) for 6. can be greatly simplified; in fact if [^]<k<n,

. _ (-1)" 2k ,s, ^ ,r+k-l, s

«k - - ï ï r k(2k-n-l)! V ( ^ ) ^-k+l(^) •

(4.9)

The reduction of (4.4) to (4.9) is a fairly challenging mathematical exercise but as itisnot really germane to the main issues of this paper, it will be dealt with in an appendix.

Assuming that (4.9) holds, we are now in a position to state and prove Theorem B.

Theorem B

An ideal algorithm is always stable for Courant numbers in its

Central range.

Proof

A necessary and sufficient condition for stability is that the

polynomial G(a) defined by

G(a) =

I

6,a

k=1 ^

G(a) - 1

(4.10)

should have a value not exceeding zero for all values of a in the closed

interval [-2,0]. A sufficient condition for this to happen is that for k > 1 the

6. s are non-negative for odd values of k and non-positive for even values of k.

It has already been proved that the 6. s vanish if k

Hj].

For values of

k>[ö-]> it will now be shown that 6. is strictly positive if k is odd and

strictly negative if k is even - provided the Courant number v is in the

Central range of the algorithm.

s r+k~1

Consider the two functions v^ (v) and ^s_i(.i(^)- If the Courant number

lies within the Central range of the algorithm, the signs of ip (v) and

r+k-1 n

^

^s-k+l(^^ (k > 2") will be the same if (n- k) is odd and opposite if (n - k) is

even. Therefore the product of the two functions will be positive if (n - k)

is odd and negative if (n - k) is even and Eq. (4.9) may be rewritten as

( - i ; >k - 1

k(2k-n-l)!

.k

(-1)

n-k+1

t=(v) <:^;5(v)

, , ^ k - l 1 2 I s , \ I « ^ t k - l , A (-1) — r TTÖT T T T U' i"^) i)' 1 1 ("^ )

^ ' n k(2k-n-l) r^ ' s-k+P 'I

(4.11)

which is clearly positive if k is odd and negative if k is even. This

proves that, provided the Courant number lies in the Central range, the

polynomial

Qia)

will always be negative for negative values of a - and in

particular for -2 <a_< 0. The amplification factor g(0), defined by (1.4) must

therefore have an absolute value not exceeding one and the proof of Theorem B

is complete.

Further Remarks

It is of course natural to enquire about the stability 'prospects' of

ideal algorithms for Courant numbers outside the Central range. The

situation can be pictured more easily if we consider the case of normal

algorithms for which the Central range is [0,1] if n is odd and [-1,1] if n

is even. If vlies in the open interval (m,m+l), 0 < m < n , then if m is odd

the leading non-zero term in the polynomial expansion of G(a) as defined by (4.10)

(in other words the term 6 , if n is odd and 6 if n is even), will have

?

r^

an incorrect sign (it will be strictly negative if its suffix is odd

and strictly positive if its suffix is even). This would imply that

G(a) must be positive in some neighbourhood of a = 0 so that the

algorithm cannot be stable. If m < 0, the situation is simply a reflexion

about the midpoint of the Central range (that is, the sign of the leading

term will be incorrect in alternate intervals starting with the interval to the

immediate left of the Central range).

If m(<n) is positive and even, the leading term (5 or 6 ) will have

_'^

_{^ "}

_n+1

_n

_,

2 2"^^

the correct sign so that ïï(a) will be negative in a neighbourhood of a = 0.

Nevertheless the algorithm will still not be stable because in all intervals

at a distance of more than one unit from the Central range, the sign of 5(a)

will alteratieast once in the interval (-2,0). In fact by substituting

(4.9) in (4.10) (and recalling that 6 Q = 1 and 6|^ = 0 if l < k < [ ^ ] ) we have

G(a)

V

^ *^(-) "l „ „ R 2 I W ^-W*-)"'- f^-l^l

Since S is a function of both v and a, it would be appropriate to denote it

henceforth as G"(v,a).

Now consider the function

H(v,a) =

I uJr..^^^ €.t](^)^\

(4-13)

-n

k=[^+l]

k(2k-n-l)! ^s-k+l'

r+k-1

Each term

^ _i, -i

is a polynomial of degree (2k-n-l) in v; this attains a

maximum value when k=n and consequently H(v,a) is a polynomial of degree

(n-1) in V. In the half-plane a < 0, -00 < v < °°, there will be (n-1) curves

which form solutions to the equation H(v,a) = 0 . If n is even, one of these

curves will be the straight linev= 0. Each of the other curves will have

as one asymptote the straight line a = 0 and the other asymptote will be a straight line v = m where m is an integer (non-zero if n is even) strictly greater than -s and strictly less than -r. Fig. (3a) depicts how these curves would appear in the half-plane a < 0 if n is odd; Fig. (3b) depicts the situation

when n is even. In both cases, it is assumed that the algorithm is normal. Note that for each value of n, there are exactly (n-1) curves forming solutions to the equation H(v,a) = 0 and since the degree of H(v,a) as a polynomial in v is (n-1), it follows that in no case can one of these curves represent a multiple root of H(v,a) so that the sign of the function must change across each of these curves.

Now if n < 2 , the Central range covers the whole CFL range. For algorithms of order 3 or greater, it is easy to determine the root of Eq. (4.13) along the line V = 2 because in this case the coefficients of all the powers of a except the two lowest vanish so that the equation becomes effectively linear. That root is, in fact

a = - •^(n+3)/(n+l) if n is odd (4.14a) a = -(n+4)/(n+2) if n is even . (4.14b) Now this point must lie on the solution curve to H(v,a) = 0 whose

asymptotes a r e a = 0 and v = 1. Consequently for v > 2, H(v,a) will have a root at a point between - 2-(n+3)/(n+l) and zero if n is odd and at a point between -(n+4)/(n+2) and zero if n is even. Since n > 3, this point will always be in the interval (-2,0) thereby proving that the sign of H(v,a), and therefore of G(v,a) must alter in the interval (-2,0) whenever v> 2. For values of v in the interval (1,2), even if the sign of (a(v,a) does not alter, it has already been shown to be incorrect. By symmetry about the midpoint of the Central range, these conclusions must also be true for negative values of the Courant number. We therefore conclude that ideal algorithms can never be stable for values of the Courant number outside the Central range.

APPENDIX

It was remarked in Section 4 that the derivation of Eq. (4.9) from Eq. (4.4) is a fairly challenging mathematical exercise. Nevertheless the author of this paper admits that the chain of reasoning set forth here may not be the quickest way of establishing the equivalence between the two formulations of 6. and the reader is invited to try and find a shorter proof if possible.

The proof here will be carried through by means of a succession of lemmas. The object will be to simplify the right-hand side of (4.4) by breaking it up into smaller constituents and to this end some new quantities will be defined as we go along.

(k > J) will be defined as s-m ... j=r-m ^' ^ Defim Lemma t i o n The < 1_ I symbol ,m <. m-m r

For values of m strictly between (s-k) and (r+k), the numbers a, K,m n

vanish (note that because k > y , the second of these two numbers is always strictly greater than the first).

Proof

By definition, T . , is non-zero only if Ijl > k. For values of m within J > ^

the stated range, the minimum value of j is > -k+1 and the maximum value is < k-1. Hence the result.

Lemma 2a

n

The number a. satisfy the following symmetry property

i\, m

n n / A o \

a. = -a. , V (A.2)

k,m k,(r+s-m) . ^ ' (Note that in the case of a normal algorithm, r+s = 0 if n is even and

Proof m-r ^k,(r+s-m) ^ "S-m j ^ - s ^ ' ^ ^ ' ^ ^ J ' ^ " ^ ' " " " ^ ^ m-i" I j l m-r ., " V r .1 ( - l ) ' ' S . k " V r + j (since n-s = - r ) S-ni 111 • " - ^ j l r - m ^ " ' ^ ^ - . k " ' — J (^^•""^j,k = - - j , k ) -a" (A.3) k,m • Definition 2 Th e function $!^ (v) (k >^) w i l l be defined as 4>|; ^(v) ={'n (v + 1 ) 1 I n (v + 1)} - i - . (A.4) •^'"1 l = r l = r + k ^^^

(In other words the product is taken over all the values of 1 for which

o? / 0 ) . Note that $[1 is a polynomial in v of degree 2(n-k)+l.

K, m K J m Lemma 2b

The function ^? (v) has the following symmetry property

K J rn $" (v) = - o " , ,(-v-r-s) , (A.5) k,m^ ' k,(r+s-m)^ ' ^ ' Proof s-k -v-r-s+lH I n f-v-r-s+lH—: -m !>!^ , ,(-v-r-s) = { 11 (-v-r-s+1)} { n (-v-r-s+1)} — k,(r+s-m)^ '.j^^^ " i^^^iJ -^^-^ (A.6)

In each bracket in (A.6), there are (n-k+1) terms being multiplied; hence the total number of multiplications being carried out in the right-hand side of (A.6) is odd and itmay be rewritten as

O k,(r+s-m) s-k s , ( - v - r - s ) = - { n (v+r+s-1)} { n (v+r+s-1)} - ^ l = r l=r+k •^'" (A.7) Let t = r + s - 1 . Then r+k $ k,(r+s-m) ( - v - r - s )

.{ n (v+t)} { n (v+t)} ^

t=s t=s-k

v+m s-k s ,

i n (v+t)} { n (v+t)} ^

t=r t=r+k

v+m

= -K (v)

_{k,m^ '}

D e f i n i t i o n 3

The polynomial PI^(V) (k >^) w i l l be defined as

Pk<W = I oJ_„ . 1 ^ (V,

Note that P1^(V) has the same degree as 4)" (v), that is, 2(n-k) +1.

Lemma 2c

The polynomial P^(v) has the following symmetry property

P|;(v) = p|;(-v-r-s). Proof Pk(-v-r-s)

I a" ^

n ^n

m=r

k,m k,(r+s-m) (v) (from lemma 2b) (A.8) (A.9) (A.10)

I

a"

m=r

k,(r+s-m) k,(r+s-m)

(v) (from lemma 2a)

Let t = r+s-m. Then

which, after reversing the order of addition»is precisely pl^(v)

(A.11)

Lemma 3

For fixed values of n and k, the polynomial pll(v) reduces to a constant.

Remark

Lemma 3 may be considered the key lemma in the proof of the equivalence

of the formulations (4,4) and (4.9) for 6,^(k>J).

Proof of Lemma 3

From Lemma 1, Eq. (A.9) may be written as

m=r m=r+k

Now

p?

is a polynomial of degree 2(n-k) + 1 in v. To prove that it is a

constant, it would suffice to show that it has a constant value at each of

the 2(n-k)+2 points (-s,-s+l... .-r-k; -s+k -r-l,-r). Further, in view

of Lemma 2c, only the (n-k+1) points (-s+k....-r-l,-r) need be taken into

consideration.

Let i be an integer such that -s+k <i <-r. From Eq. (A.4), it is clear

that <!>. „(i) will vanish unless m = -i in which case

K,m

<^|;^_i(i) = [ ( i + r ) . . . ( - 2 ) ( - l ) ] [ ( l ) ( 2 ) . . . ( i + s - k ) ] [ e + r + k ) . . . ( i + s ) ]

- I i r ( ^ ' + ' " ) f i r ^ M + c ; k^l ( i + s ) ! (A.14) - (-1) ( - i - r ) ! ( i + s - k ) ! (T^.^^.|,.i)i . s u b s t i t u t i o n of which i n t o (A.13) y i e l d s

nC.,.,[(-l)-<'*^'(-l-r)!{1«-k)! j^^^] H (-""'^J.k^^-i-r^j

S+1

I

j = r + i (A.15) But • i - r ~ ( n + i + r ) ! ( - i - r ) ! ~ (i+s) ! ( - i - r ) ! n I r^ I

Therefore s+i

P^l)=n![(-ir<^-)^|ilfz^] X (-l)m.j_,nC.,.

j = r t i '"*^ S+1 - "-Hur.,-!)'.^ j i ^ ^ . f l ) ^ J , k " ' : j - ( r . l ) - (A.17)

Now T. i s non-zero only i f | j l > k ; therefore since | r + i | = - r - i <k (since i >-s+k >-r-k), Eq. (A.17) i s equivalent to

n^(i) - „ i r . ( i + s - k ) ! I H\ , J - ( r + i ) ^

\^'> - " ' ^ i + r + k - D H j i ^ ( - ^ ^ ^ j , k " ^ j - ( r + i )

- ( ^]%^ . 2'' r ( i + s - k ) ! -, y , , > ( s + i ) - j ( j + k - p ! . (1) n! ^2kyT L ( i + r + k l ) ! J .i^^^'^' ( j k ) ! " ^ ( s + i ) j

-(A.18) With t = s+i-j, Eq. (A.18) becomes

n,.. , ,.n , 2^ r (i+s-k)! -, ^V~^ , ,.X. (s-t+i+k-1)! ^ ,. , Q ,

P„(1) = (-1) n! -^2klT[(Ur+k-l)!^ J Q ("^^ (s-t+i-k)! " S

^^'^^^

and with q = s+i-k,

Now consider the summation

Each term is a polynomial in (2k) of degree (q-t); this assumes a maximum when t=0 and consequently the degree of Q as a whole is q. But each term

[(q+2k-t-l)!]/[(2k-l)!(n-t)!] (A.22) is also a polynomial in t with degree q + 2 k - ( n + l). Recalling the combinatoric

identity (4.8), it follows that Q much vanish whenever 0<(n + l) -2k <q, that is, n - ( q - l ) <2k<n. There are precisely q integers in this range and since Q is a polynomial of degree q in (2k), these integers are precisely the roots of Q. In other words

Q = a[2k-n][2k-(n-l)]....[2k-{n-(q-l)}]

= a(2k-n+q-l)! i (2k-n-l)!

= a(i+s+k-n-l)! i (2k-n-l)!

= a(i+r+k-l)! f (2k-n-l)!

where a is a constant. Now substitution of (A.23) in (A.20) gives

p;;(i) = (-i)"(n!)^|^

-^2k4::Tyi-(A.23)

(A.24)

which is independent of i. This shows that the polynomial pP(v) has a constant value at each of the (n-k+1) points (-s+k...-r-l,-r). By Lemma 2c, it will have the same value at the (n-k+1) points (-s,-s+l.,.-r-k). Since it is a polynomial of degree 2(n-k)+l, it must therefore be a constant everywhere and the proof of Lemma 3 is complete.

Remark

n

The constant value of p. (for a fixed n & k) may be easily found by considering the case when i = -s+k. Then q = 0 and Eq. (A.20) reduces to

(A.25) n

n , ,.n 2'^'^n! ^k = ("^^ k(2k-n-l)!

Note that regardless of the value of k, p. is negative if n is odd and positive if n is even.

Lemma 4

The formulations (4.4) and (4.9) for 6, are equivalent if k >5-. Proof

Consider 6, as defined by (4.4). In view of the definition of a. , one may write

K,m

\ =

2l

:(v)Y

_-ÏÏT;

m=r

k,m

v+m (A.26)

or, from lemma 1, 2

r-n

s-k o

m=r

k,m v+m

I

m=r+k

v+m

s-k

I

m=r k,m k,m

<

JW + I

< J^)

m-

r+l^ k,m k,m r+k-1 n (v+i l=s-k+l ^ ^ ^ ( v ) \ r+k-1

2lw)^k(^)^^^.,^/-i)

(A.27)

and substitution of (A.25) in (A.27) yields (4.9). Note that in the special case when n is odd and k = -j-, the integer (s-k+1) will exceed

(by 1) the integer (r+k-1). This means that there are no terms at all in the product '^j'^"'^ (v+1), but Eq. (4.9) would still hold with ^s!|j+J(v)

l=s-k+l replaced by unity.

ACKNOWLEDGEMENTS

Many of the ideas which went into this research stemmed from spirited

discussions with Professor P.L. Roe. The author would also like to thank

the Association of Commonwealth Universities and the British Council for

financial assistance during his stay at Cranfield.

REFERENCES

ENGQUIST, B. & OSHER, S. 1981 One sided difference approximations for

nonlinear conservation laws.

Maths Comput., 36, 321-352.

GODUNOV, S.K. 1959 A difference method for the numerical calculation of

discontinuous solutions of hydrodynamic equations.

Mat Sbornik, 47, No 3, 271-306.

(Translated by US Joint Publications Research Service as JPRS 7226,

November 1960).

HALL, M.G. 1985 Private Correspondence.

ISERLES, A. & WILLIAMSON, R.A. 1984 Stability and accuracy of semi-discretized

finite difference methods.

IMA J. Num. Analysis, 4 289-307.

PIKE, J. 1984 Grid-adapted algorithms for solution of the Euler equations

on irregular grids.

RAE TR 84104.

ROE, P.L. 1981 Numerical algorithms for the linear wave equation.

RAE TR 81047.

A

B

C

X

X

X

= V

4 * 1 - 2v

Figure l.a. Regions of the ( X , v ) plane

wi^ere conditions (3.8.a) hold.

B

Equations of the curves

A :

B:

C :

D:

Figure l.b. Regions of the ( X , v ) plane

where conditions (3.8.b) hold.

X = v^

X = | ( 3 - v )

X = l ( 3 + v ' )

A : \ = ^ ( 3 + V^)

B: X = ± v '

4

C : X = V ' . 3 V ' - 8 T ; » 1

4 1-2V

Figure 1.c. Stability regime for the family of

second order algorithms spanning the four

A: X =

i-(3

+ v0

B : X = ± v'

4

C: X = l v H 5 - v ' )

Figure 2.Q. Regions of the ( X v ) plane

^ ..,K.r„ ronditions (3.15) hold.

_{where condition}

A : X

B: X

C : X

= i - ( 3 * V ^ )

Figure 2.b. The subset of So in which

condition (3.17) does not hold.

A :

B :

C :

X = ±(3*V^)

X = 1 {5v'-v')

4

^ = l{2*Sv'-v')

Figure Zc. The subset of So in which

condition (3.18) holds.

?.•'.'.•'..•'..". • l'.U...'*?

-1

B^

Equations of the curves

A: X = 4 - ( 3 + v')

4

B : X = l ( 5 v ' - V ^ )

4

Figure 2.d. Stability regime for the family of

third order algorithms spanning the five

Central range

v = 1

v = 0

Figure 3.a. Solution curves to the equation

H ( v , a ) = : 0 for odd values of n.

Central range

v = l

V = : - l

Figure 3.b. Solution curves to the equation

H(v,ot)=o for even values of n.