A frequency analysis of finite difference and finite element methods for initial value problems

R [/ithnevelQ~ i

.~ltllle dtJjctence andJ'tntle elemenl melhod~

179

A FREQUENCY ANALYSIS OF FINITE DIFFERENCE AND FINITE ELEMENT METHODS FOR INITIAL VALUE PROBL~4S* R. V~chnevetsky and F, De S c h u t t e r Rutger s U n i v e r s i t y , New Brunswick, N.J. 1. INTRODUCTION The popularity of finite element and spline methods has brought to the fore the need for criteria to compare the accuracy of computing algorithms. Several recent papers have addressed this problem, investigating the use of new analytical tools to achieve this goal (see e.g. Swartz and Nendroff (1974 a and b) Vlchnevetsky (1973), Vichnevetsky, Tu mnd Steen (1974) Fix and String

(1969)). The approach used h e r e c o n s i s t s in c o n s i d e r i n g s i m p l e partial d~fferential equations as typical models of more complex ones, and in showing that the error introduced by numerical approximations may be characterized, for those slmple models by criterla which can be related to physical parameters of those equations. The mathematics are based on the expresslon of solutions as a sum of space-slnuso~dal components for which errors can be expressed as a function of the frequency.

Finite element and B-splzne methods consist in choosing a family of basis functlons ~Om (~) and in expressing an approximate solution as the l~near combinatlon:

The c o e f f i c i e n t s 9/m. a r e i m p l i c i t l y collocation conditlons

x,e D 5 ~ ~ o

The e q u a t i o n r e s i d u a l

(1)

,A,,~ ~u~a__g_ = X'(a~) ~ a~(~)~ _ =,,~, ~(~.a~,t) c~) o b t a i n e d by some form o f f i n i t e e l e m e n t , s p l l n e or f i n i t e d i f f e r e n c e p r o c e d u r e , and then r e p l a c i n g t h e t i m e d e r i v a t i v e i n (2) bv a d i s c r e t e a p p r o x i m a t i o n i n t , e . g . t h e g e n e r a l ~ m p l l c l t method:

(,~.. _ a,,'~= e x (a~)+(,-e),x(~,,) ,6£

/

cl.i~ Ltn(j'a£)5 o ( e 4 ~

(a)

(The left hand s~de may be expressed in operator notation.

is d e f i n e d as:

a~

a s s o c z a t e d t o an a p p r o p r i a t e s e t of ~ n ~ t l a l and bound a r y c o n d i t i o n s , where ~ ( . ) i s an o p e r a t o r c o n t a i n ln g d e r i v a t i v e s w i t h r e s p e c t t o ;t only . The synt h e s i s o f an i m p o r t a n t c l a s s o f a l g o r i t h m s f o r t h e num e r i c a l s o l u t a o n o f t h i s e q u a t i o n may be viewed as a t w o - s t e p p r o c e s s , c o n s i s t i n g f ~ r s t i n r e p l a c i n g (1) by the semi-discrete approximation

~"

d e t e r m i n e d by t h e

(~ )

2. DERIVATION OF ALGORITHMS Consider the equation

;~u(z,~:)=X(u(~,',:l)~ a+_

(4)

rr'l.

See (9} below).

l?e shall assume that ~ is in some sense small with respect to z~%. This will allow us to consider the errors due to the approximation contained in the spatlal dlscretization (2) alone, and thereby compare the relative accuracy of the several competing processes used in the derlvatlon of (2). The formalism of derivation of algorithms of the form (2) whlch we intend to compare for accuracy is well known :

and the semi-discrete approximation (2) is obtained by forclng it to be orthogonM to a set of weighting functlons ~Yr~ (~t) ' :

"

d~

~- ,. j - Q

=0 --I)O~

IjQ)

.

.,

This is referred to as a welghted-residual condition, and is given Galerkln's name when the ~/,~(~)=--hOoz(pGJ. It also defines ti~7~,~)as a weak solution of (i), wlth respect to the given weighting functions 9~rn%(X) In these methods, t[w(~6,~) is unlquely defined for all ~£ . By contrast, finite difference methods are obtained by assumlng a low-order polynomlal approximation of t~ (X, h ) , different around each pelnt X ~ (say ~ ( X , ~ ) ), which is determined by simple collocatlon conditions. One then Wrltes, point by point:

,d.t which Is of the form (2) with A , . ~ = ~ ( ~ ) ( t h e identity operator). Whenever A, is not-diag~onal, the semidiscretizatlon (2) is called implicit, and explicit otherwise. 3. ACCURACY ANALYSIS One may perceive zn the theoretlcal work which surrounds the numerical approximation of partial dlfferential equations two different viewpoints. The first, whlch we call anMyti£, ms prlmarily concerned with the study of the order of accuracy of methods, or accuracy in the small, on the general premise that if errors are ~nown to decrease as some power of the mesh size, then any accuracy may theoretically be obtained by mesh re-

* Reprinted from the proceedings of the AICA International Symposlum on Computer Methods for Partial Differential Equations held at Lehigh UnlVerslty - Bethlehem, Pennsylvania, U.S.A., June 17-19, 1975.

l~lna/,t,s de l'As.sucu~tlo~l mtc~l(ttlon(de t m t u le C~d~_ul analog~que - N ° 3 - Jutlle¢ 1975

180

fznement. How fine the mesh must be chosen i s not predzcted in thzs approach. The second, whzch we shall call applied, zs the reverse viewpoint and may be stated as follows: gzven a problem, what is the maximum mesh slze which will ensure that an a-priori accuracy requirement will be met. Answering this question leads to analyzing accuracy z__nthe large. Whereas classical analyses based on Taylor series expansions are restrzcted to analyses zn the small, the "frequency method" described in the next section lends itself well to analyses in the large as well as in the small. Illustrating thls shall be one of the objectives of thls paper.

4. THE FREQUENCY METHOD We describe In thzs section the principles of a method of accuracy analysis which shall be used in evaluating the relatxve merits of the competing methods of approximation (2). It zs based on a systematic use of Fourier series, and is referred to as the frequency method. The mathematics of this method are not new s - ~ e . g . Fix and Strang [D ,l~] and Thom~e [ 14 ]). But as just remarked, the frequency method contains more information than has been exploited in previous studies. Consider the case where X (°) is a linear, constant coefficlent operator. The approxlmationX(u,)is t h e n also linear, and we may rewrzte (2) as

fl,, . cl~t~ de

where A , '

=

Az

andA2'

u.~

(s)

A,~A,,~

(E

is

the

~

by

)

~

We assume that formulae such as (2) have been multiplied by an appropriate normalizing factor such that:

_____~A I , ~

=

I

(10)

~ e normalized o p e r a t o r ~ l i s then an averaging operat o r , and the a p p r o p r i a t e c o n s t a n t f a c t o r i s i n general ~ / ~ X .To the p e r i o d i c i n i t i a l conditzon

d. (: ;~, O)

=

(14)

LL~({) = e

where ~ ( o A ) i s t h e s p e c t r a l f u n c t i o n of the operator

A ~ AT'. A ~

That is:

The approximate solution (14) also remains periodic for ~ > O , but differs from (129 in its tlme evolution: the difference between the exponents2(~) ) and ~'(aO) which describe these time evolutions is a measure of the error introduced by the semi-discrete approxzmation (S). We thus define a spectral truncation error function as:

=

I

c2o)

I t wzll be seen t h a t in the case of simple p a r a b o l i c and h y p e r b o l i c e q u a t i o n s , ~ ' < c O ) has ; h e i n t e r p r e -

classical space shift operator defined

~.~

(13)

We note that (12) remains periodic for ~ > O . The semz-dzscrete approximation (2) has, for the same initial condition (ii] the analogous solution;

~(w)

a r e the o p e r a t o r s :

X (e

X (u~) -

e

~cO~

tatzon of a relative error on the physical parameters of these equations. Strictly speaking, it is implied that the domain in~6 is the whole real axls ( - ~ , ~ ) . But the intent is of course to consider this periedzc problem as a local model o£ equations of finite ~ -domain, much in the manner zn which Von Neumann-type solutions are classically used to analyze numerical stability of difference methods. Non-linear problems are also included in the analysis, since one may consider linear equations as the model of locally iznearized forms of non-linear equations.

5. PARABOLIC CASE The simple heat equation:

"~ (11)

DX2

/k for which the spectral function >(CoO)is:

A

corresponds, f o r (1) assumed l i n e a r , t h e s o l u t i o n

(22)

(23)

J%

&oz ÷ X C~),/:

(12)

where ~ 0 0 ) i s t h e s p e c t r a l f u n c t i o n (sometimes c a l l ed symbol) of t h e o p e r a t o r X ( ' ) , defaned by:

has the exact sinusoidal solution

u.(z,~: ) =

e

(24)

R. Vtchnevetsk), Ftmtedt~feremeand fimte element methods By

analogy,

express

we

tz,(~)= =

(14)

as

!

e ;~jz~ + ~ ( ~ ) £

,'(~z.~ ~(~,c~)G- ~-~)'~'~ (2s)

with

a3z (26)

having the interpretation of a "numerical" dzffusion constant, and I~%(~).~/~) being a phase error of the approximation. It is easily shown that all symmetrical approximatlons in ~ lead to ] ~ ( ~ ( ~ ) ) = (and we limit our analysls to these cases). Thus,

