No title

Fundamentals of Computational Fluid Dynamics Harvard Lomax and Thomas H. Pulliam NASA Ames Research Center David W. Zingg University of Toronto Instit

Author: wilber david huamani maihuire

3 downloads 117 Views 4MB Size

Report

DOWNLOAD PDF

Recommend Stories

No title

Story Transcript

Fundamentals of Computational Fluid Dynamics Harvard Lomax and Thomas H. Pulliam NASA Ames Research Center David W. Zingg University of Toronto Institute for Aerospace Studies August 26, 1999

Contents 1 INTRODUCTION

1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Problem Speci cation and Geometry Preparation . . . . . . 1.2.2 Selection of Governing Equations and Boundary Conditions 1.2.3 Selection of Gridding Strategy and Numerical Method . . . 1.2.4 Assessment and Interpretation of Results . . . . . . . . . . . 1.3 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 CONSERVATION LAWS AND THE MODEL EQUATIONS 2.1 Conservation Laws . . . . . . . . . . . . 2.2 The Navier-Stokes and Euler Equations . 2.3 The Linear Convection Equation . . . . 2.3.1 Dierential Form . . . . . . . . . 2.3.2 Solution in Wave Space . . . . . . 2.4 The Diusion Equation . . . . . . . . . . 2.4.1 Dierential Form . . . . . . . . . 2.4.2 Solution in Wave Space . . . . . . 2.5 Linear Hyperbolic Systems . . . . . . . . 2.6 Problems . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

3 FINITE-DIFFERENCE APPROXIMATIONS 3.1 Meshes and Finite-Dierence Notation 3.2 Space Derivative Approximations . . . 3.3 Finite-Dierence Operators . . . . . . 3.3.1 Point Dierence Operators . . . 3.3.2 Matrix Dierence Operators . . 3.3.3 Periodic Matrices . . . . . . . . 3.3.4 Circulant Matrices . . . . . . . iii

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

1

1 2 2 3 3 4 4 5

7

7 8 12 12 13 14 14 15 16 17

21

21 24 25 25 25 29 30

3.4 Constructing Dierencing Schemes of Any Order . . . . . 3.4.1 Taylor Tables . . . . . . . . . . . . . . . . . . . . 3.4.2 Generalization of Dierence Formulas . . . . . . . 3.4.3 Lagrange and Hermite Interpolation Polynomials 3.4.4 Practical Application of Pade Formulas . . . . . . 3.4.5 Other Higher-Order Schemes . . . . . . . . . . . . 3.5 Fourier Error Analysis . . . . . . . . . . . . . . . . . . . 3.5.1 Application to a Spatial Operator . . . . . . . . . 3.6 Dierence Operators at Boundaries . . . . . . . . . . . . 3.6.1 The Linear Convection Equation . . . . . . . . . 3.6.2 The Diusion Equation . . . . . . . . . . . . . . . 3.7 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

31 31 34 35 37 38 39 39 43 44 46 47

4 THE SEMI-DISCRETE APPROACH

51

5 FINITE-VOLUME METHODS

71

4.1 Reduction of PDE's to ODE's . . . . . . . . . . . . . . . . . . . . . . 4.1.1 The Model ODE's . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 The Generic Matrix Form . . . . . . . . . . . . . . . . . . . . 4.2 Exact Solutions of Linear ODE's . . . . . . . . . . . . . . . . . . . . 4.2.1 Eigensystems of Semi-Discrete Linear Forms . . . . . . . . . . 4.2.2 Single ODE's of First- and Second-Order . . . . . . . . . . . . 4.2.3 Coupled First-Order ODE's . . . . . . . . . . . . . . . . . . . 4.2.4 General Solution of Coupled ODE's with Complete Eigensystems 4.3 Real Space and Eigenspace . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 De nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Eigenvalue Spectrums for Model ODE's . . . . . . . . . . . . . 4.3.3 Eigenvectors of the Model Equations . . . . . . . . . . . . . . 4.3.4 Solutions of the Model ODE's . . . . . . . . . . . . . . . . . . 4.4 The Representative Equation . . . . . . . . . . . . . . . . . . . . . . 4.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Model Equations in Integral Form . . . . . . . . . . . . . . . . . . . 5.2.1 The Linear Convection Equation . . . . . . . . . . . . . . . 5.2.2 The Diusion Equation . . . . . . . . . . . . . . . . . . . . . 5.3 One-Dimensional Examples . . . . . . . . . . . . . . . . . . . . . . 5.3.1 A Second-Order Approximation to the Convection Equation 5.3.2 A Fourth-Order Approximation to the Convection Equation 5.3.3 A Second-Order Approximation to the Diusion Equation . 5.4 A Two-Dimensional Example . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

52 52 53 54 54 55 57 59 61 61 62 63 65 67 68 72 73 73 74 74 75 77 78 80

5.5 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 TIME-MARCHING METHODS FOR ODE'S

6.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Converting Time-Marching Methods to OE's . . . . . . . . . . . . 6.3 Solution of Linear OE's With Constant Coecients . . . . . . . . 6.3.1 First- and Second-Order Dierence Equations . . . . . . . . 6.3.2 Special Cases of Coupled First-Order Equations . . . . . . . 6.4 Solution of the Representative OE's . . . . . . . . . . . . . . . . 6.4.1 The Operational Form and its Solution . . . . . . . . . . . . 6.4.2 Examples of Solutions to Time-Marching OE's . . . . . . . 6.5 The ; Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.1 Establishing the Relation . . . . . . . . . . . . . . . . . . . . 6.5.2 The Principal -Root . . . . . . . . . . . . . . . . . . . . . . 6.5.3 Spurious -Roots . . . . . . . . . . . . . . . . . . . . . . . . 6.5.4 One-Root Time-Marching Methods . . . . . . . . . . . . . . 6.6 Accuracy Measures of Time-Marching Methods . . . . . . . . . . . 6.6.1 Local and Global Error Measures . . . . . . . . . . . . . . . 6.6.2 Local Accuracy of the Transient Solution (er ; jj ; er! ) . . . 6.6.3 Local Accuracy of the Particular Solution (er ) . . . . . . . 6.6.4 Time Accuracy For Nonlinear Applications . . . . . . . . . . 6.6.5 Global Accuracy . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Linear Multistep Methods . . . . . . . . . . . . . . . . . . . . . . . 6.7.1 The General Formulation . . . . . . . . . . . . . . . . . . . . 6.7.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7.3 Two-Step Linear Multistep Methods . . . . . . . . . . . . . 6.8 Predictor-Corrector Methods . . . . . . . . . . . . . . . . . . . . . . 6.9 Runge-Kutta Methods . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Implementation of Implicit Methods . . . . . . . . . . . . . . . . . . 6.10.1 Application to Systems of Equations . . . . . . . . . . . . . 6.10.2 Application to Nonlinear Equations . . . . . . . . . . . . . . 6.10.3 Local Linearization for Scalar Equations . . . . . . . . . . . 6.10.4 Local Linearization for Coupled Sets of Nonlinear Equations 6.11 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 STABILITY OF LINEAR SYSTEMS 7.1 Dependence on the Eigensystem 7.2 Inherent Stability of ODE's . . 7.2.1 The Criterion . . . . . . 7.2.2 Complete Eigensystems .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83

85

86 87 88 89 90 91 91 92 93 93 95 95 96 97 97 98 99 100 101 102 102 103 105 106 107 110 110 111 112 115 117

121

122 123 123 123

7.3

7.4 7.5 7.6

7.7

7.8 7.9

7.2.3 Defective Eigensystems . . . . . . . . . . . . . . Numerical Stability of OE 's . . . . . . . . . . . . . . 7.3.1 The Criterion . . . . . . . . . . . . . . . . . . . 7.3.2 Complete Eigensystems . . . . . . . . . . . . . . 7.3.3 Defective Eigensystems . . . . . . . . . . . . . . Time-Space Stability and Convergence of OE's . . . . Numerical Stability Concepts in the Complex -Plane . 7.5.1 -Root Traces Relative to the Unit Circle . . . 7.5.2 Stability for Small t . . . . . . . . . . . . . . . Numerical Stability Concepts in the Complex h Plane 7.6.1 Stability for Large h. . . . . . . . . . . . . . . . 7.6.2 Unconditional Stability, A-Stable Methods . . . 7.6.3 Stability Contours in the Complex h Plane. . . Fourier Stability Analysis . . . . . . . . . . . . . . . . 7.7.1 The Basic Procedure . . . . . . . . . . . . . . . 7.7.2 Some Examples . . . . . . . . . . . . . . . . . . 7.7.3 Relation to Circulant Matrices . . . . . . . . . . Consistency . . . . . . . . . . . . . . . . . . . . . . . . Problems . . . . . . . . . . . . . . . . . . . . . . . . . .

8 CHOICE OF TIME-MARCHING METHODS

8.1 Stiness De nition for ODE's . . . . . . . . . . . . 8.1.1 Relation to -Eigenvalues . . . . . . . . . . 8.1.2 Driving and Parasitic Eigenvalues . . . . . . 8.1.3 Stiness Classi cations . . . . . . . . . . . . 8.2 Relation of Stiness to Space Mesh Size . . . . . . 8.3 Practical Considerations for Comparing Methods . 8.4 Comparing the Eciency of Explicit Methods . . . 8.4.1 Imposed Constraints . . . . . . . . . . . . . 8.4.2 An Example Involving Diusion . . . . . . . 8.4.3 An Example Involving Periodic Convection . 8.5 Coping With Stiness . . . . . . . . . . . . . . . . 8.5.1 Explicit Methods . . . . . . . . . . . . . . . 8.5.2 Implicit Methods . . . . . . . . . . . . . . . 8.5.3 A Perspective . . . . . . . . . . . . . . . . . 8.6 Steady Problems . . . . . . . . . . . . . . . . . . . 8.7 Problems . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

123 124 124 125 125 125 128 128 132 135 135 136 137 141 141 142 143 143 146

149

149 149 151 151 152 153 154 154 154 155 158 158 159 160 160 161

9 RELAXATION METHODS

9.1 Formulation of the Model Problem . . . . . . . . . . . . . . 9.1.1 Preconditioning the Basic Matrix . . . . . . . . . . . 9.1.2 The Model Equations . . . . . . . . . . . . . . . . . . 9.2 Classical Relaxation . . . . . . . . . . . . . . . . . . . . . . 9.2.1 The Delta Form of an Iterative Scheme . . . . . . . . 9.2.2 The Converged Solution, the Residual, and the Error 9.2.3 The Classical Methods . . . . . . . . . . . . . . . . . 9.3 The ODE Approach to Classical Relaxation . . . . . . . . . 9.3.1 The Ordinary Dierential Equation Formulation . . . 9.3.2 ODE Form of the Classical Methods . . . . . . . . . 9.4 Eigensystems of the Classical Methods . . . . . . . . . . . . 9.4.1 The Point-Jacobi System . . . . . . . . . . . . . . . . 9.4.2 The Gauss-Seidel System . . . . . . . . . . . . . . . . 9.4.3 The SOR System . . . . . . . . . . . . . . . . . . . . 9.5 Nonstationary Processes . . . . . . . . . . . . . . . . . . . . 9.6 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

163

164 164 166 168 168 168 169 170 170 172 173 174 176 180 182 187

10 MULTIGRID

191

11 NUMERICAL DISSIPATION

203

10.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 10.1.1 Eigenvector and Eigenvalue Identi cation with Space Frequencies191 10.1.2 Properties of the Iterative Method . . . . . . . . . . . . . . . 192 10.2 The Basic Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 10.3 A Two-Grid Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 10.4 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202 11.1 11.2 11.3 11.4

One-Sided First-Derivative Space Dierencing The Modi ed Partial Dierential Equation . . The Lax-Wendro Method . . . . . . . . . . . Upwind Schemes . . . . . . . . . . . . . . . . 11.4.1 Flux-Vector Splitting . . . . . . . . . . 11.4.2 Flux-Dierence Splitting . . . . . . . . 11.5 Arti cial Dissipation . . . . . . . . . . . . . . 11.6 Problems . . . . . . . . . . . . . . . . . . . . .

12 SPLIT AND FACTORED FORMS

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

204 205 207 209 210 212 213 214

217

12.1 The Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217 12.2 Factoring Physical Representations | Time Splitting . . . . . . . . . 218 12.3 Factoring Space Matrix Operators in 2{D . . . . . . . . . . . . . . . . 220

12.4 12.5 12.6 12.7

12.3.1 Mesh Indexing Convention . . . . . . . . . . . 12.3.2 Data Bases and Space Vectors . . . . . . . . . 12.3.3 Data Base Permutations . . . . . . . . . . . . 12.3.4 Space Splitting and Factoring . . . . . . . . . Second-Order Factored Implicit Methods . . . . . . . Importance of Factored Forms in 2 and 3 Dimensions The Delta Form . . . . . . . . . . . . . . . . . . . . . Problems . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

13 LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS

13.1 The Representative Equation for Circulant Operators . . . . . . . 13.2 Example Analysis of Circulant Systems . . . . . . . . . . . . . . . 13.2.1 Stability Comparisons of Time-Split Methods . . . . . . . 13.2.2 Analysis of a Second-Order Time-Split Method . . . . . . 13.3 The Representative Equation for Space-Split Operators . . . . . . 13.4 Example Analysis of 2-D Model Equations . . . . . . . . . . . . . 13.4.1 The Unfactored Implicit Euler Method . . . . . . . . . . . 13.4.2 The Factored Nondelta Form of the Implicit Euler Method 13.4.3 The Factored Delta Form of the Implicit Euler Method . . 13.4.4 The Factored Delta Form of the Trapezoidal Method . . . 13.5 Example Analysis of the 3-D Model Equation . . . . . . . . . . . 13.6 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . .

220 221 221 223 226 226 228 229

233

233 234 234 237 238 242 242 243 243 244 245 247

A USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA 249 A.1 A.2 A.3 A.4 A.5

Notation . . . . . . . . . . De nitions . . . . . . . . . Algebra . . . . . . . . . . Eigensystems . . . . . . . Vector and Matrix Norms

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

B SOME PROPERTIES OF TRIDIAGONAL MATRICES B.1 B.2 B.3 B.4

Standard Eigensystem for Simple Tridiagonals . . Generalized Eigensystem for Simple Tridiagonals . The Inverse of a Simple Tridiagonal . . . . . . . . Eigensystems of Circulant Matrices . . . . . . . . B.4.1 Standard Tridiagonals . . . . . . . . . . . B.4.2 General Circulant Systems . . . . . . . . . B.5 Special Cases Found From Symmetries . . . . . . B.6 Special Cases Involving Boundary Conditions . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

249 250 251 251 254

257

257 258 259 260 260 261 262 263

C THE HOMOGENEOUS PROPERTY OF THE EULER EQUATIONS265

Chapter 1 INTRODUCTION 1.1

Motivation

The material in this book originated from attempts to understand and systemize numerical solution techniques for the partial dierential equations governing the physics of uid ow. As time went on and these attempts began to crystallize, underlying constraints on the nature of the material began to form. The principal such constraint was the demand for uni cation. Was there one mathematical structure which could be used to describe the behavior and results of most numerical methods in common use in the eld of uid dynamics? Perhaps the answer is arguable, but the authors believe the answer is armative and present this book as justi cation for that belief. The mathematical structure is the theory of linear algebra and the attendant eigenanalysis of linear systems. The ultimate goal of the eld of computational uid dynamics (CFD) is to understand the physical events that occur in the ow of uids around and within designated objects. These events are related to the action and interaction of phenomena such as dissipation, diusion, convection, shock waves, slip surfaces, boundary layers, and turbulence. In the eld of aerodynamics, all of these phenomena are governed by the compressible Navier-Stokes equations. Many of the most important aspects of these relations are nonlinear and, as a consequence, often have no analytic solution. This, of course, motivates the numerical solution of the associated partial dierential equations. At the same time it would seem to invalidate the use of linear algebra for the classi cation of the numerical methods. Experience has shown that such is not the case. As we shall see in a later chapter, the use of numerical methods to solve partial dierential equations introduces an approximation that, in eect, can change the form of the basic partial dierential equations themselves. The new equations, which 1

CHAPTER 1. INTRODUCTION

2

are the ones actually being solved by the numerical process, are often referred to as the modi ed partial dierential equations. Since they are not precisely the same as the original equations, they can, and probably will, simulate the physical phenomena listed above in ways that are not exactly the same as an exact solution to the basic partial dierential equation. Mathematically, these dierences are usually referred to as truncation errors. However, the theory associated with the numerical analysis of

uid mechanics was developed predominantly by scientists deeply interested in the physics of uid ow and, as a consequence, these errors are often identi ed with a particular physical phenomenon on which they have a strong eect. Thus methods are said to have a lot of \arti cial viscosity" or said to be highly dispersive. This means that the errors caused by the numerical approximation result in a modi ed partial dierential equation having additional terms that can be identi ed with the physics of dissipation in the rst case and dispersion in the second. There is nothing wrong, of course, with identifying an error with a physical process, nor with deliberately directing an error to a speci c physical process, as long as the error remains in some engineering sense \small". It is safe to say, for example, that most numerical methods in practical use for solving the nondissipative Euler equations create a modi ed partial dierential equation that produces some form of dissipation. However, if used and interpreted properly, these methods give very useful information. Regardless of what the numerical errors are called, if their eects are not thoroughly understood and controlled, they can lead to serious diculties, producing answers that represent little, if any, physical reality. This motivates studying the concepts of stability, convergence, and consistency. On the other hand, even if the errors are kept small enough that they can be neglected (for engineering purposes), the resulting simulation can still be of little practical use if inecient or inappropriate algorithms are used. This motivates studying the concepts of stiness, factorization, and algorithm development in general. All of these concepts we hope to clarify in this book. 1.2

Background

The eld of computational uid dynamics has a broad range of applicability. Independent of the speci c application under study, the following sequence of steps generally must be followed in order to obtain a satisfactory solution.

1.2.1 Problem Speci cation and Geometry Preparation

The rst step involves the speci cation of the problem, including the geometry, ow conditions, and the requirements of the simulation. The geometry may result from

1.2. BACKGROUND

3

measurements of an existing con guration or may be associated with a design study. Alternatively, in a design context, no geometry need be supplied. Instead, a set of objectives and constraints must be speci ed. Flow conditions might include, for example, the Reynolds number and Mach number for the ow over an airfoil. The requirements of the simulation include issues such as the level of accuracy needed, the turnaround time required, and the solution parameters of interest. The rst two of these requirements are often in con ict and compromise is necessary. As an example of solution parameters of interest in computing the ow eld about an airfoil, one may be interested in i) the lift and pitching moment only, ii) the drag as well as the lift and pitching moment, or iii) the details of the ow at some speci c location.

1.2.2 Selection of Governing Equations and Boundary Conditions

Once the problem has been speci ed, an appropriate set of governing equations and boundary conditions must be selected. It is generally accepted that the phenomena of importance to the eld of continuum uid dynamics are governed by the conservation of mass, momentum, and energy. The partial dierential equations resulting from these conservation laws are referred to as the Navier-Stokes equations. However, in the interest of eciency, it is always prudent to consider solving simpli ed forms of the Navier-Stokes equations when the simpli cations retain the physics which are essential to the goals of the simulation. Possible simpli ed governing equations include the potential- ow equations, the Euler equations, and the thin-layer Navier-Stokes equations. These may be steady or unsteady and compressible or incompressible. Boundary types which may be encountered include solid walls, in ow and out ow boundaries, periodic boundaries, symmetry boundaries, etc. The boundary conditions which must be speci ed depend upon the governing equations. For example, at a solid wall, the Euler equations require ow tangency to be enforced, while the Navier-Stokes equations require the no-slip condition. If necessary, physical models must be chosen for processes which cannot be simulated within the speci ed constraints. Turbulence is an example of a physical process which is rarely simulated in a practical context (at the time of writing) and thus is often modelled. The success of a simulation depends greatly on the engineering insight involved in selecting the governing equations and physical models based on the problem speci cation.

1.2.3 Selection of Gridding Strategy and Numerical Method

Next a numerical method and a strategy for dividing the ow domain into cells, or elements, must be selected. We concern ourselves here only with numerical methods requiring such a tessellation of the domain, which is known as a grid, or mesh.

CHAPTER 1. INTRODUCTION

4

Many dierent gridding strategies exist, including structured, unstructured, hybrid, composite, and overlapping grids. Furthermore, the grid can be altered based on the solution in an approach known as solution-adaptive gridding. The numerical methods generally used in CFD can be classi ed as nite-dierence, nite-volume, nite-element, or spectral methods. The choices of a numerical method and a gridding strategy are strongly interdependent. For example, the use of nite-dierence methods is typically restricted to structured grids. Here again, the success of a simulation can depend on appropriate choices for the problem or class of problems of interest.

1.2.4 Assessment and Interpretation of Results

Finally, the results of the simulation must be assessed and interpreted. This step can require post-processing of the data, for example calculation of forces and moments, and can be aided by sophisticated ow visualization tools and error estimation techniques. It is critical that the magnitude of both numerical and physical-model errors be well understood. 1.3

Overview

It should be clear that successful simulation of uid ows can involve a wide range of issues from grid generation to turbulence modelling to the applicability of various simpli ed forms of the Navier-Stokes equations. Many of these issues are not addressed in this book. Some of them are presented in the books by Anderson, Tannehill, and Pletcher [1] and Hirsch [2]. Instead we focus on numerical methods, with emphasis on nite-dierence and nite-volume methods for the Euler and Navier-Stokes equations. Rather than presenting the details of the most advanced methods, which are still evolving, we present a foundation for developing, analyzing, and understanding such methods. Fortunately, to develop, analyze, and understand most numerical methods used to nd solutions for the complete compressible Navier-Stokes equations, we can make use of much simpler expressions, the so-called \model" equations. These model equations isolate certain aspects of the physics contained in the complete set of equations. Hence their numerical solution can illustrate the properties of a given numerical method when applied to a more complicated system of equations which governs similar physical phenomena. Although the model equations are extremely simple and easy to solve, they have been carefully selected to be representative, when used intelligently, of diculties and complexities that arise in realistic two- and three-dimensional uid

ow simulations. We believe that a thorough understanding of what happens when

1.4. NOTATION

5

numerical approximations are applied to the model equations is a major rst step in making con dent and competent use of numerical approximations to the Euler and Navier-Stokes equations. As a word of caution, however, it should be noted that, although we can learn a great deal by studying numerical methods as applied to the model equations and can use that information in the design and application of numerical methods to practical problems, there are many aspects of practical problems which can only be understood in the context of the complete physical systems. 1.4

Notation

The notation is generally explained as it is introduced. Bold type is reserved for real physical vectors, such as velocity. The vector symbol ~ is used for the vectors (or column matrices) which contain the values of the dependent variable at the nodes of a grid. Otherwise, the use of a vector consisting of a collection of scalars should be apparent from the context and is not identi ed by any special notation. For example, the variable u can denote a scalar Cartesian velocity component in the Euler and Navier-Stokes equations, a scalar quantity in the linear convection and diusion equations, and a vector consisting of a collection of scalars in our presentation of hyperbolic systems. Some of the abbreviations used throughout the text are listed and de ned below. PDE ODE OE RHS P.S. S.S. k-D ~ bc O()

Partial dierential equation Ordinary dierential equation Ordinary dierence equation Right-hand side Particular solution of an ODE or system of ODE's Fixed (time-invariant) steady-state solution k-dimensional space Boundary conditions, usually a vector A term of order (i.e., proportional to)

6

CHAPTER 1. INTRODUCTION

