Numerical Solution of Parabolic Equations by the Box Scheme Citation Fong, Kirby William (1973) Numerical Solution of Parabolic Equations by the Box Scheme. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/SV2E-JZ49. https://resolver.caltech.edu/CaltechTHESIS:04012013-103558891 Abstract The box scheme proposed by H. B. Keller is a numerical method for solving parabolic partial differential equations. We give a convergence proof of this scheme for the heat equation, for a linear parabolic system, and for a class of nonlinear parabolic equations. Von Neumann stability is shown to hold for the box scheme combined with the method of fractional steps to solve the two-dimensional heat equation. Computations were performed on Burgers' equation with three different initial conditions, and Richardson extrapolation is shown to be effective. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: (Applied Mathematics) Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Applied Mathematics Thesis Availability: Public (worldwide access) Research Advisor(s): Keller, Herbert Bishop Thesis Committee: Unknown, Unknown Defense Date: 22 September 1972 Record Number: CaltechTHESIS:04012013-103558891 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:04012013-103558891 DOI: 10.7907/SV2E-JZ49 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 7571 Collection: CaltechTHESIS Deposited By: Dan Anguka Deposited On: 01 Apr 2013 17:49 Last Modified: 16 Jul 2024 21:21 Thesis Files Preview PDF - Final Version See Usage Policy. 19MB Repository Staff Only: item control page