o- *CW) -

2 (~)

Fiaure I: Linear and quadratic B-Splines.

(26-a)

(c) The Stormer-Numerov maximal order three point formula (see Birkhoff and Gulati, zef [ I]] d£

The spectral error function (20) becomes: =

=

|

__

C27]

(r.( tLn_'

--d% _ ~

/

+t(~+0

(34)

(d)

Quadratic and Hermite cubic fznite elements (see .~, The corresponding expressions obtained for ~((O)are listed in the following table, and are illustrated in Figure 2.

refD~,16D.

O--

and has the physical interpretation of a relative error on the di~fus~en constant of the equation introduced by its semi-discrete approximation[8). Compar~sen of Several Schemes Several schemes for' the semiqdiscretlzation of the parabolic equation ~22) are compared. These are: Ca] Standard 5 point and 5 point central difference schemes:

MBTHOD 3-Point

Linear

finxte

finite

d~fferences

elements

S~ormer-Numerov method

and S-PolnZ flni~e d~farences ~ [ ~ - - - -~.

-

-

!

_,

12 . Z ~ ~

quadratic B-Spline~

(b) Galerkin approximations w i t h lin.ear and q u a d r a t i c B-SpI ines :

aA

(32) Cubic Hermite elements

A~z

and( ~-o ~

ell/. +

z~

......

+

~

_~

finite

- ~'e,('=l'~ * ,o~ e,,~ C - , ~ z ) ) ( ' ~ , C== 4 = ) - 4

-/ d IXlri-~ ~{

/

Quadra~Tc finzte

elements

TABLE I: Spectral error function ~qr semi-di~eretizatv-~f the equation - ~ = cr- ~

182

Annale~ de I'A~so~tatton mtertmnanale pour le C a h u l analogtque - N ' 3 - Jtnllet 1975

/ 0

LL~

/

,02

'el

/

(44)

_

With a minor modiflcatzon of definition, the spectral truncatlon error is the same as (20). 7. A PROPERTY OF THE ERROR PUNCTION NEAR cO.At = 0 An important property of the error function~(cO ) which is satisfied by all the algorithms considered here is as follows:*If the seml-dzscretization ~8). is of. order of accuracy** p , then its s~ectral error function (20) and its (p-l) first derivatives, with respect to ( o A . m ~ 5 ) are zero at the origin. Example. The three point central-dlfference approximation (30) has an order of accuracy 2. The corresponding spectral error function

<

",el

~(~)=

~ (CoS~W,'X)--')

2.(~,az)z

_

[

(4s)

-.I)2

4! verifies : A

e(o):o

FiRure 27 Speotral error funotions ef different semidiseretiaation8 for the pomabolie eq~tion ~ = 0 - ~ 1: Linear finite elements. ~6 D~ 2: ~uadratie finite elements. ~: Quadratic B-Svlines. 4: 6: 8: 7:

Cubie Hermite finite elements. Stormer-Numerov approzimation. 5-point finite differences. 3-point finite differences.

6, RELATION TO THE APPROXIMATION OF AN ELLIPTIC PROBLEM The preceding analysis applies with minor modification to the elliptic problem

~'(~)

=

~zU_

d~ z

C

(4o)

approximated with the same finite difference and flnlte element technlques by

Al' ~P,-,. = A,z' u.~

(41)

where A , ' a n d A ~ ' have the same definition as before. (See e.g. Birkhoff and Gulati, ref. [ [ ] for a recent discussion of this problem.) The parallel with the initial value parabolic problem discussed here is that (41] is used from left to right (given~ , find [% ] whereas (8) is used from right to left (given ti , find ~[i/~4:). lqith/? taken sinusoidal as e i~um , the exact solution of [40) is etw ~

d (~0,~)

(46)

8. ACCURACY IN THE LARGE A zero-error algorithm would have the zero-line as its spectral error function ~(c~(at least tip to the folding frequency ~ m~=IT). The accuracy of non-perfect algorithms may thus be compared in a global sense by comparing how well their respective spectral error functions approximate that line. We know, as a corollary of the property stated in the preceding section, that an "optimal order" algorithm is one that approxlmates ~(~A)=O as well as possible near 60.~X = o , or in the "maximally flat" sense (borrowing this term from filter theory). But solutions of actual problems are superpositions of sinusoidal components of different frequencies which all contribute to the global error. In order to have a measure whlch describes the quality of the approximation not only near~OOAX)=ohut beyond, we define, somewhat loosely and with due apology, 5/ /-accuracy of the semz-dlscretzzation (8) the length of the segment of the ( ~ a X > axis between zer---o -~here ~ ( ~ ) i s always zero) ani the polntf~ax)~ beyond wh'iohl~(oJ>l becomes larger than the pOSltlve number For all practlcal purposes, we shall use ~=.ol (or one percent) in the remainder. Taking the one-percent point is somewhat arbitrary and other values would produce different results. It is however believed that the one-percent value i~ s reasonable choice for practical applications, and that other choices not drastically different would result in a similar ordering of the algorithms in terms of ~h#zr relative merits. In the fol]nwing table, the algorithms of Table | are listed in order of increasing .Ol-accuracy.

(42)

with A

X By c o n t r a ~ t h e same

Go)

=

- w:

the approximate solution ~) is

(43) g i v e n by (41) f o r

*See r e f . [ I~ ] f o r a d e t a i l e d p r o o f . ** The o r d e r o f a c c u r a c y i s d e f i n e d as t h e o r d e r i n of the first term which does not cancel out in Taylor series developemnts of left and right hand sides of (8) in powers o f ~ when an exact solution is substituted for ~ (with A, normalized)

R Vtdmevet~hy

Ftnne d~j/eteme and/istte elemenlmethud~

Eauation Method I~ Text Number 3-Point Finite Differences (30) .... Linear Finite Elements (32) (B-Spline) Quadratic Finite Elementa 5-Point Yinit~ Differences (32) Stormer-Numerev (34) Quadratic B-SpZine~ (33) Hermite Cubic Finite Elements

I~3

(o_),~X),o,

ii. HYPERBOLIC CASE We consider the simple hyperbolic equation

.35 .35

-~

%-- =

o

(s0)

.85

1.00 1.2~ 1,48

for which

3,15

It has the exact sinusoldal solution

TABLE 2: Orderin~ of semi-disoretizationa of the equation D ~ / ~ = O - ~ z ~ / ~ ~ in order of increasing "accuracy in the large"

tL(;<~)=

e ~°3

(;o-c/z>

(s2)

the parameter c has the physical interpretation of a

velocity (the solution is displaced without d i s t o r t i o n 9. DISCUSSION A comparison of the several three point discretizatiens is interesting: Somewhat unexpectedly, the linear finite-element algorithm (32) has about the same .01accuracy as the simple 3 point central difference formula (30). The Stormer-Numerov algorithm (36) which has maximum order of accuracy is, as expected, reasonably better than either of the preceding two.

at the velocity C in the ~t directlen). Solutions of the semi-discrete approximation of (8) may be written, separatlng real and imaginary parts o£

~(~3

=

e

e (s3) :s:~

where

C ~ (cO) =

(~ ( o , ) ) CO

10. A "BEST" THREE-POINT DISCRETIZATION The three-poin~ semx-discretizatlon

(,$4)

F d u - , . - ~ s,~6 d d ~ . t. , ~ . a . , 7

= 0-- LI"tZn-~ -~z2U'~ + Ut ~.,i,] obtained entirely empirically by seeking the coefficients of A 1 which maximize ~o~),ot gives:

( a b ~ X ) . o I =. I . B 4

C48)

which is about 50 percent more than the .01-accuracy of the Stormer-Numerov method. While the latter approximates the line ~{00)=O optimally in the Mac Laurzn or "maximally flat" sense, the "best" method (47) approximates this llne optimally in the Tchebychev sense, and results in a better distribution of errors. The two ~(6o) curves are compared below:

I ~'(
\

/ I.

(53) zs the expression o£ sinusoidal solution with phase velocity C ~ J a n d wmth (non-conservative) amplltude exp(.~e(~(uJ))~A~. For a l l symmetric approximations, however, ~ e (~'(dJ))= 0, and

tO ,,~X.

2..

ti n

(~) =

e_.

(s s )

Then,

- -

--

c

(s6)

becomes simply the relative error en the phase velocity introduced by the semi-discretizatzon. (Again, we llmit ourselves te these syrmnetr~c cases here). Within a ratio of ~ T , the definition (56) of the eTro~" is identical to that used for hyperbolic equation by Krelss and Oliger [4] and by Swartz and Wendroff [IS] . ThefT[factor stems from the fact that we use velocity" error, rather than the phase-error-per-period which they use. Comparison of Several S~hemes Ne compare several semi-discret~zations of the hyperbolic equation [SO). These are (a) Standard 3 point and 5 point central difference approximations

d'--~" = - C \--

z.Ex/

>~

(s7)

and -~-C " = -

le_ , , ~

(ss)

OI

(b] Galerkzn approximation with linear and q u a d r a t i c B-sphr~¢~ (sg) =_

Fi~: Speotral er~er functions of the StormerIVume~ov approximation(1) and the "best" approximation Q for the pemabo~e equation ~ 7 ~ = d - ~ / ~ i , ~Tt

C ('-uc--,

+ "~+,)

2 z~JC

I~4

4 nnule~ de l', s ~ o ¢

qt =_c(.-u.-~_-,ou~-,

~t

. . . .

e4 (c]

l(HIon

~nlermmomde

lm{U

le C(~&ld anato~tql~e . ~

J.

JuHl¢~

1075

-~7--/

. , o ~ . ,

~

u~.2)_

~X

Quadratxc and H e r ~ i t e Cubic finite-elements (see

The correspondin~ expressions obtained £ o r ~ ( u O ) are l~sted in the fol!u~ing table, and are illustrated zn F~gure 4.

@Cw] =,(:~C~ HZTtlOD ~-Poin~

f~ni~e

C.. d~ffercacea

~6~--

,

-- --

) --

Linear £1nlt~ =lemen~s

Box merhod

(~

S-Point finite d~ffe~en~es

•

~,~Z

~ : Spectral error-functions of d~fferent semiscretizations for the hyperbolic equation B_u +C ~ _ O

- ~

I: 2: 3: 4: 5: 6:

Quadrate= B-Sphne3

Cubic Harmlte £1nlte elements

3-point finite differences, 5-point finite differences. Linear finite elements. Quadratic B-Spline~ Quadratic finite elements. C~tbio Hermite finite elements.

~ a d r s t l c £Z~Z~e ~lements

TABLE ~o~f

3:

Spectral error function for semi-discretizsthe equatio~ ~LL ~t£ .~.~ + c

.01

/

1 2. COMPARISON IN TERMS OP ACCURACY IN THE LARGE As in the parabolic case, one may characterize accuracy in the large for the semi-discre~izatx0n~o£ table by the numbers ((O.~X).O~ The resulting comparison and ordering of the competing schemes analysed is shown in table 4

Me vn'~o~ ~" I~Eauati°nTex~Number Z-Point F i n i t e ~ r e n c e s (6P) 5-Point Finite Differences (58) Linear Finite Elements (59) Quadratic B-Spline~ (60) Hermite ~abi~ Finite Elements Quadratic Finite Elements

\

/

~-~-~o

/

(~3,~),o~' .75 I. 12 1. ? 2.8

-'04

\ I

I. 98

TABLE 4: Ordering of semi-discretization8 of the equatzon 70 ./~£ + C ~cL/~2C = 0 zn order of zncreasing "aecuracy in the lar-geu

Figure 6: Spectral error functions of the lineom fin~-~-~ment ~ and the "best" ~ppro~imation @ for the hyperboli~equation B~/~, +C ~'a/~ = 0 / ot

{ dTC

R l, tahneret~kv

185

Fmtte difference and/intIe element method~

13. "BEST" THREE-POINTS DISCRETIZATION Maximizing (~)~x).m , as in section I0 , produces in this case the three point seml-dlscretzzatlon'

for which

16.

REFERLNChS

[I]

BIRKHOFF, G. and GULATI, S. (1974) "Optimal fewpoint dlseretizations of linear source problems" sIAM J. Numer. Anal., September.

[2]

COLLATZ, L. 61960) "The numerlcal treatment of differential equations" Sprlnger-Verlag.

[3]

FIX, G. and STRANG, G. "Fourler analysis of the finlte element method in Ritz-Galerkin theory" Studies in Appl. Math. Vol. 48, pp. 265-273.

[4]

KR[ISS, H. and OLIGER, J. (1972) "ComparJ~on of accurate methods for the integration of hyperbelle equations" Tellus, 24 (1972).

[5]

LAX, P.D. (1957) "Hyperbollc systems of conservatlon laws If" Comm. Pure and Appl. Math., Vol. I0, pp. 537-566.

[6]

LAX, P.D. and NENDROFF, B. (1960) "Systems of conservations laws" Comm. Pure and Appl. Math., Vol. XIII, pp. 217-237.

[7]

LAX, P.D. and NENDROFF, B. (1964) "Difference schemes for hyperbolic equatlons wlth hlgh order of accuracy" Comm. Pure and Appl. Math., Vol. 17, p. 381.

[8]

RICHT~ER, R. D. and MORTON, K. W. Difference Methods for Initial Value Problems Interscience Publishers, New York, 1967.

[9]

SCHOENBERG, I. J. (1946) "Contributions to the problem of approxlmatlon of equidistant data by analytlc functions, Parts A and B" Quart. Appl. Math., !(1946) pp. 45-99, and pp. 112141.

(w ~x),o , = I, 6~ whlch, again, gives an accuracy In the large which is about 50 percent larger than for the 3-point linear finite element discretization(59) which is of maximal order in thls case. The spectral error functions of those two 3-point algorithms are shown for comparison in flg.

14. A NUMERICAL ILLUSTRATION A concrete comparison of the three semi-discretizations (57), (59) and (62) has been obtained by computing numerically step responses uslng tbese schemes. The results illustrated in fig6 clearly show the improved quality of these numerical results as (~O~X).el increases (From(CO'~X).olm'?~inthe 3-point finite difference case to 1.66 in the "best" case.)

ILL .

..

-' ",,

[10] SCHOENBERG, I. J. (1973) "Cardlnal spline interpolatlon" Regional Conference Serles in Applied Mathematics SIAM, Vol. 12, 1973.

'6-'\

[ii] SCHULTZ, M. H. (1973) Splino Analysis, Prentice Hall, Inc., Englewood Cliffs, N. J. [12] STRANG, G. and FIX, G.J. (1973) An Analysis of the Finite Element Method, Prentice-Hall, Inc., Englewood Cliffs, N. J.

2,

4,

~.

Fi~ 6: Solution of the equation 3 ~ / ~ +~ ~Cll/~ =O by different 3-Potnt semi-discretizations. I.

Finite dlfferenees.

2. Linear finite elements. 3.

'Zest" approximation.

4.

Exact solution.

5.

Initial condition.

15, ACKNOWLEDGEMENTS : Part of thls work has been supported by the National Science Foundation of Belgzum and by the U.S. Dept. of Interlor, Office of Water Resources Research, Grant No. OWRR A-045 N.J.

[13] SWARTZ, B. and WENDROPF, B. (1974) "The relative efficiency of finlte difference and finite element methods. I: hyperbolic problems and splines" SIAM J. Num. Anal., October. [14] THOM£E, V. and WENDROFF, B. [1974) "Convergence estimates for Galerkln methods for variable coefficient initial value problems" SlAM J. Num. Anal., Vol. ii, No. 5, October. [15] VICHNEVETSKY, R. [1973) "Physical criteria in computer methods for partial differential equations" Proceedings, 7th AICA International Congress, Prague i~73. Reprinted in Prec. of AICA, Vol. XVI, No~ i, Jan. 1974, Brussels. [16] VICHNEVETSKY, R. (1975) "Introduction to Finite Element Methods for Initial Value Problems" Course Notes, Rutgers University, Dept. of Computer Science, New Brunswick, N. J. [17] VICHNEVETSKY, R., TU, K. W. and STEEN, J. A. (1974) .,Quantltative error analysis of numerical methods for partial differential equations" Proceedings, Eighth Annual Princetor Conference on Information Science and Systems. Princeton University, March 1974. [18] WENDROFF, B. (1960) "On central difference equations for hyperbolic systems" J. Soc. Indust. Appl. Math., Vol. 8, p. 549.