Chapter 2 CONSERVATION LAWS AND THE MODEL EQUATIONS We start out by casting our equations in the most general form, the integral conservation-law form, which is useful in understanding the concepts involved in nite-volume schemes. The equations are then recast into divergence form, which is natural for nite-dierence schemes. The Euler and Navier-Stokes equations are brie y discussed in this Chapter. The main focus, though, will be on representative model equations, in particular, the convection and diusion equations. These equations contain many of the salient mathematical and physical features of the full Navier-Stokes equations. The concepts of convection and diusion are prevalent in our development of numerical methods for computational uid dynamics, and the recurring use of these model equations allows us to develop a consistent framework of analysis for consistency, accuracy, stability, and convergence. The model equations we study have two properties in common. They are linear partial dierential equations (PDE's) with coecients that are constant in both space and time, and they represent phenomena of importance to the analysis of certain aspects of uid dynamic problems.

2.1 Conservation Laws Conservation laws, such as the Euler and Navier-Stokes equations and our model equations, can be written in the following integral form: Z

V (t2 )

QdV ;

Z

V (t1 )

QdV +

Z t2 I

t1

S (t)

n:FdSdt =

Z t2 Z

t1

V (t)

PdV dt

(2.1)

In this equation, Q is a vector containing the set of variables which are conserved, e.g., mass, momentum, and energy, per unit volume. The equation is a statement of 7

8

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

the conservation of these quantities in a nite region of space with volume V (t) and surface area S (t) over a nite interval of time t2 ; t1 . In two dimensions, the region of space, or cell, is an area A(t) bounded by a closed contour C (t). The vector n is a unit vector normal to the surface pointing outward, F is a set of vectors, or tensor, containing the ux of Q per unit area per unit time, and P is the rate of production of Q per unit volume per unit time. If all variables are continuous in time, then Eq. 2.1 can be rewritten as d Z QdV + I n:FdS = Z PdV (2.2) dt V (t) S (t) V (t) Those methods which make various numerical approximations of the integrals in Eqs. 2.1 and 2.2 and nd a solution for Q on that basis are referred to as nite-volume methods. Many of the advanced codes written for CFD applications are based on the nite-volume concept. On the other hand, a partial derivative form of a conservation law can also be derived. The divergence form of Eq. 2.2 is obtained by applying Gauss's theorem to the ux integral, leading to @Q + r:F = P (2.3) @t where r: is the well-known divergence operator given, in Cartesian coordinates, by ! @ @ @ (2.4) r: i @x + j @y + k @z : and i; j, and k are unit vectors in the x; y, and z coordinate directions, respectively. Those methods which make various approximations of the derivatives in Eq. 2.3 and nd a solution for Q on that basis are referred to as nite-dierence methods.

2.2 The Navier-Stokes and Euler Equations The Navier-Stokes equations form a coupled system of nonlinear PDE's describing the conservation of mass, momentum and energy for a uid. For a Newtonian uid in one dimension, they can be written as @Q + @E = 0 (2.5) @t @x with 2 3 2 3 u 3 2 0 6 7 6 7 6 Q = 66 u 777 ; 4 5

e

E

6 6 = 666 4

7 6 7 6 7 2 u + p 77 ; 666 5 4

u(e + p)

4 @u 3 @x 4 @u 3 u @x

+ @T @x

7 7 7 7 7 5

(2.6)

2.2. THE NAVIER-STOKES AND EULER EQUATIONS

9

where is the uid density, u is the velocity, e is the total energy per unit volume, p is the pressure, T is the temperature, is the coecient of viscosity, and is the thermal conductivity. The total energy e includes internal energy per unit volume (where is the internal energy per unit mass) and kinetic energy per unit volume u2=2. These equations must be supplemented by relations between and and the uid state as well as an equation of state, such as the ideal gas law. Details can be found in Anderson, Tannehill, and Pletcher [1] and Hirsch [2]. Note that the convective

uxes lead to rst derivatives in space, while the viscous and heat conduction terms involve second derivatives. This form of the equations is called conservation-law or conservative form. Non-conservative forms can be obtained by expanding derivatives of products using the product rule or by introducing dierent dependent variables, such as u and p. Although non-conservative forms of the equations are analytically the same as the above form, they can lead to quite dierent numerical solutions in terms of shock strength and shock speed, for example. Thus the conservative form is appropriate for solving ows with features such as shock waves. Many ows of engineering interest are steady (time-invariant), or at least may be treated as such. For such ows, we are often interested in the steady-state solution of the Navier-Stokes equations, with no interest in the transient portion of the solution. The steady solution to the one-dimensional Navier-Stokes equations must satisfy @E = 0 (2.7) @x If we neglect viscosity and heat conduction, the Euler equations are obtained. In two-dimensional Cartesian coordinates, these can be written as @Q + @E + @F = 0 (2.8) @t @x @y with 2 3 2 3 2 2 q1 u 3 v 3 6 7 6 7 6 u2 + p 7 6 uv 7 7 7 7 Q = 664 qq2 775 = 664 u 7; E = 6 6 7; F = 6 6 2 7 (2.9) 5 4 5 4 5 v uv v + p 3 q4 e u(e + p) v ( e + p) where u and v are the Cartesian velocity components. Later on we will make use of the following form of the Euler equations as well: @Q + A @Q + B @Q = 0 (2.10) @t @x @y @E and B = @F are known as the ux Jacobians. The ux vectors The matrices A = @Q @Q given above are written in terms of the primitive variables, , u, v, and p. In order

10

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

to derive the ux Jacobian matrices, we must rst write the ux vectors E and F in terms of the conservative variables, q1 , q2, q3 , and q4 , as follows:

E

F

2 6 6 6 6 6 = 666 6 6 6 4

E1

2 3 6 6 7 6 7 6 7 6 7 6 7 6 7 7=6 6 7 6 7 6 7 7 6 5 6 6 4

3 7 7 7 2 2 q2 ;1 q3 7 3 ;

( ; 1)q4 + 2 q1 ; 2 q1 77 7 7 7 q3 q2 7 7 q1 7 7 7 q3 q2 q 5

; 1 2 q q

4q12 ; 2 q212 + 3q12

(2.11)

2 6 6 6 6 6 = 666 6 6 6 4

3 2 F1 7 66 7 6 7 6 7 F2 7 66 7 7 = 666 7 F3 77 66 7 6 5 6 4 F4

3 7 7 7 7 q3 q2 7 q1 7 7 7 7 2 2 ( ; 1)q4 + 3;2 qq31 ; ;2 1 qq21 77 7 7 5 q2 q 3 q q q

; 1 3

4q13 ; 2 2q12 + q132

(2.12)

E2 E3 E4

q2

q3

We have assumed that the pressure satis es p = ( ; 1)[e ; (u2 + v2)=2] from the ideal gas law, where is the ratio of speci c heats, cp=cv . From this it follows that the ux Jacobian of E can be written in terms of the conservative variables as

A=

2 6 6 6 6 6 @Ei = 666 @qj 66 6 6 6 4

3 7 7 q q 7 a21 (3 ; ) q21 (1 ; ) q31 ; 1 777 7 7 q q q q 7 2 3 3 2 ; q1 q1 0 77 q1 q1 7 q 7 5 a41 a42 a43

q21

0

1

0

where

a21

= ;2 1 qq3 1

!2

; 3 ;2 qq2 1

!2

0

(2.13)

2.2. THE NAVIER-STOKES AND EULER EQUATIONS 2 ( ; 1)4

+ qq3 1

!2

a42 = qq4 ; ;2 1 43 qq2 1 1 ! ! a43 = ;( ; 1) qq2 qq3 1 1 and in terms of the primitive variables as

!2

a41 =

q2 q1

!3

2

!

2 6 6 6 6 6 A = 666 6 6 6 4

where

0

1

q2 q1

!3 5;

+ qq3 1

!2 3 5

11

q4 q1

!

q2 q1

!

(2.14)

0

0

a21 (3 ; )u (1 ; )v ( ; 1)

;uv

v

u

0

a41

a42

a43

u

3 7 7 7 7 7 7 7 7 7 7 7 5

(2.15)

a21 = ;2 1 v2 ; 3 ;2 u2

a41 = ( ; 1)u(u2 + v2) ; ue

a42 = e ; ;2 1 (3u2 + v2 )

a43 = (1 ; )uv (2.16) Derivation of the two forms of B = @F=@Q is similar. The eigenvalues of the ux Jacobian matrices are purely real. This is the de ning feature of hyperbolic systems of PDE's, which are further discussed in Section 2.5. The homogeneous property of the Euler equations is discussed in Appendix C. The Navier-Stokes equations include both convective and diusive uxes. This motivates the choice of our two scalar model equations associated with the physics of convection and diusion. Furthermore, aspects of convective phenomena associated with coupled systems of equations such as the Euler equations are important in developing numerical methods and boundary conditions. Thus we also study linear hyperbolic systems of PDE's.

12

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

2.3 The Linear Convection Equation 2.3.1 Dierential Form

The simplest linear model for convection and wave propagation is the linear convection equation given by the following PDE: @u + a @u = 0 (2.17) @t @x Here u(x; t) is a scalar quantity propagating with speed a, a real constant which may be positive or negative. The manner in which the boundary conditions are speci ed separates the following two phenomena for which this equation is a model: (1) In one type, the scalar quantity u is given on one boundary, corresponding to a wave entering the domain through this \in ow" boundary. No boundary condition is speci ed at the opposite side, the \out ow" boundary. This is consistent in terms of the well-posedness of a 1st-order PDE. Hence the wave leaves the domain through the out ow boundary without distortion or re ection. This type of phenomenon is referred to, simply, as the convection problem. It represents most of the \usual" situations encountered in convecting systems. Note that the left-hand boundary is the in ow boundary when a is positive, while the right-hand boundary is the in ow boundary when a is negative. (2) In the other type, the ow being simulated is periodic. At any given time, what enters on one side of the domain must be the same as that which is leaving on the other. This is referred to as the biconvection problem. It is the simplest to study and serves to illustrate many of the basic properties of numerical methods applied to problems involving convection, without special consideration of boundaries. Hence, we pay a great deal of attention to it in the initial chapters. Now let us consider a situation in which the initial condition is given by u(x; 0) = u0(x), and the domain is in nite. It is easy to show by substitution that the exact solution to the linear convection equation is then u(x; t) = u0(x ; at) (2.18) The initial waveform propagates unaltered with speed jaj to the right if a is positive and to the left if a is negative. With periodic boundary conditions, the waveform travels through one boundary and reappears at the other boundary, eventually returning to its initial position. In this case, the process continues forever without any

2.3. THE LINEAR CONVECTION EQUATION

13

change in the shape of the solution. Preserving the shape of the initial condition u0(x) can be a dicult challenge for a numerical method.

2.3.2 Solution in Wave Space

We now examine the biconvection problem in more detail. Let the domain be given by 0 x 2. We restrict our attention to initial conditions in the form u(x; 0) = f (0)eix (2.19) where f (0) is a complex constant, and is the wavenumber. In order to satisfy the periodic boundary conditions, must be an integer. It is a measure of the number of wavelengths within the domain. With such an initial condition, the solution can be written as u(x; t) = f (t)eix (2.20) where the time dependence is contained in the complex function f (t). Substituting this solution into the linear convection equation, Eq. 2.17, we nd that f (t) satis es the following ordinary dierential equation (ODE) df = ;iaf (2.21) dt which has the solution f (t) = f (0)e;iat (2.22) Substituting f (t) into Eq. 2.20 gives the following solution u(x; t) = f (0)ei(x;at) = f (0)ei(x;!t) (2.23) where the frequency, !, the wavenumber, , and the phase speed, a, are related by ! = a (2.24) The relation between the frequency and the wavenumber is known as the dispersion relation. The linear relation given by Eq. 2.24 is characteristic of wave propagation in a nondispersive medium. This means that the phase speed is the same for all wavenumbers. As we shall see later, most numerical methods introduce some dispersion; that is, in a simulation, waves with dierent wavenumbers travel at dierent speeds. An arbitrary initial waveform can be produced by summing initial conditions of the form of Eq. 2.19. For M modes, one obtains

u(x; 0) =

M X m=1

fm (0)eimx

(2.25)

14

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

where the wavenumbers are often ordered such that 1 2 M . Since the wave equation is linear, the solution is obtained by summing solutions of the form of Eq. 2.23, giving

u(x; t) =

M X

m=1

fm (0)eim(x;at)

(2.26)

Dispersion and dissipation resulting from a numerical approximation will cause the shape of the solution to change from that of the original waveform.

2.4 The Diusion Equation 2.4.1 Dierential Form

Diusive uxes are associated with molecular motion in a continuum uid. A simple linear model equation for a diusive process is

@u = @ 2 u (2.27) @t @x2 where is a positive real constant. For example, with u representing the temperature, this parabolic PDE governs the diusion of heat in one dimension. Boundary conditions can be periodic, Dirichlet (speci ed u), Neumann (speci ed @u=@x), or mixed Dirichlet/Neumann. In contrast to the linear convection equation, the diusion equation has a nontrivial steady-state solution, which is one that satis es the governing PDE with the partial derivative in time equal to zero. In the case of Eq. 2.27, the steady-state solution must satisfy @2u = 0 (2.28) @x2

Therefore, u must vary linearly with x at steady state such that the boundary conditions are satis ed. Other steady-state solutions are obtained if a source term g(x) is added to Eq. 2.27, as follows: "

#

@u = @ 2 u ; g(x) @t @x2

(2.29)

giving a steady state-solution which satis es

@ 2 u ; g(x) = 0 @x2

(2.30)

2.4. THE DIFFUSION EQUATION

15

In two dimensions, the diusion equation becomes " # @u = @ 2 u + @ 2 u ; g(x; y) (2.31) @t @x2 @y2 where g(x; y) is again a source term. The corresponding steady equation is @ 2 u + @ 2 u ; g(x; y) = 0 (2.32) @x2 @y2 While Eq. 2.31 is parabolic, Eq. 2.32 is elliptic. The latter is known as the Poisson equation for nonzero g, and as Laplace's equation for zero g.

2.4.2 Solution in Wave Space

We now consider a series solution to Eq. 2.27. Let the domain be given by 0 x with boundary conditions u(0) = ua, u() = ub. It is clear that the steady-state solution is given by a linear function which satis es the boundary conditions, i.e., h(x) = ua + (ub ; ua)x=. Let the initial condition be

u(x; 0) =

M X m=1

fm(0) sin m x + h(x)

(2.33)

where must be an integer in order to satisfy the boundary conditions. A solution of the form M X u(x; t) = fm(t) sin mx + h(x) (2.34) m=1

satis es the initial and boundary conditions. Substituting this form into Eq. 2.27 gives the following ODE for fm : dfm = ;2 f (2.35) m m dt and we nd fm(t) = fm(0)e;2m t (2.36) Substituting fm(t) into equation 2.34, we obtain

u(x; t) =

M X m=1

fm (0)e;2mt sin m x + h(x)

(2.37)

The steady-state solution (t ! 1) is simply h(x). Eq. 2.37 shows that high wavenumber components (large m) of the solution decay more rapidly than low wavenumber components, consistent with the physics of diusion.

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

16

2.5 Linear Hyperbolic Systems

The Euler equations, Eq. 2.8, form a hyperbolic system of partial dierential equations. Other systems of equations governing convection and wave propagation phenomena, such as the Maxwell equations describing the propagation of electromagnetic waves, are also of hyperbolic type. Many aspects of numerical methods for such systems can be understood by studying a one-dimensional constant-coecient linear system of the form @u + A @u = 0 (2.38) @t @x where u = u(x; t) is a vector of length m and A is a real m m matrix. For conservation laws, this equation can also be written in the form @u + @f = 0 (2.39) @t @x @f is the ux Jacobian matrix. The entries in the where f is the ux vector and A = @u

ux Jacobian are @fi aij = @u (2.40) j The ux Jacobian for the Euler equations is derived in Section 2.2. Such a system is hyperbolic if A is diagonalizable with real eigenvalues.1 Thus = X ;1AX (2.41) where is a diagonal matrix containing the eigenvalues of A, and X is the matrix of right eigenvectors. Premultiplying Eq. 2.38 by X ;1 , postmultiplying A by the product XX ;1, and noting that X and X ;1 are constants, we obtain

@X ;1u

z }| { ; 1 @ X AX X ;1u

=0 (2.42) @t + @x With w = X ;1u, this can be rewritten as @w + @w = 0 (2.43) @t @x When written in this manner, the equations have been decoupled into m scalar equations of the form @wi + @wi = 0 (2.44) i @t @x 1

See Appendix A for a brief review of some basic relations and de nitions from linear algebra.

2.6. PROBLEMS

17

The elements of w are known as characteristic variables. Each characteristic variable satis es the linear convection equation with the speed given by the corresponding eigenvalue of A. Based on the above, we see that a hyperbolic system in the form of Eq. 2.38 has a solution given by the superposition of waves which can travel in either the positive or negative directions and at varying speeds. While the scalar linear convection equation is clearly an excellent model equation for hyperbolic systems, we must ensure that our numerical methods are appropriate for wave speeds of arbitrary sign and possibly widely varying magnitudes. The one-dimensional Euler equations can also be diagonalized, leading to three equations in the form of the linear convection equation, although they remain nonlinear, of course. The eigenvalues of the ux Jacobian matrix, orqwave speeds, are u; u + c, and u ; c, where u is the local uid velocity, and c = p= is the local speed of sound. The speed u is associated with convection of the uid, while u + c and u ; c are associated with sound waves. Therefore, in a supersonic ow, where juj > c, all of the wave speeds have the same sign. In a subsonic ow, where juj < c, wave speeds of both positive and negative sign are present, corresponding to the fact that sound waves can travel upstream in a subsonic ow. The signs of the eigenvalues of the matrix A are also important in determining suitable boundary conditions. The characteristic variables each satisfy the linear convection equation with the wave speed given by the corresponding eigenvalue. Therefore, the boundary conditions can be speci ed accordingly. That is, characteristic variables associated with positive eigenvalues can be speci ed at the left boundary, which corresponds to in ow for these variables. Characteristic variables associated with negative eigenvalues can be speci ed at the right boundary, which is the in ow boundary for these variables. While other boundary condition treatments are possible, they must be consistent with this approach.

2.6 Problems 1. Show that the 1-D Euler equations can be written in terms of the primitive variables R = [; u; p]T as follows: @R + M @R = 0 @t @x where 2 3 u 0 M = 64 0 u ;1 75 0 p u

18

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS Assume an ideal gas, p = ( ; 1)(e ; u2=2). 2. Find the eigenvalues and eigenvectors of the matrix M derived in question 1. 3. Derive the ux Jacobian matrix A = @E=@Q for the 1-D Euler equations resulting from the conservative variable formulation (Eq. 2.5). Find its eigenvalues and compare with those obtained in question 2. 4. Show that the two matrices M and A derived in questions 1 and 3, respectively, are related by a similarity transform. (Hint: make use of the matrix S = @Q=@R.) 5. Write the 2-D diusion equation, Eq. 2.31, in the form of Eq. 2.2. 6. Given the initial condition u(x; 0) = sin x de ned on 0 x 2, write it in the form of Eq. 2.25, that is, nd the necessary values of fm (0). (Hint: use M = 2 with 1 = 1 and 2 = ;1.) Next consider the same initial condition de ned only at x = 2j=4, j = 0; 1; 2; 3. Find the values of fm(0) required to reproduce the initial condition at these discrete points using M = 4 with m = m ; 1. 7. Plot the rst three basis functions used in constructing the exact solution to the diusion equation in Section 2.4.2. Next consider a solution with boundary conditions ua = ub = 0, and initial conditions from Eq. 2.33 with fm(0) = 1 for 1 m 3, fm (0) = 0 for m > 3. Plot the initial condition on the domain 0 x . Plot the solution at t = 1 with = 1. 8. Write the classical wave equation @ 2 u=@t2 = c2 @ 2 u=@x2 as a rst-order system, i.e., in the form @U + A @U = 0 @t @x where U = [@u=@x; @u=@t]T . Find the eigenvalues and eigenvectors of A. 9. The Cauchy-Riemann equations are formed from the coupling of the steady compressible continuity (conservation of mass) equation @u + @v = 0 @x @y and the vorticity de nition @v + @u = 0 ! = ; @x @y

2.6. PROBLEMS

19

where ! = 0 for irrotational ow. For isentropic and homenthalpic ow, the system is closed by the relation

;1 1

; 1 2 2 = 1; 2 u +v ;1 Note that the variables have been nondimensionalized. Combining the two PDE's, we have @f (q) + @g(q) = 0 @x @y where u ; u ; v q = v ; f = v ; g = ;u One approach to solving these equations is to add a time-dependent term and nd the steady solution of the following equation: @q + @f + @g = 0 @t @x @y (a) Find the ux Jacobians of f and g with respect to q. (b) Determine the eigenvalues of the ux Jacobians. (c) Determine the conditions (in terms of and u) under which the system is hyperbolic, i.e., has real eigenvalues. (d) Are the above uxes homogeneous? (See Appendix C.)

20

CHAPTER 2. CONSERVATION LAWS AND THE MODEL EQUATIONS

!"$#&%%'#&)(+*-,/. ,/.1003254768,9*:#&1;=&0$@91*-10c,/.10f@90$;941OT,/*:F6l,/*->03;cn1@9;R,3Jh,/.F03@90$YoEQ1@R#bD14F"3*:F&0$D\tu.1*-;+*:;u,/.10=6Q1QF@9#o6".03%'Q1.76;R*:$03D.F03@908\!,/.1*-;f".768QF,/0$@3J7,/.10w"3#1"30$QF,Z#qmn71*T,/0Wr D1*Tgh0$@903F"3068Q1Q1@9#Cab*:%68,9*:#&1;Z,9#Q768@R,/*:6OjD10$@9*->868,9*->03;M*-;KQ1@90$;903o,903D\tu.103;R0"C6HYs06Q1Q1O-*:0$D 03*T,/.10$@f,9#d;9Q168,/*:6OmDF03@9*T>868,/*T>&03;f#&@+,/*:%'0wD103@R*->868,/*T>&0$;3\c]4F@c03%'Q1.768;9*:;f*:,/.1*-;_.76Qb,/03@c*:;#& ;9Q168,/*:6OhD10$@9*T>6l,/*->03;$F,/*-%0fD103@R*->868,/*T>&0$;u6@90Z,9@90368,/0$De*:_.76Qb,/03@+b\PN,/@968,/0$&0#&@'O:#&;R*:1< Q1@96"$,9*:"C68O6Q1Q1O-*:"368,/*-#&\ ®?¯

°

±²^³±²µ´M¶A·

¸º¹C¶A¹$»L±f¼h½¾¹$¿±Àj±¶AÁ±

ÂÄÃ`»L´=»^¹3Ã'¶

tu.10;R*:%'Q1O:0$;R,%'03;R.*:o>&#&OT>N*:1868,/*T>&0*:e,9.10K;R,/4FDFE#q^]Kã ;u(_0K41;R0 Ýbä

Ý

æ¥F\ç Þ

!,/.F*:;+,/0WaN,;R41Y1;R"3@9*-QF,/;#D103QV03FD103o,>86@9*:6Y1O-03;f6@90=10$>03@41;R03D,/#d0WabQ1@90$;9;KD10$@9*->868,9*->03;3\ tu.541;Ý (+*-O:OhF#,Ys041;R03D,/#`@90$Q1@90$;903o,f,/.10n7@R;R,D103@R*->868,9*->&0#8qLÝ(+*-,/.@90$;9QV03"W,,9#'ß\ tu.10cF#,68,9*:#&[q#&@XD1*Tgh0$@90$1"30c6Q1Q1@R#?ab*-%d68,9*:#&F;^q#&O:O-#I(+;^,9.10;/6%'0fQF.1*:O-#&;9#&QF.5EUY14F,Zæ(+*-,9. #&10P0Wab"$03QF,9*:#&ëç*-,§*-;§F#,L4F1*:254108 \ vE,9.1*:;(v0_%'0C6K,/.76l,m,9.10v;BEN%UYs# O u*-;§4F;903D=,/#@903QF@903;R03o, 6[D1*-gs03@R031"$0M68Q1Q1@9#Cab*:%68,9*:#&d,9#'6`DF03@9*T>868,/*T>&0K;941".,/.76l,CJFq#&@u0aF6%QFO:0J

! à

æ¥F\ &ç

Y14F,u,9.10KQ1@90$"3*-;90K768,941@90'æ61D#&@9D10$@ç_#q§,/.F0K6Q1Q1@R#?ab*-%d68,9*:#&d*:;u1#,+"36@9@R*:0$D*-,/.10K;BEb%MYs#&O 5\U]c,/.F03@Z(_6?EN;w6@R0[41;90$Dk,9#D10$,903@9%'*:F0U*-,9;ZQ1@90$"3*-;90[%'0C61*-1N*:DF0Y5Eeî[ß,/#'#&YF,/6*:

Ý ÝbíRü#"&$ÝNí äE8 ø

ßF9 í î[ß ¦

Ý

æîUßsç8

ß 9 í

ø868,9*->03;3Û

8 8

Ý

ß 9 í

Ý

ß 9 í ä

¦ æÝNíRü#"%$HÝbíP"Bçø~\æzîUß ç î ß U ä

¦ æÝbíRü#"&$ î ß U

ÝNí^ø¬ÝbíP"Bçø~\æzîUß ç

æz¥F\T¦?¥oç æz¥F\T¦C5ç

tu.10$;906@R0f,9.10cY768;9*:;Xq#&@j7ly|7xXy Z}W/}Wh}7}WRxÔl/X;R*:1"$0f,/.F0$E'0Z6d6Q1Q1@R#CaF*-%d6l,/*:#`,9# 6=DF03@9*T>868,/*T>&0f68,_#10cD1*-;9"$@90$,90cQV#&*:o,v*:6=%'03;R.*:',/03@R%;P#q;R41@9@R#&411DF*:1&03@$J 10$*-,/.F03@w#q_,9.103;R0d0WabQ1@R03;R;9*:#1;=,903O:Oj4F;=.1#?( #8,/.10$@=QV#&*:o,9;w*-),9.10'%0$;9.68@90'D1*-gs03@R031"$03DH#@ .1#?( YV#&411D768@RE"3#1D1*-,9*:#&F;c6@R0=0$Fq#&@R"303Dh\cN41".A6D1DF*-,/*-#&76O*-Fq#&@R%d68,9*:#&e@R032541*-@903;c6%'#&@90 ;9#Q1.1*:;B,/*-"C68,903Deq#&@9%M41OS68,9*:#&h\

`ab`ao

|%jrgolBgn{pRrpRiutvpxwmypRr|%je7r}

_#&1;R*:D10$@u,/.10K@90$OS68,9*:#& æ Ýsçzíä

¦ æÝNí9ü#"&$ î ß U

ÝbíLøÆÝNí)P"Bç

æz¥F\T¦G&ç

(+.1*-".e*:;u6MQV#&*-5,_D1*-gs03@R031"$06Q1QF@9#Cab*:%68,/*-#&,9#[6U;90$"3#&1DD10$@9*T>6l,/*->0\jò#?(ÿO-0$,_41;uD103@R*->06 {'xWy w#&QV03@/6l,/#&@+@90$Q1@90$;90$5,/68,/*-#&q#@+,9.10K;/6%'0=6Q1Q1@R#?ab*-%d68,9*:#&h\P_#&1;9*-D103@+,/.F0wq#&4F@QV#&*-5, %'03;9.(+*T,/.YV#&41FD76@REQV#&*:o,/;Z68 , 61D6K;9.1#?(+kYV03O:#?(K\ò#,/*-"30M,9.768,f(+.103A(v0[;9QV0C6÷#q ,9.10=541%MYs0$@#8q^QV#&*:o,/;+*-6`%0$;9 . FJë(v0w%036,9.10wN4F%UYs0$@f#8qy|7x!}WWyl_Qs#&*-o,/;+0Wab"3O-41D1*-1< ,/.F0wYV#&41FD76@9*-03;$\

ßeä ê[ä

$ ¦

¦

$ ;

$

¥

;

$

¤ë#41@uQs#*:o,u%'03;9.h\vîUßäRH1æ

øª¦Iç

ò #?( *-%QV#&;R0ãZ*:@R*:".1O-0$,=Ys#411D76@BEÅ"3#1D1*-,9*:#&F;3JjÝmæ ç[ä ÝIà[Ýmæç'ä ÝK681D¬41;90e,/.F0 f "30$o,/03@R03DAD1*Tgh0$@90$1"30=6Q1Q1@R#CaF*-%d6l,/*:#e03YoEj2ë\s¥b\-¦G`68,0$>03@BEQV#&*:o,*:,/.F0=%'03;R.\_ 0

ov@@¡ =£¢¤@¥¦%§¨G©Iª£«¦¬[ =I¢@¥§@¥¦«.®ZI¢¤)§«¦¢©¨G«¦¢§@¯¦°

Ç^ÈÉZÊcËmÌXÍÏÎ5Ð«ÑLÒÔÓÒ3ËmÌPÕzÖZÒÔÑ^Ñ^ÌXÍÌPÓMÇ^Ì«ÉZÊXÊXÍ×jØZÒBÙÉcË§ÒW×uÓZÚ

6@R@9*->0K68,+,/.10q#&4F@u03254768,9*:#&1;$Û æ Ýsç " ä

æ Ýsç

ä

æ Ýsç æ Ýsç

¦ î ß U ¦ î ß U ¦ î ß U ¦ î ß U ä

T ä

U

æÝk$

Ýv"ø¬Ý

ç

æÝ"&$ Ý

øÆÝ T ç

æÝ

$ Ý

T øÆÝ U ç

æÝ

T $ Ý

øÆÝ!ç

U

æz¥F\T¦?oç

@9*-,9*:186@9*-#&41;fQs#&*-o,/;c*-,/.F0Mn703O-D*:F"3O:4FD1*:10c68,v,/.10f@90Wq03@90$1"30QV#&*-5,LêU*:,/03@R%;P#q10$*:&0ÛR41;R,uDF03@9*T>&0$D,/.10q6%*-O:*:6@j¥lr£Qs#&*-o,_"$03o,/@96OhD1*Tgh0$@90$1"3*-1N*:D10$;,/.F0 O:036D1*-1&0D10$@9*->03D©6);R03"$#&1Dbr£#&@9D10$@`Y76"÷o(68@9D«DF*-gs03@90$1"306QFQ1@9#Cab*:%68,/*-#&#qc6n7@9;B, D10$@9*->868,9*->0Û

8 `aßõao

Ý

ß 9 í

¦ æÝNí)P $Å&ÝNíP"hø©¥ÝNíWçjäB\æîUß ç î ß U

$

æz¥F\ ç

×p&iupsr|ã£gv|%jgekiæe,lBgn{pRrpRiutvp ue7rÇâã|}

!)868,/*T>&0\ £,m"3#&41O-DMYV0_"$#&%QF4F,/0$@L64b,/#&%68,/0$DM681DM0ab,9031DF03D[,/#f.1*:&0$@3Jv*:;U,/.F0q6"$,U,9.768,`*T,U"C68ÆYV0q41@B,/.10$@[68O:410$;P#8qh,/.F0 q411"W,/*-#&AVUy|x|wN}Wy lxzy "} !z$ #u68,v03dQV#&*-5,9;_*-;RQ76"308\v]4F@v*-O:O:4F;R,/@968,/*-#&[*:;Xq#&@X,/.10f"C6;R0 *:k(+.1*-".D1*:;R"3@R0$,/0`>86O:4F03;w#q,/.F0'q41F"$,/*-#&681DH*T,/;Kn7@9;B,MD10$@9*T>6l,/*->0d68@90'41;90$DJ^Q1@R#bD14F"3*:F< ,/.F0w0abQ1@90$;9;9*-#& ÝLæ|ßsçjä

oýoæßsçÔÝëýPø

ýoæßsç 8

Ý

ß 9 ý

æz¥F\ &5ç

] Yo>b*-#&41;RO-EU.F*:868,9*->03;^"$#&41O:D`YV0u*:1"$O:41DF03D`6;^,9.10uQ1@9#Y1O:0$%;^DF*:"$,/68,/08\Li "$#&%`r Z Q1O-0$,/0MD1*:;R"341;R;9*:#A#qj,/.10$;90[QV#&O-EN1#%*:6O:;+"C68AYs0Mq#411DA*:%d6oE@90Wq03@90$1"30$;K#&541%'03@9*-"C6O %'0$,/.F#bD1;$J7Y14F,+.10$@90K(_0K10$03D#&1O-E,/.10K"$#&1"30$QF,C\ tu.10MQF@90$>N*-#&41;c0aF6%QFO:03;#q^6tm6?ENO-#&@u,6Y1O-0="$#&1;B,/@94F"$,/0$D©I} $V éyWy|xQs#*:o,D1*-gs03@R031"$0M#&Qbr 03@968,/#@9;q@9#&%pô6868,/*T>&03;+6l, QV#&*:o,/;^u ê $Æ¦K681D`ê+ø«¦&Jb(+.F*:".6O:;R#U%U41;B,uYV00WabQ761DF03D41;9*-1&0,903@R%;[1#?(.76C>&0"3#N0$¨'"3*-03o,/;æ,9.10"3#N0$¨'"3*-03o,'#&Æ,9.10'ê QV#&*:o,`*-;U,6÷803 6;`#&10,/#);9*-%QFO:*-q|Ek,9.1068O:N*:186O:0$5,u,9#

P" ² ð Ý ä 7ù¸ÁæB¦&à/FàC¦IçÔû

¦ ² Ç Áeæ£$w¦&à à3¦Iç ÝM ø J )² º O î ß U

æz¥F\@ ç

(+.1*-"."36YV0K@90$0WabQ1@90$;9;R03DYoEe,/.F0 Q1@90$D1*:"W,/#&@Ôr!"3#@9@90$"$,9#&@Q`;90$2541031"$0 Ý ²9

Ý²

ä ä

¦ ² Ç Áæ£$w¦&à àC¦?ç ÝM ø J )² º O î ß U P" fýù ÁeæÔ¦&à/1àC¦?çÔû Ý ²9

æz¥F\@b¦Iç

*-,9.Å@903;RQs0$"$,`,/#Q1@/6"W,/*-"C6Ov*:%'Q1O:0$%0$5,/68,/*-#&J§,/.10%'0C681*:1&0$@9;90w%68,/@R*Ta #&QV03@968,/#&@w*:%'Q1O-*:03;Z,9.10d;R#&O:4b,/*:# #qu6"3#&4FQ1O:0$DH;90W,[#quO-*:1036@K03254768,9*:#&F;3\tu.103;R0#&QV03@/6l,/#&@R; 6@R0+>03@BE"$#&%'%#&`*-`n11*-,90DF*-gs03@90$1"30f6Q1QFO:*:"368,/*-#&1;$\mtu.10WE6Q1QV0C6@X*:`,/.F0+q#@9%Ï#qhY761DF03D %68,/@R*:"30$;.16?>N*-1868,/*T>&0"36¬6O-;9#YV00C68;9*:OTE)D103@R*->&0$DÅYoE%'0C61;M#q6

tm6?ENO-#&@,6YFO:0\X¤7#&@u0aF6%QFO:0

P"

¦ ² ø J )² º O ÁeæÔ¦&àN$ àC¦Iç ÝU æz¥F\@ ç î ß U *:;h\æîUß U çP681Dd%6÷0$;_41;R0c#qh#&1OTE[,/@9*-D1*S68868,/*T>&03;6@R0KYs0$*:186O-4768,903D\X!6Q1Q1O-*:"C6l,/*:#JN,9.1*:;+"36 YV0c6Db>685,/6N*:"$*:1*T,!EU#qh;R,90303QdN*:#&4F;3J5#q "3#41@9;R0JF,9.768,Pn1>&0ZQV#&*-5,v;9".F03%'03;u41;R*:1&0c6K(+*-D103@X;RQ1@9036Dd#qh;9Q768"30c*-1D1*-"303;$\L!Q76@R,9*:"$41OS6@$J&,!(_#Mô6868,9*->03;Z#7 q þ a_ \ZH0M(+*-O:O "3#1"30$5,9@/68,90Jb.103@R0JF#n7@R;R,DF03@9*T>868,/*T>&0f6Q1Q1@R#?ab*-%d68,9*:#&F;3J56OT,/.1#&4F&0$N4F%UYs0$@3Ji`J86Q1QV0C6@R;m*-w,9.10v0aF6"$,L0WabQ1@90$;9;R*:#&\Ltu.N4F;m,/.F0vD10$&0$N4F%UYs0$@X*:;^6Z%'0C68;941@R0 #q§,/.10w6"$"341@96"$E#8qm,9.10w6Q1Q1@R#?ab*-%d68,9*:#&h\ ¤7#&@,/.10A;R03"$#&1Dbr£#&@9D10$@"$03o,/0$@90$D D1*Tgh0$@903F"30k#&QV03@/6l,/#&@d,9.10k%#ND1*Tn703Dº(_6C>&03541%MYs0$@*-; &03YoE

`Lhvä òf#8,/0,/.76l,j`

h

;9*-f`ëîUß î[ß

æz¥F\@L?oç

6Q1QF@9#Cab*:%68,/0$;f`,/#`;90$"3#&FDbr!#&@RD103@f6"$"341@96"$EJë6;+*:;,/#'YV00WabQs0$"$,903DJV;9*-1"30 ;R*:f`7îUß îUß

äk`Ó$

`

T î[ß

ø

W5WNW

j254768,/*-#&ÿ¥F\ Z?¬*-;Q1O:#,R,/0$Dð*- ¤m*: ¬:§¤£¢¤Ô©#«¦Ñkª@¢¤ª£ÒG@¦¬5§#Ô§¤¢r qF«¦®L£¢b§¤«¦¢nsBtN§¨G¦¬.§¤@©¯NÔÉÔk£§¢¤@ªR®[)¢§@©Ú«/AG§¤)@¬GQ Ñ¢¤«¦Ôcuvsxw-s+yCz )¬Z F§¤¨G©¯NÔÉÔÉ§¤¢ª®[)¢§Ñ¢«.Ô|uvs3}~s y z EU{ °

E/{

ÎNÐ[Z5Ð«ÑP×ÓèmÍ+ÒÔÌXÍóÌXÍfÍZ×+ÍóÉcÓfÉ]\RêXÚNÒBÚ

1¦

ª.10$w,9.10^n71*T,/0Wr£D1*-gs03@R031"$0j#&QV03@/6l,/#&@*-;6o,/*-;REN%'%0W,/@9*-"uæ"303o,903@90$DëçWJl,/.10j%'#ND1*-n70$Dc(_6?>03541%[r YV03@_*-;vQF41@90$O-E'@9036O\§ª.103,9.10c#&QV03@968,/#@v*-1"3O-41D103;_6=;BEb%'%'0$,/@R*:""3#&%'QV#&103o,3JN,9.10c%'#bDF*-n70$D (_6?>03541%UYV03@*:;`"3#%Q1O-0WasJj(+*-,9.Æ,9.10*:%6&0$"$,9*:#&0$254768,/*-#&J,/.1003@R@9#&@R;="36HYV0&0$ 6dQ1.oEN;9*-"C6O§*:o,903@9QF@90$,/68,/*-#&\c_#&1;R*:D10$@c#&1"$0[6&03o,9476O:OTE)O:#&;R03;K*-,9; #&@R*:&0e*-;=@R03;R#&O->03D*-; &03ÅYoE RHo`7îUß\ tu.10e@R03;R#&O->N*:F&0$541%UYV03@ O:036D1;',/#Å6 0$@9@R#&@d*-«,/.10Q1@9#&Q160,9.N4F;Z,/.10[O-0$q|,Rr .76FD)YV#&41FD76@RE*:;c,9.10`*:FG7#?(åYs#411D76@BE&J6FDk,9.10`@9*-03"$,9#&@#qL4F1÷51#?(+1;f*-; Ý ² äÏù Ý"à+Ý e

à W5WNW à+ÝÅ=û Æ

æz¥F\

61DÝä+*:;u;9QV03"$*-n70$D\ _#&1;R*:D10$@^n7@9;B,^,/.10_*:FG7#?(ÆYV#&411D768@RE&\m£,j*:;m"3O-0C6@L,9.768,^6;^O:#103@3Jc*-q(v0)4F;90,/.F0q#&41@R,9.br!#@9D10$@*:o,/0$@9*-#&@d#&QV03@968,/#@03 *-CX2ë\u¥b\8FJ+,9.103ª,/.F0 6Q1QF@9#Cab*:%68,/*-#&e68,PêMäå¦K@9032541*-@90$;f6U>86O-410Z#qmÝNí)P J1(+.F*:".*:;u#&4b,/;9*-D10Z,/.F0wDF#&%d68*:\031"$0 6D1*-gs03@R03o,w#&Qs0$@/68,9#&@K*:;K@90$2541*:@R03D6l,cêä ¦`(+.1*:". 0WaN,/0$1D1;=#&1O-E,/#êëÿ $ ¦&J§(+.F*:O:0[.16?>N*-1< ,/.F0Z6Q1QF@9#&Q1@R*S68,90f#@9D10$@_#q68"3"34F@/6"WE&\vb41".6#&QV03@/6l,/#&@$Jb÷51#?(+6;_6d7¢b{d}Wy9 #WÀ l¢NVoW

?7}W{d}WJP"C6©.76C>&068«#&@RD103@`#qZ6"3"$41@/6"WE¬(+.1*:".©*:;lh}AlRN}We :lè_}WM,/.768©,/.168,'#qZ,/.F0 *:o,903@9*-#&@c;R".103%'0J681DA,/.10U&0$Dd41;R*:1&0$"$,/*-#&k03254768,9*:#&ADF*:;9"$41;9;R03D)6YV#?>&0\ª*T,/.,/.F0 ;90$"3#&FDbr!#&@RD103@"30$o,/03@R03D*:o,/0$@9*:#@#&QV03@/6l,/#&@ æ s Ý çí ä

¦ æÝNí9ü#"&$ î ß U

ÝbíLøÆÝNí)P"Bç

æz¥F\@'¥oç

1#f%#ND1*-n1"C68,9*:#&1;6@90P@90$2N4F*:@90$DU1036@mYV#&41FD76@9*-03;$J8O-0C6D1*-1N*:#&41;RO-E (+.1*-"."36YV0K#&10K#&@RD103@uO-#I(v03@u,9.76,/.10K*:o,903@9*-#&@;R".103%'0\ H0`"36H68O:;9#e#YF,6*-A,/.10`#&QV03@968,/#&@*-6j2V\§¥F@\ '41;R*:1&0$"$,9#&@Zq#&@9%'03DY5E41;9*-1868,/*T>&0*:e,9.10q#&@9% æ ÝVçíuä

¦ æoÝNí)P ø~ÝNíP"ø~ºWÝNíLø î ß T U

¤L*-1De,/.10KO-0C6D1*-1868,9*->0w#Qs0$@/68,9#&@,/#[,/.F0Kq4F1"$,9*:#&þ [_ 03;

þ [_ ä $ ` þ [_

ß ifQ1Q1O:*-"C68,9*:#&d#q^6`D1*Tgh0$@903F"30K#&QV03@/6l,/#&@q#&@,9.10K;90$"3#&1DD103@R*->868,9*->&0&0$; æ þ

a_ í e

çíuäX$]` h þ

[_

,9.N4F;wDF0$n71*-1&03541%MYs0$@`q#@w,9.10d;R03"3#1Dbr!#@9D10$@U"30$5,903@R03D¬D1*-gs03@R031"$0 #&QV03@968,/#@`q#@6k;90$"3#&FDºD103@R*->868,9*->&08Jv,9.10F#&1"3#%Q768"$,'q#&41@R,9.br!#@9D10$@'#Qs0$@/68,9#&@æ j2ë\ ¥F@\ &çJ16FD,/.10"$#&%Q16"$,q#41@R,9.br!#&@RD103@_#&Qs0$@/68,9#&@_DF03@9*T>&0$D*-e254103;B,/*:#F\XÜjO:#f , ` h îUß >N;3S \ `7îUßq#&@ m l `7îUß l j\

DF\¤L*-1D©,9.1003O:0$1?!#@6A.2BDC"':E"1#=F"".0G!#86&%('H%(?)I)"$CJ$")"6K1(/&>? C-%(#L/&.:M!F"'2%(#NC">"$"1BD!"1O9J$.2"PM$1"6A.:)"!#8!)RQS?M!%(6A31(4K/A>"3"1"':.2"!%(#UTWVWXZY[6\1(4 .2]/A!#A$68/^.:P_"F".2)`)?;a-%5Bb.:M!6!c*6L/d%(9".:'2.0/">"%5!/v68W)"%(1(#ABq~%(.2.:BDRc-Cl6A11(#8Bb/&%(w]/"1%5")#ABxq.:15/4yA!.:#8/A!6m6868/A1.2"'2F?K/AM!.211U."M$2!6mC?%(/A'6!u cj%z9";aF?6W/ %5/8v/&#&)"%51wM$/&"1$)o(//A)1b/&!>?'2K'l16&%5q/BbA>"@!+?Bij!.:) /&>?.26 u 1#AsQHTvF"#y9-%(68.2M}M$1"M!$#AD.26 u .0/&>E/&.:Bb^/&>? M!1a$DM$.2!]/@Bb%5/.2M!$6=%(#A^.:")"!Cl!?)"!]/w1(49l1(/A>I%(")|6A$pa!M$/A.21"QjQE u .:'2'#8$4t$# /&1^68F"Md>Bb%5/.2M$!6Z%(6,(l(Q>-%(C?/A!#A6~D%(")o E)"$G(!'21(CJ$)U/A>"=#A!C?#A!68!]/d%/&.:G(w4t1#8Bb6 1(4TWVWX}Y6e(!"$#&%5/A!)`4t#A1B¡/A>"N9-%(68.2MN¢*VWX}Y6K9;/&>"N68!BD.2)".:6AM$#A$/AN%(C"C?#A1]%(Md>£c*%(")h/&>?! /&>?vTW^X}Y6!"$#&%5/A!)D4t#A1B¤/A>"m#8!C"#8!68!]/d%5/A.:G(WTWVWX}Y6H9];q%5C"C"'2.:Mz%5/A.21E154£/&.:Bb?.2" BD$/&>?1j)"6$Qy{}>"!68K$¥aF"%5/&.:1"6~%(#AW#8!C"#8!68!]/&$)9]; ¦y § b§ « ¬|§ 8® ¯ Q:°±® ¦ ©¨Oª %(") -§ ²z³J´ ¨Oµ -§ ²W« ¶ § ² ¯ Q·(® #A$6ACl!M,/&.:G(!'0;Q}¸-1#~%^1",¹68/A!CBD$/A>"1j)£c"/&>"@'º%5/8/&$#74t1#ABx.:6~1(9?/d%(.:"!)U9];U%(C"C"'0;a.2"b%^/A.2BD, Bb%(#AMd>".:"pBD$/A>"1j)/A1p/&>"E(!"$#A.2MET=VWX»4t1#AB¼.2%N4½%5.2#A'0;\68/A#&%(.:>]/A4t1# u %(#8)©Bq%5""!#$QK¸-1# ,i?%(BDC"':(c-/&>"@,ijC"':.2M$.:/mXyF"'2$#~BD$/A>"1a)'2!%()"6}/&1q¾ ¨À¿=ÁÃÂ-ª cl%(") ¶ § ² ¨ Ä§¬ ½Å Â ®Q~Æo$/A>"1a)"6 .2]G(1':Ga.2?Z/ u 1~1#ÇBb1#8y68/A!C"6ÈMz%(w%(' u %z;a6Ç9J u #A.0/A/A[email protected]/&>"y4t1#8BÀ1(4?Xy¥lQ ¯ Q}9;v.2]/Qj)"F?M!.2? " u )?!Cl!")"$]/WG(%5#A.º%59"'2$6!Q~Ém1(/Ae%(':6A1D/&>-%5/m1"'0;UBb,/&>"1a)"6.: u >".:Md>/&>?w/&.:Bbw%(")6AC-%(M$ )".:6AM!#8$/A.2Êz%/&.21("67%(#8m/A#Az%/&!)N68!C-%5#&%5/A!':;bMz%5N9J u #8.:/8/&!q.2N%(q.2]/&$#ABD!)".2%5/&m6A!BD.0¹)".268M!#8$/& 4t1#8BË68F"Md>P%(6EXy¥sQ ¯ Q:°QU{}>"q4tF?'2':;]¹)".268M!#A,/&^4t1#ABUcXy¥sQ ¯ Q·%5")I/&>?N%(6A681jM$.º%5/A!)h6L/d%(9?.2'2.0/];aCJ$#A9l1'2.:Me1"$6d®,cÇ%p6L/d%(9?.2'2.0//&>?e/&.2BDe%5")©6AC-%5M!^)".:øJ!#8!"M$.2"?Q@{}>?.26v.:6v)".268M!F"686A$).: j$M$/&.:1 ¯ Q?Qy-%('0;a6A.26Z1(4È/&>?!6A@!.:!?68;a68/A!BD6m>-%56~/A>"=.2BDCJ1#L/d%(]/Z%5)")"!)o%()?G5%(]/&%(W/&>-%/ .:/e.:G(!6D%(P!6L/&.2Bb%5/AN1(4m/A>"8(rùe%5/ u >".:Md>Ã%©6A1('2F?/A.21P%(C?C"#A1]%5Md>"!6D%68/&!%()?;]?.2Md>hM!%(6AbC"#&%(M,/&.:Mz%('H%(C"C?'2.2M!%5/&.:1"6WM!%( Bb%((W/&>?@C?#A1Cl!#L/&.2$6}1(4È/&>"='º%/A/&$#Z%(C"Clz%(#7/A1b)?1Bb.:-%5/A(Q S¹4 ª %(") µ %(#A=68/&%5/&.:1-%(#L;c u KM!%(Rcs.2N/&>?!1#L;U%5/~'2z%568/zc-$68/A.2Bb%5/&W/&>?!.2#}4tF"?)-%(BD!]/d%(' C"#81CJ$#8/A.2!6$Q¸-1#@,i?%(BDC"'25c|.:`j$M$/&.:1h"QÍ?Q q4t1F"")4t#A1Bè1F"#@BD1j)"$'7TWVWX}Y[6w4t1#@)".:4º 4tF"68.21o%(")oCJ$#A.21a)".:MWM!1]G$M$/A.21 u >"%5/M$1F"':)u 9lwi?Cl!M,/&$)4t1#Z/&>"@!.:!]G5%('2F?@68CJ$M$/A#AF"BD6 1(4C"#&%5M$/&.:Mz%('|C">];j68.2M!%('*C"#A19?'2!BD6WM!1(/&%(.2?.2"p/A>"!68bC">?!"1BD!"%?cÄ68!D¸|.:"QR"Q0°QE{}>"$6Aqi Cl!M$/&%5/&.:1"6%(#A~#A,4t!#8#A!)b/&1@Bq%(];E/A.2BD!6H.2^/&>"Z4t1'2':1 u .2"@%(-%5':;a6A.:6H154Û68/&%(9".2':.:/"W4t1':'21 u .2?"û ü #8-¸ z 1%(#v'Ç)?%.:ijøJ.:F"6!Q 6A.:1©)"1BD.2-%/&!)_"1 u 6W/A>"bý"¹!.2(!]G5%('2F"$6/&$")/&1\'2.:E%('21("q/&>?D?!]%5/A.:G( ü .:-¸ qB 1#H%((CJ.2$-%#A(.:1j#L;b)".:%MmijM$.26!1Q]G$M$/&.:1j¹)"1BD.2-%5/A!)D_-1 u 6/&>"Wý"¹!.2(!]G5%('2F"$6H/A!")p/A1K':.2m%('21?w/&>? S"q6L/d%(9".:'2.0//A>"qM$1BDC"'2iIýP%5") M!1(BbC"':,i C"'2%("$6!Q !"#$&%'($)+*-,/. 01%2. 3!4%657&84%693 3: % &4*&;

±Õ Ö × ãpä)A BÇê©ðçDC*çvêNïFEpïGIH©ò'Jzò$ïRõ íLKNMnæ ñPO8ô

°Ì1@

QSRUTSR&V WYXDZN[]\_^U`9Z3\_^acb

m÷ $#A u 36L/d%5/A/&>"6L/d%(")"%(#A)»68/&%(9".:'2.:/"wT=VWX}Y[6!Q QSRUTSRUT [YaciPjlkZ`9Znmo^pqZrbDs tusv`9ZriPs

¹S 4-%ZBb%5/.0iW>"%(6|%ZM!1(BbC"':$/AH$.2$"68;a6L/&!BUc%5'2'1(4?.:/A6È$.2$]G!M,/&1#86%(#8*'2.2?z%(#8':;v.:")"$CJ$")"!]/!c %(")\/&>?wBb%5/.0iqMz%(U9lw)".2%(1"%('2.:Ê!!)q9;U%^6A.:Bb.:'º%(#8.:/%DM!%(6A@.:/ 4t1':'21 u 6=%5/@1"M$D4t#81B¡Xy¥sQÄ jQ· ¯ cÇ4t1(#@i?%(BbC?'2(cÛ/A>-%5/=/&>"qTWVWX}Y6@%(#8b.:">"!#8!]/&'0;©68/&%(9"':D.04 %(")U1"'0;p.04 w ý"ÿ}y® x{z 4t1#Z%(':} ' | ¯ Q]® {}>".:6b6L/d%5/A!6b/A>-%5/!c4t1#D.2">"$#A$/q68/&%(9".2':.:/"ý $.2$G5%(':F"!6DBeF"6L/p':.2N1Rc1#D/A1 /&>?q'2,4f/K154rcÄ/&>?q.2Bb%(.:-%(#L;©%[email protected]/&>"bM!1BDC"':,iIý`C"'º%(?(Q\S"%5/w/A>".26@M!#8.:/A!#A.:13.26@6&%5/A.268+"!)I4t1#=/&>?qBD1j)"$'HTWVWXZY[6K#8!C"#8!68!]/&.:"9l1(/&>h)?.:øJF"6A.:1I%5") 9".:M!1]G$M$/A.21RQ}S¹/m6A>?1F"'2)o9J@!BDC">-%56A.2Ê$!) %(6~.:/~.26m%(U.2BDCJ1(#8/d%5/ZC?#&%(M,/&.2M!%('ÇM!1?6A.2)?!#&%/&.21( .2UM!1(G(!M,/&.21(j?e%(9l1±GKM$1")".0/&.:1o.:6mBD$/!QZ¸|.2"%('2'0; u K"1(/Aw/A>-%5/~4t1#vTWVWX}Y[6 uM!#8.:.:/A/A>Ã!#A.:M$11RQBbC?'2$/A\!.2(!"6L;j6L/&$Bb6D/&>"o!.2(!]G$M$/&1(#A6qC?'º%z;P"1I#81'2\.:g/&>?o.:">"!#8!]/p6L/d%(9?.2'2.0/"$.2#H6A1':F?/&.:1"6H.2p$.2$"6AC-%5M!(Q*¸-1#*/&>?.26 u v)"#A% u 1(q/&>"m#A$6AF"'0/&67.: j$M$/&.:1"6="QjQÍU%5")$6ACl!M$.º%(':':;1(IXy¥6!QÇ"Q0° p/&1U"Q0° N.2©/&>-%/w68!M,/&.21(RQ^S"1Fj/9l1F"") ':.2"!%(#A'0;=1#¥F-%()?#&%5/L .2M!%('2'0;?®y.04Û z®¨ ¬ zK1#7R´ z]®l¨ ¬ zjQy{}>"!1#8$/A.2Mz%5'2':;E/A>".26M!1?)".:/A.21b.2668F?DM!.:!]/}4t1#6L/d%(9?.2'2.0/"6A!?6Av1(4Ä/d%5/A!BD!]/ ¯ QÍK68.2"M$m ®"¡°¯²± ³± £ g z@%(6y g h 4t1#7%5'2's"1j¹Ê!$#Aµ1 ´!Q*÷~1 u ,G!#$c .2`C?#&%(M,/&.2M!%('7%5C"C"'2.:Mz%5/A.21?6w/&>?NM!#A.0/&$#A.21(hBq%z;I9J u 1#8/A>"'2$6A6E68.2"M$p/&>"$#A\Bq%z;9J\%G(!#L; 'º%5#Am#A1 u /&>N154È/A>"vCl1':;a"1(Bb.2%('s9J,4t1#A/A>"vi?Cl1"$]/&.º%5·' ¶L/d%5(!6}1±G(!# ¸D%(?)N9"#A.:"67%59J1Fj/ /&>?N)"!M!%±;(QP¸-F"#8/A>"!#8Bb1#8(c*1k%©M!1BDC"F?/A!#E68F"Md>Ã%©#81 u /&>PBb.:>]/E)"$68/Q±;h/&>"N681'2Fj/&.21( C"#81jM$!6A6~9l$4t1#8W.:/ZM!1(F"'2)\9l=/A!#8Bb.:-%5/&$)RQ Ém1(/A/&>"%5//&>?768/&%(9".2':.:/"}.:Bq%((.2-%(#L;=%ij.26 u >".2Md>E/A!")"6Ä/&1W9l 1aM!M!F?C".2$)b9];e/A>"Z!.2(!]G5%('2F"$6*#A!'2%5/&$)E/&1=9".2M$1]G$M$/&.:1DC"#A19?'2!BD6!Q|÷m1 u $G$#!cjM!1")?.:/&.:1 ¯ Q .26*1(4£'2.0/A/&':71#*"1@C"#&%5M$/&.:Mz%('-.:BbCl1#L/d%("M$}.:4£6A.2(".:+-M!%(]/%(Bb1(F"]/&6*1(4R)?.26A68.2C-%/&.21(b%5#AZC"#A$6A$/!Q ãpä º »«¼©ö lç C*òz·ì GIJ½Epï G?H©"ò Jzò$ïRõ Lí KnM¿¾ ñ O8ô QSR$R&V WYXDZN[]\_^U`9Z3\_^acb

{}>"KTWDXgM!1BDC-%(".:1p/&1ba/&%5/&$Bb$]/ ¯ QÍE.26 ¸s1(#w%©,¹(½l(^Bq%5/A#A. i cXy¥sQ ¯ Q·\.26^ Àqù,A1eÁe N,¹f'e2ùW.:4rc u >"! ¶ § .26 ¯ Q ¯ ® M!1"6L/d%(]/!c "§ ²^#8!Bb%(.2?679J1(F"µ ")"$)%(6 Å g h Q ~68!Z/&>-%/F"BD!#8.2M!%('"68/&%(9".:'2.:/".:6~)?$+-".0/&.:1N1(4Ä68/d%59".2':.:/"%(C?/&$#!c"6L/d%(9?.2'2.0/-%/o)?1g"15/o?!M!$6A6A%(#A.:':; 4t1':'21 u /&>"36A$Bb.:)".268M!#8$/&#81F?/&5Q S< 6AF?Md>ÀMz%(68!6U.:/\.:6 %(C"C?#A1C"#8.º%5/AN/&13M!1"68.2)"$#b68.2BeF?':/d%5"!1F?6A':;/&>"U$øJ!M,/&6b1(4v9l1(/&>P/&>?U/A.2BD\%(")Ã68C-%(M$%(Cj C"#81±ij.:Bq%5/A.21?6!QÀ/A.2BD,-%/ >-%(C?CJ$"6/&1K6A1BD"1(#AB1(4R/&>?v6A1':F?/&.:1 u .0/&>".:b/A>".26)"1Bb%(.2b%(6/A>"vBD!6A>q.2]/&$#8G5%('261w/A1 Ê!$#A1b%5/Z6A1(Bb=M!1("68/&%(]/m#A%5/&.:1"Q|=)".268M!F"686~/A>".26}Cl1.2]/Z1(4ÈGa.2 u .2oj$M$/A.21 ¯ Q["Q

Õ±ÖÂaÖáÒÈÙLâoÓqÃ×?ÑJÐ@ÎÓ ×ÇÒ?ÐvØÙrÚÇÙ!ÒZÜgÐvà/=Î}Ý}àÄÓ*ÔlÅZÓHàKÎÓ ÝZß Ý}^Ó?> ×