NUMERICAL SOLUTION OF PARABOLIC EQUATIONS BY THE BOX SCHEME Thesis by Kirby William Fong In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute Pasadena, of Technology California 197 3 (Submitted September 22, 1972) -ii- ··~ ACKNOWLEDGMENTS The author wishes to thank Professor Herbert B. Keller for suggesting the problems considered in this thesis and for valuable discussions during the course of the research. He wishes to thank also the National Science Foundation for fellowship support and the California Institute of Technology for graduate research and tea-c hing assistantships. He is also thankful to Mrs. Alrae Tingley and Mrs. Virginia Conner for their expeditious and meticulous typing of his manuscript. -111- ABSTRACT The box scheme proposed by H. B. Keller is a numerical method for solving parabolic partial differential equations. We give a convergence proof of this scheme for the heat equation, for a linear parabolic system, and for a class of nonlinear parabolic equations. Von Neumann stability is shown to hold for the box scheme combined with the method of fractional steps to solve the two-dimensional heat equation. Computations were performed on Burgers' equation with three different initial conditions, and Richardson extrapolation is shown to be effective. -iv- TABLE OF CONTENTS Introduction Chapter 1 1 The Box Scheme in One Space Dimension I.l A Basic Description of the Box Scheme 4 !.2 An Energy Inequality for the Heat Equation 10 I.3 Convergence of the Box Scheme for the Heat Equation 17 I.4 An Energy Inequality for a Linear Parabolic System 28 I.S Convergence of the Box Scheme for a Linear Parabolic System 34 I.6 Solving the Linear Finite Difference Equations 44 I.7 Nonlinear Parabolic Equations 54 Chapter 2 The Two-Dimensional Heat Equation II. 1 The Method of Fractional Steps 69 II. 2 The Box Scheme and the Method of Fractional Steps 73 Von Neumann Stability 75 II. 3 -v- TABLE Chapter 3 OF CONTENTS (Continued) Burgers' Equation: Computational Examples III. 1 Burgers' Equation 80 III. 2 Example 1: Smoothing of a Sharp Front 84 III. 3 Example 2: Formation of a Shock 94 III. 4 Example 3: Interaction of Two Shocks 116 References 128 - 1- INTRODUCTION This the sis deals with the numerical solution of parabolic equations by the box scheme. Chapter I is devoted to the analysis of prob- lems in one space dimension. We begin with a de scription of the box scheme and list some situations in which it would be preferable to other methods of computation. for the box scheme. We then give three convergence proofs In each proof we use discrete analogues of energy inequalities to show that the finite difference solutions are accurate approximations of the continuous solutions. Energy inequalities are generally used to prove uniqueness of solutions of initial value problems; however, in our work we have used modified forms for initial boundary value problems with inhomogeneous terms. Section I. 2 gives the derivation of such an energy inequality for the heat equation. In Section I. 3 we show how the energy inequality can be used as a model for. finite difference equations. We then generalize the convergence proof for the heat equation to a linear parabolic system. Section I. 4 gives a derivation of an energy inequality, and Section I. 5 shows how the energy inequality may be discretized with some complications to - 2 - yield a convergence proof of the box scheme for parabolic systems. The emphasis of Sections I. 6 and I. 7 is on the computation of the finite diffe renee solution. In Section I. 6 we prove that an upper and lower block triangular matrix factorization may be used to solve the finite difference equations for a special linear equation. Using an argument involving principal error functions, we also show how to resolve a problem about the "smoothness" of the finite difference solution that arises in Section I. 5. Finally in Section I. 7 we give a constructive proof that the nonlinear difference equations resulting from applying the box scheme to a particular class of nonlinear parabolic equations have a unique solution. The mean value theorem enables us to adapt the convergence proof for linear systems to this class of nonlinear equations. In Chapter II we give an example of how the box scheme may be extended to solving the heat equation in two space dimensions by using the method of fractional steps. We show that the initial value problem with periodic data leads to a numerical scheme which is stable in the sense of von Neumann. With this type of problem our numerical scheme is consistent to second order accuracy so that the numerical solution will converge to the continuous solution with second order accuracy. - 3 - In Chapter III we present the results of computations on Burgers' equation with three different initial conditions. We sought to identify the formation of a shock numerically -that is, without looking at a graph of the solution- but have not obtained a conclusive result. We have also performed Richardson extrapolations with solutions from successively refined nets and have found empirical conditions under which the extrapolations appear to be most effective for producing more accurate solutions. - 4- CHAPTER I THE BOX SCHEME IN ONE SPACE DIMENSION I. 1 A Basic Description of the Box Scheme The box scheme for the numerical solution of parabolic equations was originally proposed by Keller [1971]. Since this scheme is of fundamental importance in this thesis, we present here a brief review of the method and indicate some of the ways in which it is superior to other numerical methods for solving parabolic problems. We consider the following special problem defined in the rectangle 0 ~x ~ 1 and 0 ~ t ~~ = a: ~ T: (a (x) ;~) + c (x)U + S (x, t), U(x, 0) = g(x), (l.la) (1. lb) (l.lc) a 1 U(l, t) + [3 1 a(l )U (1, t) = g 1 (t). X (1. 1 d) An important step in the method is to reformulate the problem as a first order system of equations: - 5- au ox = V, a(x) (1. 2a) au - c (x) U - S (x, t ), = at av ox (1. 2b) U (x, 0) = g (x), (1. 2c) = a(x) d~~) , (1. 2d) a 0 U(O, t) - {3 0 V(O, t) = g 0 (t), (1. 2e) a 1 U(l, t) + {3 1 V(l, t) = g 1 (t). ( 1. 2£) V(x, 0) We define a mesh over the rectangle: (1. 3a) < tN = T. (1. 3b) The mesh spacings are then defined by h. J - X. - X. J' JJ (1. 3c) k n - t n - t n-1' (1. 3d) for j = 2, • J and n = 2, • • • N. For net functions { cp~} and coordinates of the net we use the following notation: x.±1 J 2 (1. 4a) t n ±1 - (1. 4b) 2 n cp. ± 1 J 2 (1. 4c) -6n ±.l cp . 2 t( 2 cp.n + cp.n±l) 1 J . J - J n D X n cp. J (1. 4d) n cp . - cp. 1 1 ]- - (l.4e) h. J n J n-1 J cp . - cp . n Dt cp. J - For functions k (1. 4f) n \jt{x, t) defined everywhere in the rectangle we use the notation t~ - \jt(x.1 t n J 2 - t (x. ± 1 t ) 1 J 21 n (1. 4h) t (x.1 t ± 1 ). J n 2 ( l. 4i) J \jr. ± 1 . J n (1. 4g) )1 1 t~±2 J The box scheme for the numerical approximation of problem (1. 2) is [u~} and [ v~} with all the difJ J ference approximations centered in the middle of the box given in terms of the net functions [x. • x.] J- 1 J X [t n- 1 •t n ] or on an appropriate edge of the box when coefficients do not depend on the time variable. We have n a . 1 D u. J-2 x J D for 2 ~ j X (1. 5a) 1 n-2 v. (1. 5b) J ~ J and 2 ~ n ~ N. The initial data are taken as - 7u l. J 1 v. J for 1 ~j • = g (x.) J = a.J ~ J, n d g (x . ) J dx (1. Sd) and the boundary conditions become n for (1. Sc) ao ul - i3o vl n n = go • (1. Se) n n Q'l UJ + f}l VJ n = gl • ( l. 5£) ~ 2 ~ N. We immediately see two advantages of this method over other numerical methods which have been proposed. First the mesh spacings need not be uniform so that we may use a finer net in regions where we expect rapid changes in the solution. Second, the scheme is well adapted for problems in which a(x) is not continuous. For instance in a diffusion problem where a(x), the diffusivity of the medium, is dis continuous, flux ~~ .will also be dis continuous, but the a(x) aa~ which is one of the dependent variables in the box scheme will be continuous so that we need not make any special modifications for discontinuous coefficients other than to pick the points of discontinuity to be mesh points. There are other desirable features which are not so apparent. Being implicit, the method will be unconditionally stable. has second order accuracy even with nonuniform nets. The method Richardson extrapolation is valid if the continuous solutions are sufficiently - 8 differentiable and yields an improvement of two orders of accuracy for each extrapolation. Both to the same order of accuracy. U(x, t) and oU(x, t) ox are approximated Although the box scheme requires a little more computation than the Crank-Nicolson scheme, it will nevertheless be preferable for the types of problems described in the preceding paragraph. All of the virtues cited for the box scheme are discussed by Keller [1971]; however, he does not give the complete convergence proof. Subsequently Varah [ 1971] presented a general stability result for difference approximations to mixed initial boundary value problems for parabolic systems and included as an example the box scheme. This result uses Fourier transforms in x net spacings to be uniform. and t and requires the Once a finite difference scheme has been shown to be stable, we need only check its consistency with a particular problem to show its convergence to the solution of that problem. That is to say, stability is a property of a difference scheme only and does not refer to any particular parabolic problem. Con- sistency on the other hand refers to a specific problem and tells whether or not the difference equations accurately approximate the differential equation and boundary conditions as the mesh is refined. Stability and consistency together imply convergence of the finite difference solution to the continuous solution as the mesh is refined. Obviously c onvergence is the behavior we seek when we c ompute approximate solutions to differential equations, and we would like to have an a priori guarantee of c onvergence whenever possible. Our - 9goal is to enlarge the class of parabolic problems for which we can guarantee convergence of the box scheme. We shall give a different proof for the convergence of the box scheme for the heat equation and then show how it can be extended to a parabolic system and to a class of nonlinear parabolic equations. There will be some mild restraints on the net spacing, but basically we will be allowing nonuniform time and space steps. Implementation of the box scheme will entail the solution of linear systems of algebraic equations. While stability or convergence will imply that the systems are nonsingular and have unique solutions, it still remains for us to chaos e an appropriate algorithm for obtaining these solutions. In general an algorithm for solving a linear system will require that further conditions be satisfied in addition to nonsingularity before we can prove that it will produce the desired solution. We shall use the method of factorization of block tridiagonal matrices recommended by Keller [1971]; however, we shall give an alternative proof based on an analysis suggested by Varah [ 1972] that this algorithm will work. - 1 0- I. 2 An Energy Inequality for the Heat Equation The convergence proofs which we shall present are based on energy inequalities. Energy inequalities are often used to prove well- pos ednes s of initial value problems or to show that a solution depends continuously on the initial data of an initial value problem. It is frequently possible to construct a dis crete analogue of an energy inequality which can then be used to prove convergence of a finite difference scheme. Indeed, this is how we shall obtain our conver- gence proofs, and to this end we wish to study thoroughly a simple parabolic equation - namely the heat equation - beginning with an energy inequality. We consider the following problem: v = u vX = ut • (2. 1 b) U(x, 0) = g (x) • (2. 1 c) V (x, 0) = d g(x) dx (2. 1d) a 0 U(O, t) - ~ 0 V(O, t) = go (t) , (2. 1 e) a 1 U(1, t) + ~ 1 V(1, t) = g 1 (t), (2. 1 f) 0 > x' (2. l a) (2. 1 g) - 1 1- o. (2. lh) The conditions (2. lg) and (2. lh) are physically natural requirements for a mixed boundary condition. If these quantities are for some reason less than zero, our energy inequality would contain integrals along the boundaries x = 0 and x = l. With the present assumptions these integrals will have signs such that they can be dropped from the inequalities we will derive. where {3 0 or {3 1 The c ase of Dirichlet boundary conditions is zero is actually simpler than the mixed case and would require only minor modifications of the proof. We therefore consider only the mixed cas e. Typically, an energy inequality argument is used when one wishes to show problem (2. 1) has a unique solution. U and V are one solution pair and u and v If we assume are another solution pair and define the difference functions e (x, t) - U (x, t) - u (x, t) , (2. 2a) f(x, t) - V(x, t) - v(x, t), (2. 2b) and £ are solutions of we find that e = e = et ' f f X X e (x, 0) = 0 ' (2. 3a) (2. 3b) (2. 3 c) - 12 f(x, 0) = 0, (2. 3d) ao e(O, t) - f3 0 f(O, t) = 0, (2. 3e) al e (l, t) + f3 1 f(l, t) = 0, (2. 3f) ao ~ ;e: 0 , (2. 3g) :2: 0. (2. 3h) al ~ One then notes that the integrals with respect to x of the squares of e and f are zero at time zero and must be non-increasing functions of time. Since these square integrals are non-negative, they must be zero; therefore, e and pairs must be identical. f are zero for all time, and the solution This, however, is not the manner in which we wish to use the energy inequality. u and v satisfy (2. I) but that conditions (2. le) and (2. 1£). tions governing U and V satisfy only the boundary Then we would have the following rela- and f: X = f + p , (2. 4a) X = et + a , (2. 4b) e f e Let us suppose instead that ao e (0, t) - f3o f(O, t) = 0, (2. 4c) al e(l, t) + f3 1 f(l, t) = 0, (2. 4d) ao ro >. 0 , (2. 4e) - 13 al ;;;:: 0 p - u - vI (2. 4g) (J - v X - ut. (2. 4h) p and a account for the fact that 13; The terms (2. 4f) I X U necessarily satisfy the differential equations. and V do not Before proceeding with the derivation we should like to furnish some motivation by saying and v that in the finite difference problem u finite difference equations while U and the continuous system we are modeling. be the computed approximations to V will be solutions of the will be the solutions of That is, U and V. u and v are to Since the difference equations are only approximations to the differential equations, we cannot expect U and V, the solutions of the differential equations, to satisfy the difference equations exactly. The extent to which they fail to satisfy the difference equations is embodied in the truncation errors p and o. The terms p .and a are generally small. ically, for the box scheme, they are 0 (h 2 ). Usually zero at time zero, but we shall not assume this. e and Specif- f are With this interpre- tation of the various terms in (2. 4 ), we see that the results we need are bounds on e and £ at times greater than zero in terms of and the initial values of e result because p and and f. p, a, This is called a convergence a can be made arbitrarily small by taking a sufficiently fine mesh and because the initial error also can presumably be made arbitrarily small. Then sin c e e and f go to zero as - 1 4the mesh is refined, the finite difference solution must converge to the continuous solution. Our plan is to derive an energy inequality from (2. 4) and then try to duplicate the derivation for a discrete system. It turns out that it is convenient to take the time derivative of (2. 4a): (2. 4i) If we multiply (2. 4i) by f and (2. 4b) by e. add the products, inte- grate with respect to x from 0 to 1, and make a substitution using l (2. 4a), we obtain, with (cp, 'f) = / cp(x)t(x)dx and 0 1 d 2 2 dt lie d 1 2 == (rp, rp ) , ./ 2 II + 2 dt llfll 1 II cp II (e. a) J (2. 5) l + [f et]o - (fx• et) + [e f]o - (f, f) - (p, f). Further substitutions, a careful examination of the boundary conditions, a time integration from 0 to t, and the Schwarz inequality lead us to t ~ c + 2 I t II f II ·II Pt II + II e II • IIa II + II f 0 X II ·llall (2. 6) +llfll·ll Pll }ds, where (2. 7) - l 5We recognize that on the left side of (2. 6) we may c hange t 0 ~ T ~ t. without destroying the inequality provided grate both sides with respect to T from 0 to t. to T Then we inte- Also the four t erms on the right side of (2. 6) may be separated using the generalized ab ~ ·hea 2 + .!_ b 2 } where arithmetic-geometric inequality arbitrary positive number. the other two we take € t t 0 0 € In two cases we take = (4t)- 1• e is an -1 e = (2t) , and in We then find t f f lle ll 2 ds - if II£ ll2 ds ~ C t + 2 t 2 f llol l2 ds 0 X + 4t (2. 8) t t 2 2 2 2 f II PI I ds + 4t f lloll ds. 0 0 Starting again with (2. 6) we find a judicious use of the generalized arithmetic -geometric inequality gives t ll e(t) ll 2 + ll £(t) ll ~ C + {i / ll e ll 2 ds 2 . 0 t t t - i f II£ ll2 ds} + f 11 Ptl l2 ds + -38 f l! o ll 2 ds 0 0 X (2. 9 ) 0 The final energy inequality results from substituting (2. 8) into the braces in (2. 9): lle(t)l l 2 + ll £(t) ll t ~ C(l+t) + f 11 Pt ll ds 2 2 0 2 t +(~+6t )/ ll ol l ds + (l+4t ) / ·o 2 2 0 (2. I 0) t II P II 2 ds. - 16 If we restrict t to lie between 0 and T, we see that the three terms involving the truncation errors can be bounded independently of t. Referring to (2. 7) we note that and c an be made small. C is determined by the initial errors Hence the energy inequality (2. I 0) is in the form we need for a convergence proof of the difference equations. The derivation given in this s ection was c; uggested by Lees [1960] and Lees [1961]. - 17 I. 3 Convergence of the Box Scheme for the Heat Equation Let [u~} and [ v~} be net functions which we shall use to J J approximate U and V respectively. The box scheme for the heat equation then takes the form n v. D - u.n = 1 J-z (3. la) J X 1 D- v~-z = Dt- u.n 1, X J J-z (3. 1 b) u.l J = g (x. ), (3. 1 c) v l. J = dg (x.) J dx (3. ld) J n n n "'o ul - l3o vl = go , (3. 1 e) n n Q'l UJ + j3l VJ = gln (3. If) "'o ~ ~ 0, (3. 1 g) ~ 0 . (3. lh) Q'l 13; Let U and V be the solutions of (2. 1) and define the error net functions n e. J - n U(x., t ) - u. J n J ' (3. Za) f.l - n V(x ., t ) - v. J n J (3. Zb) J - 18 - [e~ } and J t£~} are then solutions of J D- e~ = J X f~ I + p~ 1 , (3. 3a) J -z- J -z 1 = Dt- e.Jn-z-1+ 0 Jn-z. r, -z- (3. 3 b) el: = 0 , (3. 3 c) fl: = 0 , (3. 3d) J J (3. 3e) 0:'1 n eJ + f31 fnJ = 0 ' (3. 3f) ;;:: 0 ' (3. 3 g) ;;:: 0 (3. 3h) O:'o ~ 0:'1 K where the local truncation errors are defined by p~ 1 J -z- _ = { D U(x.j' t n ) X + {v(x. 1, t J -z- 1 n-zI J -z- a. n 8U(x . _.!_, t ) } J z n 8x ) - t[V(x., t ) + V(x. 1 , t J n Jn ={ t D-x [V(x.,J t n ) + V(x.,J t n- 1 )] - t Dt- [U(x ., t ) + U(x . l' t J n Jn )J} . (3. 3i) >J}, (3. 3 j) - 1 9If we apply the operator n; to (3. 3a), we obtain (3.3k) 1 where n-z( . is defined by 1 J- z 1 ( n- z1 J. -z- - [ vt (x.J --;-1, t n- -;-1 ) c. Dt- V.n 1 } J-z c. (3. 3 J,) + [ Dt- D- U~ - U t (x . X J X 1, t J -z- n-z1 ) } • It is possible to give alternate expressions for the truncation errors if we use Taylor expansions in the above definitions. however, we wish to introduce some new notation. and for some fixed r > 0 we assume max k n n = rh. Let h Let First, = max h. , j J 8 (x) and cp (t) be piecewise continuous functions such that for some fixed 5 > 0 we have h. J = ~ k n 6 9 (x . 1 )h, J-z- e(x) ~ 1 , = cp (tn-z-!)h, ~ cr:(t) ~ r , O ~ x ~ l, } 2 ~ n ~ N.} (3. 4b) O ~ t ~ T. Proceeding as in Keller [19 71] we assume continuous derivatives of order net points. (3. 4a) Then if 2m + 2 ~ M, M U and V have piecewise where any jumps must occur at the we c an show that - 2 0- ~ (h2 ) = 'J=l I n--z a. 1 J-2 2\J 2 2 [U,V;x. 1,t }+O(h m+ ), \! J--z n (3. Sa) m (h f ~ - ) \J S 1..U,V;x. 1,t 1 } +O(h 2m+2 ), = \J=l 2 \! J --z n--z (3. Sb) 2 m (h ) \J 2m+ 2 ~ Z [ U, V ;x. 1, t 1 } + 0 (h ), 2 \!=I \! J --z n-z- (3. s c) R 2 = where 2 2 R [U, V;x, t} \! \J( ) { I a \J+ I U(x, t) = e__ x • ( \J)! 2\J+I • ax2'J+I (3. Sd) 2 2\J _a vJ~· t)} . ax S [ U,V;x,t)} \! I 2/l 2\)- 21-L \! cp (t) .8 (x) = ~ 1-L=O (2/l)! (2\J- 2J.L)! 2 a \J+I U(x, t) ax2'J-2/l at2/l+I \) Z\! [V, V ;x, t} _ ~ 2\)- 2/l cp (t) (3. Se) ~ ) ' 8 2/l (x) (3. Sf) 1-L=O (2/l)! (2'J-2/l+l)! For the purpose of this section the most important feature of the truncation errors is that all three errors are O(h 2 ). We next introduce an inner product for net functions and [ \fr~}: J - 2 l - J :6 j =2 n CfJ . I J -z n I h .. J -z J \jr. (3. 6a) If a net function is differenced, we will have J :6 cp.n 1. D - \jr.n h .. j =2 J- 2 X J (3. 6 b) J The norm associated with this inner product is n n (c:p ' cp )h. (3. 6 c) We note that (3. 6c) is actually a seminorm rather than a norm since a net function which oscillates along the mesh can have norm zero without itself being zero. convergence proof. We shall say more about this after the Finally, with our inner product the following identities hold: (3. ?a) (3. ?b) As mentioned earlier our plan is to construct a dis crete analogue of the energy inequality derived for the differential equations. The first step then is to construct the energy quantity exactly as was done for (2. 5): - 22 1 1 II e nll h + 2 D t ll 1...nll h = - (1...n.- 2 , Cn- 2 )h 2 1 2 D t- 1 1 2 - (3. 8) 1 n-2 n-2) - (e , a h Beyond equation (3. 8) the discrete nature of the variables causes some difficulties which did not occur before. We therefore intro- duce new quantities which will help us notationally: - k2' (3. 9a) llf2 llh - 0 , (3. 9b) kl 1 1 IIC 2 11h - 0 , (3. 9c) 1 li n~ f2 llh =o, (3. 9d) 1 llo 2 llh - 0 , (3. 9e) 1 2 IIP IIh s - 0 , (3. 9f) (3. 9g) m (3. 9h) - 2 3- This notation plus additional substitutions, careful examination of the boundary conditions, a time summation from t 1 to t , and the Schwarz n ·inequality lead to the analogue of (2. 6): (3. 1 0) 1 n + ll fm- 2 11 h ) ~ C + 2 :6 2 m=l k m • where (3. II) We see that on the left side of (3. 10) the index n may be changed to where 1 ~ i ~ n yielding a set of valid inequalities. both sides by k. and sum from i = 1 to i = n. 1 We then multiply We again apply the arithmetic -geometric inequality to each of the four products on the right side of (3. 10). With the notation D i = t n + k 1 we now obtain -24(3. 12) n + 2 D 2 .6 k m=l m s 2 m (3. 12) corresponds with (2. 8) but has additional terms in .!. 2 ll fm- 2 ll h ll fml l ~ and because a c ancellation which occurred in the continuous case does not occur for discrete equations. Returning to (3. 10) we use the arithmetic -geometric inequality with different parameters to deduce (3. 13) n Inequality (3. 12) is still valid if we omit the term from the left side. 2 .6m=l km l l~ l lh We then substitute the remaining inequality into the braces in (3 . 13 ). With the notation - 2 5- (3. 14a) max l ~m ~n T 2 m , (3. 14b) we conclude the result (3. 15) + [(t +k1 ) n 3 2 + t +k 1 ] T (n) . n (3. 15) is the convergence result we sought. It says that the errors at a given fixed time may be made small if the initial errors are small and if the truncation errors are small. The latter error we noticed earlier was O{h 2 ) and can be made smaller by refining the mesh. We note also that km+l /km is bounded by r I 6; hence, no further conditions are needed to guarantee that are approximated to O{h 2 ) T (n) is 0 {h 2 ). or better, then we see If the initial data llen llh and llfl llh are also O{h 2 ). There remain two points which must be clarified. The first is to show that there exists a unique solution of the finite difference problem (3. 1 ). that We shall defer this to a later section. 11·11 his a seminorm rather than a norm. The second is This latter problem is fully discussed by Keller [1971 ], but we will reproduce the - 2 6explanation here since possible oscillations in the finite difference solutions are of concern in practical computation. Two net functions fcpj} and t tj} jjcp- tJJh = 0 if and satisfy only if cp. = t. + ( -1 )Jp for some constant p. Thus if lien -~Jjh = 0 J J and jjfn-Tnjjh = 0, there eXist p and q such that~= e~ + (-l)j p J and f.n = f.n + (-l)j q. In order that fe.n}, tf.n}, J J J J satisfy the boundary conditions we must have ao P - 13o q J r;;·_n}, and ff.n} J J ::}. Four cases can occur. (3. 16) First, if a 0 13 1 + a 1 13 0 f: 0, then p = q = 0, and the seminorm is actually a norm for net functions satisfying (3. 3e) and (3. 3f). Second, if 13 0 13 1 = 0, then p = 0 so that fv.n} but not J = 0, then q = 0, f u.n} may have oscillations. Third, if a 0 a 1 and J [u~} may have oscillations. Finally, none of the above may happen J so that both fu~} and f v~ } could have oscillations. In the latter J J cases oscillations are eliminated by averaging neighboring values. Define -n u. 1 - -n v. 1 - J-z J-z for 2 ~ j ~ J. 1 n n z-{u. + u. 1 ), J J- .!. n 2 (v. Then J (3. I 7a) + v.n 1 ) , J Jjli"njjh = Jjunjjh and (3. I 7b) Jj-;; njjh = JJ vnjjh, but now JJ • JJ h is a norm for net functions defined on fx. 1: 2 ~ j ~ J}. J-z Any - 2 7- oscillations are now ren10ved, and {3. 15) therefore tells us the errors are not worse than O(h 2 ) if U and times continuously differentiable. V are piecewise four - 2 8I. 4 An Energy Inequality for a Linear Parabolic System We now wish to extend our analysis to a larger class of parabelie equations. Consider the following problem: A 2 (x) U (x, t) -x = Y._(x, t), (4. la) Yx(x, t) = !:\(x, t) - C(x) Q(x, t) (4.lb) - E(x) y_ (x, t) - ~(x, t) , = g_(x)' Q(x, 0) (4. lc) y_ (x, 0) (4. ld) dx a 0 Q(O, t) - f3o Y._(O, t) = _[o (t), (4. l e) 0'1 !:!(1, t) + !31 Y._(l, t) = [1 (t)' (4. 1 f) -1 Here -1 f3 0 and !31 2 -1 (4.lg) exist , -1 A (O)f3 0 a 0 and ~(1)!3 1 a 1 are positive semidefinite and symmetric, (4. lh) A(x) is symmetric and positive definite uniformly in x (4. li) Q, ~ ~(x, t), g_ 0 , and [ 1 are vectors of dimension p A(x), C(x), E(x), a 0 , f3 0 , a 1 , and !3 1 domain of the problem is 0 ~ x ~ are p Xp matrices. 1 and 0 ~ t ~ T. and The Condition (4. lh) has been imposed so that boundary integrals which will arise will have signs such that they may be dropped from inequalities in which - 2 9they would otherwise have to be retained. (4. lh) is a convenient assumption but not an essential one. We now seek an energy inequality which we can use as the As before we suppose £ and v basis for a convergence proof. are functions which satisfy the differential equations, the initial conditions, and the boundary conditions. Let U and y_ be another pair of func- tions which satisfy only the boundary conditions. ~(x, t) - :Q(x, t) f_(x, t) - We find that e y_(x, t) and f = .f+_e, f = ~t -f t £(X, t) 1 (4. Za) ~(x, t) (4. Zb) satisfy A2 e -x - x We define (4. 3a) C e - E f + £_, (4. 3c) + ~- .' a 0 ~ (0, t) - [3 0 .£ (0, t) (4. 3b) = Q_, (4. 3d) (4.3e) (4. 3£) 0 (4. 3g) (4. 3h) - 3 0- p_, £_, and .G_ are error terms resulting from the fact that do not necessarily satisfy the differential equations. U and V If we take the dot products of ..[ with (4. 3c) and A 2 ~ with (4. 3b), add the results, integrate over 0 .,;; x .,;; 1, and integrate by parts, we arrive at (4. 4) 1 - (A2 Lx' ~t) + [A2e • f Jo - (AxA ~ ..[) - (AAx e, f ) - (A 2 e , f ) • - x We introduce the notation \ A\ \ for the maximum for x E[ 0, 1] of the Euclidean norm of the matrix A(x). Let E: be the positive arbitrary parameter in the generalized arithmetic-geometric inequality. constants K Define K, C 1 and C 2 as follows: 2 max { 1 + 1- \\ ACA - - \\ + 1- \\ Ax \\ 2 2 \\ AE \\ 2 1 + \\E \\ + 1 + 2 \\ Ax \\ 2 2 + ~ \I AA) \ (4. Sa) 2 + (3;e) \\AxA \\2}, 2 cl C2 ~ - ~ \\A -1 \\ - (4.5b) ' 2 - 2 + 4\ A\\ . (4. 5 c) - 3 1We next make substitutions using (4. 3a) and (4. 3b), simplify (4. 4) using the boundary conditions, the Schwarz inequality, and the generalized geometric inequality, multiply both sides of (4. 4) by the integrating factor Ze II A~(t) ll -Kt , and integrate both sides from 2 2 + ll .£(t) ll 2 2 to t: t 2 + f eK(t-s) (Z I[i(s) ll 0 + Cl i! Al_x(s) ii + il£(0) 11 0 (4. 6) 2 )ds ~ eKt( ii A ~(O) il t ) + /eK(t-s) 0 { C 2 11£.(s) li 2 + Z II A~(s) li · II A £(S) II + Zll i_(s) ll ·< li ~(s) ll s)·~(x, s)] 0 } ds 1 2 + ll .eJs) li) + 2[A (x)_£(x, We have assumed € has been chosen so that C 1 is positive. I In further I simplifying the boundary terms it will turn out that C 3 , another constant, arises naturally: - A 2 -1 (0) !3 0 a0 + II A~(O) II ~(0, 0)·~(0, 0) 2 (4. 7) 2 + lli_(O) il . Inequality (4. 6) is also valid if we replace t by if 0 ~ T t, simplify the boundary ~ t. We integrate T from 0 to T on the left side and terms, and use the arithmetic -geometric inequality to obtain - 3 2 t 2 2 I ( II Ae(s} ll +II f(s) ll )ds s; 2 c 3 t eKt 0 - (4. 8) - t 2 + 4 t eKt I e -Ks c2 II £_(s> ll ds 0 + 2 t 2 2 t 2 e Kt I e- Ks CIIAo(s}ll 0 - + 2 II P (s) II 2 2 + 2 11 C (s) - II 2 Jds . Inequality (4. 6) may also be reduced using the geometric inequality to (4. 9) t + eKt / e-Ks {c211£( s) ll2 + e-Ks ii A £ (s) \12 0 + 2e-Ks llf._(s )\12 + 2 e-Ks ll _e_(s) \12} ds + e Kt { / ( II 0 A~ (s ) II + \1 !_( s ) II G)d s} . Inequality (4. 8) is now substituted into the braces in (4. 9): (4. 10) t +2eKt(l + 2t 2 e 2Kt) I e -2Ks ll<: (s) I!Gd s 0 t + 2eKt(l + 2e e2Kt) I e -2Ks ll_e(s) 11 2 ds 0 - 3 3 (4. 1 0) is the desired form of energy inequality. and f can be bounded for a fixed time t It tells us that in terms of 0 _Q, ~ e and c 3 which depends on the initial conditions. In c onclusion we give a brief summary of the technique of energy inequalities as used in our work. One first derives a differ- ential inequality for suitable variables such as ~ ~ ~(t) ll 2 + ll _!.(t) ll 2 • The differential inequality is solved in the manner of Gronwall's inequality. This process c an be used both for c ontinuous and dis c rete equations, but since the dis crete case tends to be more c omplicated, we first derive the c ontinuous inequality to use as a model. - 3 4I. 5 Convergence of the Box Scheme for a Linear Parabolic System In this section we wish to analyze the convergence of the box scheme applied to problem (4. 1) using a dis crete analogue of the energy inequality derived in Section 4. functions approximating Q and Let {u~} and J Y. the solutions of (4. 1 ). We make the same assumptions on the matrices as before -namely (4. lg), (4. lh), and (4. li). The finite difference equations are (5. la) 1 = D- v~- 2 X - J (5. 1 b) 1 1 n-2 - E. 1 v 1 J -2 - j -2 ul . - J l v. - J = g_(x.)' J = A 2 (x.) J n-2 s- J-2 . 1 ' (5. 1 c) d g_ (x . ) n a o ~~ (5.ld) J dx - f3 o ::0n = .[o (tn) ' (5. 1 e) g (t ) . (5. lf) :::J. n We define the errors U(x ., t - J V (x . , t - J n n (5. Za) ) n -J. ) - v. (5. Zb) - 3 51.r fn. } The errors sa t"1s f y -J f~ I - J -2 + 0~ I (5. 3a) J::..J -2 I D - f n-2 . X -J = I - E. (5. 3b) I n-2 n-2 + 0. 1 J-2 J-2 f . 1 J-2 I 1 n n-2 = Dt- -f J. -21 + £. J. -21 (5. 3 c) = Q_, (5. 3d) (5.3e) where {.e.~]-2.!.} • { )a~-}}. -2 and C ~-t} are truncation errors. { -J-2 in the case of the heat equation, they will all be 0 (h 2 ). As In place of the function exp(-Kt) we will use its discrete analogue: gl g where k n m - - 1 • g (t ) n ( 2-K k 2 +K k : =fr m=2 ) for n ~ 2 • } (5. 4) < 2/K for all m. 1 We begin by taking the dot products of l.J· .!. with (5. 3c) and -2 1 A~ n-2 n-- .!. ~-J -2~ with (5. 3b), adding the products, multiplying by h.,J and J- 2 summing from j = 2 to j = J. We would like to sum by parts 1 D- f n- 2 ) xh' but now we are faced 1 with the problem that, for instance, A2 f n -2 means evaluating A 2 - 36I at x . .!. and multiplying by the averaged value of ~-zJ-2. I -t n-n-) f. 2 and A 2 (x.) f. 2 • The summation rather than averaging A 2 (x. J- 1 -J- 1 J -J by parts formula (3. 7a) requires the latter quantity. In order to n-.!. I proceed we shall reinterpret A 2 1. 2 and A 2 en-z- to fit the formula, I but we shall then have to accept new terms I wn-z- and y_n-z- which account for the difference between the terms we actually have and the terms we need for summation by parts. An analysis of Y':! n- -21 n--I and y_ 2 is not needed now so that it will be postponed to a later I section. I At the present time we simply accept :y;:n-z- and .:yn-z- as net functions which make the following modified summation by parts formulas true: (5. 5a) (5. 5b) With (5. 5), the analogue of (4. 4) becomes - 37 (5. 6) I n-z-) 2 1 h + [A2 fn-z-. D- n J t ~ - J 1 As a result of the modified summation by parts formulas, the A2 in the next to last term on the right side of (5. 6) means the average of A 2 (x.) A 2 (x. and J J- 1 2 ) rather than A (x . I 1 ). J -z- Thus we cannot use A 2 D- en-z- in this term. (5. 3a) to substitute for x- As a matter of fact, we really should have a different notation for this A2 since it has a different meaning in this term than in the other terms of (5. 6 ). Define A.2 I - 2 I z-(A (x.) + A 2 (x. 1) ), (5. ?a) X. - ~ J (5. ?b) J-z I J-z J 2 ,._,n We define _£. I J-z A~ I J-z J- D X by means of n e. -J -n . I = -~J-zI + £. J-z (5. 8) - 3 8We claim now that ---n .e_. 1 J-2 n = -p.J-2 1 + O(h 2 ), (5.9) The validity of (5. 9) depends although the proof will be postponed. 2 both on the smoothness of A (x) and the (yet to be defined) ''smoothness'' of the ~n- 1 2 and [e~}. -t We add at this time that the terms c ontaining :tn-2; will also require "smoothness" on the part of the finite difference solution. Define two constants K K These points will be fully discussed later. and by C1 = max { [2 + \I ACA - 11 + II ACA- 11 1 1 2 2 (5. lOa) + llcA-1 11 2 + II (D~A 2 )A- II J. [ II A- D~ A2 ll 1 2 1 + II AE II + IIX Ell + 3 11 D~ A2 ll +l+IIEIFJ} . 2 2 2 cl = l+ II XII • (5. lOb) Substitution of (5. 3b) and (5. 8) into (5. 6) followed by applications of the Schwarz, arithmetic -geometric, and triangle inequalities and 1 multiplication by 2 gn- 2 yield (5. ll) -39Let 2 -1 A (O)f3 o 1 1 1 2 -1 ao!:1·~1 · + A (1)[31 1 , al~J·~J (5. 12) + I ~ ~ II ~ + 1 1.£ II ~ . 1 Then a time summation of (5. 11) and an analysis of the boundary terms lead to 1 n (5. 13) C2 g Let k 1 - k2 and define (5.14) On the left side of (5. 13) we change to i = n. n to i then sum from i =2 Also we use the geometric inequality and the triangle inequality in order to get norms of time averages: -40- (5. 15) C2 (l +C3) t n g +C n 3 1 + (ym--z1 D- fm--z) - • x- h J + (5.15) corresponds to (4. 8). To find the discrete version of (4. 9) we return to (5. 13) where we estimate the terms involving truncation errors with the arithmetic-geometric inequality: :;;; _1_ n g c2 (5. 16) -41- The terms in the last set of braces are estimated with (5. 15 ). Suppose we define (5.17a) = max 2 ::; m ::; n 0 (5. 1 7b) -f , m max~ :n (I (5.17c) Then the final inequality may be written as (5. 18) 2 -1 + A (I) f31 1 a1 1 l ~J • ~J J -42(5. 18) tells us that the error at time t n can be bounded by an ex- pression involving three types of quantities. The first type consists {~j} and {Lj} made in the initial values. of the errors type consists of {wr:-}} and -J-z The second {v~-f} combined with {e~ IJ_ and {f~L ..cJ-z -J -JJ which arose when we summed by parts two terms for which the summation by parts identity was not really valid. truncation errors { The third type involves the Q_T-i-~J, {~i} and {£.~t} . It is clear that the first type of term can be made 0 (h 2 ) by a sufficiently accurate approximation of the initial condition. if the { £:j~t} do not differ from the has previously been claimed. The third type will be O(h 2 ) m I 2 { .e_j-tf bymorethanO(h )as We must investigate further the net {£.j~t } , {Y:jm-·n -t J , and {y_ m-tl j -t J . 1 functions We start with { £.~ .!.} J-z pansions and (5. 3a), = where (5. 7a), It can be shown by Taylor series ex- (5. 8), and (5. 9) that h.a m _.L_ 01 + 8 D J:;__J -z X n e. ' (5. 1 9) J S· 1 is some point between x. 1 and x.. J-z JJ If the second derivative of A 2 (x) is uniformly bounded and if D- e~ is 1 the {£:.m.!.} J-z will be O(h 2 ). As for the {w:ni 2 } -J-z treatments are similar so we shall consider using Taylor series expansions, we find X-J and 0(1 ), then .!. , their {y_.m-t} J-z m-z} alone. Again {w. -J-2 1 1 -431 ..n_ -z- A 2 (x. .!.1. 1 ) • - J-z: 1 n-z- + i]· - 1 2 = - (5. 20) 2 h. + _l. • 4 8 A 2 (x . I) J -z- ox 1 n-z: w. 1 -J -2 is defined as all the terms on the right side of (5. 20) except .!. the first one. The same definition holds for 1 placed by r:-} y_:J-2 ~- 2 • J If D-e~ and X "J will be 0 (h 2 ). y_~- i n-1. 2 with f. J -z -J 1 D- f~ are 0(1), then X- J -J-z re- wr:-1 and ~ii} and Furthermore, the :n1:i} {y J-2 - n occur in inner products and time summations with D Dt e . and X -J - n D f . . These latter must be 0(1) to guarantee that our terms of x- J the second type be 0 (h 2 ). If U and in x and the cross derivative of conditions are that D for { ~j} and { ~j} - n- V {w -J-2 are continuously differentiable U is continuous, then our three n - -n u., D v . , and D Dt u . be 0 (1) as X -J X -J X -J solutions of the difference equations. h ...... 0 We shall show in the next section that these conditions can be satisfied. -44I. 6 Solvin~ the Linear Finite Difference Equations Problem (5. 1) may be written as a block t ridiagonal linear algebraic system of equations with Zp X Zp blocks. We would like to show that factorization into upper and lower block triangular matrices can be used to solve this have a proof only for the case p syst~m; however, = 1. we currently In this section we restrict ours elves to an equation with constant coefficients and a net with uniform space steps h and uniform time steps k. In the next section on nonlinear parabolic equations we will generalize the proof to include variable coefficients and nonuniform meshes. After we show that the finite difference solution exists and is unique, we can show that it is "smooth" in the sense described at the end of Section 5. In this section we consider the equation = au XX + b u X + c u, where (6. 1) a, b, and c are constants with a > 0. One can show that the matrix of the system of linear difference equations is 0 .. B3 A= 0 c3 ·•···...·•....··•··•···...·.•·...··c J -I •• •• ••A. ••••. B J J (6. Za) -45where the A., B., and C. are given by J J J ao - ( ~- c~ (6. 2b) 0 = (!;_ ch k 2 (6. 2c) A. J (6.2d) B. (6.2e) c. (6. 2£) J J (6. 2g) (6. 2h) -46We wish to show that A can be facto red in upper and lower block triangular matrices with 2 X 2 blocks in the following for j = 2 to J - 1. manner: (6. 3a) I !J!= L~·········... ••• 0 0 ••• •• ••• •••• •••• •• •• ••••• LJ ul cl (6. 3 b) •• I ..... 0 ••••• U2 •• ••••• ••••••• OM = 0 (6. 3 c) ••• •••• ••• ••• ••C2 •• •••• UJ If some right hand side vector f is given, then we could solve Ax = f !Z!Y:! = f for w and then solving OM?!::_= '!:!· fZ! is clearly nonsingular so that Y:! can be found by working recursively down through !Z!. The back substitution to find ?:::_ will require each by first solving of the U . to be nonsingular. As a matter of fact, 6l/ cannot be conJ structed unless ea ch U . is non singular for j = 1 to J - 1. What we J shall do is ve rify the block tridiagonal facto riz;ation and then ch e ck -47U J for invertibility. Once this is done we will know that Ax = f has a unique solution for each _£or that the box scheme advances the· solution uniquely for each time step. If we multiply!/! and uti. we see that the following relations · must hold: U1 = L. = U. = J J B1 , (6. 4a) -1 ~ A.J U J. (6.4b) . . . • J. J = 2, B.-L . C. I J (6.4c) J J- Define (6. Sa) It can then be shown by induction, which we omit, that (e~ )(~ - ~h) h 2a + (-I - b; )(~. ) 1 1 = 1 (6. 5 b) U. h k ch l 2 bh ---y det U. 1 = (6. 5 c) These recursions were suggested by Varah [1972]. to show that all of the the first of the e. : 1 e. 1 are negative. The first step is We start with e 2 which is -48(~o_ = e2 · Po !:_) - h (~0 ·zb + f) + k ao 1 + Po e2 Po (6. 6) h • 2a will be negative under our hypotheses if h Let us suppose that e2 ~ -M where M is sufficiently small. is a positive number. We can show under a mild restriction that all of the e.1 will. be less than or equal to -M. From (6. 5c) we see that this would imply all of the U. J at least through J -1 would be nonsingular. Eliminating U. between 1 (6. 5b) and (6. 5c) results in (6. 7) ~~) • ei} • { 1 - ~ • ei + ~ • ( b h + ak We assume to ei ~ -M. -1 - ~~ ) } . If we ask that ei+l also be less than or equal then we are imposing a condition on the right side of (6. 7). -M, We find that this condition implies k ~ 2 (6. 8) M2 + Mb + c a A similar examination of det U J shows that it is negative if we take into account the boundary conditions and if h is sufficiently small. Hence for a sufficiently fine net Ais nonsingular. -49Now that the { ut} and { vt} are known to exist uniquely, we return to the question of whether D - u.n , J X all 0 (1 ). X and n. X are J An examination of the difference equations shows that it would be sufficient to show D j n Dt D u. and n v. J - n and Dt v. J were 0(1) for all We will do this in an inductive manner using an argument similar to one given by Strang [1960]. The essence of the argument is to interpolate the finite difference solution { ut} and { vt} at time t n with functions U and V. sufficiently smooth at time t n We insist that U and V be and that the coefficients and boundary conditions of the differential equation be sufficiently smooth so that U and V will have five continuous derivatives at time tn+l . di££erencebetween at time tn+l {~+l} and U attime tn+l 1 and {vf+ } and V will then be equal to the first principal error terms which are O(h 2 ) plus some residual terms. The point of the argu- ment is to show that the residual terms are at worst 0 (h 2 ). since The U and V are smooth, Then {u~+l} and { vj+l} will be "smooth" also. We begin by introducing additional notation. Let (~:) be a vector consisting of the {ut} and {vt}. Let ( ~:) be a vector consisting of continuous functions U (x, t) and V (x, t) evaluated at x. and t n . J error terms for U and (1, n)) ~I, n) ( are the first principal V in the Richardson extrapolation of the -50r finite difference solutions at time t n [Keller, 1971 ]. An+z- is the matrix multiplying the vector of unknowns I Bn+z- is the matrix multiplying the vector of knowns I f n -t 2 is a vector of inhomogeneous and boundary data. We define w n by (1, = ( ;(1, n)) (6. 9) n) (6. 9) says that the difference between the net functions Jl u.n l and J ) ·tr vjn lJ and the functions U(x, t) and -V(x, t) evaluated at the net points is equal to the first principal error terms which involve U plus some residual vector ~n· V and The system representing all of the finite difference equations and boundary conditions in the box scheme n to tn+l = B n+l.z for advancing from time t An + ! ( n+l) ~ v n+l can be written as (~:) 1 + fn +z- (6.10) For the single parabolic equation considered earlier in this section 1 An+z- would be function of x { vt}. A for all n. At time which interpolates t n we construct a smooth { ut} and has derivatives matching Let this fun c tion be an initial condition for an initial boundary value problem starting at time t n and having the same boundary conditions as the continuous problem we are dis cretizing. V Let U and be the solutions of this problem, and let them have five continuous derivatives. This will in general require the initial condition, bound- ary data, and inhomogeneous terms in the differential equation to -51- satisfy some differentiability conditions. not the same as U and V. Notice that V are U and The latter are solutions of an initial boundary value problem starting at time zero and which we are trying to approximate by { ur} and { vr} while the former are solutions of a problem starting at time tn with initial data based on { ur} { vr}. The principal error terms are but here we are substituting U and and g~ne;ally functions of U and V, V. Combining (6. 9) and (6. 10), we c an show that w n+l (6. 11) (1, n+l )) ( ;l,n+l) + By definition of the principal error terms, the quantities in the braces must add up to a result which is O(h 4 ). Our choice of U and n guarantees that ':!:!_ n+l = Q. Therefore w multiplying a vector whose terms are O(h of ., A n+.l2 is at worst O(h - 2 ) so that ':!:!_ I pri._ncipal error terms are also O(h 2 ). In particular since O(h 2 ) at the point n+l is equal to (A 4 ). V n+.l2 -1 ) The norm of the inverse = O(h 2 ). Since the first 0 (h 2 ), the left side of (6. 9) must be V is smooth and v.n+l differs by only J must be 0 (1 ); hence, the desired smoothness c onditions on the finite difference solutions can he s atisfied, and the c onvergen c e proof is es .s ent:i a lly -52For the linear parabolic system where p > 1 and for which complete. we do not have a proof of nonsingularity based on block factorization, we will include nonsingularity as an assumption. The proof of smooth- ness will then be formally identical to the one we have just given for p=l. A further remark is that while U and V may be taken to have an arbitrary number of derivatives, the magnitude of the derivatives need not be 0(1 ). In particular if at time zero there is a sharp change in the initial data over an interval of length h, then our present analysis is not adequate to show that the finite difference solutions will be smooth. U and V On the other hand if the derivatives of the solutions are small compared to the inverse of the mesh spacings, then the preceding argument when applied at each time step for 0 (h time steps shows that no oscillations greater than O(h) c an form. recall from an earlier discussion that for certain boundary conditions O(h) unless -1 II • ll h is a s eminorm and that {uf} or . U and V {vt} might have oscilla- have derivatives which are large c ompared Finally, it should be noted that such small amplitude os cilla- tions are allowed under our definition of smoothness for net functions ; that is, smoothness and freedom from oscillations are not equivalent. We summarize the results of Sections 4, 5, and 6 in the following theorem. ) We It is now clear that these oscillations will not be worse than tions. to h -1 -53- Theorem 1: Assume (1) that the box scheme formulation (5. 1) of the linear parabolic system (4. 1) ~~unique solution and (2) that the coefficients ~ (4. 1) ~sufficiently smooth~ that ~initial boundary value problem posed at any non-negative~ with piecewise five times continuously differentiable initial functions will have solutions which ~also five times piecewise continuously differentiable. l£ points of discontinuity of the derivatives ~always taken~ be mesh points, then the box scheme solution converges to the continuous solution of 2 (4. 1) as the mesh is refined, and the errors ~ O(h ). -54I. 7 Nonlinear Parabolic Equations In this section we wish to study the problem v = a(x, U)U vX = U t - S (x, t, U, V) , (7. 1a) , X (7. 1 b) U(x, 0) = g(x), (7. 1 c) V(x, 0) = a(x, g(x)) dg(x) dx (7. 1d) = g 0 (t), ao p(O, t) - 13o V(O, t) a 1 U(1, t) + 13 1 V (1, t) = g 1 (t) , (7. 1 e) (7. If) ao ~ ~ 0 , (7. 1 g) ~ 0 , (7. 1h) al K 0 <a* ~ a (x, U) ~ a * <oo, (7. li) lau (x, U) l ~ a** < oo , (7. 1j) lS (x, t, U, V) l ~ s* < oo , (7. 1 k) * <oo, (7. 1 £,) s* , (7. 1m) ~ s ~ where in (7. li) through (7. 1m} the inequalities hold uniformly in x, t, U, and V. The box scheme applied to this problem yields the following finite difference equations: -55- (7. 2a) ' D X r n--z V. J I 1 u. = g (x.J ) , v.1 = a (x.,J g (x.J ) ) J J I - n n--z n--z = Dtu.J --zr-S(x.J --zr,t n--zr,u.J --zr,v.J --z1 ) , (7.2b) (7.2c) dg (x.) (7. 2d) dx n ao u1 - f3o v1n = go(tn)' n Q"l UJ + f31 VJ n (7. 2e) = gl (tn ) . The domain of this problem is (7. 2£) 0 .:;; x .:;; 1 and 0 .:;; t .:;; T. The net is the same as that used earlier. We would first like to discover under what conditions (7. 2) will have a unique solution. to construct the solution. Furthermore we would like to know how Let us then consider the matter of advancing the finite difference solution from time t n- 1 to time t . n have a nonlinear system of equations in the form r(y) is a vector consisting of the unknown { uj} and n nt , uJ, v J ) . Basically we = Q where y_ { vj} arranged in the The equations are ordered in the following way: (7. 3a) -56I -2 h.~­ v.n-z J ~X J - n Dt u. I J-2 (7.3b) F 2j -1 (~) - n-f - U.n-z - V.n-~I'" G 1 -2 h. a(x. I , u. 1 )D J J -z J -z x J -z J J ' (7. 3c) (7. 3d) where J ranges from 2 to J. We wish to solve this system itera- tively using the chord method so we next calculate ofF. $, the Jacobian The order of the unknowns and the equations was chosen so that $would be block tridiagonal: 0 $ = .··... ..•• •••• ...... •••••••• ••....•• ..··. •••••• ••··c J -1 0 •••• •••• •••• AJ .. (7. 4a) •• BJ Elements of the blocks are labeled in the following manner: A. J (7. 4b) -57B~ J ( )· B: J = B. J (7. 4c) B:J B~ J c. (7. 4d) J If we use the notation n-z-I a. n-z-1 a(x . I,U. 1 ) J -z- J -z- - I J-z - - au J-z J-z a S (x. _!_ t I - av J-z J-z (7. 4£) J-z a S (x. 1, t I sn-?. v. 1 u. J -z-. au J-z sn-?. u. I I n-z1 ) a a(x. 1 I a n-zu. I (7.4e) 1 n-z-I n-z-1 n-z-1, u.J-zI, v.J-z1 ) 1 n-z- _!_, u. n- z I , J-z n-z-11 ) v. I (7. 4g) I (7. 4h) J-z the matrix elements take the form Bll - Cl'o B/ - - ~0 B3 - Ql - ~~ - - __.]_ a n-z- D - u.n-z- + a.n - f- , J B4 J l A. J (7. 4i) I (7. 4j) I (7. 4k) I (7. 4£,) h. 1 2 u. 1 J-2 1 X J 1 J-z (7. 4m) -58A~ J h. _J_ (7. 4n) 2 h. _j_ D 2 u~ J X n-z-I a. (7. 4o) 1 J-2 h. B~ J + B~ J -~ sn-z- B~ J (7. 4p) I h.+l 2 h . +l + ..J..!..:._ k u.+l n J -z.!. hj+l 2 S n-zv.+ .!. J I -z (7. 4q) 1 c~ J B~ c~ J -~2 sn-z- + I (7. 4r) (7. 4s) J h.+l 1 - 1 (7. 4t) v.J + l -z .!. Since the matrix$ comes from linearizing a system of equations, it has the same general form as the matrix for the box scheme solution of a linear parabolic equation with variable c oefficients. We have already examined the special case of a linear equation with constant coefficients and uniform net spacing, and we shall use the previous study as a model for the current nonlinear problem. wish to show that $ can be factored asPCfUwhere As before we P and CfU are lower and upper block triangular (in fact bidiagonal) matrices. shall use the same notation as in (6. 3) and (6. 4). As before w e define a sequence { ej :j = 2, • • • , J} by We -59det B 1 (7. Sa) - e2 !3o A~+ ll'o B: J U. J = + 1 e. J A(# 1 B~ + e. J c~J- 1 J (7. Sb) B~ J B~ J det U. (7.5c) We wish to show that all the J will be negative. This is equivalent is nonsingular for j = 1, • • • , J -1. We start with J In the case where !3 0 = 0, e 2 will be negative for h 2 sufficiently to showing e2• e. U. If !3 0 f. 0, it will in general also be required that small. n-t n-t n-t - z1 a n-t uz 1. ( u 2 - u 1 ) + a 2 -21. > 0 (7. 6) -z as h 2 .... 0 so that e2 < 0 for sufficiently small h 2 . more about this requirement later. e2 ~ - M where M > O. Let us assume now that ~ - M for j = 2, • • ·, J. J and see what condition is neces- We wish to show that e. Let us then assume it is true for sary to insure it for We shall say j + I. on the right side of (7. Sc) by j Multiply the numerator and denominator 1 4 (Aj+l Cj-l /ej) -l . We introduce new notation for several important groupings of terms: -60B .1 B.4 (A1"+1 C.4 1 ) -1 J J J J- B.:a J 3 (Al Bj 1 c4. j+l 2 J -1 1 :a J J > -1 (7. 7a) h-.1 J (7. 7b) A.+l B. - B. A.+l J J 1 (7. 7 c) 4 Aj+l Gj -I hj k C .3 1 B4. - B3. C .4 1 ]J 1 J- n h."" • l J (7. 7d) 4 Aj+I cj-I Aa"+ 1 c.3 1 J J1 (7.7e) 4 A.+ 1 c. 1 h. J JJ n--aI If uj hj n--a1 - uj-l and hj+l = O(h) as h ... 0, then Dj+l ~ 0 and sufficiently small. can estimate ej+l Ej+l ~ 0 for If we assume this is the case, we by (7. 8) 3 1 4 + C. I (A.+ I C. I ) JJ J- -1 4 Under our current assumptions Gj+l small. 3 l B. - B. (A. + I ) J J J ~0 If we also assume Hj+l < 0 for for hj -1 } • and hj+l sufficiently hj, hj+l' and kn sufficiently small and if we ask for the right side of (7. 8) to be less than or equal to - M, we arrive at the condition -611 k n Ej+l (7. 9) and Ij+l kn -+0 if and only if hj+l/kn is remain bounded as bounded which we have already assumed. Hj+l 0(1) as is is O(h.) as as kn .... 0. J k n kn .... 0. Gj+l is independent of kn We are thus left with showing that Dj+l - 1 .... 0 in order for the right side of (7. 9) to be bounded If we write out Dj+l - 1, . n-21 n-21 1f uj+l - uj = O(h). we will discover that it is We find therefore that if the net function is smooth, then for sufficiently fine net spacing all of the e. 1 less than or equal to -M. n-n- The requirement that uj+l - uj is a result of the discretization of U(x, t). with the a u J1 2 O(hj) { n} uj will be = O(h) It arises only in conjunction terms and corresponds to a difference approximation of the derivative. This discrete condition is analogous to asking to have a c ontinuous x derivative. factorization o f $ is possible. U(x, t) At any rate we now know that the It further turns out that det UJ < 0 without any additional assumptions so that $ i s nonsingular. We wish to find the solution of .[(y:_) = Q_ by the chord method. This is an iterative method of the form v+l Y.. = Y.. v - A -1 :r (y_ v ) • (7. 1 0) 0 We take y_ , the initial guess for the solution at time same as the s<?lution already known at time A is chosen to be the Jacobian of ~ tn-l" tn, to be the In the chord method evaluated at the initial guess -62o y__ • f n-11f If l uj is smooth, then ciently fine net. where A will be nonsingular for a suffi- Let 'L - A -1 !: (y__) ' y__ is any vector of length (7. ll) If we c an show that 2J . tracting map in some neighborhood containing 0 y__ and that g_ is a cong_ maps the neighborhood into its elf, then we will know there exists a unique solution of the nonlinear difference equations in that neighborhood. Let s and t be vectors. (7. 12) We apply the mean value theorem with vector r lying between s and t : -1 f:A _ g ( s ) - g (t ) = A L: We define the matrix Let the vector o~E_) u A... J (~ - •.!) to be the matrix in the brackets in (7. 13 ). M E_ have the components to evaluate the matrix elements = 0' M 2J = 0, 'm Ml, m J _ h. _j_ M2. . J- 2 , 2 J- 3 2 (7. 13) M. [ru. } and J We proceed .: 1, J m = I, • • • , 2J , ~ n- -1 S [rv.}. J (7. l4a) n -1 (x . r,t 1,u . f,v.-f) u J- 2 n- 2 J- 2 J- 2 r)J, -S (x . 1,t r,ru . r,rv. u J-2 n- 2 J- 2 J- 2 (7. 14b) -63- h. ~ ~ = M2j-2, 2j-2 2 .!. n-z 1 n-- S (x. 1, t 1, u. f , v. 1 v J -2 n-2 J -2 J -2 (7. 14c) for j M2j-2, 2j-l = M2j -2, 2j = M2j -2, 2j -2 = 2, • • •, J. All other M h. ~ = M2j-l, 2j-3 (7. 14d) M2j -2, 2j -3 2 (7.14e) . . are zero. 2 J- 2 > 1 r. n-2I n- l - (x. 1, u . 1 ) D u. ~..u J-2 J-2 x J a (7. 14f) -a (x. z, ru. 1 )D- ru . l U J-2 J-2 X ]...; + ra(x. l M2j-l, 2j-l for j M = 2, • • ·, J. I 1 J-2 U~-1) - a(x.J-2I ruJ-2 . 1 )] J-2 1 1 (7. 14 g ) = M2j-l, 2j-3 All other is block tridiagonal with 2 M . . 2 J- 1 > 1 are zero. X 2 blocks. Thus the matrix The terms in the even numbered rows can be made arbitrarily small by taking enough. h The odd numbered rows require in addition that be smooth. If s smooth, and if h and ..!. are restricted to be near _y_ 0 , is sufficiently small, we find that 0 in a neighborhood of Y...· The fact that M small £ and r then r will g_ is contracting is block tridiagonal guaran- tees that its maximum absolute row sum c an be made small even as h goes to zero. Our previous investigation of the linear algebraic system shows that A -1 exists as h goes to zero so that the norm -64of A -l (not necessarily the maximum absolute row sum) must be bounded in the limit. The size of the neighborhood in which con- traction occurs for any pair of vectors portional to (a tions on s ** ). -l and t of h. that d is and !_ ie essentially pro- However, we wish to avoid smoothness assumpso we shall restrict them to a smaller spherical neighborhood of radius r1 s 0 about y_ R where is some fixed multiple d = Il l - y_ ll . 0 We will specify r1 shortly. 0 (h) g0 and g1 are continuous. we can c hoose r1 so that provided ficiently small h an even smaller h Let R R > d. We note Thus for suf- We then take II A-l I\ • II M II < 1 - (d/R). so that we also have Then (7. 15) for s 1 v ll ,;._ - and !_ in the sphere. P'(t_)_l\ R- - - · J.S If s is taken to be }!_ l essth anorequa l toR-. d within a distance R 0 of y_ , Hence if t , we find is any vector g_(!J will also be within R of y_0 • In 0 other words there is a neighborhood of y_ its elf by g_ and in which 0 whi ch is mapped into g_ contracts. Therefore s_ has a unique fixed point in this neighborhood, and the nonlinear difference equations have a solution. We further note that as with the linear parabolic equation we may ask the linearized initial boundary value problem to have very smooth solutions given sufficiently smooth initial data when posed at times greater than or equal to zero. As before, we can then show that the iterates are smooth by studying a series of linear problems one for each iterate. Indeed this study is ne c essary for an inductive -65proof of the existence of a solution of the nonlinear difference equations for succeeding time steps, but it has the further implication that the matrix A might be re-evaluated at each iteration rather than at each time step as in the chord method. In fact the Jacobian of F can be shown to be nonsingular in a neighborhood of y_0 under certain conditions of smallness on h and smoothness of the "point" of evaluation. The previous methods c an again be used to show contraction of the series of maps g_v and of the mapping into themselves of successive neighbor- hoods. Of c ourse we may have to start with a smaller initial neighbor- hood or equivalently a finer net spacing to insure that successive Jacobians remain nonsingular, but otherwise the use of Newton's method is justified. Returning now to the question of convergence, we find that nearly all of the analysis we exhibited for linear systems can also be adapted to the present nonlinear problem with the use of the mean value theorem. As before u is the c ontinuous solution. is the finite difference solution and U The mean value theorem must be applied several times, but since it is not necessary to keep track of each application, we will use u as a generic symbol to indicate some function value intermediate to u and U. The three basic equations involving the truncation errors may then be written as I -n- - - 1 n-2 a (x. 1, u . r 2 )D e. J-2 J-2 X J (7. 16a) 1 - [ a -n-2 1 - -n-2 (x . 1, u. 1 )D u . U J -2 J -2 X J J f , 1 n- - e. J -2 -66I D- f.n-zX J I n-zI J--z = + 0. (7. 16 b) I - S u (X . I1 t J-z n--z 1 f J-z I J-z J --z 1 -n-z- - V. I n- - f J-z ) e. I - [ a (x. 1, u. 1 U J --z J -z [ Dt -n-z1 + -,rn_-zI = D t- l..n . 1 -n--z - - n a(x. 1, u . I )Dt D e . J -z J -z X J - I -n- - I1 U . 1 ) D (7. 16c) J -z- - -n-zu . X J JDt- e.Jn-z 1 1 -n - n-za(x . .!.1 u . .!. )] D e . , J-z J-z x J where n e. J - n U (x., t ) - u. • J n J (7. 17a) f.n J - n V (x., t ) - V. J n J (7. 17b) 1 We multiply (7. 16c) by 1 f.l-1 h., multiply (7. 16b) by a(x . .!.. u.nl_2 ). J-z J I e .n-z I h., add the results, and sum from J --z J J-z j = 2 to J: J- z -67- (7. 18) Equation (7. 18) has the same form as (5. 6) except for the variable coefficients depending on the solution. This means additional appli- cations of the mean value theorem will be necessary in order to carry out the analysis of Section 5. As before, the requirement of smooth- ness on the finite difference solution arises and is handled as in the linear case by an examination of the system of linear algebraic equations obtained from the Jacobian of the nonlinear system and by requiring U and V and solutions of the linearized initial boundary value problem to have x derivatives which are small compared to the inverse of the net spacing. in the following theorem. The results of this section are summarized -68Theorem 2: Assume (l) ~~nonlinear problem. G.-.!> has~ unique solution and (2) that the linearization ~problem (?_·_!) ':! ':lEY piecewise five times continuously differentiable .!:! and '!_ results in a differen- tial equation with coefficients sufficiently smooth so that an :!:._nitial boundary value problem posed~ any non-negative time with initial functions both piecewise five times continuously differentiable will have solutions which are also piecewise five times continuously differ entiable. ll_ the points of dis continuity of the derivatives ~ always taken t.Q_be mesh points, then for sufficiently small .h ~~scheme forum- lation (7.. ~) of the nonlinear parabolic problem (1.1) ~~unique solu2 tion which approximates the solution Qi (1.1) 12. Q(h ). -69- CHAPTER II THE TWO DIMENSIONAL HEAT EQUATION II. 1 The Method of Fractional Steps The numerical solution of multidimensional parabolic problems is of concern to scientists and engineers because there are many physical processes in which boundary conditions or properties of materials prevent realistic modeling by one-dimensional equations. In particular we are interested in two space dimensions. We could of course write various sets of difference equations coupling net points in both space directions, but this means we would have to solve a large algebraic system for all the net points in the domain at each time step. Instead we restrict ourselves to a rectangular domain and ask if the two-dimensional problem can be reduced to a series of one-dimensional problems. Just such a reduction is accomplished by the method of fractional steps. This method has not yet received a complete theo- retical justification for all the problems in which we would like to use it and is still undergoing active investigation. The reader is referred to Yanenko [ 1971] for an exposition of the techniques and applications -70of the method of fractional steps. For our present purposes we shall present a simple example using this method. In this case the method will be theoretically justified. Consider the equation u XX +U yy (1. 1) on the domain D = [ 0 ,1] X [ 0 ,1] to be zero on the boundary of D U(x,y,t) in the (x,y) for all time plane. t 2': 0. is equal to some given function g(x,y). We require U At time t = 0, As is well-known, this problem can be solved by separation of variables and Fourier series. The solution is then represented as a sum of terms of the form where m in time t and n 1 are positive integers. and t close together. 2 such that t 2 > t 1 > Now pick two different points 0. t 1 We pose two more problems. and t 2 need not be First ( 1 . 3 a) V(x,y,t) = 0 on an , (1. 3b) (1. 3c) -71This is a heat equation in one space dimension with y V is chosen to coincide with U at time t . 1 as a parameter. Now consider the second problem: (1. 4a) W(x,y,t) = 0 on an, (1.4b) ( 1 • 4c) This is another one-dimensional heat equation but with x parameter. t • 2 Notice that the initial value of W as a is the value of V at We now assert that (1. 5) This is an example of how the method of fractional steps reduces a two-dimensional problem to two one-dimensional problems. We can easily verify (1. 5) for this simple case when separation of variables is legitimate. A component of the form (1. 2) when used as an initial condition for the 2 2 -m w t e 2 where the factor constant. V problem evolves to - n 2 2 'IT' t 1 s in m wx s in n 'IT'Y , 2 2 exp (-n w t ) 1 ( 1 • 6) acts as though it were a multiplicative When (1. 6) is used as an initial condition for the W prob- -72- lem, the other factor, so that at time -m e 2 2 11' t 2 t2- n 2 2 exp ( -m w t ), acts as a multiplicative constant 2 we get 2 2 'If t2 sin mwx sin n"Y • ( 1 • 7) Thus we see that Fourier components of the two-dimensional problem evolve to the same extent as Fourier components of two successive one-dimensional problems. -73II. 2 The Box Scheme and the Method of Fractional Steps Having reduced the two- dimensional heat equation to two onedimensional heat equations with parameters, we easily see how to apply the box scheme. square D. We place a rectangular grid over the unit For each horizontal grid line we solve an equation in x, then we change directions to solve an equation in y grid line. for each vertical The procedure is then repeated for the next time step. The box scheme requires however that we give the derivatives of the initial data as well as the data themselves. x direction we deal with u sweep in th e dition. y and y X direction, we must give u x u as part of the initial con- y after an x using the computed values of u. interpolation as a way to do this. u u , but when we wish to perform a What we must do is ignore struct u For a sweep in the sweep and con- We have chosen Lagrange At any net point take the value of and combine it with the values of u at the next two nearest net points on the grid line to form the Lagrange quadratic polynomial which interpolates all three function values. of the quadratic. We then use the derivative Furthermore we do not change directions at every half step since it is not really necessary. in solving for u at t followed by a y sweep. 2 we start at t 1 Suppose for instance that and perform an x It will be necessary to fabricate when changing directions; however, in moving from perform the y sweep first followed by the x t sweep. 2 sweep y to t derivatives 3 , we Proceeding in this way, we need create derivatives only once per time step rather than twice. -74We have performed computations 1on the two-dimensional heat heat equation with the following initial condition: 3 5 g(x,y) = 2m +3n sin mwx sin nW\T ··; (2. 1) m,n= 1 The continuous solution is U(x,y,t) = 3 l ( 2. 2) m,n= 1 Using several different mesh spacings we have found that the error of the computed solution u is 2 O(h ) where h was the size of the steps in both space directions and in time. In the present problem there is no doubt about the consistency of the numerical scheme with the continuous problem since the fractional steps are theoretically correct and since the box scheme is cons is tent with one-dimensional heat equations to second order; that is, the truncation errors are be shown. 2 O(h ). Convergence however is yet to For this problem we choose to show that the numerical scheme is stable in the sense of von Neumann. This stability analysis is appropriate for difference schemes with constant coefficients and which correspond to pure initial value problems with periodic initial data [Isaacson&. Keller, 1966] • -75II. 3 Von Neumann Stability The basic idea of a von Neumann stability analysis is to decompose the finite difference solution into Fourier components and then to show that none of the components can grow in amplitude as time increases. In this section we will use complex harmonics and coefficients rather than sines and cosines with real coefficients. shall study the evolution of the general harmonic i = .{::1 and a eiaxeif3y where and f3 are fixed, arbitrary real numbers. that all steps in the x direction are of size direction are of size - n+1 u and - n +1 w approximate t n We assume hx, all steps in the hy, and all time steps are be the net functions at time We k. which approximate Let U un and y and vn U • x Let be the net functions after a half time step which Finally let u n+1 V and V • y functios at time tn +1 which approximate and U w and . iax if3y functions must each be some multlple of e e n+1 U • y be the net These net evaluated at the net points: u n = c n1 e ia x e if3y n = c 2e (3. 1 a) n iax if3y e , (3. 1 b) -n+1 u = -c n1 +1 e iax e il3y , (3.1c) -n+1 w = -c n3 +1 e iax e if3y (3. 1 d) n+1 =c1 n +1 iax if3y e e , (3. 1 e) n+1 =c3 n +1 iax if3y e e {3. 1£) v u w -76- The above net functions are involved in going from time with an x sweep followed by a time to tn+i tn+Z with a y y sweep. n 2 are given. sweep followed by an x that c n+Z 2 does not exceed c n+Z does not exceed 1 c~ in amplitude. sweep. Now suppose The difference scheme will determine What we must show is that tn+ 1 We in turn proceed from is, one complete cycle covers two time steps. c tn to c n 1 c n+Z 1 That c n 1 and n+Z and c2 in amplitude and This would then imply that no harmonic can grow in amplitude; hence, the numerical scheme is stable. We introduce two symbols which we shall use in simplifying notation: 1 () = l ahx , (3. Za) ({' = 2 f3hy • (3. 2b) 1 Since u y derivatives are constructed from the rnost recent values of and do not involve any information about it is clear that = -n+1 c 3 - n+1 c1 ~ U X -n+1 can depend only on c 1 or its approximation, In fact i sin f3hy • (3. 3) and the substitution of {3. 1a) through (3. 1d) into the finite difference equations for an x between the coefficients: sweep yield the following relations (3. 3) -77- (3. 4a) where a= 2 . 2 stn hx2 cj A= A A 2 4 2 ) ' (3. 4b) (3. 4c) 1 2 k cos () i s ln "ll = hx o cos i A 3 = hy k A k A4 A3 A1 = e + 1 cos 2 () ' (3. 4d) () , (3.4e) 2 sin l3hy cos () , (3. 4f) 1 =- "hy1iX sin l3hy sin () cos () • (3. 4g) Similarly the substitution of (3. 1 c) through (3. if) into the difference equations for a y sweep yields (3. Sa) where 2 2 1 2 b = - - sin cp + k cos cp , 2 hy (3. Sb) -78- (3. Sc) = k1 cos 2 cp B (3. Sd) i . 2 = hy s 1n cp cos cp • (3. Se) 2i . = hyT s 1n cp cos cp 2 (3. Sf) 2 B 4 = - :--z sin cp • hy (3. Sg) These two transformations must be composed and followed with two more similar transformations representing the progression from time tn+i to tn+ • 2 We shall not give the details. but we write the com- position of the four transformations symbolically as (3. 6) where G is a 2 X 2 matrix. One of the eigenvalues of G The other is ~( 1 k cos 2 8 - ~ 2 s1n . 28 )~ • ~(1k cos 2cp - "1;2 2 sin 2 ({') ~ . .2 1 2) 2 ( -:--2 stn 8 + k cos 8 hx n (:--z 2 sin 2 cp + 1< 1 cos 2 cp ) hy is zero. -79- ( ~ cos e - ~ ~ 2 tz e sin cos f) sin ahx )~ • 2 .2 1 2) ( hx2 Sln f) + K cos f) (~ cos qJ -~sin qJ cos qJ sin ~hy)~ 2 .2 1 2) 2 ( -:-z sm qJ + k cos qJ (3. 7) hy It will be found upon study that each of the four terms in braces must lie between -1 and +1; hence, the eigenvalues of G are less than or equal to one in magnitude. Since a are real and and ~ were arbitrary, we have shown that no harmonic can increase in amplitude, and the numerical scheme is stable • . -. -80- CHAPTER III BURGERS' EQUATION: COMPUTATIONAL EXAMPLES III. 1 Burgers 1 Equation In this chapter we shall describe the results of computing on Burgers 1 equation with three different initial conditions. Burgers 1 equation is au+ u au at ax where U = U(x,t) ( 1. 1) and v is a positive number. This equation is discussed by Cole [ 1951] who describes its applications and its general solution. His result is that if B(x, t) is any solution to the heat equation (1. 2a) then e X U(x, t) =- Zv T ( 1. 2b) -81- is a solution of (1.1). We shall use this result in our third computa- tional example. We are particularly interested in the situation when v compared to unity. is small In fact let us first examine the case when v = 0: au+ u au = 0 ax · ( 1 • 3) Tt The characteristic ordinary differential equations are ciS= u, dx (1. 4a) dt (1. 4b) d£ = 1 ' ( 1 • 4c) If we take the initial condition x = 11 at time t = 0 and if U(x,O) = u 0 (x), we find (1. 5a) t (1. 5b) =~ ' (1. 5c) ~ is a parameter along the characteristic curve in the (x,t) plane on -82- which U is a constant. An interesting feature of (1. 5) is that characteristics originating at two different points the x t =0 axis at time can intersect at TJ 1 TJ and x \x,t) where on 2 and t are given by TJ1 uo<TJz> - TlzUo(TJ1) x= {1. 6a) uo<TJz)- uo(TJ1) (TJz- TJ1 > (1. 6b) t = - 0 o<T'lz) - uo<TJ1) Since we are interested only in positive time, we restrict our attention u 0 (TJ) to initial functions This means T). u0 such that intersection times should be a decreasing function of T) t are positive. in some range of The first instant in time at which an intersection occurs is deter- mined by the greatest negative slope of allowing T) 1 and T] 2 u0 • We may see this by to approach each other in (1. 6b). At this point we say a shock is forming in the solution; that is, the solution will become discontinuous. The initial function u0 determines which characteristics will lead into the shock and what the magnitude of the discontinuity will be. We note that the value of U must always lie in the range between the minimum and maximum values of u 0 (T)) u0 so that if has a bounded range, the jump in the solution must also be bounded. The speed with which the shock propagates is determined by the values of U write (1. 3) as just ahead of and just behind the shock. We re- -83- au + 1 a [ uz] zax at =0 . ( 1 • 7) If we seek a steadily propagating solution of the form U= U(s) where s=x-vt, (1.7)becomes -v 8 ~ u + _!_ ~ u 2 = o z as which we integrate from s1 to ( 1 • 8a) Z:z: (1. 8b) (1. 8b) may be solved for v: (1. 9) If s 1 and sz are chosen on opposite s.ides of the shock, then (1. 9) tells us that the velocity of a steadily moving shock is equal to the average of the function values just ahead of and just behind the shock. If v > 0, we no longer have shocks in the sense of intersecting characteristics. Nevertheless, if v see shocks trying to form. is small, we should be able to The small amount of diffusion will prevent the complete formation of a shock. On the other hand if u0 is a step function, the diffusion will smooth the step into a continuous function. We shall demonstrate both of these situations in the computational examples~ -84- III. 2 Example 1: Smoothing of a Sharp Front The domain for all three examples will be example the initial function 0. 52::5 x S 1. u0 1 for 0 ~ x :S 0. 48 In the interval 0.48 ::S x::::: 0. 52 cubic polynomial with value 1 zero derivatives at both ends. and is U{1 ,t) = 0. We take v 0 < x:::::: 1. U 0 and In this 0 for is the unique at the left end, 0 at the right end, and The boundary conditions are = 3 X 10 -3 • U(O, t) =1 With the exception of the initial condition, this example is the same as that given by Swartz & Wendroff [ 1969]. Their initial condition was a step function. We use a cubic transition function instead because we need to specify x derivatives as initial data for the box scheme. We have performed the computation using both uniform and nonuniform net spacings in the x always taken to be uniform. direction. The time spacing was Newton's method was used for solving the nonlinear difference equations in all of the examples. All compu- tations were performed in double precision on an IBM 370/155. first four figures 1 show the solution for uniform meshes. The The curves were plotted at intervals of 0. 1 time units so that the last curve corresponds to time t = 0. 5. The oscillations which appear most promi- nently in Figure 1 are actually a part of the numerical solution since according to our boundary conditions we cannot have any oscillations of constant amplitude over the entire net which must be averaged out. 1 Tables and figures are at the ends of the sections in which they are discussed with tables (if any) preceding the figures. -85- In Figures 5 through 8 we show a nonuniform mesh which is successively refined. These curves are also given at time intervals of 0. 1. In Figure 8 we see most clearly the behavior of the solution. The "corners" of the initial function are quickly rounded off, and a shock-like profile is moving to the right with speed one half. For an economical computation, one should probably change the space net as the shock propagates by deleting net points behind the shock and interpolating additional points in the neighborhood of the shock. . :::::J 0 . I I I I 0.1 I I I I I 0.2 I I I I I I 0.3 I EACH VERTICAL BAR BEL~W REPRESENTS A NET P~INT . I I I I 0.4 I NUMBER ~F SPACE STEPS = 50 MAXIMUM STEP SIZE= 2.000E-2 MINIMUM STEP SIZE = 2.000E-2 TIME STEP SIZE= 2.000E-2 NU = 3.0E-3 0.0 o1 ~a /'. xLD. ....... 0 I I I X 0.5 I Figure l I 0.6 0.7 0.8 I I I I 0.9 I I I I 1.0 I I 00 0' . ::J . EACH VERTICAL BAR BELeW REPRESENTS A NET PeiNT . NUMBER eF SPACE STEPS = 100 MAXIMUM STEP SIZE = 1 .OOOE-2 MINIMUM STEP SIZE= 1.000E-2 TIME STEP SIZE= 1.000E-2 NU = 3.0E-3 0.0 0.1 0.2 0.3 0.4 Figure 2 X 0.5 0.6 0.7 0.8 0.9 1.0 O.lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll o1 '-/ 0 "' xLD . ........ 0 I I 00 -.J . . 1 0.1 0.2 0.3 0.4 Figure 3 X 0.5 0.6 0.7 0.8 0.9 1.0 111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 I 1111111111111111111111111111111111111111111111111111111 EACH VERTICAL BAR BELQW REPRESENTS A NET PQINT . NUMBER QF SPACE STEPS = 150 MAXIMUM STEP SIZE = 0.667E-2 MINIMUM STEP SIZE = 0.667E-2 TIME STEP SIZE = 0.667E-2 NU = 3.0E-3 0.0 0. o '-../ 0 ::) xlD -" . ....... 0 I I 00 00 ~ '-J 1 I I I 0.0 0.1 0.2 0.3 0.4 Figure 4 X 0.5 0.6 0.7 0.8 0.9 1.0 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 J - I J - - - - - - - I ~ \\\\\\. \ '\ \ \ \ , I - I r- r- t--EACH VERTICAL BAR BELOW ~~REPRESENTS A NET POINT. 0 x~I- /'"'"'.. I OF SPACE STEPS = 200 ~MAXIMUM STEP SIZE = O.SOOE-2 MINIMUM STEP SIZE = O.SOOE-2 1-TIME STEP SIZE = O.SOOE-2 ~NU = 3.0E-3 r- NUMBER ~L ...... I I 00 I -.() . I I [EACH VERTICAL BAR BELaW I- 0.0 I I 0.1 I I I 0.2 I I I 0.3 I I I I \_ 0.4 Figure 5 X 0.5 0.6 I ~ \ \ I l I 0.7 I\\\\ \ I \ I I 0.9 I 0.8 - - - - - - I - - - 1.0 I I I J I I I I I 1111111111111111 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I NUMBER ~F SPACE STEPS = 60 r- MAXIMUM STEP SIZE = 3 .500E-2 MINIMUM STEP SIZE = 1 .OOOE-2 r TIME STEP SIZE = 1 .OOOE-2 ~ I ~ ~ REPRESENTS A NET P~ INT. ........ 0 I I I -.D 0 ::J 0 ~ 1 I \\\\\\ I I I I 0.0 l 1---- I 0.4 Figure 6 X 0.5 0.6 0.7 0.8 0.3 0.1 0.2 L I 0.9 1.0 I I I I I I I I I I I I I I I I 1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 J - I J - - - - - - I I- EACH VERTICAL BAR BELeW REPRESENTS A NET PeiNT. 1 \\ \ \ \ - 1 - 1 I- I- "' x~l- '-/ 0 I NUMBER eF SPACE STEPS = 120 f-- MAX I MUM STEP SIZE = 1 . 750E-2 MINIMUM STEP SIZE = O.SOOE-2 rTIME STEP SIZE = O.SOOE-2 ~ NU = 3.0E-3 ~L I I I ...... -..o ::::J I \ ~\_\_\_\ I I I - 0.0 0 .1 0.2 0.3 0.4 Fi gure 7 X 0.5 0.6 0.7 0.8 0.9 1.0 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 ~ I J - - - - - - - \ \ \ I J- I I - I J- f-- EACH VERTICAL BAR BEL~W ~~REPRESENTS A NET P~INT. '-/ 0 I NUMBER ~F SPACE STEPS = 180 ~MAXIMUM STEP SIZE = 1 .167E-2 MINIMUM STEP SIZE = 0.333E-2 ITIME STEP SIZE = 0.333E-2 ~ NU = 3.0E-3 f-- X~J- " . ....... 0 I I I N -.D ::::J '\,J 0 EACH VERTICAL BAR BELeW REPRESENTS A NET PeiNT. NUMBER eF SPACE STEPS = 240 MAXIMUM STEP SIZE = 0.875E-2 MINIMUM STEP SIZE = 0.250E-2 TIME STEP SIZE = 0.250E-2 NU = 3.0E-3 0.0 0 .1 0.2 0.3 O.L! Figure 8 X 0.5 0.6 0.7 0.8 0.9 1.0 . IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIAIIIIIIHHHIIIII-IUIIIIIHHIIIHIIIHIIIHIHHIHHIIIIIIHIIIIIHIUIIIIIIIIIIIIIIIHIHHIHHIIIHIHIIIHHIIIIIIIIIIIIIIIIIIIIIIIII ~~ 0 X~ .!'"' . ....... 0 I I (j.l -..o -94III. 3 Example 2: Formation of a Shock u0 In this example the initial function is for O.O:S:x:!S0.1~ 1 )) ( x-0. 0.4 11' for 0.1.SxS0.5 (3. 1) for 0, 5!:: X!:: 1.0 The boundary conditions are take v = 3 X 10 -3 • U(O,t) = 1 and U(1 ,t) = 0. We again If we examine the characteristics when v = 0, we learn that a shock will start to form at will be fully developed at (x,t) (x, t) = (0. 42 73,0. 2546) and = (0.5,0.4). For this example we have used uniform nets consisting of 100, 200, 300, and 400 space steps. equal to the space step h. The time step k was always Figures 9 through 12 show plots of the solution for the four meshes. The curves were plotted at time inter- vals of 0.1 so that the rightmost curve corresponds to time Visually it appears that by time t t = 0. 8. = 0. 4 the solution has reached its final shape and is traveling to the right at speed one half without further change. This leads us to ask whether or not it is possible to detect numerically when a shock has formed; that is, can we tell when the solution has reached a steady profile. case We decided to recompute the h = k = 0. 01, but now at each time we performed an inverse Neville interpolation on the solution for five different values of U: -950.10, 0.25, 0.50, 0.75, and 0.90. values of x In other words we ask for what does the solution take the five prescribed values. We then took time differences over one time step to approximate the speed at which each of the five given values of point out that when U moves to the right. We v = 0, a point of the solution will propagate to the right at a speed equal to its amplitude. For example, the point U = 0. 6 will appear to move at speed 0. 6 to the right until it is absorbed into the shock after which it moves at the speed of the shock. When v =f. 0, we would hope to see the speed decrease gradually (or increase gradually when 1 U < Z) to speed one half, and we would want it to reach this ultimate velocity at or before time t = 0. 4 shock for v =0 would be fully formed. of these computations. when the In Table 1 we give the results In the first column are the time values midway between the time mesh points. In the succeeding colwnns are the velocities found by inverse interpolation and time differencing. numbers are presented graphically in Figures 13 and 14. These We have plotted also a line segment from the first velocity point to the value on the velocity axis to which it corresponds. U We add that the values = 0. 25 and U = 0. 75 should be absorbed into the shock at (x, t) = (0.4333,0. 2667). Since it appeared that the velocities were oscillating, we decided to average them over two time steps or equivalently to difference the locations over two time steps. results. Table 2 gives these The time in column one is the midpoint of two time steps. We see that on the average these five points of the solution do approach speed one half though it is difficult to say when. One further compu- -96tation was done with v = 3 X 10 -2 and three prescribed values of U. These results are presented in Table 3 and Figure 17. We see that with greater diffusivity the velocities no longer oscillate. We are also interested in achieving more accurate solutions through Richardson extrapolation. Suppose u 1 solution and is related to the continuous solution u (x,t) 1 is a finite difference U by = U(x,t) + h 2f 2 (x,t) + h 4 f 4 (x,t) + h 6f 6 (x,t) + O(h 8 ) • u 1 of course is defined only for errors f , f , and f 4 6 2 (3. 2) (x,t) in the mesh. are defined for all x and t The principal and are independent of h, but in (3.2) they are used only at mesh points. k = h. If u 2 We are assuming denotes the finite difference solution for the same prob- lem but with the mesh refined by cutting all intervals in half, then u 2 must satisfy (3. 3) We can now take a combination of u of u 1 1 and in such a way as to eliminate the 4 an O(h ) approximation to U u at each of the net points 2 2 O(h ) term leaving us with at those mesh points. The combination -97is If we divide the basic mesh into thirds and quarters, {4u - u )/3. 1 2 we can compute u (3. 2) but with h replaced by h/3 3 and u 4 which will have expansions similar to or h/4. With four finite difference 4 O(h ), solutions we have six possible pairs which we can extrapolate to or we might take three solutions and extrapolate to we do not know what U 6 O(h ). Even though is, we can perform the extrapolations and look for an agreement to a greater number of significant figures in the extrapolated solutions than in the original finite difference solutions. Let us introduce a notation for extrapolated solutions where subscripts indicate the and u . 4 u ui 124 used. For example, is extrapolated from u 14 is extrapolated from u , u , and 1 2 1 u • 4 In this example we have selected the solution at (0. 35, 0. 1) for extrapolation. u (x,t) = The results are given in Table 4. In Table 5 we give the equivalent extrapolations for the flux V. We see that pairwise extrapolations give very good agreement, but for our particular data an extrapolation of three solutions appears not to yield much further improvement. Iterations in Newton 1 s method were stopped when the relative change in two succeeding iterates was less than 2 X 10 -7 except we did not compute the relative change for com- ponents less than 8 X 10 -9 in magnitude. Since the maximum value of the solution is very nearly unity, the maximum absolute error in any component should be 2 X 10 -7 , and for most components it should be even less. This assumes that round-off error can be neglected. We also performed extrapolations in Example 1, but the -98increase in agreement was not as great as in Example 2. Another computation in which the transition from 1 to 0 was spread out over a wider interval yielded better extrapolations. On the basis of a limited number of computations involving different initial functions and different values of v, we believe, although we cannot prove, that the deciding factor in whether or not Richardson extrapolation will yield a significant improvement is the ratio of the magnitudes of the derivatives of the solutions (including their initial values) to the inverse of the mesh spacing h. Examination of seven computations indicates that this ratio should not exceed 1/10. We therefore pro- pose that in the numerical computations of shocks in the presence of a small amount of diffusion, an examination of the solution should be made every few time steps to detect regions of rapid change and that new net points should be interpolated in accordance with the empirical ratio given above. -99Table 1, Part 1 PROPAGATION VELOCITIES OF POINTS OF THE SOLUTION TO BURGERS' EQUATION OBTAINED BY DIFFERENCING OVER ONE TIME STEP T u = .10 0.005 O.C15 0.025 0.035 0.045 0.055 0.065 0.075 J.085 0.095 0.105 0.115 0.125 0.135 0.145 0.155 .0 .165 0.175 0.185 0.195 0.205 0.7.15 0.225 0.235 0.245 0.255 0.265 0.275 0.285 0.295 J.305 0.315 0.":\25 0.335 0.345 0.355 0.3t5 0.375 0.385 0.395 0.131896 0.133366 0.134943 0.136955 0.139206 J.141434 0.144154 0.146672 3.149018 0.150915 0.154249 0.157861 0.161544 0.165826 0.169317 0.170938 J.175938 0.181671 0.188251 0.195237 0.197095 0.202681 0.213199 0.226338 J.240473 0.232511 0.251871 0.284354 0.295886 0.265001 0.352190 0.348102 0.276€67 0.514069 0.1822'32 0.543654 0.273807 0.477959 0.529732 0.248122 = .25 u = • 50 u = • 75 u = .90 0. 26 39 50 0.500229 0.500385 0.500418 0.500454 0.500494 J.5u0539 0.500589 0.500644 0.500704 0.500771 0.5J0844 0.500927 0.501012 0.501101 0.501209 0.501313 J.5J1428 0.501507 0.501691 0.501701 0.501983 0.501764 o. 50238 7 0.501641 0.502931 0.501104 0.50":\774 0.500032 0.504988 0.498274 0.506678 0.495878 J.508778 0.492992 0.511175 0.489895 0.513646 0.486864 J.516J05 0.484124 0.735932 0.735190 0.734038 0.732675 o. THl 70 0.729910 0.728056 o. 12 5949 0.723863 0.721798 0.718828 0.715464 0.712926 0.709069 0.704125 0.699780 o.6S5773 0.688387 0.680944 0.678189 0.666758 0.655675 0.656417 0.637038 0.629535 0.628673 0.594255 0.620296 0.579894 0.583154 o. 583511 0.520129 0.615186 0.486143 0.620607 0.464828 0.617004 0.449309 0.611680 0.438972 0.867148 0.865906 0.864276 0.862196 0.859946 0.857525 0.854938 0.852194 0.850064 o. 847822 0.844465 0.840790 0.836811 0.832552 0.828908 u o. 264944 0.266:J91 0.267448 0.268947 0.270197 0.272047 0.274142 0.276230 o. 278279 0.281238 0.284647 0.287128 0.291021 o. 296058 0.300l79 0.304443 0.312055 0.319649 0.322264 0.334374 o. 345836 0.344235 0.366175 0.371462 o. 3 7 36 49 0.4141<;6 0.374014 0.432928 0.405740 0.433709 0.464705 0.403774 0.533543 O."l60~12 0.564241 0.351342 0.591057 0.342Hl8 0.613707 0.8268~2 0.821377 0.815165 0.808207 0.80055 8 o. 799516 o. 79140 8 0.779761 0.765761 0.757152 0.7~6301 0 .11~844 0.703383 0.714271 0.699956 0.637966 0.679412 0.648747 0.538052 0.720912 0.514237 0.660511 0.5l4815 0.550754 0.606499 -100- Table 1, Part 2 VEL OCITIES OF POINTS CF ThE SOLUTION TO BURGERS' EQUATION OBTAINED BY DIFFERENCING OVER ONE TIME STEP PROPAGATIC~ T u :: .10 u = .25 u = .so 0.405 0.415 0.425 J.435 0.445 0.455 0.465 0.475 0.485 0.495 0.505 0.515 0.525 o. 535 0.545 0.555 0.565 0.575 0.585 0.595 0.(:05 0.(:15 0.625 0.635 0.877381 -0.072516 0.991424 -0.159028 1.096824 0.335297 0.632246 0.328991 0.647034 0.323949 0.658496 0.320022 0.667245 o. 317076 0.673741 0.314878 0.678557 0.313301 C.682015 0.312150 0.684542 0.518095 0.481803 0.519848 0.479933 0.521249 0.478493 J.522329 0.477415 0.523140 0.476632 0.523732 0.476071 0.524162 0.475679 0.524465 0.475406 0.524681 0.475220 0.524830 0.475092 0.524936 0.475007 0.525006 0.474949 0.525051 0.474911 0.525089 0.474885 0.525114 0.474869 0.525129 0.474857 0.525141 0.474851 0.525147 0.474846 0.525153 0.474843 0.525155 0.474841 0.~45 0.655 0.665 0.675 0.685 0.695 0.705 0.715 0.725 0.735 0.745 0.155 0.765 0.775 0.785 0.795 -0.2~6592 1.189696 -0.303239 1.26849 7 -0.357689 1.332260 -0.400985 1.382645 -0.433887 1.420 732 -0.458746 1.449442 -0.47667S 1.470C86 -0.489877 1.485286 -0.499025 1.495790 -0.505121 1.503500 -0.510187 1.508(:23 -0.513504 1.512446 -0.515613 1.514862 -0.517240 1.516745 -0.518203 1.517f45 -0.519004 1.518776 -0.519426 1.519256 -0.519826 0.311'\~2 0.686304 0.310772 0.687591 0.310384 c. 68 8460 0.310097 0.689104 0.309915 0.689521 o. 30977 5 0.689841 0.309692 0.690036 O.JOS623 0.690195 0.309587 0.690282 0.309553 C.690363 0.309538 0.690399 o. 109520 0.69Q441 = .75 u :: .90 0.61)5468 0.432635 0.599165 0.429077 0.593422 0.427235 0.588570 0.426374 0.5€4678 0.426016 0.581711 0.425894 0.579480 0.425874 a. 577886 0.425880 0.576715 0.425900 0.575922 0.425907 0.575334 0.425919 0.574961 0.425916 o.57467C 0.425920 0.574504 0.425912 0.574359 0.425914 0.574291 0.425906 0.574216 0.425908 0.574192 0.425901 0.574151 0.425903 0.574147 0.4258<i8 0.398020 0.724722 0.289208 0.804172 0.255958 0.813541 0.231211 0.819653 0.213138 0.823654 0.200089 0.826185 0.190796 0.827933 0.184145 0.828987 0.179499 0.829616 0.176159 0.830251 0.173886 u 0.8"1068~ 0.172218 o. 830845 0.171133 o. 831093 0.170299 0.8311~2 0.169798 0.831228 0.169376 0.831273 0.169157 0.831379 0.168937 0.831343 0.168851 0.831421 0.1687'\1 0.831380 -101Table 2, Part 1 PROPAGATION VELOCITIES OF POINTS OF THE SOLUTION TO BURGERS' EQUATION OBTAINEO BY DIFFERENCING OVER TkO TIME STEPS T u = .10 u = .25 u = .50 u = • 75 u = .90 0.010 0.020 0.030 0.040 o.u5o 0.060 0.010 0.080 0.090 0.100 0.110 0.120 0.130 0.140 0.150 0.160 0.170 0.180 0.190 0.200 0.210 0.220 0.230 0.240 0.250 0.260 0.27J 0.280 0.290 0.300 0.310 0.320 0.330 0.340 0.350 0.360 O.l70 0.380 0.390 0.400 0.132631 0.134154 0.135949 0.138080 0.140320 0.142794 0.145388 0.147820 0.149<:i66 0.152582 0.156055 0.159702 0.163685 0.167601 0.170157 0.173438 0.1788J4 0.184961 0.191744 0.196166 J.199888 0.2(7940 0.219768 0.2334J5 0.236492 0.242191 0.268ll2 0.290120 0.280443 0.308595 0.350146 0.312484 0.395468 0.348150 0.362943 0.408730 0.375883 0. 5 03 84 5 0.388927 0.562152 0.264447 0.265517 0.266769 0.268197 0.269572 0.271122 0.273094 0.275186 0.277254 0.279758 0.282942 0.285887 0.289074 0.293539 0.298218 0.302411 0.308249 0.315852 0.320956 0.128319 0.340105 0.345035 0.355205 0.368818 0.372555 0.393922 0.394105 0.403471 0.419334 0.419724 0.449207 0.434239 o. 46 8658 0.446927 o. 462276 0.457791 0.471199 0.466937 0.478262 0.474507 0.500307 0.500401 0.500436 0.500473 0.500516 0.500564 0.500616 0.500674 0.500737 0.5008J7 0.500885 0.500969 0.501056 0.501155 0.501261 0.501370 0.501467 .) .501599 0.501696 J.501842 0.5J1873 o. 502075 0.502014 0.502286 0.502017 0.502439 0.501903 0.502510 0.501631 o. 502476 0.501278 0.502328 0.500885 0.502083 0.500535 0.501770 0.500255 J.501434 0.500064 0.501109 0.735561 0.734614 0.733356 0.731922 0.730540 o. 728983 0.121002 0.724906 o. 12 2830 0.720313 0.717146 0.714195 0.710997 0.706597 0.701952 0.6c;7776 0.692080 0.684665 0.67Cj566 0.672473 0.661216 0.656046 0.646727 0.633286 0.629104 0.611464 0.607275 0.600095 0.581524 0.583332 0.551820 0.567657 0.550664 0.553375 0.542717 0.540916 0.533156 0.530494 0.525326 0.522220 0.866527 0.865091 0.863236 0.861071 0.858735 0.856231 0.853566 0.851129 0.848943 0.846143 0.842627 0.838800 0.834681 0.830730 0.827870 0.824104 0.818271 0.811686 0.804382 0.800037 0.795462 0.785584 0.772761 0.761456 0.756726 0.745072 0.718613 0.708827 J.7J7113 0.668961 0.65868~ 0.664079 0.593399 0.629492 0.617584 0.587374 0.597661 0.542784 0.578626 J.50225CJ -102Table 2, Part 2 PROPAGATICN VELOCITIES OF POINTS OF THE SOLUTION TO BURGERS' EQUATION CHTAINEO BY DIFFERENCING OVER TkO TIMF. STEPS T u = .10 0.410 0.420 0.43J 0.440 0.450 0.460 0.470 G.480 0.490 0.500 ;).510 0.520 0.402434 0.459454 0.416198 0.4688<;8 0.430116 0.476552 0.443229 0.482629 0.455404 0.487286 0.465638 0.490830 0.5~0 0.474~7<1 0.540 0.550 0.560 0.57J 0.580 0.590 0.600 0.610 J.620 0.630 0.640 0.650 0.660 0.670 0.680 C.690 J.7ou 0.710 0.720 0. 730 0.740 0.150 0.760 o.77u 0.780 0.790 0.493422 0.480993 0.495":448 0.486331 0.496704 0.490105 0.497704 0.493130 0.498383 0.495035 0.498890 J.4<i6656 0.499218 0.497560 0.499471 0.498417 .).499625 0.498811 0.49<1752 0.499271 0.499821 0.499421 0.499886 0.499675 0.499915 0.49<1715 = .25 u = • 50 u = • 75 u = .90 0.483771 0.480618 0.488J12 0.485491 0.491222 0.489259 0.493633 0.492160 0.49540d 0.494109 0.496717 0.495929 0.497658 0.497082 0.498346 0.497947 0.498828 0.498538 0.499181 0.498987 0.499422 J.499278 0.499600 0.499509 0.499718 0.499648 0.499808 0.499766 0.499864 0.499829 0.499909 0.499391 0.499934 0.499917 1.).499958 0.499950 0.499968 0.499959 0.499980 0.499949 0.500825 J.499890 0.500591 0.4998 71 0.500411 o. 499872 0.500277 0.499886 o. 500182 J.4999J1 0.500116 0.499920 J.500J72 0.499935 0.500041 0.499950 0.500025 0.499961 0.500014 o. 499971 J.5JOJJ6 0.499977 0.500001 0.499984 0.500000 0.499987 0.499999 0.499991 J.499999 0.499993 0.4999'19 0.499996 0.499999 0.499996 0.499999 0.499998 1).499199 0.499998 0.519051 0.515900 0.514121 0.511249 0.510328 0.507902 0.507472 0.505526 0.505347 0.503863 0.503802 0.502687 0.502677 0.501880 0.501383 0.501297 0.501307 0.500911 J.5J0914 0.500620 0.500626 0.500440 0.500438 o. 50 0293 J.500295 0.500212 0.5002C8 0.500135 0.5001":\6 0.500102 0.500098 0.500061 0.500062 0.500050 0.500046 0.500026 o. 50002 7 J. 5JJQ25 0.500022 0.561371 0.506965 0.546690 0.530065 o. 534 749 0.522376 L 0.5254~2 0.516395 o. 518 396 0.511871 0.513137 0.508490 0.509364 0.506039 0.506566 0.504243 0.504657 0.502987 0.503205 0.502068 0.502284 J.5~1450 0.501531 0.500989 0.501113 0.500696 o.sJo715 0.500465 0.50051':\ 0.500302 0.500324 0.500215 0.500268 0.500158 D.5JJ14J o. 50009 7 0.500136 0.500076 0.500055 -103- Table 3, Part I PROPAGATION VELOCITIES OF PCINTS OF THE SCLUTION TO BURGERS' eQUATION OBTAINED BY DIFFERENCING OVER CNE TIME STEP T u = .25 u :: .so u 0.005 0.015 0.025 0.035 J.045 0.055 O.C65 0.075 o.C85 0.095 0.105 o. 115 0.125 0.135 0.145 0.155 0.16 5 0.175 0.185 0.195 0.205 0.215 0.225 0.215 0.245 0.255 o. 26 5 0.275 0.285 0.295 0.305 o. 315 0.325 0.335 0.345 0.355 0.365 0.375 0.385 0.395 0.~92306 0.4998l6 ).499951 0.499956 0.607234 0.594442 0.5ti0939 0.568132 0.555963 0.545359 0.536250 0.528484 0.522167 0.516598 0.512309 0.508316 J.505455 o. 502517 0.500703 0.498617 0.497421 0.495904 0.495173 0.494071 0.493655 0.492861 0.492657 0.4<12096 0.492032 0.4<11649 0.491676 0.491432 0.491512 0.491379 0.491488 0.491444 0.491563 0.4915<13 0.491708 0.491802 0.491908 0.492027 0.492153 0.492357 J.405520 0.419353 0.432555 0.444864 0.455285 0.464216 0.471848 0.478C82 0.483595 0.487853 0.4<H822 0.494684 0.497534 0.499445 0.501464 0.502738 0.504203 0.504805 0.506253 0.5065~0 0.506812 0.507502 o. 50 84 72 0.507854 0.508197 0.508744 0.508641 0.508358 0.508976 0.508747 0.508325 0.508544 0.508810 G.508348 0.508091 0.508277 0.508178 0.507919 0.507722 0.49998~ J.500059 0.500159 0.500256 0.500321 0.500354 0.500363 0.500357 0.500342 J.500322 0.500299 0.500217 0.500255 0.500234 0.500214 0.500197 0.500180 0.500167 ·) .500151 0.500143 0.500127 0.500123 0.500116 0.500091 0.500127 0.500114 J.499986 0.500102 0.500266 J.499980 0.499854 0.500254 0.50J238 0.499854 0.500008 0.500234 0.499996 = • 75 -104Table 3, Part 2 PROPAGATICN VELOCITIES OF POINTS OF THE SOLUTION TO BURGERS' EQUATION OBTAINED BY DifFERENCING OVER ONE TIME STEP T u = .25 u = .50 0.405 0.415 0.425 0.507758 0.507514 0.507436 0.507186 0.507187 0.506891 0.506900 0.506590 0.506626 C.506292 0.506353 0.5060C2 0.506085 0.505718 0.505825 0.505444 0.505571 0.505178 0.505326 0.504921 0.50508':1 0.504674 0.50485':1 0.499991 0.500209 0.500019 0.499881 0.500178 0.500212 0.499904 0.499943 0.500198 0.500119 0.499942 0.500020 0.500127 0.500621 0.500007 0.500047 0.500072 0.500047 0.500039 0.500048 0.500051 0.500044 0.500044 0.500044 0.500044 0.500041 0.500041 0.500038 0.500037 0.5JJ033 0.500031 0.500027 0.5J0023 0.500017 0.500011 0.5JOJJ2 0.499994 0.499980 0.499967 0.499':148 0.4~5 0.445 0.455 0.465 0.475 0.485 0.495 0.505 0.515 0.525 0.535 0.545 0.555 0.565 0.575 o. 585 0.595 0.605 0.615 o. 62 5 0.615 0.645 0.655 0.665 0.675 0.685 0.695 0.705 a. 715 0.725 0.73 5 0.745 0.755 o. 765 0.775 0.785 0.1':15 0.5044~4 0.504635 0.504201 0.504416 0.503972 0.504200 0.503745 0.503982 0.50".\514 0.503756 0.503274 0.503516 0.503015 0.503251 0.502725 0.502943 0.502384 u = • 75 0.492~25 0.492563 0.492791 0.492921 0.492731 0.493249 0.493423 0.4933':12 0.493308 0.493952 0.493831 0.49394~ 0.494048 0.494456 0.494156 0.494633 0.4946':11 0.4948~4 0.4':14599 0.495273 0.495140 0.495264 0.495124 0.495752 0.4954':18 0.495754 0.495611 0.496131 0.4':15851 J.496222 0.496021 o. 496481 0.496202 0.4':16627 0.496372 0.496817 0.496529 0.496972 0.496680 0.4':17127 -105- Table 4 Richardson Extrapolations of Example 2 at (x, t) = (0.35, O.l) ul 0.50032235050 u2 0.50008107259 u3 0.50003613622 u4 0.50002035943 ul2 0.50000064662 ul3 o.5ooooo35944 ul4 0.50000022669 u23 0.50000018712 u24 o.5ooooo12171 u34 0.50000007499 ul23 0.50000012969 ul24 0. 5 00000086 72 u 0.50000005602 134 u234 u 1234 0.50000003761 0.50000003147 -106- Table 5 Richardson Extrapolations of Example 2 at (x, t) = (0.35, 0.1) vl -0.018540499230 v2 -0.018499265130 v3 -0.018491656850 v4 -0.018488994690 vl2 = -0.018485520430 vl3 -0.018485551553 vl4 -0.018485561054 v23 -0.018485570226 v24 -0.018485571210 v34 = -0.018485571913 vl23 -0.018485576451 vl24 -0.018485574595 vl34 v 234 -0.018485573270 vl234 -0.018485572210 -0.018485572475 -107- . m 0 . co 0 . [""'- 0 . CD 0 0' lJ) 0 • X <1) 1-< ;:l b.O .. . ~ . ::::t' 0 . ('(") 0 N N ('(") I w 0 0 . I w 0 . 00 ........ ('(") ........ II :::J z 0. 1 s·o cx )n 0 II I a· a 0 0 • -108- . m 0 . CD 0 . r0 . (D 0 0 ...... Ln ll> 0 ..tlD. . X '-< ;:::l ~ . ::t' 0 . (Y"') 0 N (Y"') I w 0 . N I 0 . w 0 0 Ln ........ 0 0 (Y"') II II :::::) z 0. l s·o cx )n I o·o 0 0 . -10 9 - . m 0 . (X) 0 . r-0 . (.D 0 lf) 0 • X ...... ...... a) '"' ;::l bD .... . =:1' 0 . (Y) 0 N N I (Y) . w I w 0 0 . (Y) (Y) (Y) ....... 0 0 (Y) II :::::> z o· t s·o cx )n II I o·o 0 0 • ~ -110- . m 0 . CD 0 . r.-0 . CD 0 N ...... Ln 0 • X Cl) ';:J""' .....bD ~ . :::1' 0 . (Y1 0 N (Y1 I w 0 N I 0 . w 0 . Ln N ........ (Y1 0 II :::> z 0. t s·o cx )n 0 II :r: o·o 0 0 . -111- . CD 1, l 1 0 l I -= ) ( ) ) f.- . r- - '' 0 ) ) --== ) 0 ) f- --== . CD - ( ( ) ) --=== 0 ) ~ > --== ==-('() ::::-- --== -( .....:::::::: ....... ::rw . 2: ::::::- - ) :> <s: ~ . l[) - 0 1--i t-- 0 f- N - 0 1- f- . (Y1 - . 0 0 I o·c I o·1 I s·o AllJIJl::lA o·o I s·o- . ...-4 0 . 0"1- (I) 1-4 ;:I tlD ...... ~ -112- . CD 0 . ('. 0 . CD 0 . Ln 0 :::t'w • L: 0~ f- "" ,_, C!) '"' b.O ;:j .. . ~ . (Y') 0 . 0J 0 . ......... 0 0 . oo·~----------~--------~----------~----------~0 r SL"O os·o sc·o 0'0 ,.\JJJOl~l\ -113- 0 0 . CD 0 . Lf) 0 1.[") :::::1'w • :L: 0 .,..._._. I- Q) ';:l"' .....0.0 ~ . (Y) 0 . 0J 0 0 0 . . L----L-----J----------~-----------L----~----~0 oo· t SL" 0 os·o A1IJOl3t\ sc·o o·o -11 4 - 0 . ["'-- 0 . CD 0 . Lf) 0 ::rw • L: a~ 1- -.o ....... Q) H ;:l ....OlJ ~ . (Y") 0 . (\J 0 0 0 . . ~--------~~--------~----------~--------~0 oo· r sco os·o J..liJOl~A o·o -115- . [" 0 . CD 0 . t.n 0 r:...... ::f'w 0 . ~ ...._. 1-- . (Y") 0 . 0J 0 . ....... 0 0 SL" os· A1IJOl3A . v J..< ;:l .....0.0 ~ -116- III. 4 Example 3: Interaction of Two Shocks In (1. 2) we gave a transformation which converts solutions of the heat equation into solutions of Burgers 1 equation. This transfor- mation has been used to advantage by Whitham [ 197 2] to construct an exact solution of Burgers' equation representing the overtaking of one shock by another. This solution is 1 (x-7) W(x, t) = {:o e- 3 (x- "S") - z-v-+e 1 (x-2) 4V +e In the limit as w0 (x) = For v =0 v 20v } 99t - 400v 1 (x-2) 1 + ze 1 (x-2) { e- +e 3t - Tbv 99t - 400V 3 (x- -g-) 3t - Tbv 20v 4v -~ r ( 4. 1) goes to zero we have the initial condition : ~ {~ for X < 0. 25 for 0 • 25 < X < 0. 5 for 0. 5 <X (4. 2) the solution would be a shock moving at speed 0. 75 and starting at x = 0. 25 overtaking a shock moving at speed 0. 3 and starting at x = 0. 5. The shocks would merge at (x, t) = ( 2/3, 5 /9) and continue as a single shock of strength 0. 9 and speed 0. 55. For a -117- small, positive v the initial condition would resemble a staircase function with the corners rounded. v = 3 X 10 In Example 3 we take -3 . The initial and boundary conditions are U(x,O) = W(x,O), (4. 3a) = W(O,t), (4.3b) U(1 ,t) = W(1 ,t). (4. 3c) U(O,t) We use uniform meshes consisting of 100, 200, 300, and 400 space steps over the interval 0 !S x !S 1. We always take k = h. Graphs of the solutions for each of the four nets are given in Figures 18 through 21. The curves are plotted at time intervals of 0. 1 so that the right- most curve, which goes off the right edge at about responds to U = 0. 96, cor- t = 1.2. As in Example 2 we have selected a point at which to perform all of the possible Richardson extrapolations with four solutions u 1 , u that 2 u , u 12 3 , and u • 4 These are given in Tables 6 and 7. is more accurate than u . 4 We notice In order to find if this might be true in general, we select twelve different points, and at each tenth of a time unit, for comparison. In Table 8, Part 1 we give u and , U, and u u3" In Part 2 we list u 4 u 12 . We see that best three digit accuracy as is the case with u 12 • 1 4 , u 2 has at We conclude that -118- u 12 gives roughly the same number of correct digits as u . 4 This is a significant result when we recall that the amount of computation required by various meshes is proportional to the square of the ratio of steps. as u , 1 In particular and u u 2 requires four times as much computing requires sixteen times as much as 4 is comparable to computing as u 4 u 4 Thus if u 12 in accuracy but requires only 5/16 as much , it will clearly be preferable to compute and to extrapolate rather than to compute We note incidentally that more accurate. u • 1 u 4 u , u , u and u 1 2 3 4 u 1 and with a very fine mesh. become progressively This means we have not yet refined the mesh to such an extent that the increased amount of computation causes significant round-off errors. tions involving u • 3 For reference, Part 3 gives two more extrapola- -119- Table 6 Richardson Extrapolations of Example 3 at (x, t) = (0.56, 0.2) u 1 u2 • 0.30225365244 = 0.30059271979 0.30026632734 0.3CXJ15060111 u12 = 0.30003907557 u 13 = 0.30001791170 ~4 0.30001039769 u23 0.30000521338 u24 0.30000322822 u34 0.30000181024 u123 0.30000098061 u124 0.30000083839 = 0.30000073681 C.3 0000067586 u 1234 Exact So1t:tion: 0.30000065555 u c 0.30000056686 -120- Table 7 Richardson Extrapolations of Exanple 3 at (x, t) = ( 0 .56, 0.2) vl -0.020013464?31 v2 -0.020008700290 V3 -0.020004435706 v4 = -0.020002664316 -0.020007112076 vl2 c -0.020003307053 vl4 = -0.020001944275 v "' -0.020001024039 v 13 23 -0.020000652325 v24 v34 = -0.020000263034 v123 v = -0.020000221675 = -0.020000192132 124 vl34 -0.020000386815 v2J4 vl234 Exact Solution: -0.020000174406 -0.020000168498 v - 0 . 0200001822 0h -121- Table 8, Part 1 Example 3 t X ul u2 U3 0.1 0.50 0.45162932809 0.45213909765 0.!!52238136!13 0.2 0.50 o.L9260003939 0.49284055156 0.49288720790 o.; 0.55 o.L!737h691548 0.47464422617 o.47h8186oo4o o.L 0.60 0.42107292706 o.L230l439573 0.42335725720 0.5 0.65 0.33964540345 0.34073709968 0 •.34090349587 0.6 0.70 0.33h73215370 0.32641934827 0.32523826859 0.7 0.75 0.!,4725589532 o.L3ll4939356 0.42808375484 0.8 0.80 0.60225117856 o.591J!tl49643 0.58882496729 0.9 0.85 0.?4645749560 0.74395011947 0.74323437130 1.0 0.90 0.85541977691 0.85650835749 0.85674654929 1.1 0.97 0.44512538419 0.42816.511.566 0.42491123354 1.2 0.99 0.99248461873 0.98966.51810.5 0.98983381566 -122- Table 8, Part 2 Example 3 t X u4 u u O.l 0.)0 0.45227329626 0.45231905503 0.!.6230902084 0.? 0.)0 0.49290374042 o.h92925209o6 0.49292072228 0.3 0.)5 0.47488025263 0.47496004739 0.47494332873 0.4 0.60 0.42347622629 0.42362844757 0.42366155195 0.5 0.65 0 •.34095825381 0.34102590553 0.34110099842 0.6 0.70 0.32485439735 0.32438289934 0.32364841313 0.7 0.75 0.42702022924 0.42566448866 0.42578055964 o.e 0.80 0.)8790276994 0.58668558549 0.)8770493572 0.9 o.e5 0.74296136413 0.74259164202 0.74311432743 1.0 0.90 0.85683280384 0.85694563756 0.85687121768 1.1 0.97 o.l.o2378437932 0.42235376627 0.!.,2251169282 1.2 0.99 0.99015457747 0. 99161484 77 9 0.988?2536849 12 -123- Table 8, Part 3 Example 3 t X u23 u 123 0.1 o.so 0.45231736745 0.45231841078 0.2 0.50 0.49292453297 0.1.:9292500931 0.3 0.55 0.47495809978 0.1.:7495994616 O.h 0.60 0.42363154638 o.h2362779568 0.5 0.65 0.34103661282 0.34102856462 0.6 o. 70 0.32429340485 0.32437402881 0.7 0.75 0.42563124386 0.425612)7939 o.e 0.80 0.58681174398 0.)8670009501 0.9 0.85 0.74266177276 0.74260520343 1.0 0.90 0.85693710273 0.85694533836 1.1 0.97 0.1.:2230812784 0.1.:2228268222 1.2 0.99 0.98996872335 0.99012414271 -124- . c..o 0 00 ...... lJ) 0 • X <1) "-< ;j .....t:lD ~ . ::1' 0 N (Y') I LL.J 0 . N I LL.J 0 0 C) C) ~ ~ 0 (Y') II 0. 1 s·o cx )n • II ::) z . :r: o·o 0 0 • -125- . ([) 0 0' ...... lf) 0 • X Q) !-< ;::l tlD ...... ~ . :::r 0 N ('f") I w 0 N 0 . I w 0 . 0 lf) ........ 0 0 ('f") II :::::> :z: 0. l s·o cx )n II :r:: o·o 0 0 • -126- . (i) 0 0 N IJ) 0 • X v '"' ;:J .....OJ) ~ . ::1' 0 N (Y) I w 0 . N I 0 . w (Y) (Y) (Y) ........ 0 0 (Y) II II :::::J z j_ 0. l s·o cx )n I o·o 0 0 . -127- . CD 0 ....... N lJ) 0 • X <l.l J..; ::l 00 ...... ~ . :::t' 0 N (r) I w 0 . N I 0 . w 0 lJ) N ,_.., 0 0 (r) II ~ z 0. T s·o cx )n II ::r:: o·o 0 0 . -128- REFERENCES Cole, Julian D. [ 1951]: "On a Quasi-Linear Parabolic Equation Occurring in Aerodynamics." Quarterly of Applied Mathematics, Vol. 9, No. 3, pp. 225-236. Isaacson, Eugene, and Herbert Bishop K~ller [ 1966]: Analysis of Numerical Methods. John Wiley & Sons, Inc., New York. Keller, Herbert B. [ 1971 ]: "A New Difference Scheme for Parabolic Problems. " Numerical Solution of Partial Differential Equations -- II: S"Nnspade 1970. Edited by Bert Hubbard. Academic Press, ew York, pp. 327-350. Lees, Milton [ 1960]: ''Energy Inequalities for the Solution of Differential Equations." American Mathematical Society Transactions, Vol. 94, pp. 58-73. Lees, Milton [ 1961]: "Energy Inequalities for Partial Difference Equations." Lecture notes, Mathematics Department, California Institute of Technology, Pasadena. (Mimeographed.) Strang, W. Gilbert [ 1960]: "Difference Methods for Mixed BoundaryValue Problems." Duke Mathematical Journal, Vol. 27, No. 2, pp. 221-231. Swartz, Blair, and Burton Wendroff [ 1969] : ''Generalized FiniteDifference Schemes." Mathematics of Computation, Vol. 23, No. 105, pp. 37-49. Varah, James M. [ 1971]: "Stability of Difference Approximations to the Mixed Initial Boundary Value Problems for Parabolic Systems." SIAM Journal of Numerical Analysis, Vol. 8, No. 3, pp. 598-615. Varah, James M. [ 1972] : On the Solution of Block Tridiagonal Systems Arising from Certain Finite-Difference Equations. Technical Report No. 72-02, Department of Computer Science, University of British Columbia, Vancouver. Whitham, Gerald B. [ 1972] : private communication. -129Yanenko, N. N. [ 1971]: The Method of Fractional Steps: The Solution of Problems of Mathematical Ph sics in Several Variables. Trans ated by orthrop orporate a oratortes. rans edited by Maurice Holt. Springer- Verlag, New York.