QSR$RUT [YaciPjlkZ`9Znmo^pqZrbDs tusv`9ZriPs

°Ì ¯

1"68.2)"$#y%68$/y1(4JTW^XPY6(1ÌG(!#8"!)D9];e%M$1BDC"'2,/&$.2$"68;a6L/&!BUQ{}>"68/d%59".2':.:/"bM!1")?.:/&.:136A$/e.2X*¥sQ ¯ Q ¯ cÄ4t1':'21 u 6@%5/K1("M!q4t#81Bè%o6L/&F")j;I1(4}X*¥sQ| ?Q·° %(").0/&6WM!1(BbC-%5".214t1(#vBeF"'0/&.2C?'2?l".:6UM$1")".0/&.:1»68/d%/&!6N/A>-%5/zcm4t1#NF"Bb$#A.:Mz%('W6L/d%(9".:'2.0/"È $.2$G5%(':F"!6 9J1(/A> C"#8.2"M$.2C-%('%(")36AC"F"#8.21F?6!cÇ.:4H/&>"$#Ab%(#A^%5;?®WBeF"6L/w':.2e1(1#=.:"6A.:)"e/&>"^F"?.:/=M!.:#AM$'2E.2/&>? M!1(BbC"':,i J".26*)",+-".0/&.21(E1(4£6L/d%(9?.2'2.0/"*CJ$#A.:1j)".:M,¹u M!1]G$M$/A.21=BD1j)"$' u >".:Md>KC"'2%(M!$6 /&>?É¶8M!1#8#A$M$/Ê¸\':1jM!%5/&.:11(4*/&>?EC"#A.:"M!.:C-%('Ë_Ìa#A1a1(/WC"#8!M$.26A$':;1©/&>"EF"?.:/WM$.2#AM$'2 u >?!#A^/&>? 6A1('2F?/A.21=.26È1"'0;W"$F?/&#A%('2'0;W68/&%(9"'25Q|¸-F"#L/&>"$#!c4t1#È%mM!1(BbC"':$/AH$.2$"68;a6L/&!BUc5/A>"H!.2(!]G$M$/&1(#A6 C"'2%±;q"1D#A1':W.2p/&>?@F"BD!#8.2M!%('R68/&%(9".:'2.:/"*)".:6AM!F?6A6A.:1=4t1#£/&>?!6AH68;a6L/&!BD6ÈC-%(#&%5'2'2$'26l/&>"y)".:6AM$F"6A68.21=4t1#Û)",4t!M$/A.:G(TWVWX}Y[6!QÄXi?%(BD.2? Xy¥sQÄ jQ:° U%(")3"1(/Ab.0/&6=6A.:Bb.:'º%(#8.:/IXy¥sQ ¯ Q @jQDD6A!D/&>-%/@4t1(#@)?$4t!M,/&.0GqT=^X}Y[6w/&>? #A$¥F".2#8!)UBb1a)".:+"Mz%5/A.21q/A1 ¯ Q E.:6 Æ "ÿ7® ® Æ © ° 4t1#Z%(':' |%5")Ç ¯ Q ® 6A.:"M!)"$4t$M$/&.0Gv6L;a68/&$Bb6})?1e "1(/F-%(#A%(]/&!9J1F?")"!)?"!686Z4t1(# Æ Æ ¨ °ca4t1#i"%5BbC"':m.:NXy¥sQ ¯ Q @e.:4 Æ Æ ¨ °w%5")o$.:/A>"!#Z z]®71#}R´ z]®D¨ ¬ z u @$/Z':.2"!%(#71#}¥F-%()?#&%5/A.2M=#A1 u /A>RQ ãpä Í Î òzö /ç Ï E\é GeìZÐ ç Epï G?H©"ò Jzò$ïRÑ õ Geê©ë ÒíDÔê ÓDlç Cyó^çWê©ìZçáí KoM¿¾ Pñ O8ô Õ $/UF"6U"1 u i"%5Bb.:"©/&>"M!1("M!$C?/154K6L/d%(9".:'2.0/-%5// w ýly® x{z?Q Í?QZ{}.:BbpBD$/A>"1a)"6Z%(#AW)?$G$'21Cl!) u >".:Md>oF"%(#&%(]/A!v/&>"%5/ Æ ý Â ® Æ xÞ°w%5")N/&>".:6 .:67/d%(5!o/A1^9JW/&>?@M$1")".0/&.:1p4t1#}F"BD!#A.:Mz%('£68/d%59".2':.:/-%5/W%p68/d%/&.21(-%(#8;U6L;j6L/&$BocÛ$"!#A%5/&$)©4t#A1B %N¢*VWX 1©681BD6Ör¹jù 6AC"%(M!bBb$6A>Rc u .2'2'|>-%zGp%UF"Bb$#A.:Mz%('y681'2F?/A.21/&>"%5/@.26@9J1F?")"!)I%(6= ¨ Å Â g h Qp{}>".:6 a±ùdw£Ä(F-%(#&%5/A!K/A>-%5/)"!68.2#&%59"'2@6A1':F?/&.:1"6m%(#AK!"$#&%5/A!).:/&>?w/&.:Bb@Bb%(#AMd>C?#A1aM!!686 %(6Z9l1(/&>N/&>?=/A.2BDW%(")U6AC"%(M!@Bb$6A>\.2]/A!#8G5%(':6Z%(C"C"#81]%(Md>UÊ!$#A1"Q Ém1 '2,/F"6m)"$+""e68/&%(9".:'2.:/"w/A.2BD,¹6AC-%(M$w68!"68(Q¸|.2#868/mM!1"6L/F"M$/=%^+-".0/&w/A.2BD, 6AC"%(M!@)"u 1Bb%(.2N'0;a.2" u .:/A>".2×zLxªØx½Ùþ%(")Úzxþyx{Û@Q1±G$#m/&>?.26})"1Bb%(.: u .0/&>U%^ #A.:) /&>"%5/Z.26}$¥aF?.26AC"%(M!$).2\9l1(/&>\/&.:Bbv%(")\6AC-%5M!w%(")\+?iN/&>?@BD!68>8(½v9];N/A>"=!¥F-%5/A.21 Ü ² ¨ µ EØ Émia/*#A!)?F"M!}1F?#*TWDX3%(C"C?#A1zij.2Bb%5/&.:1K1(4-/A>"}¢*VWX©/&1=%/ u 15".2674t1(#ABeF"'2%e.26 -§ ²z³J´ ¨ µ -§ ² ¯ Q0° z]® Xy¥sQ ¯ Q0° z\.26=6&%5.2)/A1o9lq6L/d%(9?'2^.:47%(];9l1F"?)"!).2".0/&.º%5'ÄG$M$/A1#!c _ ÝzcÄC"#A1a)"F"M$!6K%U9l1F"")?!) 6A1('2F?/A.21G(!M$/A1#!c -§ ²?cÈ%(6W/A>"^Bb$6A>36A>"#8.2"6W/&1UÊ!$#A1\4t1#=%\+ji?$) Ü§ ²jQD{}>?.26=.26v/A>"^M!'º%56A6A.:Mz%(' )",+-".:/A.21@1(4-6L/d%(9?.2'2.0/]/&BK;(!#Ä68/d%59".2':.:/?@BD!68>o.:]/&!#LG5%('26}1^/&1^Ê$!#A1?cs/A>"=F"Be9l!#Z1(4È/&.2BDW68/A!C"6$²c ßcsBeF"6L/m(1^/A1^.2?+-?.:/".26}6L/d%(9".:'2.0/"~6A.2(".:+-M!%("M$}1(4J/A>".26y)?$+-".0/&.:1^1(4£6L/d%(9?.2'2.0/"#A1F?> àÛ" ¹âá L ãd-ùd&4ù ÀEc >?.2Md> 68/&%5/&$6~/A>-%5/!c".:4|%EaF?Bb$#A.2M!%('RBD$/A>"1j)\.:6,¹ f'e2ù .:\/A>"=6A$"6Aw1(4 Õ %i"®}%(")Yd?,2,?! .:/.2L6 d ä(ù, åaù,-½Q Bb,/&>"1a).2L6 d?,2,rù,-È.:4*.0/C?#A1a)"F"M!$6W"1N!#8#A1# .:/&>"K{|%z;a'21#m6A$#A.2$6 6A$"6Az®.2D/&>"m'2.:Bb.0/*%(6*/&>"mBb$6A>p68C-%(M!.:"e%5")q/&>?m/&.:BbZ6L/&!Cq1w/A1KÊ!$#A1 u .:/A> Ü ²@+ji?$)Rc".: /&>?~>];aCl!#A9l1':.2M}Mz%56A±®Qy{}>".26*.:64tF"#8/A>"!#*)".:6AM!F?6A6A$)p.2bj$M$/A.21 ¯ Q ?Q»BD$/&>?1j)^.266 d ä(ù,+ åaù,- .:4|.:/~M!1(G(!#8!6/A1D/A>"Ki?%(M$/m6A1':F?/&.:1o%(6~/A>"@Bb$6A>68C-%(M$.2"q%(")o/A.2BD=68/&$C1D/&1DÊ!$#A1p.: /&>?.26}Bb%(""$#!Q '2!%(#A'0;c?/A>".26}.:6Z%(\.2BDCJ1(#8/d%5/7C"#A1Cl!#L/-%zG( Æ¦Æ -§ ² Æ¦Æ ¨ Æ¦Æ µ ² _ § Ý Æ¦Æ x ÆçÆ µ ² ÆçÆèÆ¦Æ ² § Ý Æ¦Æ x ÆçÆ µ Æ¦Æ ² èÆ¦Æ ² § Ý Æ¦Æ ¯ Q0°°Ì® 0 *!4éñì ò: %6 ($+ q* +)ï2ó+ ô!cDõö4 : )-+0 202**+!(y+ï20 !*+*+.÷"&:v+r .ê!î1ø*ï°!4' * ë"!4.ú01%Dù*&ì°ð'°!ï2ûÊ:)üË+ýçþ' !40 ü%* !4; ì ï°: ($.u &í"ø S+.6îÿø*!c%'& %*)*+ : !cë·($.Êî1ë$&*+ì1&%'3% (y.)&% .u %ê!4' *;0234*&4*+ì°Ë!ø: *%"ï ê(S**!S &)c:%" +ï !4(S: ï+* !4 &%µ):ù +!4°: ûÊïü+ýç þ !4%I'þËù)!·2ûyüý¦ þ +"4þâ.+ë/!y ÿ+ + +,!4: Uï + !4!%v:ï ; + !4%u!4v ² õ 3; v)+*øâø : :*+ê*+!y+2

±Õ ÖÂaÖáÒÈÙLâoÓqÃ×?ÑJÐ@ÎÓ ×ÇÒ?ÐvØÙrÚÇÙ!ÒZÜgÐvà/=Î}Ý}àÄÓ*ÔlÅZÓHàKÎÓ ÝZß Ý}^Ó?> × °Ì° j.:"M!W/&>?@.:".:/A.º%('l)-%5/&%eG!M,/&1#Z.:6}9J1F?")"!)£cs/A>"=6A1':F?/&.:1pG$M$/A1#m.:6}9J1(F"")"$)o.04 Æ¦Æ µ Æ¦Æ xÞ° ¯ Q0°Ì® u4t1>"#}$6L#A/d %(9"Æ¦Æ.:µ '2.0/".26W.:6v1(4f/&$F"6A$)I%(6=%o4 ª,ù,-|M$1")".0/&.21( Ém1 u»u "!$)p/&1K#8!'2%5/&~/&>"m68/d%59".2':.:/D/&>-%/ .0G!I.2hXy¥sQ ¯ Q:°ÌaQoS"U-ù,½81e8(Á]=1(4 cÈ.Q5Q:cÈ.:/A6W!.:!]G5%('2F?D154Bb%ij.2BeF?BnBb%(".0/&F")"5Q^S"$o/A>"wBb%5/A#A.0ip.26Z"1#8Bq%5'c?.Q5Q:c .0/}M!1BDBeF?/A!6 u .:/A>\.:/&67/A#&%(?6ACl16A5Q jQZ{}>"=68CJ$M$/&#A%('Ç#&%()?.2F"6}.:67/&>"L e2 ù, f,2 a² =1(4%('2'J"1#8Bb6$Q ¸-F"#8/A>"!#8Bb1(#A(c u >"$N¾ .:6?1#ABb%('c/A>"68!M$1")N.2?!¥F-%(':.:/?DF"BD!#8.2M!%('y68/&%(9".2':.:/?1j)"6}F?6A!)U/&1^6A1(':G@¢*VWX}Y6$Q ü {}>">"^$`/6L/d&%(>"9"p.:'2(.0/"$6A3M!#8.:/&$#A.2%h%(#A9l1(/&>þ?!M!$6A6A%(#8;O%(")»68F?DM!.2$]/U4t1(# BD$/A>"1j)?6*/A>-%5/!?!#&%/&~6AF"Md>pBb%5/A#A.2M$!6*%(?)q)"$CJ$")p6A1('2!'0;eF"Cl1D/&>"~!.:!]G5%('2F?!6*1(4 /A>"=Bq%5/A#A.:M!!6$Q ü S¹.:6q4/AF">""6KL/d68%(CJ9?$'2M$/A#&9]%(;þ'È#&%5%(;)?g.2F"6~M!#81(.:/A4~!#A(.:-1=RQÀ1±G{}$>#AF"".6p:"4tp1(#qBb%5$/&"#8.0!i\#A%(.:6m'(Bb#Az%5%//&.2!M$#Z!6!/&c}>"%(/A1>""6ACl5cl!M,/A/&>"#AK%('BD#A$%(/A)">".21aF?) 6 M$1")".0/&.21(N.26}"$M!!686&%(#L; W9"F?/~"1(/~68F?DM!.2$]/~4t1(#}68/d%59".2':.:/",/&>"$#1#*"1(/*/A>"m68!BD.0¹)".268M!#A,/&~%(C"C?#A1]%(Md> u %56H/&%(($q/&1=+-")D/&>"~)".:øJ!#8!"M$.2"w%(C"C?#A1zi .2Bb%5/A.21p154%E6A,/~1(4¢*VWX}Y[6!c-/&>"W+-"%('Û)".:øJ!#8!"M$=!¥F-%5/A.21"6ZM!%(o9lW#A$C"#A$6A$/A!)9]; -§ ²z³J´ ¨ µ -§ ²W« ¶ § ² ¸-F"#8/A>"!#8Bb1(#A~.:4 µ >-%(6%wM!1(BbC"':$/A: 97$.2$"68;a68/A!BUcj/A>"681'2F?/A.21D/&1@/&>"m>"1BD1!?!1F"6C-%(#L/ Mz%5%(' u %z;a6Z9J=i?C?#A!686A$)©%(6 -§ ² ¨ Ü ´ ²´ ØR§ ´ Á èèè Á Ü qÿ ÿ² Ø-§ ÿ Á èèè Á Ü"'2 %5/&"ÿ.:1»q9%(l#8$W/ /&>"!W$N/$.2&>"$DG5%%(':5F"")!q6Z1(/&>?4 =µ ý\QyS¹$4|.2/&$>?]G@(68%5!'2F"BD$.06!¹Qy)".2{}68M!>"#A.:,67/&6AW$%5#8GC"$C"6Z#A1%(%(6}Md%w>G(2!#LF";p68M!!1()RGc (u !w".:M!!%(]/ F".:)"W%(67/A1 u >"!#8 u u =Bb.:>]/7i?Cl!M,/Z/&>"J¹#A1a1(/&6/A1E'2.:#A$'º%5/A.:G(v/&1e/&>?WF"".:/7M!.2#8M!':v.2N/&>? M!1(BbC"':,Éi J¹C"'º%5"(Q\¸-1#W/&>?.26@#Az%56A1 u u .:'2'yC"#81jM$!$)I/&1U/&#&%5M!b/A>"D'21aM!F"6@1(47/A>" J¹#A1a1(/&6 %(6%D4tF""M,/&.:11(4/A>"KC"%(#&%(BD$/A!#ý Â 4t1#Z/&>?K!¥F-%/&.21("6BD1j)?!'2.:"^)".:øJF"68.21%(")Cl!#8.21a)".2M M!1(G(!M,/&.21(>=ÌQ ?A@CBEDGFH@JILKMONPNQRBSKTKVUR>BWN

|¸ .:F"#A ¯ Q:°E68>"1 u 6W/&>"Ei"%5M$/W/&#A%(M!E154H/A>"µJ¹#A1a1(/v.04.0/v.:6v!?!#&%/&!)39];ö¡ ¢X #A!C?#A!68!]/&.:" !.0/&>"$#^)".:øJF"6A.:1k1#^9".2M$1]G!M,/&.:1RQ SgM!%(6A$6b/A>" ü #A!C?#A!68!]/&6b/&>"U68/&%(#8/A.2"G5%('2F" uBb>"%($#A". 0/&Â`F")"¨ *1(z4Jý %5")].2?M!#A!¨ %(6A°($6!Q\c¸s/A>"1(#@/A/A#&>"%(qM$)"#A.0ø£$C"F?#A6A$.21(6A$B/AD.2"1jm)?/&!>"'c|ý )"Â [email protected]"#A%5z/&%5.:1'wB%(?)DI1a)""!$']>"%5!/&%(.0G)"56|Qo/A1 m6w%(/>?)"6 /&>?@1(#A.2(.2p%56ý ÂpÂ ¨ g «Dh QHTvN/&>"=1(/A>"!#}>"%(")Rc"4t1#/&>"=9".:M!1]G$M$/A.21\BD1j)"$'c"ý Âp{¨ u 'YÂ .26K%(' u %z;a6e.2Bb%(.:-%(#8;(Q\6K/A>"pBb%(".0/&F")"D1(4 YÂ .:"M!#8z%(68!6!cy/A>"b/A#&%(M$p#A$C"#A$6A!]/A.2"©/&>? 9".:M!1]G$M$/A.21bBb1a)"!'"/A#&%zG$'26%(#A1(F"")^/&>"~M!.:#AM!F?BD4t$#A!?M!~1(4£/A>"mF?".:/HM!.:#AM$'2(c u >".:Md>q.0/*"$G(!# '2!%zG!6$QD^BeF"6L/w9l^Mz%(#8$4tF"'y.:.2]/&$#AC"#8$/&.:" u >"!.:/=.26W#8!C"#8!68!]/&.:"o¡:Z\[ X Q^{}>"E4½%(M,/ /&>"%5/.0?/ e ùd^p/A>"KF"?.:/mM!.2#8M!':@BDz%5"61?':;N/&>"%5//A>"N1 À]²e f aù}1(4/&>?K#A$C"#A$6A$/&%5/&.:1.:6 M!1(#A#A$M$/!c£.:/Z/A!'2':6mF?6"1(/A>".2?q1(4/A>"^ d?±ù7!#8#A1# 6A$eXy¥sQR ?Q["°Ì®QZ{}>"KC">"%(6Ae$#A#81#m#8!'2%5/&$6 /&1E/A>"_ -±,f^/&>?@F?".:/ZM$.2#AM$'25Q

. +`ceï2d *& ÷·!4%6 + ï ë4!)ì !4aÿïµ&+ : b 3f4!)4 +0°÷ . +2êu)&:v02+& 4484* +S44(y; b 4+&&%2*h +ë"gI +ê&!($*A i ,2&4.I Ë8ì14:ï %µ . !4.k*j& & ð'.ï2;c4:_ô2*+!m!({ln·% 4!%2øÉ.!oâ%úê!4ø*y ø /: :_. -!4î(Sï2 + !4êï %L*4 %2. * ì !%'8&+ !4%/&4 +& &÷*+ 01&+ 84&:ë;

Õ±ÖqpÖÃàsr|âÓ*ÔZÙ,Î|ÐvÚP×ÇÒ?ÐvØÙrÚÇÙ!ÒZÜÀÎ}Ý}àeÎÓ*ÑÒ|×oÙrà ÒÄÏÓáÎ}ÝZâÑ*ÚÇÓutNÃÑ*ÚÛÐvàÓ °±Íj° v]QCR>wyx z{N|F}@ILF:@wyN~wyNEKMO@GF

mÉ 1 u '2,/vF"6WM$1BbC"%(#AK/A>"E,i?%(M,/J©6A1BDK/&>"%5/=%(#8EC"#A1a)"F"M$!)9];©%(M,/&F-%(' /&.:Bb?.2"EBD$/A>"1j)?6!Q7{Ä%(9"': ¯ Q0°@6A>"1 u 6/A>"eý?U3#A$'º%5/A.21?6Z4t1(#m%^G5%(#A.:$/"1a)"6!Q ¸|.:F"#A$6 ¯ Q·K%5") ¯ QÍK.2':'2F"6L/&#A%5/&}/&>?#A$6AF"'0/&67C?#A1a)"F"M!$)N9];bG5%(#8.21F"6HBD$/&>?1j)"6 u >"!p/A>"$;q%(#8 %(C"C?'2.2$)/&1b/&>"eBD1a)"!'|TWVWXZY[64t1#)".:øJF"6A.:1%5")©CJ$#A.:1j)".:M,¹M!1]G$M$/A.21Rc£Xy¥a6$QR"Q[N%(")©?Q@jQ S¹/W.:6v.2BDC"'2.:!)U/&>-%5//&>"E9l!>-%zGa.21(#v6A>"1 u .26m/"q.2BDCJ1#L/d%(]/ Cl16A68.29".:'2.:/A.2$67%(#8WM!1±G$#A!)9;q/&>"=.2':'2F"6L/&#A%5/&.:1"6!Q ° ©« þ°Z«Pý Â¨ z |X ijC"'2.:M!.0/~X*F?'2!# « ký Â ©« þ° ¨ z ´ Õ z%(C?4t#81 Í « ° Á U ý Â ® Á ý ´ Â¨ z } _ « ° Á ´ ý Â ® Á ´ = ý Â ©« ´ 9 ý Âo¨ z 7Í _ @ G°Z«ký ´ Â ®|«þ° ¨ z ´ S z ÞLW °±Í °~« ´ ý Â Â ®|« ý Â Â « = ° ÁÂ ´ ý Â ® ¨ Â z ¨ Æo.:'2"W(/A> {Ä%(9"': ¯ Q0°QHj1BD@ýb« ÉÞ~$'º%5/A.21?6

Rv]QCxGz'BEK}vmDOz'NUTNKMO@C

|¸ .:F"#A ¯ QN6A>"1 u 6K#8!68F"':/A6@4t1#W/A>"bijC"'2.:M!.:/@XyF"':!#@Bb,/&>"1a)RQD»>"!F"6A$)I4t1#=)".2686A.:C-%5/&.:1j )"1BD.2"%5/&$)KMz%56A!6y.0/.:6Ä6L/d%(9"':4t1#È/&>"7#A%("*L xþý Â x/z?Q ' 68F-%('2'0;W/&>"Bq%((".:/AF")"H1(4JýD>-%56 /&1K9Jm!68/A.2Bb%5/&$)p%(?)p1(4f/&$p.:/.26*4t1F"?)p9];b/A#A.º%5'J%5")p!#8#A1#&®,Q|»>?!NF"68!)p4t1#9".2M$1]G!M,/&.:1 /&>?l J*/&#A%(M!v4½%(':'261F?/&68.2)"v/A>"=F"".0/}M!.:#AM!':4t1#Z%(':'J+-?.:/& Â cs%5")N/&>"WBD$/A>"1j)\>-%56}"1^#&%("( 1(4Ä68/&%(9".2':.:/?.26}Mz%56A(Q ?ANR>xOIU@NEKMk@C

{}>".:6N.26\%I/ u 15"1a)Rc~6A.2?M!©/&>"$#A%(#8/ u 1 Y[6NC"#A1a)"F"M$!) 9];þ$G(!#8;OýÛQ »>?! %(C"C?'2.2$)h/&1©)".:6A68.2C-%5/A.21I)"1BD.2-%5/A!)hC?#A19"':!BD6 u U68!N4t#A1(B ¸|.:"Q ¯ Q·/A>-%5/e/A>"NC"#A.:"M!.:C-%('

°±Í]

ÎÏmÐÑvÒÄÓ*ÔxÕ±ÖÃ×ÇÒ"ÐØÙrÚÛÙ!ÒZÜÞÝ}ßPÚÛÙràÓyÐvÔá×ÇÜ*×ÇÒÄÓHâ× Ι(σ)

Ι(σ)

λ h= 0

R (σ)

λ h= οο

σ = e λ h , λh

ωh= 0

- oo

a) Diffusion

R (σ)

σ=e

i ωh , ωh

oo

b) Convection

|¸ .2F?#A ¯ Q0°û*X|i?%(M$/Z/&#A%(M!$6m154ËJ¹#A1a1(/&64t1#}Bb1a)"$'R!¥F-%5/A.21?6!Q #A1a1(/@.26=6L/d%(9"':E4t1#@%o#&%5"^1(4Zý Â c|9"Fj/@/A>"D6AC"F"#8.21F?6w#81j1(/@.:6@?1(/zQpS"D6AC"F"#8.21F?6 #A1a1(/W6L/d%(#L/&6W1/A>"EF"".0/WM!.:#AM$'2E%(?)4½%('2':61F?/A6A.:)"E1(4H.:/m4t1#^1eÁe7w ý Â ®l©z?Qe÷~1 u $G(!#$cÇ4t1(# 9".:M!1]G$M$/A.21Mz%56A!6$c u >"!Iý3.:*6 ²j&ù^ÁÀDå(fl(cJ/A>"EBb,/&>"1a)©.26"1(/=1"'0;6L/d%(9?'2(c£9"F?/=.:/ %(':6A1©C"#81j)"F?M!!6^%þ/A>-%5/e4½%(':'26wC"#A$M!.268!'0;h1h/A>"NF"".0/EM!.:#AM!':q.:h/&>"N#A%(" z x x¡°Q 6 u %(6^CJ1(.2]/&$)k1(F?/D%(9l1ÌG((c*/&>".:6E)"1a!6^"1(/^BDz%(£c*/A>-%5/E/A>"\Bb,/&>"1a)`.26 u .0/&>"1FjY/^Â !#A#81#!Q '0/&>"1(F">o/A>"=+-F"#8K6A>?1 u 6m/&>-%5/~/&>?!#Aw.26m%b#A%("@1(4 YÂ .: u >".:Md>/&>?K'2!%(C?4t#81bBD$/A>"1a) C"#81j)"F?M!!6H"1=!#8#A1#*.:b%5BbC"':.:/AF")"(c5.:/*6A%±;a6H"1(/A>".2"@%(9l1F?/y/A>"Z!#8#A1#*.2d?±ù,QÆo1#8~.:6y6&%(.:) %(9l1F?/}/A>".26}.:>"%(C?/&$u# ?Q B>N|BW@OC HUONUkRwFRF:MOI@JUKMNEKMO@

{}>".:6E.26^%(':6A13%/ u 15¹#A1a1(/EBD$/A>"1j)P9"Fj/zcHF""'2.:(p/A>"\'2z%5C?4t#A1©68Md>"!BD(cH/&>"N68C"F"#8.21F"6^#81j1(/ 68/&%(#8/A6@%5/W/&>"D1#8.2.:Rc£#A%5/&>"$#W/&>-%(31©/&>?DF?".:/=M$.2#AM$'25cÈ68!D¸|.:"Q ¯ QjQE{}>"$#A$4t1(#A(cÇ/&>?!#A^.:6 %o#A%("D1(4}#8z%('H"$]%5/&.0Gqý Â 4t1(# u >".:Md>I/&>?qBD$/&>?1j) u .2':'y9Jb68/d%59"'25QN{}>"q+-F?#AD6A>"1 u 6 /&>"%5/}/&>"=#A%(" W$")"6 u >"!ý Â © «@W° z^6A.:"M!@%5/}/A>-%5/ZCl1.2]/7/A>"=6AC"F?#A.21(F"6Z#A1a1(/}'2!%zG!6Z/&>? M!.:#AM$'2^%(") Æ Æ 9l!M$1Bb$6=#A!%5/&$#=/A>-%(31"5Qq{}>"D6A.0/&F-%5/A.21©.:6=¥aF?.:/&^)".0ø£$#A$/ u >"$I/&>? ý?? C"#8.2"M$.2C-%('H#A1a1(/@4½%('2':6w1(F?/&68.2)"q/A>"pF"?.:/eM$.2#8M!'2D4t1#E%(':' YÂ z?cy%(")h4t1(#w/&>?N9".2M$1]G!M,/&.:1

Õ±ÖqpÖÃàsr|âÓ*ÔZÙ,Î|ÐvÚP×ÇÒ?ÐvØÙrÚÇÙ!ÒZÜÀÎ}Ý}àeÎÓ*ÑÒ|×oÙrà ÒÄÏÓáÎ}ÝZâÑ*ÚÇÓutNÃÑ*ÚÛÐvàÓ °±Í(Í σ1 σ1 λ h=-2

a) Euler Explicit ωh = 1

σ2

σ1

σ2

σ1

b) Leapfrog

σ2

σ1

σ2

σ1

λ h=-1 c) AB2

Diffusion

|¸ .:F"#A ¯ Q·jûy{È#&%5M!!6Z1(4âJ¹#A1a1(/&674t1#7G5%(#8.21F?6}Bb,/&>"1a)"6!Q BD1j)"$'R!¥F-%/&.21(\/A>"=Bb,/&>"1a)\.26}F"?68/d%59"'2W4t1#Z%(':' Â Q A¡hUR¢xN|£E@'OR>z#NKMO@C

Convection

}{ >"7/A#&%(Cl!Ê$1.2)-%5'jBD$/A>"1j)K.26%G(!#L;eCl1C"F"'2%(#|1?4t1#|#8z%(681"6/A>-%5/*%(#87C-%(#L/&.2%('2'0;W.2':'2F"6L/&#&%/&!) .23¸|.:"Q ¯ QÍ?QES¹/A6J"^F"".0/=M!.2#8M!':e4t1#W9J15/&>©/&>"^)".:6A68.2C-%5/A.2"\%5") /&>?qCl!#A.:1j)?.2MEM!1(G(!M,/&.2?Mz%(68N%(")£c|.24½%(M,/zcfm2^,¹f'e2ùL¤zb1eÁeä1esùdN6¤Ký Â ¤zËdad ýÀft±4ù e ¤q2\f9 d-ù,&ù,-+ e , f'e2ù,Q (F"6L/D'2.25q/&>?U':z%(Cj4t#A1©BD$/A>"1j)P.0/E>-%(6e/A>"\Mz%(C-%59".2':.:/-%56A$#A#81#b4t1#E/A>"Cl!#8.21a)".2MNM$1]G!M,/&.:"`Mz%(68(c}9"Fj/b/A>"!#8.26b%Bb% ¥81# )".0ø£$#A!?M!~9J,/ $!p/&>?Z/ 1e68.2"M$Z/A>"Z/&#&%5CJ$Ê!1.:)-%('-BD$/&>?1j)DC"#A1a)"F"M$!6H"1e%5BbC"':.:/AF")"7$#A#A1(# 4t1#@(- YsÂ ¦ u "1( / ¥8F"6L/%Eu ':.2BD.:/&$)p#&%5"W9l$/ u !!zLx YÂ xÞ°Q

°±Í(

ÎÏmÐÑvÒÄÓ*ÔxÕ±ÖÃ×ÇÒ"ÐØÙrÚÛÙ!ÒZÜÞÝ}ßPÚÛÙràÓyÐvÔá×ÇÜ*×ÇÒÄÓHâ×

NJ§R>£kR>NEKMO@

{}>"*=%(Ê!)-%5WBb,/&>"1a) u %56H)?!6A.:"$)D/A1WC"#A1a)"F"M$~':1 u C?>-%(6AZ$#A#A1(#!Q*j.:"M!}.0/&6Md>"%(#&%(M,/&$#A.26L/&.:M Cl1':;a"1BD.º%5'|4t1# .26@%UM!F"9".:M {|%(9?'2 ¯ Q:°(cÄ"1?Q°z®,cÈ.:/@BeF"6L/e>-%zGD/ u 16AC"F"#8.21F?6K#A1a1(/A6w.: %()")?.:/&.:1/A1\/A>"EC"#8.2"M$.2C-%('Ä1"5Qe{}>"!68q%5#AE6A>?1 u 3.2¸Ä.2"Q ¯ QÍjQKS/&>"^)".:6A68.2C-%5/A.21 %(")o9".2M$1]G$M$/&.:1M!%(6A$6!c£%b6AC?F"#A.:1F"6~#A1a1(/~'2.2BD.:/A6/&>?w6L/d´ %(9".:'2.0/"K)".:6A68.2C-%5/A.2"^M!%(6A5c %p6AC?F"#A.:1F"6v#81j 1(/'2!%±G(!6v/A>"eF"".0/WM!.:#AM$'2 u >?!ý Â ©¤« cÛ%5")4t1#m/A>"E9".:M!1]G$M$/A.2"NM!%(6A5c u >"$ Y¨Â Q*Ém1(/A=/A>-%5/Z9l1(/&>\68C"F"#A.:1F"6Z#A1a1(/A6Z%(#A='21aMz%/&!)U%5/}/A>"=1#A.:.2 u >"!ý ¨ z?Q Ie©ªkN|BW@GRO¬«k@DkUKMC HUkN|UsDGkJNW®PDkKVKVRyNEKMk@CGF3©u®~¯°RO±®P²

{È#&%(M$!6\1(4=/&>"öl"6A!M$1")jb%5")þ4t1F"#8/A>j"1a)"6\%(#8 6A>?1 u Ã.2g¸|.2"Q ¯ QÍjQ {}>"o+"F"#A$6b68>"1 u /&>-%/q9l1(/&>kBD$/&>?1j)"6b%(#8o6L/d%(9"':N4t1#D%#&%5"\1(4 ý >"$3ý .26m#A!%('%(?)"!%5/&.0G(c£9"F?/v/A>-%5//A>"e#A%("e154*68/d%59".2':.:/? u %±;^/A1@LaQ ?c u >"$#A!%(6SÞL@@.26H'2.:Bb.0/&$)e/&1@rjQyTvD/&>"m1(/&>"$#>"%(")b4t1(# 9".:M!1]G$M$/A.21öÞL@^.:6mF?"68/&%(9"'2=4t1#m%('2' YÂ c u >"$#Az%5l6 Þ*Wb#A!Bb%(.:"6Z.2"68.2)"=/&>"@F"".0/M$.2#8M!'2 4t1c# zx YÂ x ³ aQ*Tv"vMz%5p6A>"1 u /&>"%5//A>"l ÞLWK68/&%(9".:'2.:/".2MdÚ> w ýJS® x{z?Q QSR.SRUT $µ `>4A¶l^k^U`1tF)ac\·µ/i¸4kkT¹»º

¹S /H.26*.2]/A!#A$68/A.2"=/&1@C"F"#86AF"~/&>"~¥F"!6L/&.21(q1(4£68/d%59".2':.:/?!b/A>"Z/&.2BD}68/A!Cb6A.2Ê$(c Â ca.:6*6ABb%('2' 6A1b/&>-%5/%(M$M!F"#A%(M$;o1(4*%(':'Û/&>"eý"¹#A1a1(/&6~.26~1(4y.:BbCl1#8/&%("M$(Qmj.:/AF-%5/A.21"6}4t1# u >".:Md>/&>".:6m.:6m"15/ /&>?@M!%(6A@%(#8=M!1"68.2)"$#A$).2U>"%(C?/&$#6jQ 8z'5'OF"1a)"6vC"#81j)"F?M!^%NC"#A.:"M!.:C-%('È#A1a1(/v/A>-%5/W.:6G$#8;©M$'2168 /&1 ¡ ¢X 4t1(#p6ABb%('2'~G5%('2F"$6\1(4Eý Â Q {}>"!#8$4t1#8(cv1(»/A>"9-%(68.26N1(4@/&>"C"#8.2"M$.2C-%5'#A1a1(/!cm/&>? 68/&%(9".:'2.:/"1a)N/&>-%5/Z.:6}#A$¥aF?.2#A$)U/A1^#A!681':G(K%e/A#&%("68.2$/Z681'2Fj/&.21(N1ÌG(!#m%E#A!'2%5/&.0G$':; 6A>?1#8/m/&.2BD@6AC-%(Bq%z;U9JK%qBD1j1(/~.:6A6AF?(QWaF"Md>©Mz%(68!6W%5#Aw/?1j)"6 >"$/&>",;%5#Aw%(C"C"':.2!)\/&1q%b9".:M!1]G(!M$/A.21C"#A19?'2!BUQ}¸|.:6!Q ¯ Q5Mw%(") ¯ QÍ4y6A>"1 /&>"%5/m4t1(#~9lu 1(/A>BD$/&>?1j)"6}/A>"KC"#8.2"M$.2C-%5'R#A1a1(/Z4½%('2':6}1F?/A6A.2)?@/A>"wF?".:/~M!.:#AM$'2@%(")o.26~F""68/&%(9"u': 4t1#m%('2' YÂ QZ÷~1 u $G(!#$cR.:4Ä/&>"=/&#A%("6A.:!]/v681'2Fj/&.21(U154y.2]/&$#A$68/vM!%(9l@#A!681':G(!)©.:%D'2.2BD.:/A!) F"Be9l!#*1(4l/&.:Bb76L/&$C"6H/A>-%5/H%(#8}6ABb%('2'j.2E/&>?~68!"68m154J/A>"}+-F"#8(c]/&>"Z$#A#A1(#yMz%(F"68!)b9;K/&>".:6 .2?68/d%59".2':.:/?w#81j1(/7>-%()\4½%('2':!N.2?6A.2)?v/&>"=M!.:#AM$'2v/&>? BD$/&>?1j) u 1F"':)>-%zGE9l!$)"!M$'º%(#8!)6L/d%(9"':e9"F?/W%(!#8#A1#1(4y/&>"K6&%(BD ÀD å(-f aù u 1F"':) >-%zG~9l!!bM!1BDBD.:/A/A!)Rc¥8F"68/H.2E/A>"Z1C"Cl16A.0/&})".:#A!M,/&.:1RQ¸-1#/A>".26y#8z%(681b/A>"Z_}=%(")^/&>? ÞLweBD$/A>"1j)?6~>"%±G(w9l1(/A>U9l!$F"68!)o.2\6A$#A.21(F"6Z¥aF"%(]/&.:/&%5/&.0G=68/AF")".:!6Z.2]G1(':Ga.2"^Cl!#8.21a)".2M M!1(G(!M,/&.21(RQþ{}>".26D.2")k1(4.2"6L/d%(9".:'2.0/"1a)"6 M$1]/zY[)s®,Q

2

°±Í

ÎÏmÐÑvÒÄÓ*ÔxÕ±ÖÃ×ÇÒ"ÐØÙrÚÛÙ!ÒZÜÞÝ}ßPÚÛÙràÓyÐvÔá×ÇÜ*×ÇÒÄÓHâ×

½ R¢KRF KU@xGMOB¨OF/H%59J1Fj/ 9];K/A>"7,ij.:68/&$"M!Z1(4sM$!#L/d%(.:e/"wBb16L/m%(M!M$F"#&%5/A(cJ,ijC"'2.:M!.0/zc"/ u 15"1a) 6A$w{Ä%(9"': ¯ Q0°Ì®,û -²±³J´ ¨ «~-² Á @s² ¯ ´ Á AÂ ¾ 5 ² Á ² ¯ 6´ ¿ ¯ Q0°±Í]® Tv"NMz%5g68>"1 u cF"68.2"©/A>"NBb,/&>"1a)"6E.0G$k.:gj!M,/&.21(g jQ ?cy/A>-%5/e/A>".26EBD$/A>"1j)P.:6K/&>".:#A)j 1#8)"!#%(M!M$F"#&%5/AK9l1(/&>o.2U/&!#8Bb6~1(4 ¡VÀ ¢ %(")ö¡VÀ:Áacl6A1q4t#81B %(©%5M!M!F?#&%(M,;oCl1.:/m1(4Ga.: u .:/~.:6 %5/8/&#&%5M$/&.0G5Qy÷m1 u $G$#!c-':$/F"67.2?6ACl!M$/.:/A66L/d%(9".:'2.0/".:6 Mz%53!%(6A.:':;U9l^%(M!M$1BbC?'2.268>"!)9;68/AF")?;a.2"\.0/&6vMd>-%5#&%(M,/&!#8 .26L/&.2MKCJ1'0;a"1BD.º%(' u >"$Iý Â g z?Q ¸-#A1B X*¥sQ ¯ Q0°±Íb.0/~4t1('2'21 u 6}/&>"%5/Z4t1#vý Â\¨ z?cÂ 'Ã ® ¨ Ã Á Ã « @jQ7¸l%5M$/&1(#A.2?sÂ Ç® ¨ z u!¥qF-%5+-'R")/&»1UÂ °c- Û®%(")×¨ «c-%^À6AC"°Ì®F? #A.21(Á F"6}@(1® "¨ 5c-!z?¥QUF-{}%('Û>"/&!1e#8N @´%(Ä%#AÄ D/ u o1 J"qC?#A.2?M!.2C"%('*1"(c S?.26e#8!6AF?':/zcH1"NBKF"68/DF"")?!#A6L/d%(?)g>"1 u BD$/&>?1j)"6 u .0/&>\6AC?F"#A.:1F"6Z#A1a1(/&6 u 1#8p.2\C"#A%(M$/A.2M!5QyIW"1 u /A>-%5/}/&>?$;U%(#8W"1(/Z6A$':4Ä68/&%(#8/L .2?"cl%(")/A>"K68CJ$M!.º%5'ÈC?#A1aM!!)?F"#A$6vMd>"168!/&1b68/&%(#8/m/&>"$B¼.:".:/A.º%(':.2Ê$/&>"wM!1a$DM!.:!]/&6m1(4y/&>? 6AC?F"#A.:1F"6Z#A1a1(/&6$c-/&>" Ü ÿ ® 4t1/# Ç °w.:Xy¥sQl ?QÍ( ?QS¹4Ä/&>"=6L/d%(#L/&.2?bC"#81jM$!686m.:6 u !':'Û)"!68.2?!) /&>?!6AbM!1a$DM!.:!]/&6K%(#A^4t1#8M!!)3/&1U9J^G$#8;6ABb%('2'cÇ%(").:4H/&>"bBb,/&>"1a).:6=68/d%59"'25cÈ/A>"$;3,/ 6ABb%(':'2!# u .0/&>g.2"M$#Az%56A.2? Å Q»÷m1 u ,G$#!cZ.:4m/&>"UBq%5".:/AF")"\´ 1(ÝL4Å 1"U1(4v/&>?o68C"F"#A.:1F"6 .:6 !¥F-%5'Ä/&1 @jcl1"KM!%(68!e)".:6&%(6L/&$#.26~.2BDBD.2"$/~9l!M!%(F"6A @® ° zWÆ!QmXG$%^G$#8;o6ABb%('2' .2?.:/&.2%('"G5%('2F?m1(4 Ü ÿ ® .26¥aF?.2Md':;b1±G!# u >"!':Bb$)RQHjF"Md>\BD$/A>"1j)?67%(#8vMz%(':'2$)qMz%/d%(6L/&#A1(C">".2M!%('2'0; F""6L/d%(9?'2@%(")U%(#8 u 1#L/&>"':!6A6Z4t1#}BD168/!c".:4È"15/~%('2'c?M!1(BbC"Fj/d%5/A.21"6$Q 8z'kNPROOR>wFK ¼xNÇwyNEKMk@CGF

¹S 4 u w.:"6ACl!M,/m/A>" J"wBKF"':/A.2C"':v#A1a1(/~Bb,/&>"1a)"6~.2o¸|.2(6!Q ¯ Qb%5") ¯ QÍ?c u +-")/&>?!B /A1\9JE1(4H/ u 1U/"p5")j¹1#A)"$#Km)-%(BD6L6 %(6A>?4t1(#8/&>I%(")@%(Ê$)-%(oBb,/&>"1a)"6!QUS".26@Mz%(68p%(':'*6AC"F"#8.21F?6 #A1a1(/A674½%('2'*©ÿ d-ùE å(f u >"! Â g z?Q {}>"m4t1(#ABD!#y/"1j)%(' u %z;j6W6L/d%(#L/&6v1(/&>"eF?".:/WM$.2#8M!'25cJ/A>"eBD$/&>?1j).:6v'2.:($':;N/&1N9l!M$1BD F""6L/d%(9?'2v4t1#7681BDvM!1BDC"':,iNýo%(6 Â C"#A1aM!$!)"6~% u %±;q4t#A1BÊ!$#A1"QHTvN/&>"W9"%(6A.:671(4Ä%e{|%z;a'21(# 6A$#A.:!6vi?C"%("6A.:1RcR>?1 u $G(!#$cÇ/A>"!68EBb,/&>"1a)"6%(#Ae(!"$#&%(':':;\/&>"KBb16L/W%(M$M!F"#A%5/&K.2"681(4½%(#v%56 /&>?$;NBD.2".:Bb.:Ê!wÿ d-µù d±{ù ª,ù,-£.:N/&>"='2!%()".:"e/&$#AB4t1·# ¡VÀ £ Q {}>"3'º%5/8/&!#b/"Z9-%(68.26y1(4s/&>?~Bb%(?.:/&F?)"1(4l/&>?}M!1a$DM!.:!]/ .2e/A>"}'2!%()".:"v{Ä%±;a':1#6A$#A.:!6y!#8#A1#/A!#8Boc(/A>"$;K6AF?øJ!#$c#A!'2%5/&.0G$':;w6AClz%(.2?"c(4t#A1(B¤%(M$M!F"#A%(M$;(Q ¸-1#\%`.:G(!O%(Bb1(F"]/\1(4KM$1BDC"F?/d%/&.21(-%(' 1#8lcm/&>"1#8)"!#\1(4K%(M!M$F"#&%5M$;þ1(4@/&>?/ 1 /"K6ABb%('2':"!686v1(4y/&>?e/&.2BD,¹68/A!CRc uM!1(q68/~6A$1(!4È©/&>"/&1oWM$Bq1Bb%5(C?^F?/&/d>"%5D/A.216AN.:Ê!D(!1("47$/&#&>?%(':.2':6=;b68)"/A!!C`Cl!%(")6@?6m6ABb1\/%(':'*&>?%(W6waFCl?1Be686A9J.29?$'2#Z(1(Qp4|Tv6L/&/$C"&6~>?q/&%(1((/A$>"U/!#@A1^>-M!%(1")BD£cÄC"/&F?>?/A %p6A1('2F?/A.21RcR%(")©/&1NBD.2".:Bb.:Ê!W/&>?.26 .26A>©/A1UBb%(5K/&>"^6L/&!C68.2Ê!^%56W'º%(#8E%(6WCl16A68.29"':(Q S"@M!1BDC"#A1(Bb.:6A(c-6L/d%(9?.2'2.0/"1a)"6$cv>"1 u ,G$#!c=/&>"6A.0/&F-%/&.21( .: u >".2Md> %('2'v154w/A>"3/A#&%(?6A.2$]/U/A!#8Bb6U%(#8 #A$6A1'0G$)oM$1"68/A.:/AF?/&$6%E#&%/&>"$#~BD.2?1##81'2W.:U6L/d%(9?.2'2.0/h9ý Ìa$.2$]G(%5'2F"$6K>-%zGa.2" u .2)"$':;

ü "F %(6A$]/!&g$)1/&"'#A:%(;P"6A/A.:1I!]/+-z")Q »%h6L/&z%5)?;]? u S`1(4m/A>"!68\Mz%(68!6E/&>"$#A\,ij.:68/^.2h/A>"\!.:!"6L;a68/&$Bb6^#8!'º%/&.:G(!'0;I'2%(#AbG(%5'2F"$6e1(4 Æ ý Â Æ %(6A681jM$.º%5/A!) u .0/&>U!.:!]G$M$/A1#A6~/&>-%5/ u u .268>U/A1q)"#8.:G(@/A>"#A1F?>U/A>"K681'2Fj/&.21(UC?#A1aM!!686 u/&>?.:e/A>"Bq1F?% /N¥81%(#];BD1(#A/&$.0G5]%5%(/.:)»1N44t1(t#q1#m/&>"/A>"$.2e#b68.2/A")F"")?.0;Ga.2)"1(4*F-%5F'"Bb%(M$$M!#AF".:Mz#A%(%('ÈM$; 68/&%(.2þ9".2':!.:.:/N/&>?@)?$+-".0/&.:1Rû ÃF"BD!#8.2M!%('Bb,/&>"1a)K.:c6 a d² (fl1 eÁe ~,¹ f'e2ùÇ.04-.:/È.:6|68/&%(9"':H4t1#Ä%(':'?TWVWX}Y[6 /&>"%5/~%(#A=.:">"!#8!]/&'0;N68/&%(9"':(Q nBD$/A>"1j) .:/A>g/A>".26DC"#A1(CJ$#8/?D#8!.:11(468/du %59".2':.:/?1j)"6}.:6}6A>"1 u o.2\{|%59"'2 ¯ QjQ Æo,/&>"1a) Tv#8)"!# Þ à â ° zz zz S-%/Z"1"=1(4È/&>"$6A=BD$/&>?1j)"6}>"%(6~%(U%(M!M$F"#&%(M,;\>".2(>"!#7/A>-%(\6A$M!1")a"^1#8)"!#=1(4%(3~"1j)?6\/A>"/&#&%5CJ$Ê!1.:)-%('vBD$/A>"1j)»>-%56U/A>"6ABb%('2':!68/ /F""M!%5/&.:1\!#A#81#!Q

Õ±ÖªÌ]ÖkàÍr|âoÓHÔZÙ,Î|ÐvÚP×ÛÒ"ÐØÙrÚÛÙ!ÒZÜÞÎ}Ý}àeÎÓ*ÑÒ|×UÙràÒÄÏmÓÞÎ}ÝZâÑ*ÚÈÓut¤ýÎèÑ*ÚÛÐvàmÓZ°zÍ Þ~,/&F"#8".2?p/&1p/&>?E68/d%59".2':.:/? M!1(")".:/A.21 ° à â « Þ ¨ Á %(")\/&>?=/ u 1D6A$/A6m154|.:"!¥F-%(':.:/A.2!6#A$)"F"M!=/&1^/&>?@6A%(BDW6A$/ u >".:Md>o.:6 ¯ Q1 z]® â x Þ «þ° ° ¯ Qj°Ì® â ß « ÷m$"M!5cj/ u 15-%(#8 /&>?q6&%5Bb n àhå â åÞ ®=C-%5#&%(BD$/A!#@6AC-%(M$(QU'0/&>"1F?>/A>"q1#8)"!#@1(4Z%5M!M!F?#&%(M,;I1(4}%(h~¹68/&%(9"'2 BD$/&>?1j)DMz%(?"1(/7ijM!!$)p/ u 1?c *ª Ü ".:>\1#A)?!#!Q S¹/N>-%(6q9J$!»6A>"1 þ/&>"%5/p4t1#q%hBb,/&>"1a)Ã/&1I9JS¹68/&%(9"'2.:/pBKF"68/p%('2681h9J~"1a)"6E>-%zGa.2"3%©Md>-%(#A%(M$/A!#A.:68/A.2MNCl1':;a"1BD.º%5'*4t1(# u/&.:1>"-%.:Md(>':':;/&>"68E/d%5M$9"1j'2,5bQNM${}.2$>"/.2v6w1(.24*"M$/&'2>"F"e)"$>?6K.2%5>"'2$'y68,/Wij1(C"#A':)".2M$$.:#/@/&Bb$#A,B¼/&>"1a.2©X)"6K»%(.2")6mIF"C?".#A:/\BD$/A>"1j)?6}%(#A#A$4t$#A#8!)N/&1"ca/&>"$#A,4t1#A5c-%(ú6 d² (f½ 1 eÁe E,¹ f'e2ùBb,/&>"1a)"6$Q QSRÒ$R $µ `>4A¶l^k^U`1tÑ[]ayby`9a1æD\s ^&b¿`XDZ [Yaci jlkZçéè1ÔÕê k64bDZcR

OG(!#L;pM!1(G(!".:!]/ u %±;b/&1eC?#A!68!]/}/&>?68/&%(9".:'2.:/".:"wBD$/A>"1a) .26E/A1IC"':1(/^/&>"U'21aM!F"6^1(4/&>"UM!1BDC"':,igý Â 4t1# u >".2Md> Æ Æ ¨ °c768F"Md>Ã/A>-%5/D/&>"U#A$6AF"'0/&.2? M!1(/A1F"#e(1j$6K/&>"#81F">/&>"pCl1.:/Eý Âk¨ zjQ©÷m$#A Æ Æ #8$4t!#86e/&1/A>"pBb%ij.2BeF?Bè%(9"6A1('2F?/A G5%('2F?1(4l%5;úc5C"#A.:"M!.:C-%('1#|68C"F"#8.21F"6$c5/A>-%5/.:6%~#A1a1(/Ä/A1/&>?Md>-%(#A%(M$/A!#8.268/A.2MCJ1'0;a"1BD.º%('(4t1(#

ë d *7&)%$ì1r($4.r&ð'ï2):2+!ï2% -ÿëyì"ë$Ë+*8"):2%2!4*(S4:bì&)+ !4% " . 8" !4%$ì"ëS!%2 ÿ)%'² %2.&0°&%2.%'!) !(q;3ì %2 )+ !40 %2*!" !)!49 ê0 !4*: :!!.ø3ï2 Uê* !4!)( )+: 2:- 4*r!'!4+° &; 3 q!"#$&%' 3!4_ ï2l$01!4: ë"% !4(S):v4*+ +ï ($ 3!48Ê4*+ !4ï°

jSg3&

°zz

ÎÏmÐÑvÒÄÓ*ÔxÕ±ÖÃ×ÇÒ"ÐØÙrÚÛÙ!ÒZÜÞÝ}ßPÚÛÙràÓyÐvÔá×ÇÜ*×ÇÒÄÓHâ× ï I( λ h)

í I(λ h)

î I(λ h)

Stable Unstable

1

Stable

ï R(λ h)

Stable

a) Euler Explicit

1

Unstable

î R(λ h)

ð θ=0

b) Trapezoid Implicit

ñò θ = 1/2

Unstable

c) Euler Implicit

í R( λ h) θ=1

ó ôöõW÷¢øùú´ûýü>þhÿ>ô%ô S÷>ø Wø¢ù ù ! û

õWô#"Wù%$'&û)(* +,{ø ÿ´ù-ô ú|û/.01 2 >ù3ô¢ù3 4>ô5-6 S÷>ø¢ù3÷ù:øô-7 ù8 !ô9:% ù +> ô% ù ;>ù2 >ù:ø ùøù@{ù:øA Çô6ù ôö ô# 9-6 X W÷>ø:TW ø,Y ù@Z!U% ôûú´ûýüù:ø ôS" ù6%{ø ^>ù_ > ù` ø: ! a ù ! õWô#W" ù¨ ôP ÿ´ù6-@eô c´ b ûe´ d ûýü>û

, f -

fgYh'igYjlkOmCn*gVk2op4qsrJtvutGixwyoi{z|g4}sm ó ô%õ>ûú´ûýü~: +3>ù:>ôöô 2-6 S÷>ø Wø¢ùù@Z!U%ô-õ

+5 [+7 !*ô ô U>ô -64 BW=¢ô ö#ô [-6 S÷>ø Wø@ù Z!U%ô -ù_- W÷>øÛù-6 ù390>ô#eù_UY Wø:ô L l¢ù%ùGW 0#|-6 [U%ù@Z$Y&! U>ù3û

?

ûS]>ùø ù:õWô % BW=¢ôöô# sô EGN @Ev >ù_' S÷0SøW ù:ø ù Wø ùôô nû +hù"Sù:ø6ôTù@ ~ô-6%÷>ù6 2U03ø: 1)>ùô~SõWôSø: %{Z¢ôC{ù@Z!-:ùU8 Wø>ùA SøôöõSôV ?:

ômô_ ÷W=% ù9W ø3>ùC~ >ùB> ô-6 W" ù-ô ?> U ø: ö ùP û_8W ÷>õ ùW" ùô>ù0ûô5-6 ùLô,: +Í ôT ó ôöõ¢û|ú´û\ü û\S]> ùa'S ÷0SøWH ô)¢ù]ô~S õWô0S øW O´ Z ô>ù~| ÷8T ùù % ùG ¢ôO´ Z ô: û S+5 ~ >ù:ø ù !80107S" ù>ô]> U ø: U ù:ø: LS øù¢ù_ +5 1*:eùUcø U ù6W ô01T ù@ cCTù@ %>÷>ù [âCù:ù

Ú0ÛÜ × Ö Ú0Û1Ý ×Þ &,ß Ú0Ûà Ü × Þ Ú0Û1à Ý ×Pá

& ß Ú àÛ7Ü × Þ Ú àÛ Þ Ú àÛOÝ ×Pá . 1eô-ù6 ùù !8Søù_ l>ùãPô>ù UùSû Ú0ÛÜ × Ö Ú0Û1Ý ×Þ

fgYh'igYjlkOmCn*gVkLugYhl}stGi{tg'h|0rvrämi¡srvotJwqBrvtJutGiAwyoi{zlg}sm å ÷WC-ù-7S ÷8 ù%¢ S ù8 !ô ô~U% ôôSøù9: +°ôÇó ô%õ>û ú´û/æ´ûÎ¢ù9 )¢ù6ù}ô3¢ù9301~PWãc W÷8 . ø:! W ø:>ù:ø3T S ù@ Åv>

ûAæ8ù:ømô3>ù}ü1!* Wø>ùù9ù !ÇõSô"Wù 2¢ùUYW ôX U ù

û ¢ . û3(*mô3W=% ù MON0FHI W ø$Öç vè +> ù6?é%ê è ê ¢ . û (*W=>ôöô# ~Y W÷80Sø: sô"Sù:ø: 2ôôSø, C0O Wø¢ùöù71U8{ø WõCù £vùû¢ú´û/d1V û

ü.

¦§¦§ÀÁT¹x¬s¯°P®,¯¸´|µ8²¶°P³4°6µ·Å²ªa²³l·)´!°W´ I(λh)

I(λh)

3.0 3.0

2.0

Stable only on Imag Axis

Stable

2.0

1.0 1.0

-3.0

-2.0

R(λ h)

-1.0

R(λ h)

b) Milne 4th Order

a) AM3

ó ôöõS÷>øù^ú´ûæ¢þØÿ>ôöô W÷>ø Wø>ù - >ô#eô 0 [:öùmô~U8öô-ù87 {+¿>ù6À ÷ù >ô 1~- : ¢ø W÷>õ W÷8~eôùcÒU-:ù1-6-6 Sø>ô>õy >ù1÷8~U8ô 8ÍõWô"Sù6Ò1' "Wù3û¸S]>ù6Õ ¢ù ô[UY ù9¢U1eô5S ø[ > ô- x1 ô> ô#eô\"1ö ÷>ù ;>ù T ù 19 1X9>ù;:|÷¢ù6:ô þ >ôk ô[Uö ô#e÷¢ùõWø: +ÂW ø3>ù-77~ ôcô ù= ùCW1 + ùùù:ø:T ô> ù6¢ c0>ô> õ ¢ù -6 >ôô 8 ÷>ù.?9@:)ÖBADCDEGFHAJILKM

nú´û

N @R@NÓK= ùJ:÷01ô ÓÅöô%ù6ô;¢ùø>õSù ô9¢ %÷8eô y ¢>ù E éCêROST > êVUûhÿ´ô -ø: !¢÷-:ù õWô#"Wù6ÅTù@ 1ù:- ù6:S øW £-6 ¢ôeô ¢ WøW=>ôöô# øùC17ùC+ Wø::UY ô%ùC-6 9¢ô01ô

sU0S øT ù@eùùxU>ø: !-÷>ø ùx-7£-ù:ù}ù@Z!UÔ3ô>ù£ Çù`Z8~U8öù6ùô#eù5>ô¢-ù:ø ù6-ø: 7Z¢ô Oeô _ 3¢ù~ >ù¢ô¢-÷ô ù{:|÷1eô T1X +x\Ûô-=Sø: ÐÑ ù8 !; 1 {S" ù:ø:ÔU¢ U ô> õ¢WeùU: û£S]> ô+AT ùXô ¢ ù6»ô» ÿ´ù6-@eô ü>û »õWô#W" ù C! :Jû ü>û þ Ú ]W

T@ ÛÜ × [ Ö ÛOÝ [ ÛD[ ß Û=[ Û=[ Ú ]W × ¤ Þ £ T> Ø Ú ] W Ü ×c¥ Ú ] W Þ Ú ] W Ý × á

nú´û üX

ÿ´÷8:eô#e÷8ô % B!:Jûú´û ô [!:û-ú´û ü}õWô#"Wù>ùøùÔ1ô

Ö Ý ×

T@ ß A ILK ` M Þ A Ý ILK ` M ÞP£ T> Ø ¥ á

Wø Ø

ü

£

T@

Þ§¦ Q> Ø ©

¥

- O*T>'%¨ ª%«

¬

¥

Ö

é

nú´û d

S]÷y:JûGú´û ô_% %÷8eô y:JûOú´û ü ô~ ô_ ø y:JûOú´û d´ûxS]>ùx +5 Pø: ! ú´û dCS øù Ø Ö Ø × ¥® ® Þ {ø +> ô-= ô#ôC-6ö ù6S øA01C ¢ ù Côx+7 ¯ û£>yù2ù ß Ú ] W Ü ×³¥ Ú ] W Ý × á

( %>ô]-7 ù cÖ

¥

nú´û b

T@ F F ôO*T> Ë T>

{ø +>ô-= ô^ô_-6%ù7Sø_1Ì ³¯ Wø) :ùùL÷>ùõ ÷[U8eô Hô[ó0 W÷¢øôöùôöô A01 ôØô\0O)>ù× ~1eø ô#Z'ùù:øP ô> ù+> ù¢ù>ô ùõ :-=>ù6ùôaU¢÷8]ôc>ù Sø¿ \y:Jû ú´û U ø: !÷>ø ù3W øT0>ô> õA¢ù @E#Ì@N0ÃOËFHQ0=¥ ¢ùù6Ó> U ø ö ù~ ùx-=S ô:- ùC 1~ äW ø>ùC:> ôö ô U03 ø ùù:ømô¢>ùHó0W ÷>ø ôöù:ø0# ! ôùùT ù?h + ù_07W" ù_ ù:ù6 ÷ ô> õHô%' lW ÷>ø> U ø ù|" ôS ÷ >ô:- ÷8 ô ÷8Søø ôS" ù6c1]{ø »[¢ô¢-ù:øùU-ùô::- øù@eùAU> U ø: -=ùø ù6 ÷83> U øù ù6ø1>ùùÛù õ^ô5 ù õA-7S ø øôöùL W÷8X ¢ùT% ù7U!{ øW õ] ù !8ûVú´û/d!ô4T ù@ ôC ÷W=8ö ù\W ø ËFGF ù:ô%õWù6"1ö ÷¢ù6Û + ô ùø: U' ù6¢ T ù@ cS ø*ôXù:õWø1eô> õùJ:÷01ô a )>ô U-ù3ûa> ùAU>ø ù6 ù6eù£¢ô3ù ô! :Jûü>û 110 ù^ó0W ÷>øô%ù:ø] ù ! ô¨ÿ´ù6-@eô ~ú´ûbú´û ( ?L ô-=0S ø: ÓÐaeô ù ô% ô ~+8 1X + û_a + ùS " ù:øù5-6 ôÔ1ù6# 31B¢ô\1U> U ø -= +5W ÷~-ùL÷W=8ö ùSû) ù6 ô#*>ô::- ø ùeù ù !TôT-71 ù*ù`Z8¢ U øù6: ù¨ô ~1eø ô#Z9 =Oeô ]¢ùW :eù¿ \Î_Ï]ÐÑ: þáã!ä ã

ä

ä

Ú Ö @

£ ÛL ? ? Ú Þ ¥ T> Ø

®å

nú´û æ

+ÛôÅ>ùc%ù7U!{ø Wõ?ù8 ! ÷ùÀ WøC>ùeôTùL~Sø-= ûÕÎ_÷>ø 10 ôHô«>ô U8eùù> U ÷¢øôS ÷ Y ø Ûô%>ùöù71U8{ø Wõ[ù8 ! ôÛ÷8:=1öù_ WøaÓ÷-= -71 ù6ûú´û/d1 û S]> ùT ù@ ä+H ÷ ù6ÅX «Û ô-=S ø: äW ø9h + ù71>ù:ø[> U øù>ô-@eô ô9v-[-6 + ùa:W ÷>ø::- ùa Ó: ù]ö ù|" ô# Wû\( ~8> U ø > ô% ô ô#:ø> õSù U> ù ù60´ ûs> ù]-6W ÷ø: ù1|< ÷8 ù>ù W ø:>ù:ø)Û ÷> õSù@Pm ÷8:=3 ù A _ôeù:õSø1ù :Jûú´û æc ô ! ù:õOeôS" ùL$4Ð {Ø ç û £ S]> ù:ø ùS ø ù3~ [+7 ¢ ô> U ÷1eù]>ù3 ÷ ù:ø ô-71YW=¢ ôö ô# C l% õW ø ô8~ ù ¢ù6 ô5 ôø >÷:- ùT ô©´ Z ù6 ô ù38LU0ô ùù %÷8eô L l>ù0:÷0¢ø1ô-mô Ö

ð - GO*Q>C

ð Ø ô Ø O*T> ¥

Þ ð S]>ùô-=ySø ù2ê Sø Çøù6 ð ô¢>ù}ø>õWùxéê ð êöõ·û S]>ôaù710Oa>ù_Tù@ ô QNÓK=MONYEG¤JEÍMONYËFGFHI9¥@¤ÍËD6F ÷ S]>ù5' "SùØø ù6 ÷| ù:ù6[4 Ûõ ! ù\eø ÷>ùù ] K=MON*ú ¥@E¥@¤ @NÓK@I +Ûô%>ù_ Wøô%õWôX,Ïû Ss x> U ø: W" ù¢ô6ùù:ø:~#ô%!:ûú´û 1 ôL9Ss7 ! Wø ù:øô%ù68¨ø ù6- :ø÷-@ ¢ùaUSø:ôÔ!>ô¢-ù:ø ù6eô>ùJ:÷0Oeô x+hùSøù-@e÷01 "|ô>õ_1,>ùTù[ô:ù-ù-6 ù6\"Sù:øW ~nûuó0 Wø>ùeôTù¢ù:øô#"11eô#"Wù+hù07"Sù

T@ ñ Ú ] W

ÛÜ × [

¥ Ú ]W

Û1Ý × [ ç Ö â Ú W Û=[ E ] Þ ò b

T@ Ø â

Û=[ ELELE Ú ] W Þ

F=F=F

% Wø>ù_Tô©Z¢ùLeôù8c:U0-:ù_>ô ù Ø

Ûù6U-:ù£>ùyù:ø:~~ô¸!:Jûmú´û +Ûô# T@?~T>Ìý ¢ é û\y > ù08 â Ú Ö â @

T@ Ø W Û=[ â Ú Ú M M ] ¥öû T ELE >yü T> Ø â W Û=[ Ú ] ¥ M M M M T@ Ø T@ Ø â Û=[ ELELELE Ú ] W û T>!ü

Ö

>ùäY "Wùù`Z8Uô ?8

â Ø Ú £ â > Ø ¥

+>ù

Ø

â Ø Ú â @Ø

Þ

]

Þ

F=F=F

nú´û.é

=213ù;>ùyöôô1 nú´û.

üæ

±\¼²3ºµs®,¯¿¦§À´|µ²3¶°P³4°6µ·¸¹]Á ³4°Pª3®)²¯Â´|·,´|µs®T¢´

:û-ú´û úô ù8ËR:ËD@MOFHEÍK û,!:Jû ú|û/. ô ÆI ù0@RD@MOFHEÍK û,S]|÷8*ô# T@mý é~8T>ý éHô%÷-= ! c+7 y1 ` ` E ø ù6~Sô-6 :7ùsùJ:÷01ô y+hù%-@e÷0 ?: "Sù2 y¢ù2Tù@ ô M 7"Wùù{:|÷1eô ô:>ô¢õ ~ :ù6U ô: ù6 ù_ ô÷01ô Í ô] ÷[~S øô: ù6ôLSB1ö ù^ú´û¢ . û)(*b T@yý é[1Q>ý é ôL ÷-= c+7 y1 ` ` E øù6~S ôA- :7ù6 Ø ô 2a Q@:`ù2 ù ô-6 ôWeù6:û þ h + ùS" ù:ø6ôM ö ù6 L-6 :øS ôX >ùCeô ù9Weù6U?+> ô-= ô~ Ù ÷:_ ùW" ùùx:> ôö ô c-6 >ôô 1T>ùx ù-6 W ø:>ù:øÛ ÷> õWù`Pm ÷8LT ù@ ? + ôL>ù=% ùSû !* Wø>ùõWù@ m÷8 Q@ê

Ï_÷ó0 Wø:´óø_1Sù6

ó0 WøLÿ>ôöô

` M Ø

T@êõ

>ô#eô 0 sÿ=1öù

8-6 ¢ôeô 0 Tÿ%ù

ó0 Wøa ô:ù6-@

a>ô Sø[ ô:ù6

¢ôeô 0 ~ 8ôWeùX6<

+Ûô# Ö E

3U8U>ø 7Z´ô~1eù Ö E

s ô Ø ß ` E £ ` M á S]¢ ù:øù@W øù

£ M

T@

Q>

£ M

SB1öù^ú´û.¢þØÿ´÷8~~Sø: B-6-ø-@ -6 8ôWeù- %- >ô#eô Wøa 1 Ï ÷Pó0W øW:óJø# _ 13ù6 ù !8:û ÖçS ø¢ ôeøS ø:s ù:ø øW ø]YW ÷ û ë2ì!

"çð)î[÷ø6ò$#

ý ã

û 8ô¢ù:ø]>ùÎÏ

ã

Ú0à Ö +Ûô# %ÙÖ

é

è)

¥

é é

¥

é # ¥

Ú Ö&% Ú ( Þ ' @

¥

é # ¥

+î * ï ?

'

Ö é

è)

¥

é +î * é ï

¦§-,X§Àº,¯¹]¶³|®T¢´

ü

J G ó ô2>ùmù:ô%õWù6"1ö÷¢ù6] .%÷ ô>õx9÷ù:ø ô-71U0-1OSõWù3û\>Õ01Ûô>ùWeù68 W=1ù_ %÷8eô SùÎÏÅ: ö÷!eô L-ù07"Wù^ôLeôTùD< v02> ø ô#eùC~-6 >ùA ôeù:õSø1ù{ø >ùô¢ôeô4-6 >ôô Ú Jé]Ö0/ ? ? 2143 {ø ô ù~@)Ö é]÷ ô>õ>ù#ù@Z!U%ô-ø7+ >ù*÷>ô#5-ù_ ùU 1Vû v-7Ûù6Uù7Oa>ùY {"Sù Sø>ùa ù û

.¢û]>Õ>ùA1UUöô%ù63 >ù,%ô>ù6SøÓ- "Wù6-@eô ù{:÷01eô ù\+Ûô>ù 18 {+â|OZ6X> ù6>ø: 2 ù ! õWô#W" ù: þ Ú ]

ÛÜ × Ö Û Ú ] ¥

× Û Û Ú ]Ü × ¥

Û Ú ]Ý × Þ

× ÛØ Ú Û Ü ] × ¥

Ú ]Û Þ Ú ]Û Ý ×

+>ùø2 Wø#ó\âl÷9ù:ø6kû7 ô>õ¨ó0 W÷>ø ôöù:ø W=> ôö ô# 2 ô6ùø¢õWù !× Û Wø+>ô-=%¢ùùË ¨ô]:öù3û ü>ûaÏ_ù@eù:ø:Tô>ùLU >ù:>ôöô 9-6 S÷>ø WøT>ù330[W :8 Wø:2ù8 !5

W:ø >ù:ø >:ø W÷>õ Tü>, û [U0Sø ù +Û#ô [> ù L8 ÷ >õW@ù m8 ÷ ù !8T 14 Wø ¢ù:ø >ø W÷¢ õ ü>û

d´ûaÛù6-6Ó>ù_-6 [U0-@-:ùXù:ø ù6.O*U' SôUU¢ø 7Z´ô~1eô Wø9ø::a>ù:ø ô"11ô"Wù3þ 98 Ú ] Ý üV98 Ú ] 9 8 Ú ] Ü Ö × Þ Ó Þ × M M M

. Q>

Ú ]Ü × ¥ Ú ]Ý ×

5 ]øùUÔ-õ>ù)U1eôSô>ù@Z;: >ù\eù~UY Wø1Sô8>ù@Z= õ3¢ôsS ø ÷¢ û4>Õ01B Wø>ùùùc[ >ù6¢ô¢-÷ô ùJ:÷01ô ä+Ûô U ù:ø ô >ô-YW ÷80Sø: -6 8>ôô õx3U8U ù>ô#Z%] û\ü8ùùù$8 øù61eô W ø]>ù !*W ø>ùøôeù:ø:ù->ù2¢ù:õ1ô"Sùsøù65OZ´ô6ù+Ûô $Y&cÖ õ~>ùùø öù æ¢ûl3U8U C_ Ï ô%øô-=8ö ù'S ÷0SøW 2- >ô#eô L1]¢ùö ùG' W÷80SøW Wû, xY W÷Sø: -6 ¢ôeô Íô U ù:ø:T ô#eùx1\>ùÛøôöõX\'W ÷80SøWW ûB> øô#eù>ù: :ù6 ÓÎÏ]Ð\+> ô-=T ø ù6 ÷\{ ø ø::W W ø:>ù:øs0- 1+3 øxU0OeôÔ1!>ô ù õ]ô9 Oeøô©ZX*W" ù- WølS ø¨ ûFa ô¢ õ33U8U ù>ô#Z û ù$ ù:ôöõSù6"1ö ÷>ù: û)> øôù>ù ] Î xÅ+> ô-=Ë øù ÷#]{ ø ¢ùAU8Uö ô-71ô

G ù@Z!U% ô ô> ä õ}ô% 1øä ô©ZÍW" ùä - SøW ø¢ < ôû\ù3û< Ö × Ú0Û ä Ú0ÛÜ × R ¥ GOÛ > øô#eù× ô«0>ùÀ~1eø ô#Z 8 =1ô «85õWô#"Wù%>ù¨ù6eø ôöù~ G û ô>õ?>ù¢$ øù61eô Ò WøL>ù ù@Z!U%ô- > 66 ;01 10 1 > 66 64 > ;1 0 > : ;2

39 > 77> > 77= 77> 5> > 0 ;

3 2 77 66 ;0ua 77 ~ 66 77 u + 66 0 15 4 0 2

(9.9)

The matrix in Eq. 9.9 has eigenvalues whose imaginary parts are much larger than their real parts. It can rst be conditioned so that the modulus of each element is 1. This is accomplished using a diagonal preconditioning matrix

2 66 10 6 D = 2x 66 0 64 0

0 1 0 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1 2

3 77 77 77 5

(9.10)

which scales each row. We then further condition with multiplication by the negative transpose. The result is

A = ;AT A = 2

2 66 ;10 10 1 66 ;1 0 66 ;1 4

1

1

32 77 66 ;10 10 1 77 66 ;1 0 77 66 154 ;1

1 0 ;1 ;1

2 66 ;10 ;20 10 1 6 = 66 1 0 ;2 0 64 1 0 ;2 1

3 77 7 1 77 1 75

1 ;2

1 0 ;1

3 77 77 7 1 75 1

(9.11)

CHAPTER 9. RELAXATION METHODS

166

If we de ne a permutation matrix P and carry out the process P T [;AT A ]P (which just reorders the elements of A and doesn't change the eigenvalues) we nd 1

1

1

1

2 66 00 66 66 0 40

1 0 0 0 1 0

0 0 0 1 0

0 1 0 0 0

32

32

2 66 ;21 ;21 1 6 1 ;2 1 = 66 64 1 ;2

3 77 77 7 1 75

0 7 6 ;1 0 1 77 66 0 77 66 0 ;2 0 1 76 1 77 66 1 0 ;2 0 1 77 66 0 75 64 1 0 ;2 1 75 64 0 1 1 ;2

0 1 0 0 0

0 0 0 1 0

0 0 0 0 1

0 0 1 0 0

3

17 0 77 0 77 0 75 0

(9.12)

1 ;1 which has all negative real eigenvalues, as given in Appendix B. Thus even when the basic matrix Ab has nearly imaginary eigenvalues, the conditioned matrix ;ATb Ab is nevertheless symmetric negative de nite (i.e., symmetric with negative real eigenvalues), and the classical relaxation methods can be applied. We do not necessarily recommend the use of ;ATb as a preconditioner; we simply wish to show that a broad range of matrices can be preconditioned into a form suitable for our analysis.

9.1.2 The Model Equations

Preconditioning processes such as those described in the last section allow us to prepare our algebraic equations in advance so that certain eigenstructures are guaranteed. In the remainder of this chapter, we will thoroughly investigate some simple equations which model these structures. We will consider the preconditioned system of equations having the form

A~ ; ~f = 0 (9.13) where A is symmetric negative de nite. The symbol for the dependent variable has been changed to as a reminder that the physics being modeled is no longer time 2

A permutation matrix (de ned as a matrix with exactly one 1 in each row and column and has the property that P T = P ;1 ) just rearranges the rows and columns of a matrix. 2 We use a symmetric negative de nite matrix to simplify certain aspects of our analysis. Relaxation methods are applicable to more general matrices. The classical methods will usually converge if Ab is diagonally dominant, as de ned in Appendix A. 1

9.1. FORMULATION OF THE MODEL PROBLEM

167

accurate when we later deal with ODE formulations. Note that the solution of Eq. 9.13, ~ = A; ~f , is guaranteed to exist because A is nonsingular. In the notation of Eqs. 9.5 and 9.7, 1

A = CAb and ~f = C~f b (9.14) The above was written to treat the general case. It is instructive in formulating the concepts to consider the special case given by the diusion equation in one dimension with unit diusion coecient : @u = @ u ; g(x) @t @x This has the steady-state solution 2

2

(9.15)

@ u = g(x) (9.16) @x which is the one-dimensional form of the Poisson equation. Introducing the threepoint central dierencing scheme for the second derivative with Dirichlet boundary conditions, we nd 2

2

d~u = 1 B (1; ;2; 1)~u + (b~c) ; ~g (9.17) dt x where (b~c) contains the boundary conditions and ~g contains the values of the source term at the grid nodes. In this case 2

Ab = 1 B (1; ;2; 1) x f~b = ~g ; (b~c) (9.18) Choosing C = x I , we obtain B (1; ;2; 1)~ = ~f (9.19) where ~f = x f~b. If we consider a Dirichlet boundary condition on the left side and either a Dirichlet or a Neumann condition on the right side, then A has the form A = B 1; ~b; 1 2

2

2

CHAPTER 9. RELAXATION METHODS

168

~b = [;2; ;2; ; ;2; s]T s = ;2 or ;1 (9.20) Note that s = ;1 is easily obtained from the matrix resulting from the Neumann boundary condition given in Eq. 3.24 using a diagonal conditioning matrix. A tremendous amount of insight to the basic features of relaxation is gained by an appropriate study of the one-dimensional case, and much of the remaining material is devoted to this case. We attempt to do this in such a way, however, that it is directly applicable to two- and three-dimensional problems.

9.2 Classical Relaxation

9.2.1 The Delta Form of an Iterative Scheme

We will consider relaxation methods which can be expressed in the following delta form: h i H ~n ; ~n = A~n ; ~f (9.21) where H is some nonsingular matrix which depends upon the iterative method. The matrix H is independent of n for stationary methods and is a function of n for nonstationary ones. The iteration count is designated by the subscript n or the superscript (n). The converged solution is designated ~1 so that +1

~1 = A; ~f

(9.22)

1

9.2.2 The Converged Solution, the Residual, and the Error Solving Eq. 9.21 for ~n gives ~n = [I + H ; A]~n ; H ; ~f = G~n ; H ; ~f +1

+1

1

1

1

where

(9.23)

G I + H; A (9.24) Hence it is clear that H should lead to a system of equations which is easy to solve, or at least easier to solve than the original system. The error at the nth iteration is de ned as ~en ~n ; ~1 = ~n ; A; ~f (9.25) 1

1

9.2. CLASSICAL RELAXATION

169

where ~1 was de ned in Eq. 9.22. The residual at the nth iteration is de ned as ~rn A~n ; ~f (9.26) Multiply Eq. 9.25 by A from the left, and use the de nition in Eq. 9.26. There results the relation between the error and the residual A~en ; ~rn = 0 (9.27) Finally, it is not dicult to show that ~en = G~en

(9.28)

+1

Consequently, G is referred to as the basic iteration matrix, and its eigenvalues, which we designate as m , determine the convergence rate of a method. In all of the above, we have considered only what are usually referred to as stationary processes in which H is constant throughout the iterations. Nonstationary processes in which H (and possibly C ) is varied at each iteration are discussed in Section 9.5.

9.2.3 The Classical Methods

Point Operator Schemes in One Dimension Let us consider three classical relaxation procedures for our model equation B (1; ;2; 1)~ = ~f (9.29) as given in Section 9.1.2. The Point-Jacobi method is expressed in point operator form for the one-dimensional case as h i jn = 21 jn; + jn ; fj (9.30) This operator comes about by choosing the value of jn such that together with the old values of j; and j , the j th row of Eq. 9.29 is satis ed. The Gauss-Seidel method is h i jn = 21 jn; + jn ; fj (9.31) This operator is a simple extension of the point-Jacobi method which uses the most recent update of j; . Hence the j th row of Eq. 9.29 is satis ed using the new values of j and j; and the old value of j . The method of successive overrelaxation (SOR) ( +1)

( ) 1

( ) +1

( +1)

1

+1

( +1)

( +1) 1

1

1

+1

( ) +1

CHAPTER 9. RELAXATION METHODS

170

is based on the idea that if the correction produced by the Gauss-Seidel method tends to move the solution toward ~1, then perhaps it would be better to move further in this direction. It is usually expressed in two steps as h i ~j = 21 jn; + jn ; fj h i jn = jn + ! ~j ; jn (9.32) ( +1) 1

( +1)

( ) +1

( )

( )

where ! generally lies between 1 and 2, but it can also be written in the single line jn = !2 jn; + (1 ; !)jn + !2 jn ; !2 fj (9.33) ( +1)

( +1) 1

( )

( ) +1

The General Form

The general form of the classical methods is obtained by splitting the matrix A in Eq. 9.13 into its diagonal, D, the portion of the matrix below the diagonal, L, and the portion above the diagonal, U , such that

A=L+D+U

(9.34)

Then the point-Jacobi method is obtained with H = ;D, which certainly meets the criterion that it is easy to solve. The Gauss-Seidel method is obtained with H = ;(L + D), which is also easy to solve, being lower triangular.

9.3 The ODE Approach to Classical Relaxation 9.3.1 The Ordinary Dierential Equation Formulation

The particular type of delta form given by Eq. 9.21 leads to an interpretation of relaxation methods in terms of solution techniques for coupled rst-order ODE's, about which we have already learned a great deal. One can easily see that Eq. 9.21 results from the application of the explicit Euler time-marching method (with h = 1) to the following system of ODE's: ~ H ddt = A~ ; ~f (9.35) This is equivalent to d~ = H ; C A ~ ; ~f = H ; [A~ ; ~f ] (9.36) b b dt 1

1

9.3. THE ODE APPROACH TO CLASSICAL RELAXATION

171

In the special case where H ; A depends on neither ~u nor t, H ; ~f is also independent of t, and the eigenvectors of H ; A are linearly independent, the solution can be written as 1

1

1

~ = c e1 t~x + + cM eM t~xM +~1 (9.37) | {z } error where what is referred to in time-accurate analysis as the transient solution, is now referred to in relaxation analysis as the error. It is clear that, if all of the eigenvalues of H ; A have negative real parts (which implies that H ; A is nonsingular), then the system of ODE's has a steady-state solution which is approached as t ! 1, given by 1

1

1

1

~1 = A; ~f 1

(9.38)

which is the solution of Eq. 9.13. We see that the goal of a relaxation method is to remove the transient solution from the general solution in the most ecient way possible. The eigenvalues are xed by the basic matrix in Eq. 9.36, the preconditioning matrix in 9.7, and the secondary conditioning matrix in 9.35. The eigenvalues are xed for a given h by the choice of time-marching method. Throughout the remaining discussion we will refer to the independent variable t as \time", even though no true time accuracy is involved. In a stationary method, H and C in Eq. 9.36 are independent of t, that is, they are not changed throughout the iteration process. The generalization of this in our approach is to make h, the \time" step, a constant for the entire iteration. Suppose the explicit Euler method is used for the time integration. For this method m = 1 + mh. Hence the numerical solution after n steps of a stationary relaxation method can be expressed as (see Eq. 6.28)

~n = c ~x (1 + h)n + + cm~xm (1 + mh)n + + cM ~xM (1 + M h)n +~1 (9.39) | {z } error The initial amplitudes of the eigenvectors are given by the magnitudes of the cm . These are xed by the initial guess. In general it is assumed that any or all of the eigenvectors could have been given an equally \bad" excitation by the initial guess, so that we must devise a way to remove them all from the general solution on an equal basis. Assuming that H ; A has been chosen (that is, an iteration process has been decided upon), the only free choice remaining to accelerate the removal of the error terms is the choice of h. As we shall see, the three classical methods have all been conditioned by the choice of H to have an optimum h equal to 1 for a stationary iteration process. 1

1

1

1

CHAPTER 9. RELAXATION METHODS

172

9.3.2 ODE Form of the Classical Methods

The three iterative procedures de ned by Eqs. 9.30, 9.31 and 9.32 obey no apparent pattern except that they are easy to implement in a computer code since all of the data required to update the value of one point are explicitly available at the time of the update. Now let us study these methods as subsets of ODE as formulated in Section 9.3.1. Insert the model equation 9.29 into the ODE form 9.35. Then ~ H ddt = B (1; ;2; 1)~ ; ~f (9.40) As a start, let us use for the numerical integration the explicit Euler method n = n + h0n (9.41) with a step size, h, equal to 1. We arrive at H (~n ; ~n) = B (1; ;2; 1)~n ; ~f (9.42) +1

+1

It is clear that the best choice of H from the point of view of matrix algebra is ;B (1; ;2; 1) since then multiplication from the left by ;B ; (1; ;2; 1) gives the correct answer in one step. However, this is not in the spirit of our study, since multiplication by the inverse amounts to solving the problem by a direct method without iteration. The constraint on H that is in keeping with the formulation of the three methods described in Section 9.2.3 is that all the elements above the diagonal (or below the diagonal if the sweeps are from right to left) are zero. If we impose this constraint and further restrict ourselves to banded tridiagonals with a single constant in each band, we are led to B (; ; !2 ; 0)(~n ; ~n) = B (1; ;2; 1)~n ; ~f (9.43) where and ! are arbitrary. With this choice of notation the three methods presented in Section 9.2.3 can be identi ed using the entries in Table 9.1. 1

+1

TABLE 9.1: VALUES OF and ! IN EQ. 9.43 THAT LEAD TO CLASSICAL RELAXATION METHODS

0 1 1

!

Method

1 Point-Jacobi 1 h i Gauss-Seidel 2= 1 + sin M + 1 Optimum SOR

Equation 6:2:3 6:2:4 6:2:5

9.4. EIGENSYSTEMS OF THE CLASSICAL METHODS

173

The fact that the values in the tables lead to the methods indicated can be veri ed by simple algebraic manipulation. However, our purpose is to examine the whole procedure as a special subset of the theory of ordinary dierential equations. In this light, the three methods are all contained in the following set of ODE's d~ = B ; (; ; 2 ; 0)B (1; ;2; 1)~ ; ~f (9.44) dt ! and appear from it in the special case when the explicit Euler method is used for its numerical integration. The point operator that results from the use of the explicit Euler scheme is ! !h n ! !h ! ! n n n n j = 2 j; + 2 (h ; )j; ; (!h ; 1)j + 2 j ; 2 fj(9.45) 1

( +1)

( +1) 1

( ) 1

( )

( ) +1

This represents a generalization of the classical relaxation techniques.

9.4 Eigensystems of the Classical Methods The ODE approach to relaxation can be summarized as follows. The basic equation to be solved came from some time-accurate derivation Ab~u ; ~f b = 0 (9.46) This equation is preconditioned in some manner which has the eect of multiplication by a conditioning matrix C giving A~ ; ~f = 0 (9.47) An iterative scheme is developed to nd the converged, or steady-state, solution of the set of ODE's ~ H ddt = A~ ; ~f (9.48) This solution has the analytical form ~n = ~en + ~1 (9.49) where ~en is the transient, or error, and ~1 A; ~f is the steady-state solution. The three classical methods, Point-Jacobi, Gauss-Seidel, and SOR, are identi ed for the one-dimensional case by Eq. 9.44 and Table 9.1. 1

CHAPTER 9. RELAXATION METHODS

174

Given our assumption that the component of the error associated with each eigenvector is equally likely to be excited, the asymptotic convergence rate is determined by the eigenvalue m of G ( I + H ; A) having maximum absolute value. Thus 1

Convergence rate jmjmax ; m = 1; 2; ; M

(9.50)

In this section, we use the ODE analysis to nd the convergence rates of the three classical methods represented by Eqs. 9.30, 9.31, and 9.32. It is also instructive to inspect the eigenvectors and eigenvalues in the H ; A matrix for the three methods. This amounts to solving the generalized eigenvalue problem A~xm = m H~xm (9.51) 1

for the special case

B (1; ;2; 1)~xm = m B (; ; !2 ; 0)~xm

(9.52)

The generalized eigensystem for simple tridigonals is given in Appendix B.2. The three special cases considered below are obtained with a = 1, b = ;2, c = 1, d = ; , e = 2=!, and f = 0. To illustrate the behavior, we take M = 5 for the matrix order. This special case makes the general result quite clear.

9.4.1 The Point-Jacobi System

If = 0 and ! = 1 in Eq. 9.44, the ODE matrix H ; A reduces to simply B ( 12 ; ;1; 21 ). The eigensystem can be determined from Appendix B.1 since both d and f are zero. The eigenvalues are given by the equation m m = ;1 + cos M + 1 ; m = 1; 2; : : : ; M (9.53) The - relation for the explicit Euler method is m = 1 + m h. This relation can be plotted for any h. The plot for h = 1, the optimum stationary case, is shown in Fig. 9.1. For h < 1, the maximum jmj is obtained with m = 1, Fig. 9.2 and for h > 1, the maximum jm j is obtained with m = M , Fig. 9.3. Note that for h > 1:0 (depending on M ) there is the possibility of instability, i.e. jm j > 1:0. To obtain the optimal scheme we wish to minimize the maximum jm j which occurs when h = 1, j j = jM j and the best possible convergence rate is achieved: jm jmax = cos M + 1 (9.54) 1

1

9.4. EIGENSYSTEMS OF THE CLASSICAL METHODS

175

For M = 40, we obtain jmjmax = 0:9971. Thus after 500 iterations the error content associated with each eigenvector is reduced to no more than 0.23 times its initial level. Again from Appendix B.1, the eigenvectors of H ; A are given by ~xj = (xj )m = sin j m ; j = 1; 2; : : : ; M (9.55) M +1 This is a very \well-behaved" eigensystem with linearly independent eigenvectors and distinct eigenvalues. The rst 5 eigenvectors are simple sine waves. For M = 5, the eigenvectors can be written as 1

2 1=2 3 2 3 2 p3=2 3 p p 66 3=2 77 66 10 77 66 3=2 77 ~x = 666 p 1 777 ; ~x = 666 p0 777 ; ~x = 666 ;1 777 ; 64 3=2 75 64 0 75 64 ; 3=2 75 p 1 1=2 ; 3=2 3 2 p 3 2 3=2 1=2 p p 66 ; 3=2 77 66 ; 3=2 77 6 7 ~x = 66 p 0 77 ; ~x = 666 p1 777 64 3=2 75 64 ; 3=2 75 p 1=2 ; 3=2 1

3

2

4

5

(9.56)

The corresponding eigenvalues are, from Eq. 9.53

p

= ;1 + 23 = ;0:134 = ;1 + 12 = ;0:5 = ;1 = ;1:0 1

2

(9.57)

3

= ;1 ; 12 = ;1:5 4

p

= ;1 ; 23 = ;1:866 From Eq. 9.39, the numerical solution written in full is p ~n ; ~1 = c [1 ; (1 ; 3 )h]n~x 2 1 + c [1 ; (1 ; 2 )h]n~x 5

1

2

1

2

CHAPTER 9. RELAXATION METHODS

176 1.0

h = 1.0 m=1

m=5

Values are equal

σm m=4

m=2

Μ=5 m=3

0.0 -2.0

0.0

-1.0

Figure 9.1: The ; relation for Point-Jacobi, h = 1; M = 5. + c [1 ; (1 3

)h]n~x

3

+ c [1 ; (1 + 12 )h]n~x p + c [1 ; (1 + 23 )h]n~x 4

4

5

5

(9.58)

9.4.2 The Gauss-Seidel System If and ! are equal to 1 in Eq. 9.44, the matrix eigensystem evolves from the relation

B (1; ;2; 1)~xm = mB (;1; 2; 0)~xm

(9.59)

which can be studied using the results in Appendix B.2. One can show that the H ; A matrix for the Gauss-Seidel method, AGS , is 1

AGS B ; (;1; 2; 0)B (1; ;2; 1) = 1

9.4. EIGENSYSTEMS OF THE CLASSICAL METHODS Approaches 1.0

1.0

h

h < 1.0

De

ng

cr

si

ea

ea

si

cr

ng

De

σm

h

0.0 -2.0

0.0

-1.0

λ mh

Figure 9.2: The ; relation for Point-Jacobi, h = 0:9; M = 5. Exceeds 1.0 at h

∼ ∼

1.072 / Ustable h > 1.072

1.0

h > 1.0

In

cr

h

ea

si

ng

h

ng si ea cr In

σm

M = 5 0.0 -2.0

-1.0

0.0

λmh

Figure 9.3: The ; relation for Point-Jacobi, h = 1:1; M = 5.

177

CHAPTER 9. RELAXATION METHODS 2 ;1 1=2 3 1 66 0 ;3=4 1=2 77 2 66 0 1=8 ;3=4 1=2 77 3 66 77 66 0 1=16 1=8 ;3=4 1=2 77 4 (9.60) 66 0 1=32 1=16 1=8 ;3=4 1=2 77 5 ... 64 ... . . . . . . . . . . . . . . . 75 ... 0 1=2M M The eigenvector structure of the Gauss-Seidel ODE matrix is quite interesting. If M is odd there are (M + 1)=2 distinct eigenvalues with corresponding linearly independent eigenvectors, and there are (M ; 1)=2 defective eigenvalues with corresponding 178

principal vectors. The equation for the nondefective eigenvalues in the ODE matrix is (for odd M ) M +1 ) ; m = 1 ; 2 ; : : : ; (9.61) m = ;1 + cos ( Mm +1 2 and the corresponding eigenvectors are given by j; m M +1 ~xm = cos m sin j ; m = 1 ; 2 ; : : : ; (9.62) M +1 M +1 2 The - relation for h = 1, the optimum stationary case, is shown in Fig. 9.4. The m with the largest amplitude is obtained with m = 1. Hence the convergence rate is (9.63) jm jmax = cos M + 1 Since this is the square of that obtained for the Point-Jacobi method, the error associated with the \worst" eigenvector is removed in half as many iterations. For M = 40, jm jmax = 0:9942. 250 iterations are required to reduce the error component of the worst eigenvector by a factor of roughly 0.23. The eigenvectors are quite unlike the Point-Jacobi set. They are no longer symmetrical, producing waves that are higher in amplitude on one side (the updated side) than they are on the other. Furthermore, they do not represent a common family for dierent values of M . The Jordan canonical form for M = 5 is 2

1

2

2 h~ i 66 h~ i 66 ; 6 X AGS X = JGS = 6 66 4 1

1

2

3 77 2~ 3 777 66 ~1 77 777 1 55 4 ~ 3

3

3

(9.64)

9.4. EIGENSYSTEMS OF THE CLASSICAL METHODS 1.0

179

h = 1.0

σm Μ=5

2 defective λ m 0.0 -2.0

0.0

-1.0

λmh

Figure 9.4: The ; relation for Gauss-Seidel, h = 1:0; M = 5. The eigenvectors and principal vectors are all real. For M = 5 they can be written

2 3 2 p 3 2 3 2 3 2 3 3 = 2 1 = 2 1 0 66 3=4 77 66 p3=4 77 66 0 77 66 2 77 66 00 77 ~x = 666 3=4 777 ; ~x = 666 p0 777 ; ~x = 666 0 777 ; ~x = 666 ;1 777 ; ~x = 666 4 777 (9.65) 64 9=16 75 64 ; 3=16 75 64 0 75 64 0 75 64 ;4 75 p 9=32 0 0 1 ; 3=32 1

2

3

4

5

The corresponding eigenvalues are

= ;1=4 = ;3=4 ) = ;1 (4) Defective, linked to (5) Jordan block The numerical solution written in full is thus ~n ; ~1 = c (1 ; h )n~x 4 3 + c (1 ; 4h )n~x 1 2 3

3

1

2

1

2

(9.66)

CHAPTER 9. RELAXATION METHODS " # n n ( n ; 1) n n ; n ; + c (1 ; h) + c h 1! (1 ; h) + c h 2! (1 ; h) ~x + c (1 ; h)n + c h n (1 ; h)n; ~x

180

3

4

4

5

+ c (1 ; h)n~x 5

1

5

1

1!

2

2

3

4

(9.67)

5

9.4.3 The SOR System

If = 1 and 2=! = x in Eq. 9.44, the ODE matrix is B ; (;1; x; 0)B (1; ;2; 1). One can show that this can be written in the form given below for M = 5. The generalization to any M is fairly clear. The H ; A matrix for the SOR method, ASOR B ; (;1; x; 0)B (1; ;2; 1), is 1

1

1

2 x 0 0 66 ;2;x2x+ x x 0 1 66 ;2x + x x x; 2;x2x+ x x ; 2 x x x 664 ;2x + x x ; 2x + x x ; 2x + x x ; 2x ;2 + x 1 ; 2x + x x ; 2x + x x ; 2x + x x 4

5

4

3

4

2

3

2

3

2

4

4

3

2

4

3

3

4

2

3

2

2

4

4

3

3

2

4

3

4

0 0 0 x

4

3

; 2x

4

3 77 77 77 (9.68) 5

Eigenvalues of the system are given by

m = ;1 +

!pm + zm 2

2

; m = 1; 2; : : : M

(9.69)

where

zm = [4(1 ; !) + ! pm] = pm = cos[m=(M + 1)] 2 2

1 2

If ! = 1, the system is Gauss-Seidel. If 4(1 ; !) + ! pm < 0; zm and m are complex. If ! is chosen such that 4(1 ; !) + ! p = 0; ! is optimum for the stationary case, and the following conditions hold: 2

2 2 1

1. Two eigenvalues are real, equal and defective. 2. If M is even, the remaining eigenvalues are complex and occur in conjugate pairs. 3. If M is odd, one of the remaining eigenvalues is real and the others are complex occurring in conjugate pairs.

9.4. EIGENSYSTEMS OF THE CLASSICAL METHODS

181

One can easily show that the optimum ! for the stationary case is !opt = 2= 1 + sin M + 1 and for ! = !opt m = m ; 1 ~xm = mj; sin j m M +1

(9.70)

2

(9.71)

1

where

q ! opt m = 2 pm + i p ; pm Using the explicit Euler method to integrate the ODE's, m = 1 ; h + hm , and if h = 1, the optimum value for the stationary case, the - relation reduces to that shown in Fig. 9.5. This illustrates the fact that for optimum stationary SOR all the jm j are identical and equal to !opt ; 1. Hence the convergence rate is jm jmax = !opt ; 1 (9.72) !opt = 2= 1 + sin M + 1 2 1

2

2

For M = 40, jmjmax = 0:8578. Hence the worst error component is reduced to less than 0.23 times its initial value in only 10 iterations, much faster than both GaussSeidel and Point-Jacobi. In practical applications, the optimum value of ! may have to be determined by trial and error, and the bene t may not be as great. For odd M , there are two real eigenvectors and one real principal vector. The remaining linearly independent eigenvectors are all complex. For M = 5 they can be written 2 3 2 3 2 p3(1 3 2 3 )=2 1 = 2 ; 6 1 p p 66 1=2 77 66 9 77 66 3(1 i 2)=6 77 66 0 77 6 7 6 7 6 7 ~x = 66 1=3 77 ; ~x = 66 16 77 ; ~x ; = 66 p 77 ; ~x = 666 1=3 777 (9.73) 0p 64 1=6 75 64 13 75 64 3(5 i 2)=54 75 64 0 75 p p 1=18 6 1=9 3(7 4i 2)=162 The corresponding eigenvalues are = ;2=3 (2) Defectiveplinked to = ;(10 ; 2p2i)=9 = ;(10 + 2 2i)=9 = ;4=3 (9.74) 1

2

34

5

1

1

3

4

5

CHAPTER 9. RELAXATION METHODS

182 1.0

h = 1.0 2 real defective λ m

σm

2 complex λ m

1 real λ m

Μ=5 0.0 -2.0

0.0

-1.0

λmh

Figure 9.5: The ; relation for optimum stationary SOR, M = 5, h = 1. The numerical solution written in full is

~n ; ~1 = + + + +

[c (1 ; 2h=3)n + c nh(1 ; 2h=3)n; ]~x c (1 ; 2h=3)n~xp c [1 ; (10 ; 2p2i)h=9]n~x c [1 ; (10 + 2 2i)h=9]n~x c (1 ; 4h=3)n~x 1

2

1

2

2

3

3

4

4

5

1

5

(9.75)

9.5 Nonstationary Processes In classical terminology a method is said to be nonstationary if the conditioning matrices, H and C , are varied at each time step. This does not change the steadystate solution A;b ~f b , but it can greatly aect the convergence rate. In our ODE approach this could also be considered and would lead to a study of equations with nonconstant coecients. It is much simpler, however, to study the case of xed H and C but variable step size, h. This process changes the Point-Jacobi method to Richardson's method in standard terminology. For the Gauss-Seidel and SOR methods it leads to processes that can be superior to the stationary methods. The nonstationary form of Eq. 9.39 is 1

9.5. NONSTATIONARY PROCESSES

183

N N ~N = c ~x Y (1 + hn ) + + cm~xm Y (1 + mhn ) 1

1

n=1

+ + cM ~xM

1

N Y n=1

n=1

(1 + M hn ) + ~1

(9.76)

where the symbol stands for product. Since hn can now be changed at each step, the error term can theoretically be completely eliminated in M steps by taking hm = ;1=m , for m = 1; 2; ; M . However, the eigenvalues m are generally unknown and costly to compute. It is therefore unnecessary and impractical to set hm = ;1=m for m = 1; 2; : : : ; M . We will see that a few well chosen h's can reduce whole clusters of eigenvectors associated with nearby 's in the m spectrum. This leads to the concept of selectively annihilating clusters of eigenvectors from the error terms as part of a total iteration process. This is the basis for the multigrid methods discussed in Chapter 10. Let us consider the very important case when all of the m are real and negative (remember that they arise from a conditioned matrix so this constraint is not unrealistic for quite practical cases). Consider one of the error terms taken from M N ~eN ~N ; ~1 = X cm~xm Y (1 + mhn) m=1

n=1

(9.77)

and write it in the form

cm~xm Pe(m ) cm~xm

N Y n=1

(1 + m hn)

(9.78)

where Pe signi es an \Euler" polynomial. Now focus attention on the polynomial (Pe)N () = (1 + h )(1 + h ) (1 + hN ) 1

2

(9.79)

treating it as a continuous function of the independent variable . In the annihilation process mentioned after Eq. 9.76, we considered making the error exactly zero by taking advantage of some knowledge about the discrete values of m for a particular case. Now we pose a less demanding problem. Let us choose the hn so that the maximum value of (Pe)N () is as small as possible for all lying between a and b such that b a 0. Mathematically stated, we seek max j(Pe)N ()j = minimum ; with(Pe)N (0) = 1

b a

(9.80)

CHAPTER 9. RELAXATION METHODS

184

This problem has a well known solution due to Markov. It is

TN 2 ; ;a ; b a b ! (Pe)N () = TN ;a;;b a b

!

(9.81)

where

TN (y) = cos(N arccos y)

(9.82)

are the Chebyshev polynomials along the interval ;1 y 1 and

q N 1 q N 1 TN (y) = 2 y + y ; 1 + 2 y ; y ; 1 (9.83) are the Chebyshev polynomials for jyj > 1. In relaxation terminology this is generally referred to as Richardson's method, and it leads to the nonstationary step size choice given by 2

(

2

"

#)

1 = 1 ; ; + ( ; ) cos (2n ; 1) ; n = 1; 2; : : : N b a hn 2 b a 2N

(9.84)

Remember that all are negative real numbers representing the magnitudes of m in an eigenvalue spectrum. The error in the relaxation process represented by Eq. 9.76 is expressed in terms Q ~ of a set of eigenvectors, xm , ampli ed by the coecients cm (1 + m hn). With each eigenvector there is a corresponding eigenvalue. Eq. 9.84 gives us the best choice of a series of hn that will minimize the amplitude of the error carried in the eigenvectors associated with the eigenvalues between b and a . As an example for the use of Eq. 9.84, let us consider the following problem: Minimize the maximum error associated with the eigenvalues in the interval ;2 ;1 using only 3 iterations. The three values of h which satisfy this problem are

"

hn = 2= 3 ; cos (2n ;6 1)

#!

(9.85)

(9.86)

9.5. NONSTATIONARY PROCESSES

185

and the amplitude of the eigenvector is reduced to (Pe) () = T (2 + 3)=T (3)

(9.87)

n p p o T (3) = [3 + 8] + [3 ; 8] =2 99

(9.88)

3

where

3

3

3

3

3

A plot of Eq. 9.87 is given in Fig. 9.6 and we see that the amplitudes of all the eigenvectors associated with the eigenvalues in the range ;2 ;1 have been reduced to less than about 1% of their initial values. The values of h used in Fig. 9.6 are

p

h = 4=(6 ; 3) h = 4=(6 ; p 0) h = 4=(6 + 3) 1

2 3

Return now to Eq. 9.76. This was derived from Eq. 9.37 on the condition that the explicit Euler method, Eq. 9.41, was used to integrate the basic ODE's. If instead the implicit trapezoidal rule

n = n + 21 h(0n + 0n)

(9.89)

0 1 1 M N 1 + hn m C ~N = X cm~xm Y B CA + ~1 B@ 2 1 m n 1 ; 2 hnm

(9.90)

+1

+1

is used, the nonstationary formula

=1

=1

would result. This calls for a study of the rational \trapezoidal" polynomial, Pt: 0 1 1 N Y B 1 + 2 hn CC (Pt )N () = B (9.91) @ 1 A n 1 ; 2 hn under the same constraints as before, namely that =1

max j(Pt)N ()j = minimum , with (Pt)N (0) = 1

b a

(9.92)

CHAPTER 9. RELAXATION METHODS

186

1 0.9 0.8 0.7

(Pe )3 (λ)

0.6 0.5 0.4 0.3 0.2 0.1 0 −2

−1.8

−1.6

−1.4

−1.2

−1

−0.8

−0.6

−0.4

−0.2

0

−1

−0.8

−0.6

−0.4

−0.2

0

λ

0.02

0.015

0.01

e 3

(P ) (λ)

0.005

0

−0.005

−0.01

−0.015

−0.02 −2

−1.8

−1.6

−1.4

−1.2

λ

Figure 9.6: Richardson's method for 3 steps, minimization over ;2 ;1.

9.6. PROBLEMS

187

The optimum values of h can also be found for this problem, but we settle here for the approximation suggested by Wachspress

!

2 = ; a n; = N ; ; n = 1; 2; ; N (9.93) b hn b This process is also applied to problem 9.85. The results for (Pt) () are shown in Fig. 9.7. The error amplitude is about 1/5 of that found for (Pe) () in the same interval of . The values of h used in Fig. 9.7 are h = 1p h = 2 h = 2 (

1) (

1)

3

3

1 2

3

9.6 Problems 1. Given a relaxation method in the form H ~n = A~n ; f~ show that

~n = Gn~ + (I ; Gn)A; f~ 1

0

where G = I + H ; A. 2. For a linear system of the form (A + A )x = b, consider the iterative method 1

1

2

(I + A )~x = (I ; A )xn + b (I + A )xn = (I ; A )~x + b where is a parameter. Show that this iterative method can be written in the form H (xk ; xk ) = (A + A )xk ; b 1

2

+1

2

+1

1

1

2

Determine the iteration matrix G if = ;1=2. 3. Using Appendix B.2, nd the eigenvalues of H ; A for the SOR method with A = B (4 : 1; ;2; 1) and ! = !opt. (You do not have to nd H ; . Recall that the eigenvalues of H ; A satisfy A~xm = mH~xm .) Find the numerical values, not just the expressions. Then nd the corresponding jm j values. 1

1

1

CHAPTER 9. RELAXATION METHODS

188

1 0.9 0.8 0.7

t 3

(P ) (λ)

0.6 0.5 0.4 0.3 0.2 0.1 0 −2

−1.8

−1.6

−1.4

−1.2

−1

λ

−0.8

−0.6

−0.4

−0.2

0

−0.8

−0.6

−0.4

−0.2

0

−3

2

x 10

1.5

1

(Pt )3 (λ)

0.5

0

−0.5

−1

−1.5

−2 −2

−1.8

−1.6

−1.4

−1.2

−1

λ

Figure 9.7: Wachspress method for 3 steps, minimization over ;2 ;1.

9.6. PROBLEMS

189

4. Solve the following equation on the domain 0 x 1 with boundary conditions u(0) = 0, u(1) = 1:

@ u ; 6x = 0 @x For the initial condition, use u(x) = 0. Use second-order centered dierences on a grid with 40 cells (M = 39). Iterate to steady state using (a) the point-Jacobi method, (b) the Gauss-Seidel method, (c) the SOR method with the optimum value of !, and (d) the 3-step Richardson method derived in Section 9.5. Plot the solution after the residual is reduced by 2, 3, and 4 orders of magnitude. Plot the logarithm of the L -norm of the residual vs. the number of iterations. Determine the asymptotic convergence rate. Compare with the theoretical asymptotic convergence rate. 2

2

2

190

CHAPTER 9. RELAXATION METHODS

Chapter 10 MULTIGRID The idea of systematically using sets of coarser grids to accelerate the convergence of iterative schemes that arise from the numerical solution to partial dierential equations was made popular by the work of Brandt. There are many variations of the process and many viewpoints of the underlying theory. The viewpoint presented here is a natural extension of the concepts discussed in Chapter 9.

10.1 Motivation 10.1.1 Eigenvector and Eigenvalue Identi cation with Space Frequencies Consider the eigensystem of the model matrix B (1; ;2; 1). The eigenvalues and

eigenvectors are given in Sections 4.3.2 and 4.3.3, respectively. Notice that as the magnitudes of the eigenvalues increase, the space-frequency (or wavenumber) of the corresponding eigenvectors also increase. That is, if the eigenvalues are ordered such that

j j j j jM j 1

(10.1)

2

then the corresponding eigenvectors are ordered from low to high space frequencies. This has a rational explanation from the origin of the banded matrix. Note that @ sin(mx) = ;m sin(mx) (10.2) @x and recall that (10.3) xx~ = 1x B (1; ;2; 1)~ = X 1x D(~) X ; ~ 191 2

2

2

2

2

1

CHAPTER 10. MULTIGRID

192

where D(~) is a diagonal matrix containing the eigenvalues. We have seen that X ; ~ represents a sine transform, and X ~, a sine synthesis. Therefore, the operation 1x D(~) represents the numerical approximation of the multiplication of the appropriate sine wave by the negative square of its wavenumber, ;m . One nds that 1 = M + 1 ;2 + 2 cos m ;m ; m (11.29) < i i R i 0 where the subscript i indicates individual components of w and g. This is achieved with (^gi)j+1=2 = 21 i [(wi)L + (wi)R ] + 21 jij [(wi)L ; (wi)R] (11.30) or g^j+1=2 = 12 (wL + wR) + 21 jj (wL ; wR ) (11.31) Now, as in Eq. 11.19, we premultiply by X to return to the original variables and insert the product X ;1X after and jj to obtain X g^j+1=2 = 21 X X ;1X (wL + wR) + 21 X jjX ;1X (wL ; wR ) (11.32) and thus f^j+1=2 = 12 (fL + fR ) + 21 jAj (uL ; uR) (11.33) where

jAj = X jjX ;1

(11.34)

11.5. ARTIFICIAL DISSIPATION

213

and we have also used the relations f = Xg, u = Xw, and A = X X ;1. In the linear, constant-coecient case, this leads to an upwind operator which is identical to that obtained using ux-vector splitting. However, in the nonlinear case, there is some ambiguity regarding the de nition of jAj at the cell interface j + 1=2. In order to resolve this, consider a situation in which the eigenvalues of A are all of the same sign. In this case, we would like our de nition of f^j+1=2 to satisfy

(

i0s > 0 f^j+1=2 = ffL ifif all (11.35) 0 R all i s < 0 giving pure upwinding. If the eigenvalues of A are all positive, jAj = A; if they are all negative, jAj = ;A. Hence satisfaction of Eq. 11.35 is obtained by the de nition f^j+1=2 = 21 (fL + fR ) + 12 jAj+1=2j (uL ; uR) (11.36) if Aj+1=2 satis es fL ; fR = Aj+1=2 (uL ; uR) (11.37) For the Euler equations for a perfect gas, Eq. 11.37 is satis ed by the ux Jacobian evaluated at the Roe-average state given by pLuL + pR uR (11.38) uj+1=2 = p + p L R p H + p H L L Hj+1=2 = (11.39) pL + pRR R where u and H = (e + p)= are the velocity and the total enthalpy per unit mass, respectively.3

11.5 Arti cial Dissipation We have seen that numerical dissipation can be introduced by using one-sided differencing schemes together with some form of ux splitting. We have also seen that such dissipation can be introduced by adding a symmetric component to an antisymmetric (dissipation-free) operator. Thus we can generalize the concept of upwinding to include any scheme in which the symmetric portion of the operator is treated in such a manner as to be truly dissipative. 3 Note that the ux Jacobian can be written in terms of u and H only; see problem 6 at the end

of this chapter.

214

CHAPTER 11. NUMERICAL DISSIPATION

For example, let

(11.40) (xa u)j = uj+1 ; uj;1 ; (xs u)j = ;uj+1 + 2uj ; uj;1 2x 2x Applying x = xa + xs to the spatial derivative in Eq. 11.15 is stable if i 0 and unstable if i < 0. Similarly, applying x = xa ; xs is stable if i 0 and unstable if i > 0. The appropriate implementation is thus ix = ixa + jijxs (11.41) Extension to a hyperbolic system by applying the above approach to the characteristic variables, as in the previous two sections, gives x (Au) = xa (Au) + xs (jAju) (11.42) or xf = xa f + xs (jAju) (11.43) where jAj is de ned in Eq. 11.34. The second spatial term is known as arti cial dissipation. It is also sometimes referred to as arti cial diusion or arti cial viscosity. With appropriate choices of xa and xs , this approach can be related to the upwind approach. This is particularly evident from a comparison of Eqs. 11.36 and 11.43. It is common to use the following operator for xs (xs u)j = (uj;2 ; 4uj;1 + 6uj ; 4uj+1 + uj+2) (11.44) x where is a problem-dependent coecient. This symmetric operator approximates x3 uxxxx and thus introduces a third-order dissipative term. With an appropriate value of , this often provides sucient damping of high frequency modes without greatly aecting the low frequency modes. For details of how this can be implemented for nonlinear hyperbolic systems, the reader should consult the literature. A more complicated treatment of the numerical dissipation is also required near shock waves and other discontinuities, but is beyond the scope of this book.

11.6 Problems 1. A second-order backward dierence approximation to a 1st derivative is given as a point operator by (xu)j = 1 (uj;2 ; 4uj;1 + 3uj ) 2x

11.6. PROBLEMS

215

(a) Express this operator in banded matrix form (for periodic boundary conditions), then derive the symmetric and skew-symmetric matrices that have the matrix operator as their sum. (See Appendix A.3 to see how to construct the symmetric and skew-symmetric components of a matrix.) (b) Using a Taylor table, nd the derivative which is approximated by the corresponding symmetric and skew-symmetric operators and the leading error term for each. 2. Find the modi ed wavenumber for the rst-order backward dierence operator. Plot the real and imaginary parts of x vs. x for 0 x . Using Fourier analysis as in Section 6.6.2, nd jj for the combination of this spatial operator with 4th-order Runge-Kutta time marching at a Courant number of unity and plot vs. x for 0 x . 3. Find the modi ed wavenumber for the operator given in Eq. 11.6. Plot the real and imaginary parts of x vs. x for 0 x . Using Fourier analysis as in Section 6.6.2, nd jj for the combination of this spatial operator with 4th-order Runge-Kutta time marching at a Courant number of unity and plot vs. x for 0 x . 4. Consider the spatial operator obtained by combining second-order centered differences with the symmetric operator given in Eq. 11.44. Find the modi ed wavenumber for this operator with = 0; 1=12; 1=24, and 1=48. Plot the real and imaginary parts of x vs. x for 0 x . Using Fourier analysis as in Section 6.6.2, nd jj for the combination of this spatial operator with 4th-order Runge-Kutta time marching at a Courant number of unity and plot vs. x for 0 x . 5. Consider the hyperbolic system derived in problem 8 of Chapter 2. Find the matrix jAj. Form the plus-minus split ux vectors as in Section 11.4.1. 6. Show that the ux Jacobian for the 1-D Euler equations can be written in terms of u and H . Show that the use of the Roe average state given in Eqs. 11.38 and 11.39 leads to satisfaction of Eq. 11.37.

216

CHAPTER 11. NUMERICAL DISSIPATION

Chapter 12 SPLIT AND FACTORED FORMS In the next two chapters, we present and analyze split and factored algorithms. This gives the reader a feel for some of the modi cations which can be made to the basic algorithms in order to obtain ecient solvers for practical multidimensional applications, and a means for analyzing such modi ed forms.

12.1 The Concept Factored forms of numerical operators are used extensively in constructing and applying numerical methods to problems in uid mechanics. They are the basis for a wide variety of methods variously known by the labels \hybrid", \time split", and \fractional step". Factored forms are especially useful for the derivation of practical algorithms that use implicit methods. When we approach numerical analysis in the light of matrix derivative operators, the concept of factoring is quite simple to present and grasp. Let us start with the following observations: 1. Matrices can be split in quite arbitrary ways. 2. Advancing to the next time level always requires some reference to a previous one. 3. Time marching methods are valid only to some order of accuracy in the step size, h. Now recall the generic ODE's produced by the semi-discrete approach

d~u = A~u ; ~f dt 217

(12.1)

CHAPTER 12. SPLIT AND FACTORED FORMS

218

and consider the above observations. From observation 1 (arbitrary splitting of A) :

d~u = [A + A ]~u ; ~f dt 1

(12.2)

2

where A = [A + A ] but A and A are not unique. For the time march let us choose the simple, rst-order, explicit Euler method. Then, from observation 2 (new data ~un in terms of old ~un): 1

2

1

1

2

+1

~un = [ I + hA + hA ]~un ; h~f + O(h ) +1

1

(12.3)

2

2

or its equivalent

~un = h[ I + hA ][ I + hA ] ; h A A i~un ; h~f + O(h ) +1

1

2

2

1

2

2

Finally, from observation 3 (allowing us to drop higher order terms ;h A A ~un): 2

1

2

~un = [ I + hA ][ I + hA ]~un ; h~f + O(h )

(12.4) Notice that Eqs. 12.3 and 12.4 have the same formal order of accuracy and, in this sense, neither one is to be preferred over the other. However, their numerical stability can be quite dierent, and techniques to carry out their numerical evaluation can have arithmetic operation counts that vary by orders of magnitude. Both of these considerations are investigated later. Here we seek only to apply to some simple cases the concept of factoring. +1

1

2

2

12.2 Factoring Physical Representations | Time Splitting Suppose we have a PDE that represents both the processes of convection and dissipation. The semi-discrete approach to its solution might be put in the form

d~u = A ~u + A ~u + (bc~ ) c d dt

(12.5)

where Ac and Ad are matrices representing the convection and dissipation terms, respectively; and their sum forms the A matrix we have considered in the previous sections. Choose again the explicit Euler time march so that

~ ) + O (h ) ~un = [ I + hAd + hAc]~un + h(bc +1

1 Second-order time-marching methods are considered later.

2

(12.6)

12.2. FACTORING PHYSICAL REPRESENTATIONS | TIME SPLITTING 219 Now consider the factored form

~un

+1

~) = [ I + hAd ] [ I + hAc]~un + h(bc ~ ) +O(h ) (12.7) ~ ) + h Ad Ac~un + (bc = [| I + hAd + h{zAc]~un + h(bc } | {z } Original Unfactored Terms Higher Order Terms 2

2

and we see that Eq. 12.7 and the original unfactored form Eq. 12.6 have identical orders of accuracy in the time approximation. Therefore, on this basis, their selection is arbitrary. In practical applications equations such as 12.7 are often applied in a predictor-corrector sequence. In this case one could write 2

u~n ~un

+1 +1

~) = [ I + hAc]~un + h(bc = [ I + hAd]~un

(12.8)

+1

Factoring can also be useful to form split combinations of implicit and explicit techniques. For example, another way to approximate Eq. 12.6 with the same order of accuracy is given by the expression

~) = [ I ; hAd]; [ I + hAc]~un + h(bc ~ ) +O(h ) = [| I + hAd + h{zAc]~un + h(bc } Original Unfactored Terms where in this approximation we have used the fact that ~un

1

+1

2

(12.9)

[ I ; hAd ]; = I + hAd + h Ad + 1

2

2

if h jjAdjj < 1, where jjAd jj is some norm of [Ad ]. This time a predictor-corrector interpretation leads to the sequence

u~n [ I ; hAd]~un

+1 +1

~) = [ I + hAc]~un + h(bc = u~n +1

(12.10)

The convection operator is applied explicitly, as before, but the diusion operator is now implicit, requiring a tridiagonal solver if the diusion term is central dierenced. Since numerical stiness is generally much more severe for the diusion process, this factored form would appear to be superior to that provided by Eq. 12.8. However, the important aspect of stability has yet to be discussed. 2 We do not suggest that this particular method is suitable for use. We have yet to determine its

stability, and a rst-order time-march method is usually unsatisfactory.

CHAPTER 12. SPLIT AND FACTORED FORMS

220

We should mention here that Eq. 12.9 can be derived for a dierent point of view by writing Eq. 12.6 in the form

un ; un = A u + A u + (bc ~ ) + O(h ) c n d n h +1

Then

2

+1

~) [ I ; hAd ]un = [ I + hAc]un + h(bc which is identical to Eq. 12.10. +1

12.3 Factoring Space Matrix Operators in 2{D

12.3.1 Mesh Indexing Convention

Factoring is widely used in codes designed for the numerical solution of equations governing unsteady two- and three-dimensional ows. Let us study the basic concept of factoring by inspecting its use on the linear 2-D scalar PDE that models diusion:

@u = @ u + @ u @t @x @y 2

2

2

(12.11)

2

We begin by reducing this PDE to a coupled set of ODE's by dierencing the space derivatives and inspecting the resulting matrix operator. A clear description of a matrix nite-dierence operator in 2- and 3-D requires some reference to a mesh. We choose the 3 4 point mesh shown in the Sketch 12.12. In this example Mx , the number of (interior) x points, is 4 and My , the number of (interior) y points is 3. The numbers 11 ; 12 ; ; 43 represent the location in the mesh of the dependent variable bearing that index. Thus u represents the value of u at j = 3 and k = 2. 3

32

My 13 23 33 43 k 12 22 32 42 1

11 21 31 41 1 j Mx

(12.12)

Mesh indexing in 2-D. 3 This could also be called a 5 6 point mesh if the boundary points (labeled

in the sketch) were included, but in these notes we describe the size of a mesh by the number of interior points.

12.3. FACTORING SPACE MATRIX OPERATORS IN 2{D

221

12.3.2 Data Bases and Space Vectors

The dimensioned array in a computer code that allots the storage locations of the dependent variable(s) is referred to as a data-base. There are many ways to lay out a data-base. Of these, we consider only two: (1), consecutively along rows that are themselves consecutive from k = 1 to My , and (2), consecutively along columns that are consecutive from j = 1 to Mx. We refer to each row or column group as a space vector (they represent data along lines that are continuous in space) and label their sum with the symbol U . In particular, (1) and (2) above are referred to as xvectors and y-vectors, respectively. The symbol U by itself is not enough to identify the structure of the data-base and is used only when the structure is immaterial or understood. To be speci c about the structure, we label a data{base composed of x-vectors with U x , and one composed of y-vectors with U y . Examples of the order of indexing for these space vectors are given in Eq. 12.16 part a and b. ( )

( )

12.3.3 Data Base Permutations

The two vectors (arrays) are related by a permutation matrix P such that

U x = Pxy U y ( )

( )

U y = PyxU x

and

( )

where

(12.13)

( )

Pyx = PTxy = P;xy Now consider the structure of a matrix nite-dierence operator representing 3point central-dierencing schemes for both space derivatives in two dimensions. When the matrix is multiplying a space vector U , the usual (but ambiguous) representation is given by Ax y . In this notation the ODE form of Eq. 12.11 can be written 1

4

+

dU = A U + (bc ~) x y dt

(12.14)

+

If it is important to be speci c about the data-base structure, we use the notation Axx y or Axy y , depending on the data{base chosen for the U it multiplies. Examples are in Eq. 12.16 part a and b. Notice that the matrices are not the same although they represent the same derivative operation. Their structures are similar, however, and they are related by the same permutation matrix that relates U x to U y . Thus ( ) +

( ) +

( )

Axx y = Pxy Axy y Pyx ( ) +

4 Notice that

( ) +

( )

(12.15)

Ax+y and U , which are notations used in the special case of space vectors, are subsets of A and ~u, used in the previous sections.

CHAPTER 12. SPLIT AND FACTORED FORMS

222 2 6 x 6 6 6 6 6 6 6 6 6 6 o 6 6 A(xx+)y U (x) = 666 6 6 6 6 6 6 6 6 6 6 6 4

j o x j o x x j o x j o

j j j j

j x o j x x o j x x o j x

j o j o j o j

x

j o j o j o j o

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 o 777 7 7 7 7 7 7 x 75

j x j x x j x j x

11 21 31 41

;;

;;

a:Elements in 2-dimensional, central-dierence, matrix operator, Ax y , for 34 mesh shown in Sketch 12.12. Data base composed of My x{vectors stored in U x . Entries for x ! x, for y ! o, for both ! . +

2 6 o 6 6 6 6 6 6 6 x 6 6 6 6 6 6 (y ) (y ) Ax+y U = 666 6 6 6 6 6 6 6 6 6 6 6 6 4

13 23 33 43

( )

j x j o j x j o j x j o

j j j

j o j x j x j o o j x j x j o j x j j x j o j x j x j o o j x j x j o j j j j

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 x 777 7 7 7 7 7 o5

j x j o j x j o j x j o

b: Elements in 2-dimensional, central-dierence, matrix operator, Ax y , for 34 mesh shown in Sketch 12.12. Data base composed of Mx y{vectors stored in U y . Entries for x ! x, for y ! o, for both ! . +

12 22 32 42

( )

11 12 13

;; 21 22 23

;; 31 32 33

;; 41 42 43

(12.16)

12.3. FACTORING SPACE MATRIX OPERATORS IN 2{D

223

12.3.4 Space Splitting and Factoring We are now prepared to discuss splitting in two dimensions. It should be clear that the matrix Axx y can be split into two matrices such that ( ) +

Axx y = Axx + Ayx ( ) +

( )

(12.17)

( )

where Axx and Ayx are shown in Eq. 12.22. Similarily ( )

( )

Axy y = Axy + Ayy ( ) +

( )

(12.18)

( )

where the split matrices are shown in Eq. 12.23. The permutation relation also holds for the split matrices so

Ayx = Pxy Ayy Pyx ( )

and

( )

Axx = Pxy Axy Pyx ( )

( )

The splittings in Eqs. 12.17 and 12.18 can be combined with factoring in the manner described in Section 12.2. As an example ( rst-order in time), applying the implicit Euler method to Eq. 12.14 gives h i Unx = Unx + h Axx + Ayx Unx + h(bc~ ) ( ) +1

( )

( )

( )

( ) +1

or h

i I ; hAxx ; hAyx Unx = Unx + h(bc~ ) + O(h ) ( )

( ) +1

( )

( )

2

(12.19)

As in Section 12.2, we retain the same rst order accuracy with the alternative h

ih i ~ ) + O(h ) I ; hAxx I ; hAyx Unx = Unx + h(bc ( )

( ) +1

( )

( )

2

(12.20)

Write this in predictor-corrector form and permute the data base of the second row. There results h

i

~) I ; hAxx U~ x = Unx + h(bc h i y I ; hAyy Un = U~ y ( )

( )

( )

( )

( ) +1

( )

(12.21)

CHAPTER 12. SPLIT AND FACTORED FORMS

224

2 x 6 x 6 6 6 6 6 6 6 6 6 6 6 6 (x) (x) Ax U = 666 6 6 6 6 6 6 6 6 6 6 6 4

x x x x x x x x

j j j j

j j j j

j x x j x x x j x x x j x x

j j j j

j j j j 2 o j o 6 o j o 6 6 6 o j o 6 6 o j o 6 6

6 6 6 o 6 6 A(yx) U (x) = 666 6 6 6 6 6 6 6 6 6 6 6 4

j o o j o o j o o j o j o j o j o j o The splitting of Axx y . ( ) +

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 U (x) 7 7 7 7 7 7 7 7 7 7 7 x 75

j x x j x x x j x x j x x 3 j 7 j 7 7 7 j 7 7 j 7 7

(12.22)

7

7 7 j o 7 7 j o 7 x j o 777 U j o 777

( )

7

7 7 j o 7 7 j o 7 j o 75 j o

12.3. FACTORING SPACE MATRIX OPERATORS IN 2{D

2 x 6 6 6 6 6 6 6 6 x 6 6 6 6 6 6 A(xy) U (y) = 666 6 6 6 6 6 6 6 6 6 6 6 6 4

j x j x j x j x j x j

j x j x j x j x j x j x j x j x j j x j x j x j x j x j x j x j x j

j j j 2 o o j 6 o o o j 6 6 o o j 6 6

6 6 6 6 6 6 6 6 6 (y ) (y ) Ay U = 666 6 6 6 6 6 6 6 6 6 6 6 6 4

j j j

225

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 U (y) 7 7 7 7 7 x 777 7 7 7 7 7 5

j x j x j x j x j x j x 3 j j 7 j j 7 7 j j 7 7

j o o j j o o o j j o o j

j j j

j j j

j o o j j o o o j j o o j

j j j

j j j

7 7 7 7 7 7 7 7 7 7 7 U (y) 7 7 7 7 7 7 7 7 7 7 7 7 7 o5

j o o j o o j o o

The splitting of Axy y . ( ) +

(12.23)

CHAPTER 12. SPLIT AND FACTORED FORMS

226

12.4 Second-Order Factored Implicit Methods

Second-order accuracy in time can be maintained in a certain factored implicit methods. For example, apply the trapezoidal method to Eq. 12.14 where the derivative operators have been split as in Eq. 12.17 or 12.18. Let the data base be immaterial ~ ) be time invariant. There results and the (bc 1 1 1 1 ~ ) + O(h ) (12.24) I ; 2 hAx ; 2 hAy Un = I + 2 hAx + 2 hAy Un + h(bc Factor both sides giving 1 1 1 I ; 2 hAx I ; 2 hAy ; 4 h AxAy Un 1 1 1 ~ ) + O(h ) (12.25) = I + hAx I + hAy ; h AxAy Un + h(bc 2 2 4 Then notice that the combination 14 h [AxAy ](Un ; Un) is proportional to h since the leading term in the expansion of (Un ; Un) is proportional to h. Therefore, we can write 1 1 1 1 I ; 2 hAx I ; 2 hAy Un = I + 2 hAx I + 2 hAy Un + h(bc~ ) + O(h )(12.26) and both the factored and unfactored form of the trapezoidal method are second-order accurate in the time march. An alternative form of this kind of factorization is the classical ADI (alternating direction implicit) method usually written 1 1 ~ I ; 2 hAx U = I + 2 hAy Un + 12 hFn 1 1 (12.27) I ; 2 hAy Un = I + 2 hAx U~ + 12 hFn + O(h ) 3

+1

2

+1

2

2

3

3

+1

+1

3

+1

5

+1

+1

3

For idealized commuting systems the methods given by Eqs. 12.26 and 12.27 dier only in their evaluation of a time-dependent forcing term.

12.5 Importance of Factored Forms in 2 and 3 Dimensions When the time-march equations are sti and implicit methods are required to permit reasonably large time steps, the use of factored forms becomes a very valuable tool 5 A form of the Douglas or Peaceman-Rachford methods.

12.5. IMPORTANCE OF FACTORED FORMS IN 2 AND 3 DIMENSIONS 227 for realistic problems. Consider, for example, the problem of computing the time advance in the unfactored form of the trapezoidal method given by Eq. 12.24 1 1 ~) I ; 2 hAx y Un = I + 2 hAx y Un + h(bc Forming the right hand side poses no problem, but nding Un requires the solution of a sparse, but very large, set of coupled simultaneous equations having the matrix form shown in Eq. 12.16 part a and b. Furthermore, in real cases involving the Euler or Navier-Stokes equations, each symbol (o; x; ) represents a 4 4 block matrix with entries that depend on the pressure, density and velocity eld. Suppose we were to solve the equations directly. The forward sweep of a simple Gaussian elimination lls all of the 4 4 blocks between the main and outermost diagonal (e.g. between and o in Eq. 12.16 part b.). This must be stored in computer memory to be used to nd the nal solution in the backward sweep. If Ne represents the order of the small block matrix (4 in the 2-D Euler case), the approximate memory requirement is (Ne My ) (Ne My ) Mx

oating point words. Here it is assumed that My < Mx. If My > Mx, My and Mx would be interchanged. A moderate mesh of 60 200 points would require over 11 million words to nd the solution. Actually current computer power is able to cope rather easily with storage requirements of this order of magnitude. With computing speeds of over one giga op, direct solvers may become useful for nding steady-state solutions of practical problems in two dimensions. However, a three-dimensional solver would require a memory of approximatly +

+1

+

+1

6

7

8

Ne My Mz Mx 2

2

2

words and, for well resolved ow elds, this probably exceeds memory availability for some time to come. On the other hand, consider computing a solution using the factored implicit equation 12.25. Again computing the right hand side poses no problem. Accumulate the result of such a computation in the array (RHS ). One can then write the remaining terms in the two-step predictor-corrector form I ; 21 hAxx U~ x = (RHS ) x I ; 12 hAyy Uny = U~ y (12.28) ( )

( )

( )

( ) +1

( )

( )

6 For matrices as small as those shown there are many gaps in this \ ll", but for meshes of

practical size the ll is mostly dense. 7 The lower band is also computed but does not have to be saved unless the solution is to be repeated for another vector. 8 One billion oating-point operations per second.

CHAPTER 12. SPLIT AND FACTORED FORMS

228

which has the same appearance as Eq. 12.21 but is second-order time accurate. The rst step would be solved using My uncoupled block tridiagonal solvers . Inspecting the top of Eq. 12.22, we see that this is equivalent to solving My one-dimensional problems, each with Mx blocks of order Ne . The temporary solution U~ x would then be permuted to U~ y and an inspection of the bottom of Eq. 12.23 shows that the nal step consists of solving Mx one-dimensional implicit problems each with dimension My . 9

( )

( )

12.6 The Delta Form Clearly many ways can be devised to split the matrices and generate factored forms. One way that is especially useful, for ensuring a correct steady-state solution in a converged time-march, is referred to as the \delta form" and we develop it next. Consider the unfactored form of the trapezoidal method given by Eq. 12.24, and ~ ) be time invariant: let the (bc 1 1 1 1 ~ ) + O(h ) I ; 2 hAx ; 2 hAy Un = I + 2 hAx + 2 hAy Un + h(bc From both sides subtract 1 1 I ; 2 hAx ; 2 hAy Un leaving the equality unchanged. Then, using the standard de nition of the dierence operator , Un = Un ; Un one nds h 1 1 ~ )i + O(h ) I ; 2 hAx ; 2 hAy Un = h Ax y Un + (bc (12.29) Notice that the right side of this equation is the product of h and a term that is identical to the right side of Eq. 12.14, our original ODE. Thus, if Eq. 12.29 converges, it is guaranteed to converge to the correct steady-state solution of the ODE. Now we can factor Eq. 12.29 and maintain O(h ) accuracy. We arrive at the expression h i 1 1 (12.30) I ; 2 hAx I ; 2 hAy Un = h Ax y Un + (bc~ ) + O(h ) This is the delta form of a factored, 2nd-order, 2-D equation. 3

+1

+1

3

+

2

+

3

9 A block tridiagonal solver is similar to a scalar solver except that small block matrix operations

replace the scalar ones, and matrix multiplications do not commute.

12.7. PROBLEMS

229

The point at which the factoring is made may not aect the order of time-accuracy, but it can have a profound eect on the stability and convergence properties of a method. For example, the unfactored form of a rst-order method derived from the implicit Euler time march is given by Eq. 12.19, and if it is immediately factored, the factored form is presented in Eq. 12.20. On the other hand, the delta form of the unfactored Eq. 12.19 is h

i

~) [I ; hAx ; hAy ]Un = h Ax y Un + (bc +

and its factored form becomes

10

h

~) [I ; hAx][I ; hAy ]Un = h Ax y Un + (bc

i

+

(12.31)

In spite of the similarities in derivation, we will see in the next chapter that the convergence properties of Eq. 12.20 and Eq. 12.31 are vastly dierent.

12.7 Problems 1. Consider the 1-D heat equation:

@u = @ u @t @x 2

0x9

2

Let u(0; t) = 0 and u(9; t) = 0, so that we can simplify the boundary conditions. Assume that second order central dierencing is used, i.e., (xxu)j = 1 (uj; ; 2uj + uj ) x The uniform grid has x = 1 and 8 interior points. 2

1

+1

(a) Space vector de nition i. What is the space vector for the natural ordering (monotonically increasing in index), u ? Only include the interior points. ii. If we reorder the points with the odd points rst and then the even points, write the space vector, u ? iii. Write down the permutation matrices,(P ,P ). (1)

(2)

12

21

10 Notice that the only dierence between the O(h2 ) method given by Eq. 12.30 and the O(h) method given by Eq. 12.31 is the appearance of the factor 1 on the left side of the O(h2 ) method.

2

CHAPTER 12. SPLIT AND FACTORED FORMS

230

iv. The generic ODE representing the discrete form of the heat equation is

du = A u + f dt Write down the matrix A . (Note f = 0, due to the boundary conditions) Next nd the matrix A such that du = A u dt Note that A can be written as 2 T3 U D 7 A = 64 5 U D De ne D and U . (1)

(1)

1

1

2

(2)

2

(2)

2

2

v. Applying implicit Euler time marching, write the delta form of the implicit algorithm. Comment on the form of the resulting implicit matrix operator. (b) System de nition In problem 1a, we de ned u ; u ; A ; A ; P , and P which partition the odd points from the even points. We can put such a partitioning to use. First de ne extraction operators 21 0 0 0 0 0 0 03 6 0 0 0 0 77 0 1 0 0 6 6 3 0 0 0 0 777 2 I 0 0 1 0 6 6 0 7 6 0 0 0 0 7 = 64 7 I o = 66 00 00 00 01 5 7 0 0 0 0 7 6 0 6 6 0 0 0 0 777 0 60 0 0 0 40 0 0 0 0 0 0 05 0 0 0 0 0 0 0 0 20 0 0 0 0 0 0 03 6 0 0 0 0 0 0 0 0 77 6 6 3 0 0 0 0 0 0 0 0 777 2 0 6 6 0 7 6 0 0 0 0 7 = 64 7 I e = 66 00 00 00 00 5 7 1 0 0 0 7 6 0 I 7 6 6 0 1 0 0 77 60 0 0 0 40 0 0 0 0 0 1 05 0 0 0 0 0 0 0 1 (1)

( )

( )

(2)

1

2

12

21

4

4

4

4

4

4

4

4

12.7. PROBLEMS

231

which extract the odd even points from u as follows: u o = I o u and ue =Ieu . i. Beginning with the ODE written in terms of u , de ne a splitting A = Ao + Ae, such that Ao operates only on the odd terms, and Ae operates only on the even terms. Write out the matrices Ao and Ae. Also, write them in terms of D and U de ned above. ii. Apply implicit Euler time marching to the split ODE. Write down the delta form of the algorithm and the factored delta form. Comment on the order of the error terms. iii. Examine the implicit operators for the factored delta form. Comment on their form. You should be able to argue that these are now trangular matrices (a lower and an upper). Comment on the solution process this gives us relative to the direct inversion of the original system. (2)

( )

( )

( )

(2)

(2)

2

( )

(2)

232

CHAPTER 12. SPLIT AND FACTORED FORMS

Chapter 13 LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS In Section 4.4 we introduced the concept of the representative equation, and used it in Chapter 7 to study the stability, accuracy, and convergence properties of timemarching schemes. The question is: Can we nd a similar equation that will allow us to evaluate the stability and convergence properties of split and factored schemes? The answer is yes | for certain forms of linear model equations. The analysis in this chapter is useful for estimating the stability and steady-state properties of a wide variety of time-marching schemes that are variously referred to as time-split, fractional-step, hybrid, and (approximately) factored. When these methods are applied to practical problems, the results found from this analysis are neither necessary nor sucient to guarantee stability. However, if the results indicate that a method has an instability, the method is probably not suitable for practical use.

13.1 The Representative Equation for Circulant Operators Consider linear PDE's with coecients that are xed in both space and time and with boundary conditions that are periodic. We have seen that under these conditions a semi-discrete approach can lead to circulant matrix dierence operators, and we discussed circulant eigensystems1 in Section 4.3. In this and the following section we assume circulant systems and our analysis depends critically on the fact that all circulant matrices commute and have a common set of eigenvectors. 1

See also the discussion on Fourier stability analysis in Section 7.7.

233

234 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS Suppose, as a result of space dierencing the PDE, we arrive at a set of ODE's that can be written d~u = A ~u + A ~u ; ~f (t) (13.1) ap bp dt where the subscript p denotes a circulant matrix. Since both matrices have the same set of eigenvectors, we can use the arguments made in Section 4.2.3 to uncouple the set and form the M set of independent equations

w10 = (a + b)1 w1 ; g1 (t) ... wm0 = (a + b)m wm ; gm(t) ... wM0 = (a + b)M wM ; gM (t)

(13.2)

The analytic solution of the m'th line is

wm(t) = cm e(a+b)mt + P:S: Note that each a pairs with one, and only one2 , b since they must share a common eigenvector. This suggests (see Section 4.4: The representative equation for split, circulant systems is du = [ + + + ]u + aet (13.3) a b c dt where a + b + c + are the sum of the eigenvalues in Aa ; Ab ; Ac ; that share the same eigenvector.

13.2 Example Analysis of Circulant Systems 13.2.1 Stability Comparisons of Time-Split Methods Consider as an example the linear convection-diusion equation: @u + a @u = @ 2 u @t @x @x2 2

(13.4)

This is to be contrasted to the developments found later in the analysis of 2-D equations.

13.2. EXAMPLE ANALYSIS OF CIRCULANT SYSTEMS

235

If the space dierencing takes the form d~u = ; a B (;1; 0; 1)~u + B (1; ;2; 1)~u (13.5) dt 2x p x2 p the convection matrix operator and the diusion matrix operator, can be represented by the eigenvalues c and d, respectively, where (see Section 4.3.2): (c)m = iax sin m (13.6) (d)m = ; 42 sin2 2m x In these equations m = 2m=M , m = 0 ; 1 ; ; M ; 1 , so that 0 m 2. Using these values and the representative equation 13.4, we can analyze the stability of the two forms of simple time-splitting discussed in Section 12.2. In this section we refer to these as 1. the explicit-implicit Euler method, Eq. 12.10. 2. the explicit-explicit Euler method, Eq. 12.8.

1. The Explicit-Implicit Method

When applied to Eq. 13.4, the characteristic polynomial of this method is P (E ) = (1 ; hd )E ; (1 + hc) This leads to the principal root 1 + i ahx sin m = 1 + 4 h2 sin2 m 2 x where we have made use of Eq. 13.6 to quantify the eigenvalues. Now introduce the dimensionless numbers Cn = ahx ; Courant number R = a x ; mesh Reynolds number and we can write for the absolute value of q 1 + Cn2 sin2 m jj = ; 0 m 2 (13.7) C m n 2 1 + 4 R sin 2

236 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS 1.2

1.2 C = 2/ R ∆ n

C = R /2 n ∆

0.8 C

C

n

0.4

0

C = 2/ R n ∆

0.8 n

0.4

0

4 R∆ Explicit-Implicit

8

0

0

4 R∆ Explicit-Explicit

8

Figure 13.1: S_ tability regions for two simple time-split methods. A simple numerical parametric study of Eq. 13.7 shows that the critical range of m for any combination of Cn and R occurs when m is near 0 (or 2). From this we nd that the condition on Cn and R that make jj 1 is 2 h i C n 2 2 2 1 + Cn sin = 1 + 4 R sin 2 As ! 0 this gives the stability region Cn < R2 which is bounded by a hyperbola and shown in Fig. 13.1.

2. The Explicit-Explicit Method An analysis similar to the one given above shows that this method produces " # q C n m 2 2 jj = 1 + Cn2 sin m 1 ; 4 R sin 2 ; 0 m 2 Again a simple numerical parametric study shows that this has two critical ranges of m , one near 0, which yields the same result as in the previous example, and the

13.2. EXAMPLE ANALYSIS OF CIRCULANT SYSTEMS

237

other near 180o, which produces the constraint that Cn < 21 R for R 2 The resulting stability boundary is also shown in Fig. 13.1. The totaly explicit, factored method has a much smaller region of stability when R is small, as we should have expected.

13.2.2 Analysis of a Second-Order Time-Split Method

Next let us analyze a more practical method that has been used in serious computational analysis of turbulent ows. This method applies to a ow in which there is a combination of diusion and periodic convection. The convection term is treated explicitly using the second-order Adams-Bashforth method. The diusion term is integrated implicitly using the trapezoidal method. Our model equation is again the linear convection-diusion equation 13.4 which we split in the fashion of Eq. 13.5. In order to evaluate the accuracy, as well as the stability, we include the forcing function in the representative equation and study the eect of our hybrid, time-marching method on the equation u0 = cu + du + aet First let us nd expressions for the two polynomials, P (E ) and Q(E ). The characteristic polynomial follows from the application of the method to the homogeneous equation, thus un+1 = un + 21 hc(3un ; un;1) + 21 hd (un+1 + un) This produces P (E ) = (1 ; 12 hd)E 2 ; (1 + 23 hc + 12 hd)E + 12 hc The form of the particular polynomial depends upon whether the forcing function is carried by the AB2 method or by the trapezoidal method. In the former case it is Q(E ) = 21 h(3E ; 1) (13.8) and in the latter (13.9) Q(E ) = 12 h(E 2 + E )

238 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS

Accuracy

From the characteristic polynomial we see that there are two {roots and they are given by the equation s

2 1 + 32 hc + 21 hd ; 2hc 1 ; 21 hd (13.10) 2 1 ; 12 hd The principal -root follows from the plus sign and one can show 1 = 1 + (c + d )h + 12 (c + d )2h2 + 41 3d + c2d ; 2c d ; 3c h3 From this equation it is clear that 61 3 = 61 (c + d)3 does not match the coecient of h3 in 1 , so er = O(h3) Using P (eh) and Q(eh ) to evaluate er in Section 6.6.3, one can show er = O(h3) using either Eq. 13.8 or Eq. 13.9. These results show that, for the model equation, the hybrid method retains the second-order accuracy of its individual components.

1 + 23 hc + 12 hd =

Stability

The stability of the method can be found from Eq. 13.10 by a parametric study of cn and R de ned in Eq. 13.7. This was carried out in a manner similar to that used to nd the stability boundary of the rst-order explicit-implicit method in Section 13.2.1. The results are plotted in Fig. 13.2. For values of R 2 this second-order method has a much greater region of stability than the rst-order explicit-implicit method given by Eq. 12.10 and shown in Fig. 13.1.

13.3 The Representative Equation for Space-Split Operators Consider the 2-D model3 equations @u = @ 2 u + @ 2 u @t @x2 @y2 3

The extension of the following to 3-D is simple and straightforward.

(13.11)

13.3. THE REPRESENTATIVE EQUATION FOR SPACE-SPLIT OPERATORS239 1.2

0.8

Cn Stable

0.4

0

8

4

R

∆

_ Figure 13.2: Stability regions for the second-order time-split method. and

@u + a @u + a @u = 0 (13.12) @t x @x y @y Reduce either of these, by means of spatial dierencing approximations, to the coupled set of ODE's: dU = [A + A ]U + (bc ~) (13.13) x y dt for the space vector U . The form of the Ax and Ay matrices for three-point central dierencing schemes are shown in Eqs. 12.22 and 12.23 for the 3 4 mesh shown in Sketch 12.12. Let us inspect the structure of these matrices closely to see how we can diagonalize [Ax + Ay ] in terms of the individual eigenvalues of the two matrices considered separately. First we write these matrices in the form 2 3 2 ~ 3 ~b1 I B b I 0 A(xx) = 64 B 75 A(yx) = 64 ~b;1 I ~b0 I ~b1 I 75 ~b;1 I ~b0 I B where B is a banded matrix of the form B (b;1 ; b0; b1 ). Now nd the block eigenvector

240 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS matrix that diagonalizes B and use it to diagonalize A(xx). Thus 2 ;1 32 32 3 2 3 X B X 6 76 7=6 X ;1 B 75 64 X 75 4 54 5 4 X ;1 B X where 2 3 1 6 7 2 7 = 664 7 5 3 4 ( x ) Notice that the matrix Ay is transparent to this transformation. That is, if we set X diag(X ) 2 ~ 3 2 ~ 3 b0 I ~b1 I b0 I ~b1 I X ;1 64 ~b;1 I ~b0 I ~b1 I 75 X = 64 ~b;1 I ~b0 I ~b1 I 75 ~b;1 I ~b0 I ~b;1 I ~b0 I One now permutes the transformed system to the y-vector data-base using the permutation matrix de ned by Eq. 12.13. There results h i Pyx X ;1 A(xx) + A(yx) X Pxy = 3 2 3 2 ~ B 1 I 6 7 6 7 ~ 6 7 I B 6 7 2 6 7 6 7+6 (13.14) 4 5 4 3 I B~ 75 4 I B~ where B~ is the banded tridiagonal matrix B (~b;1 ; ~b0 ; ~b1), see the bottom of Eq. 12.23. Next nd the eigenvectors X~ that diagonalize the B~ blocks. Let B~ diag(B~ ) and X~ diag(X~ ) and form the second transformation 2~ 3 2~ 3 1 6 7 ~ 7 X~ ;1B~ X~ = 664 ~ 775 ; ~ = 64 ~2 5 ~ 3 ~ This time, by the same argument as before, the rst matrix on the right side of Eq. 13.14 is transparent to the transformation, so the nal result is the complete diagonalization of the matrix Ax+y : h ih ih i X~ ;1 Pyx X ;1 A(xx+)y X Pxy X~ = 2 3 ~ I + 1 6 7 6 7 2 I + ~ 6 7 (13.15) 6 7 3 I + ~ 4 5 4I + ~

13.3. THE REPRESENTATIVE EQUATION FOR SPACE-SPLIT OPERATORS241 It is important to notice that: The diagonal matrix on the right side of Eq. 13.15 contains every possible combination of the individual eigenvalues of B and B~ . Now we are ready to present the representative equation for two dimensional systems. First reduce the PDE to ODE by some choice4 of space dierencing. This results in a spatially split A matrix formed from the subsets

A(xx) = diag(B ) ; A(yy) = diag(B~ )

(13.16)

where B and B~ are any two matrices that have linearly independent eigenvectors (this puts some constraints on the choice of dierencing schemes). Although Ax and Ay do commute, this fact, by itself, does not ensure the property of \all possible combinations". To obtain the latter property the structure of the matrices is important. The block matrices B and B~ can be either circulant or noncirculant; in both cases we are led to the nal result: The 2{D representative equation for model linear systems is du = [ + ]u + aet x y dt where x and y are any combination of eigenvalues from Ax and Ay , a and are (possibly complex) constants, and where Ax and Ay satisfy the conditions in 13.16. Often we are interested in nding the value of, and the convergence rate to, the steady-state solution of the representative equation. In that case we set = 0 and use the simpler form du = [ + ]u + a (13.17) x y dt which has the exact solution u(t) = ce(x+y )t ; +a (13.18) x y We have used 3-point central dierencing in our example, but this choice was for convenience only, and its use is not necessary to arrive at Eq. 13.15. 4

242 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS

13.4 Example Analysis of 2-D Model Equations

In the following we analyze four dierent methods for nding a xed, steady-state solution to the 2-D representative equation 13.16. In each case we examine 1. The stability. 2. The accuracy of the xed, steady-state solution. 3. The convergence rate to reach the steady-state.

13.4.1 The Unfactored Implicit Euler Method Consider rst this unfactored, rst-order scheme which can then be used as a reference case for comparison with the various factored ones. The form of the method is given by Eq. 12.19, and when it is applied to the representative equation, we nd (1 ; h x ; h y )un+1 = un + ha from which

P (E ) = (1 ; h x ; h y )E ; 1 Q(E ) = h giving the solution

"

(13.19)

#

n 1 un = c 1 ; h ; h ; +a x y x y Like its counterpart in the 1-D case, this method:

1. Is unconditionally stable. 2. Produces the exact (see Eq. 13.18) steady-state solution (of the ODE) for any h. 3. Converges very rapidly to the steady-state when h is large. Unfortunately, however, use of this method for 2-D problems is generally impractical for reasons discussed in Section 12.5.

13.4. EXAMPLE ANALYSIS OF 2-D MODEL EQUATIONS

243

13.4.2 The Factored Nondelta Form of the Implicit Euler Method

Now apply the factored Euler method given by Eq. 12.20 to the 2-D representative equation. There results (1 ; h x)(1 ; h y )un+1 = un + ha from which P (E ) = (1 ; h x)(1 ; h y )E ; 1 Q(E ) = h (13.20) giving the solution " #n 1 un = c (1 ; h )(1 ; h ) ; + a; h x y x y x y We see that this method: 1. Is unconditionally stable. 2. Produces a steady state solution that depends on the choice of h. 3. Converges rapidly to a steady-state for large h, but the converged solution is completely wrong. The method requires far less storage then the unfactored form. However, it is not very useful since its transient solution is only rst-order accurate and, if one tries to take advantage of its rapid convergence rate, the converged value is meaningless.

13.4.3 The Factored Delta Form of the Implicit Euler Method

Next apply Eq. 12.31 to the 2-D representative equation. One nds (1 ; h x)(1 ; h y )(un+1 ; un) = h(xun + y un + a) which reduces to (1 ; h x)(1 ; h y )un+1 = 1 + h2xy un + ha and this has the solution " #n 2 xy 1 + h un = c (1 ; h )(1 ; h ) ; +a x y x y This method:

244 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS 1. Is unconditionally stable. 2. Produces the exact steady-state solution for any choice of h. 3. Converges very slowly to the steady{state solution for large values of h, since jj ! 1 as h ! 1. Like the factored nondelta form, this method demands far less storage than the unfactored form, as discussed in Section 12.5. The correct steady solution is obtained, but convergence is not nearly as rapid as that of the unfactored form.

13.4.4 The Factored Delta Form of the Trapezoidal Method

Finally consider the delta form of a second-order time-accurate method. Apply Eq. 12.30 to the representative equation and one nds 1 1 1 ; 2 hx 1 ; 2 hy (un+1 ; un) = h(xun + y un + a) which reduces to 1 1 1 1 1 ; hx 1 ; hy un+1 = 1 + hx 1 + hy un + ha 2 2 2 2 and this has the solution 2 1 h 1 + 1 h 3n 1 + x y 7 6 un = c64 21 21 75 ; +a x y 1 ; 2 hx 1 ; 2 hy This method: 1. Is unconditionally stable. 2. Produces the exact steady{state solution for any choice of h. 3. Converges very slowly to the steady{state solution for large values of h, since jj ! 1 as h ! 1. All of these properties are identical to those found for the factored delta form of the implicit Euler method. Since it is second order in time, it can be used when time accuracy is desired, and the factored delta form of the implicit Euler method can be used when a converged steady-state is all that is required.5 A brief inspection of eqs. 12.26 and 12.27 should be enough to convince the reader that the 's produced by those methods are identical to the produced by this method. In practical codes, the value of h on the left side of the implicit equation is literally switched from h to 12 h. 5

13.5. EXAMPLE ANALYSIS OF THE 3-D MODEL EQUATION

245

13.5 Example Analysis of the 3-D Model Equation

The arguments in Section 13.3 generalize to three dimensions and, under the conditions given in 13.16 with an A(zz) included, the model 3-D cases6 have the following representative equation (with = 0): du = [ + + ]u + a (13.21) x y z dt Let us analyze a 2nd-order accurate, factored, delta form using this equation. First apply the trapezoidal method: un+1 = un + 12 h[(x + y + z )un+1 + (x + y + z )un + 2a] Rearrange terms: 1 ; 21 h(x + y + z ) un+1 = 1 + 21 h(x + y + z ) un + ha Put this in delta form: 1 ; 1 h(x + y + z ) un = h[(x + y + z )un + a] 2 Now factor the left side: 1 1 1 h 1 ; h 1 ; h 1; 2 x 2 y 2 z un = h[(x + y + z )un + a] (13.22) This preserves second order accuracy since the error terms 1 h2( + + )u and 1 h3 4 x y x z y z n 8 x y z are both O(h3). One can derive the characteristic polynomial for Eq. 13.22, nd the root, and write the solution either in the form 2 3n 1 1 1 2 3 6 1 + h(x + y + z ) + 4 h (x y + x z + y z ) ; 8 h x y z 7 7 un = c64 21 5 1 1 2 3 1 ; 2 h(x + y + z ) + 4 h (xy + xz + y z ) ; 8 h xy z

; + a + x y z 6

Eqs. 13.11 and 13.12, each with an additional term.

(13.23)

246 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS or in the form

3n 2 1 1 1 1 3 6 1 + 2 hx 1 + 2 hy 1 + 2 hz ; 4 h x y z 7 a (13.24) 7 ; un = c64 5 1 1 1 x + y + z 1 ; hx 1 ; hy 1 ; hz

2 2 2 It is interesting to notice that a Taylor series expansion of Eq. 13.24 results in = 1 + h(x + y + z ) + 21 h2(x + y + z )2 (13.25) h i + 14 h3 3z + (2y + 2x) + 22y + 3xy + 22y + 3y + 2x2y + 22xy + 3x + which veri es the second order accuracy of the factored form. Furthermore, clearly, if the method converges, it converges to the proper steady-state.7 With regards to stability, it follows from Eq. 13.23 that, if all the 's are real and negative, the method is stable for all h. This makes the method unconditionally stable for the 3-D diusion model when it is centrally dierenced in space. Now consider what happens when we apply this method to the biconvection model, the 3-D form of Eq. 13.12 with periodic boundary conditions. In this case, central dierencing causes all of the 's to be imaginary with spectrums that include both positive and negative values. Remember that in our analysis we must consider every possible combination of these eigenvalues. First write the root in Eq. 13.23 in the form + i ; + i = 11 ; i ; + i

where , and are real numbers that can have any sign. Now we can always nd one combination of the 's for which , and are both positive. In that case since the absolute value of the product is the product of the absolute values 2 + ( + )2 (1 ; ) 2 jj = (1 ; )2 + ( ; )2

>1

and the method is unconditionally unstable for the model convection problem. From the above analysis one would come to the conclusion that the method represented by Eq. 13.22 should not be used for the 3-D Euler equations. In practical cases, however, some form of dissipation is almost always added to methods that are used to solve the Euler equations and our experience to date is that, in the presence of this dissipation, the instability disclosed above is too weak to cause trouble. 7

However, we already knew this because we chose the delta form.

13.6. PROBLEMS

13.6 Problems

247

1. Starting with the generic ODE, du = Au + f dt we can split A as follows: A = A1 + A2 + A3 + A4. Applying implicit Euler time marching gives un+1 ; un = A u + A u + A u + A u + f 1 n+1 2 n+1 3 n+1 4 n+1 h (a) Write the factored delta form. What is the error term? (b) Instead of making all of the split terms implicit, leave two explicit: un+1 ; un = A u + A u + A u + A u + f 1 n+1 2 n 3 n+1 4 n h Write the resulting factored delta form and de ne the error terms. (c) The scalar representative equation is du = ( + + + )u + a 1 2 3 4 dt For the fully implicit scheme of problem 1a, nd the exact solution to the resulting scalar dierence equation and comment on the stability, convergence, and accuracy of the converged steady-state solution. (d) Repeat 1c for the explicit-implicit scheme of problem 1b.

248 CHAPTER 13. LINEAR ANALYSIS OF SPLIT AND FACTORED FORMS

Appendix A USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA A basic understanding of the fundamentals of linear algebra is crucial to our development of numerical methods and it is assumed that the reader is at least familar with this subject area. Given below is some notation and some of the important relations between matrices and vectors.

A.1 Notation 1. In the present context a vector is a vertical column or string. Thus 2 v1 3 6 7 ~v = 66 v..2 77 4 . 5 vm T and its transpose ~v is the horizontal row ~vT = [v1 ; v2; v3; : : : ; vm] ; ~v = [v1 ; v2; v3 ; : : : ; vm ]T 2. A general m m matrix A can be written 2 a11 a12 6a a A = (aij ) = 664 21 22 am1 am2 249

a1m 3 a2m 777 ...

amm

5

250APPENDIX A. USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA 3. An alternative notation for A is h i A = ~a1; ~a2; : : : ; ~am and its transpose AT is 2 ~T 3 66 ~a1T 77 AT = 666 a.2 777 4 ..T 5 ~am

4. The inverse of a matrix (if it exists) is written A;1 and has the property that A;1A = AA;1 = I , where I is the identity matrix.

A.2 De nitions 1. A is symmetric if AT = A. 2. A is skew-symmetric or antisymmetric if AT = ;A. 3. A is diagonally dominant if aii Pj6=i jaij j ; i = 1; 2; : : : ; m and aii > Pj6=i jaij j for at least one i. 4. A is orthogonal if aij are real and AT A = AAT = I 5. A is the complex conjugate of A. 6. P is a permutation matrix if P ~v is a simple reordering of ~v. 7. The trace of a matrix is Pi aii. 8. A is normal if AT A = AAT . 9. det[A] is the determinant of A. 10. AH is the conjugate transpose of A, (Hermitian). 11. If a b A= c d then det[A] = ad ; bc and A;1 = det[1 A] ;dc ;ab

A.3. ALGEBRA

251

A.3 Algebra We consider only square matrices of the same dimension. 1. A and B are equal if aij = bij for all i; j = 1; 2; : : : ; m. 2. A + (B + C ) = (C + A) + B , etc. 3. sA = (saij ) where s is a scalar. 4. In general AB 6= BA. 5. Transpose equalities: (A + B )T = AT + B T (AT )T = A (AB )T = B T AT 6. Inverse equalities (if the inverse exists): (A;1);1 = A (AB );1 = B ;1 A;1 (AT );1 = (A;1 )T

7. Any matrix A can be expressed as the sum of a symmetric and a skew-symmetric matrix. Thus: A = 21 A + AT + 12 A ; AT

A.4 Eigensystems 1. The eigenvalue problem for a matrix A is de ned as A~x = ~x or [A ; I ]~x = 0 and the generalized eigenvalue problem, including the matrix B , as A~x = B~x or [A ; B ]~x = 0 2. If a square matrix with real elements is symmetric, its eigenvalues are all real. If it is asymmetric, they are all imaginary.

252APPENDIX A. USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA 3. Gershgorin's theorem: The eigenvalues of a matrix lie in the complex plane in the union of circles having centers located by the diagonals with radii equal to the sum of the absolute values of the corresponding o-diagonal row elements. 4. In general, an m m matrix A has n~x linearly independent eigenvectors with n~x m and n distinct eigenvalues (i) with n n~x m. 5. A set of eigenvectors is said to be linearly independent if a ~xm + b ~xn 6= ~xk ; m 6= n 6= k for any complex a and b and for all combinations of vectors in the set. 6. If A posseses m linearly independent eigenvectors then A is diagonalizable, i.e.,

X ;1AX = where X is a matrix whose columns are the eigenvectors, h i X = ~x1 ; ~x2; : : : ; ~xm and is the diagonal matrix

2 0 66 01 . . . = 66 .. . .2 . . 4. . .

03 ... 77 7 0 75 0 0 m If A can be diagonalized, its eigenvectors completely span the space, and A is said to have a complete eigensystem.

7. If A has m distinct eigenvalues, then A is always diagonalizable, and with each distinct eigenvalue there is one associated eigenvector, and this eigenvector cannot be formed from a linear combination of any of the other eigenvectors. 8. In general, the eigenvalues of a matrix may not be distinct, in which case the possibility exists that it cannot be diagonalized. If the eigenvalues of a matrix are not distinct, but all of the eigenvectors are linearly independent, the matrix is said to be derogatory, but it can still be diagonalized. 9. If a matrix does not have a complete set of linearly independent eigenvectors, it cannot be diagonalized. The eigenvectors of such a matrix cannot span the space and the matrix is said to have a defective eigensystem.

A.4. EIGENSYSTEMS

253

10. Defective matrices cannot be diagonalized but they can still be put into a compact form by a similarity transform, S , such that

2J 0 0 3 66 01 J . . . ... 77 7 J = S ;1AS = 66 .. . .2 . . 4 . . . 0 75 0 0 Jk

where there are k linearly independent eigenvectors and Ji is either a Jordan subblock or i. 11. A Jordan submatrix has the form

2 1 0 66 0i 1 i 66 Ji = 66 0 0 i 64 ... ... 0 0

0 3 ... ... ... 0

... 77 7 0 777 1 75 i

12. Use of the transform S is known as putting A into its Jordan Canonical form. A repeated root in a Jordan block is referred to as a defective eigenvalue. For each Jordan submatrix with an eigenvalue i of multiplicity r, there exists one eigenvector. The other r ; 1 vectors associated with this eigenvalue are referred to as principal vectors. The complete set of principal vectors and eigenvectors are all linearly independent. 13. Note that if P is the permutation matrix

2 3 0 0 1 P = 64 0 1 0 75 1 0 0

then

;

P T = P ;1 = P

2 3 2 3 1 0 0 0 P ;1 64 0 1 75 P = 64 1 0 75 0 0 0 1

14. Some of the Jordan subblocks may have the same eigenvalue. For example, the

254APPENDIX A. USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA matrix

3 2 2 1 1 66 64 1 1 75 66 1 66 1 66 1 1 66 1 66 2 1 64 2

is both defective and derogatory, having:

3

3 77 77 77 77 77 77 75

9 eigenvalues 3 distinct eigenvalues 3 Jordan blocks 5 linearly independent eigenvectors 3 principal vectors with 1 1 principal vector with 2

A.5 Vector and Matrix Norms 1. The spectral radius of a matrix A is symbolized by (A) such that

(A) = jm jmax where m are the eigenvalues of the matrix A. 2. A p-norm of the vector ~v is de ned as

0M 11=p X jjvjjp = @ jvj jpA j =1

3. A p-norm of a matrix A is de ned as

jjAvjjp jjAjjp = max x6=0 jjv jj p

A.5. VECTOR AND MATRIX NORMS

255

4. Let A and B be square matrices of the same order. All matrix norms must have the properties

jjAjj jjc Ajj jjA + B jj jjA B jj

0; jjAjj = 0 implies A = 0 = jcj jjAjj jjAjj + jjB jj jjAjj jjB jj

5. Special p-norms are jjAjj1 = maxj=1;;M PMi=1 jaij j maximum column sum

q T jjAjj2 = (A A ) jjAjj1 = maxi=1;2;;M PMj=1 jaij j maximum row sum where jjAjjp is referred to as the Lp norm of A. 6. In general (A) does not satisfy the conditions in 4, so in general (A) is not a true norm. 7. When A is normal, (A) is a true norm, in fact, in this case it is the L2 norm. 8. The spectral radius of A, (A), is the lower bound of all the norms of A.

256APPENDIX A. USEFUL RELATIONS AND DEFINITIONS FROM LINEAR ALGEBRA

Appendix B SOME PROPERTIES OF TRIDIAGONAL MATRICES B.1 Standard Eigensystem for Simple Tridiagonals In this work tridiagonal banded matrices are prevalent. It is useful to list some of their properties. Many of these can be derived by solving the simple linear dierence equations that arise in deriving recursion relations. Let us consider a simple tridiagonal matrix, i.e., a tridiagonal with constant scalar elements a,b, and c, see Section 3.4. If we examine the conditions under which the determinant of this matrix is zero, we nd (by a recursion exercise) det[B (M : a; b; c)] = 0 if p m b + 2 ac cos M + 1 = 0 ; m = 1; 2; ; M

From this it follows at once that the eigenvalues of B (a; b; c) are p m (B.1) m = b + 2 ac cos M + 1 ; m = 1; 2; ; M The right-hand eigenvector of B (a; b; c) that is associated with the eigenvalue m satis es the equation B (a; b; c)~xm = m~xm (B.2) and is given by j ; 1 ~xm = (xj )m = a 2 sin j m (B.3) c M + 1 ; m = 1; 2; ; M 257

258

APPENDIX B. SOME PROPERTIES OF TRIDIAGONAL MATRICES

These vectors are the columns of the right-hand eigenvector matrix, the elements of which are j ; 1 j = 1; 2; ; M X = (xjm) = ac 2 sin Mjm ; (B.4) m = 1; 2; ; M +1 Notice that if a = ;1 and c = 1, j ; 1 a 2 = ei j; 2 (B.5) c (

1)

The left-hand eigenvector matrix of B (a; b; c) can be written m ; 1 c 2 mj = 1; 2; ; M 2 ; X = M +1 a sin M + 1 ; m j = 1; 2; ; M In this case notice that if a = ;1 and c = 1 m ; 1 c 2 = e;i m; 2 a 1

(

1)

(B.6)

B.2 Generalized Eigensystem for Simple Tridiagonals This system is de ned as follows 2 32 32 x 3 2e f x 3 b c 6 7 6 7 6 x 777 666 d e f x 777 a b c 6 7 6 7 6 6 7 6 7 6 6 7 6 7 6 x 77 = 66 d e x 77 a b 6 7 6 7 6 . . . 6 7 6 7 6 7 6 . . . c 5 4 .. 5 4 . f 5 4 ... 75 4 a b xM d e xM In this case one can show after some algebra that det[B (a ; d; b ; e; c ; f ] = 0 (B.7) if q m b ; m e + 2 (a ; md)(c ; m f ) cos M + 1 = 0 ; m = 1; 2; ; M (B.8) If we de ne m = Mm ; p = cos m +1 m 1

1

2

2

3

3

B.3. THE INVERSE OF A SIMPLE TRIDIAGONAL

259

q

m = eb ; 2(cd + af )pm + 2pm (ec ; fb)(ea ; bd) + [(cd ; af )pm ] e ; 4fdpm The right-hand eigenvectors are #j ; 1 " 2 a ; d m = 1; 2; ; M m ~xm = sin [ j m] ; j = 1; 2; ; M c; f 2

2

2

2

m

These relations are useful in studying relaxation methods.

B.3 The Inverse of a Simple Tridiagonal The inverse of B (a; b; c) can also be written in analytic form. Let DM represent the determinant of B (M : a; b; c)

DM det[B (M : a; b; c)] De ning D to be 1, it is simple to derive the rst few determinants, thus D = 1 D = b D = b ; ac D = b ; 2abc (B.9) One can also nd the recursion relation DM = bDM ; ; acDM ; (B.10) Eq. B.10 is a linear OE the solution of which was discussed in Section 4.2. Its characteristic polynomial P (E ) is P (E ; bE + ac) and the two roots to P () = 0 result in the solution 8" p p " #M 9 < b + b ; 4ac #M = 1 b ; b ; 4 ac DM = p ; ; 2 2 b ; 4ac : M = 0; 1; 2; (B.11) where we have made use of the initial conditions D = 1 and D = b. In the limiting case when b ; 4ac = 0, one can show that 0

0 1

2

2

3

3

1

2

2

+1

2

2

2

0

2

DM = (M + 1) 2b

!M

1

; b = 4ac 2

+1

260

APPENDIX B. SOME PROPERTIES OF TRIDIAGONAL MATRICES

Then for M = 4

B; = 1

and for M = 5

B ;1 =

2 1 666 D4 4

D ;cD cD ;c D ;aD D D ;cD D c D a D ;aD D D D ;cD ;a D a D ;aD D 3

3

2 D4 6 3 1 666 ;a2aD D5 64 ;a3DD21

aD The general element dmn is 4

0

2

2

2

2

1

1

2

0

2

1

1

1

2

1

;cD

3

1

0

2

1

1

1

2

2

3 7 7 7 5

3

3 7 17 7 2 7 7 5 3

cD ;c D c D D D ;cD D c D D ;c D ;aD D D D ;cD D c D a D D ;aD D D D ;cD ;a D aD ;aD D 2

3

1

2

3

1

1

1 3

2

2

1

2

1

3

2

2

2

1

2

2

1

2

1

4

1

3

2

1

0

3

2

1

3

4

Upper triangle:

m = 1; 2; ; M ; 1 ; n = m + 1; m + 2; ; M dmn = Dm; DM ;n(;c)n;m =DM 1

Diagonal:

n = m = 1; 2; ; M dmm = DM ; DM ;m =DM 1

Lower triangle:

m = n + 1; n + 2; ; M ; n = 1; 2; ; M ; 1 dmn = DM ;mDn; (;a)m;n=DM 1

B.4 Eigensystems of Circulant Matrices B.4.1 Standard Tridiagonals

Consider the circulant (see Section 3.4.4) tridiagonal matrix

Bp(M : a; b; c; )

(B.12)

B.4. EIGENSYSTEMS OF CIRCULANT MATRICES

261

The eigenvalues are

m ; i(a ; c) sin 2m ; m = 0; 1; 2; ; M ; 1 m = b + (a + c) cos 2M M (B.13) The right-hand eigenvector that satis es Bp(a; b; c)~xm = m~xm is ~xm = (xj )m = ei j m=M ; j = 0; 1; ; M ; 1 (B.14) p where i ;1, and the right-hand eigenvector matrix has the form ij 2m 0; 1; ; M ; 1 X = (xjm) = e M ; mj = = 0; 1; ; M ; 1 The left-hand eigenvector matrix with elements x0 is 2 j ;im M ; m = 0; 1; ; M ; 1 X ; = (x0mj ) = M1 e j = 0; 1; ; M ; 1 Note that both X and X ; are symmetric and that X ; = M1 X , where X is the conjugate transpose of X . (2

)

1

1

1

B.4.2 General Circulant Systems

Notice the remarkable fact that the elements of the eigenvector matrices X and X ; for the tridiagonal circulant matrix given by eq. B.12 do not depend on the elements a; b; c in the matrix. In fact, all circulant matrices of order M have the same set of linearly independent eigenvectors, even if they are completely dense. An example of a dense circulant matrix of order M = 4 is 2 b b b b 3 6 b b b b 777 6 6 (B.15) 4 b b b b 5 b b b b The eigenvectors are always given by eq. B.14, and further examination shows that the elements in these eigenvectors correspond to the elements in a complex harmonic analysis or complex discrete Fourier series. Although the eigenvectors of a circulant matrix are independent of its elements, the eigenvalues are not. For the element indexing shown in eq. B.15 they have the general form MX ; m = bj ei jm=M 1

0

1

2

3

3

0

1

2

2

3

0

1

1

2

3

0

1

of which eq. B.13 is a special case.

j =0

(2

)

262

APPENDIX B. SOME PROPERTIES OF TRIDIAGONAL MATRICES

B.5 Special Cases Found From Symmetries

Consider a mesh with an even number of interior points such as that shown in Fig. B.1. One can seek from the tridiagonal matrix B (2M : a; b; a; ) the eigenvector subset that has even symmetry when spanning the interval 0 x . For example, we seek the set of eigenvectors ~xm for which 2 b 6 a 6 6 6 6 6 6 6 4

a b a a ...

32x 3 1 6 7 7 x 6 2 7 7 6 7 ... 77 6 7 6 7 7 6 7 . 7 . 6 7 7 . 7 76 6 a 5 4 x2 75

... a a b a b

x

2 3 x1 6 x2 777 6 6 ... 7 6 6 = m 66 .. 777 .7 6 6 7 4 x2 5

x

1

1

This leads to the subsystem of order M which has the form 2 6 6 6 6 6 6 6 6 6 6 4

3

b a 7 7 a b a 7 7 . . 7 . a ~ ~ ~xm = m~xm 7 B (M : a; b; a)xm = 7 ... a 7 7 a b a 75 a b+a

(B.16)

By folding the known eigenvectors of B (2M : a; b; a) about the center, one can show from previous results that the eigenvalues of eq. B.16 are

m = b + 2a cos (22mM;+1)1

!

; m = 1; 2; ; M

(B.17)

B.6. SPECIAL CASES INVOLVING BOUNDARY CONDITIONS

263

and the corresponding eigenvectors are ~xm = sin j (2m ; 1) ; 2M + 1 j = 1; 2; ; M

x=

Line of Symmetry

x =0

Imposing symmetry about the same interval but for a mesh with an odd number of points, see Fig. B.1, leads to the matrix 2 3 b a 6 7 a b a 6 7 6 7 . . 6 7 . a B (M : ~a; b; a) = 66 . . . a 777 6 6 4 a b a 75 2a b

~xm = sin j (2m ; 1) 2M

!

j0 = 1 2 3 4 5 M0 j= 1 2 3 M b. An odd{numbered mesh Figure B.1 { Symmetrical folds for special cases

By folding the known eigenvalues of B (2M ; 1 : a; b; a) about the center, one can show from previous results that the eigenvalues of eq. B.17 are

; 1) m = b + 2a cos (2m2M and the corresponding eigenvectors are

j0 = 1 2 3 4 5 6 M0 j= 1 2 3 M a. An even-numbered mesh Line of Symmetry x =0 x=

!

; m = 1; 2; ; M

; j = 1; 2; ; M

B.6 Special Cases Involving Boundary Conditions We consider two special cases for the matrix operator representing the 3-point central dierence approximation for the second derivative @ =@x at all points away from the boundaries, combined with special conditions imposed at the boundaries. 2

2

APPENDIX B. SOME PROPERTIES OF TRIDIAGONAL MATRICES

264

Note: In both cases

m = 1; 2; ; M j = 1; 2; ; M ;2 + 2 cos() = ;4 sin (=2) 2

When the boundary conditions are Dirichlet on both sides, 2 6 6 6 6 6 6 4

;2

1 1 ;2 1 1 ;2 1 1 ;2 1 1 ;2

3 7 7 7 7 7 7 5

m = ;2 h+2 cos Mim +1 ~xm = sin j m M +1

(B.18)

When one boundary condition is Dirichlet and the other is Neumann (and a diagonal preconditioner is applied to scale the last equation), 2 6 6 6 6 6 6 4

;2

1 1 ;2 1 1 ;2 1 1 ;2 1 1 ;1

3 7 7 7 7 7 7 5

m ~xm

(2 m ; 1) = ;2 + 2 cos 2M + 1 (2 m ; 1) = sin j 2M + 1

(B.19)

Appendix C THE HOMOGENEOUS PROPERTY OF THE EULER EQUATIONS The Euler equations have a special property that is sometimes useful in constructing numerical methods. In order to examine this property, let us rst inspect Euler's theorem on homogeneous functions. Consider rst the scalar case. If F (u; v) satis es the identity

F (u; v) = nF (u; v) (C.1) for a xed n, F is called homogeneous of degree n. Dierentiating both sides with respect to and setting = 1 (since the identity holds for all ), we nd u @F + v @F = nF (u; v) @u @v

(C.2)

F (Q) = nF (Q)

(C.3)

Consider next the theorem as it applies to systems of equations. If the vector F (Q) satis es the identity

for a xed n, F is said to be homogeneous of degree n and we nd "

#

@F Q = nF (Q) @q 265

(C.4)

266APPENDIX C. THE HOMOGENEOUS PROPERTY OF THE EULER EQUATIONS Now it is easy to show, by direct use of eq. C.3, that both E and F in eqs. 2.11 and 2.12 are homogeneous of degree 1, and their Jacobians, A and B , are homogeneous of degree 0 (actually the latter is a direct consequence of the former).1 This being the case, we notice that the expansion of the ux vector in the vicinity of tn which, according to eq. 6.105 can be written in general as,

E = En + An(Q ; Qn) + O(h2) F = Fn + Bn(Q ; Qn) + O(h2)

(C.5)

can be written

E = An Q + O(h2) F = BnQ + O(h2)

(C.6)

since the terms En ; AnQn and Fn ; BnQn are identically zero for homogeneous vectors of degree 1, see eq. C.4. Notice also that, under this condition, the constant term drops out of eq. 6.106. As a nal remark, we notice from the chain rule that for any vectors F and Q "

#

@F (Q) = @F @Q = A @Q @x @Q @x @x We notice also that for a homogeneous F of degree 1, F = AQ and "

#

@F = A @Q + @A Q @x @x @x

(C.7)

(C.8)

Therefore, if F is homogeneous of degree 1, "

#

@A Q = 0 @x in spite of the fact that individually [@A=@x] and Q are not equal to zero.

(C.9)

Note that this depends on the form of the equation of state. The Euler equations are homogeneous if the equation of state can be written in the form p = f (), where is the internal energy per unit mass. 1

Bibliography [1] D. Anderson, J. Tannehill, and R. Pletcher. Computational Fluid Mechanics and Heat Transfer. McGraw-Hill Book Company, New York, 1984. [2] C. Hirsch. Numerical Computation of Internal and External Flows, volume 1,2. John Wiley & Sons, New York, 1988.

267

Recommend Stories

Story Transcript

Get in touch

Social