Lecture 1Wednesday, March 30, 20051 IntroductionThis course is concerned with the numerical solution of partial differentialequations. The focus is on hyperbolic PDE. We will also discuss basic meth-ods for numerically solving parabolic and e lliptic equations as they relateto hyperbolic problems. Our approach follows the historical development ofthe numerical metho ds. Specifically, we will study the numerical solution ofequations of types1. hyperbolic2. parabolic3. ellipticthrough looking at, respectively, the problems ofA. compressible flowB. heat flowC. incompressible flow2 ConvectionReference: Osher and Fedkiw, Section 3.1There are two different approaches we can take in describing fluid motion.1. In the Lagrangian formulation, we follow the motion of individualparticles as they are advected with the flow. Assuming that a velocity1~V (~x) = hu, v, wi is known for every point ~x we can compute the particlemotion by solving the ordinary differential equation (ODE)d~xdt=~V (~x) (1)This approach is simple, and the equation is easy to solve numerically,but it is hard to apply it to nonlinear phenomena such as shock waves.2. In the Eulerian formulation we use the simple convection (or advec-tion) equationφt+~V · ∇φ = 0 (2)where the t subscript denotes a temporal partial derivative in the timevariable t. Recall that ∇ is the gradient operator so that~V · ∇φ =uφx+ vφy+ wφz. This approach can be regarded as sitting still at apoint ~x and observing the changes in various quantities at that pointdue to the flow.3 Upwind DifferencingReference: Osher and Fedkiw, Section 3.23.1 NotationIn the notes that follow we will always indicate the time step in the super-script, and the spatial indices in the subscripts.At some point in time, say time tn, let φn= φ(tn) represent the currentvalues of φ. Updating φ in time consists of finding new values of φ at everygrid point after some time increment 4t. We denote these new values of φby φn+1= φ(tn+1) where tn+1= tn+ 4t.3.2 Temporal DiscretizationWe look at some possible discretizations for the time derivatives in the aboveconvection equation. The equations below are called semi-discrete becausewe have discretized only the time derivatives.• forward Eulerφn+1− φn4t+~Vn· ∇φn= 0 (3)2• backward Eulerφn+1− φn4t+~Vn+1· ∇φn+1= 0 (4)The forward Euler and backward Euler are both first order accurate, mean-ing that the error in the discretization is O(4t).3.3 Spatial DiscretizationWe begin by writing equation 3 in expanded form asφn+1− φn4t+ unφnx+ vnφny+ wnφnz= 0 (5)and address the evaluation of the unφnxterm first. The techniques used toapproximate this term can then be applied independently to the vnφnyandwnφnzterms in a dimension by dimension manner.For simplicity, consider the one dimensional version of equation 5φn+1− φn4t+ unφnx= 0 (6)where the sign of unindicates whether the values of φ are moving to theright or to the left. Since uncan be spatially varying, we focus on a specificgrid point xiwhere we writeφn+1i− φni4t+ uni(φx)ni= 0 (7)We use the forward Euler discretization in time. This means that we cansolve the equations for each grid point independently, and we will not haveto solve a linear system at each time step as in the case of backward Euler.We first introduce the difference operators Do, D+and D−.• Do, central difference operator (second-order accurate)(Doφ)i=φi+1− φi−124x(8)• D+, forward difference operator (first-order accurate)(D+φ)i=φi+1− φi4x(9)3• D−, backward difference operator (first-order accurate)(D−φ)i=φi− φi−14x(10)For the hyperbolic equations, the data is propagated in specific directions(along the characteristic curves). Therefore central differencing, which usesinformation from both directions, is not very useful for numerically solvinghyperbolic problems. Instead we prefer to use upwinding. In upwindingthe idea is to choose the s patial discretization based on the direction thatinformation is propagating. For the one-dimensional case,if ui> 0, we use D−φif ui< 0, we use D+φif ui= 0, then uiφx= 0, so we do not need to approximate φxUpwinding is first-order accurate. If we instead choose the difference op-erators so that data is taken from the direction opposite the one from whichit is propagating, then we are using downwinding which is unstable.The numerical errors resulting from upwind differencing cause dissipat ionin the numerical solution. The numerical errors resulting from central differ-encing cause dispersion in the numerical solution (see Strikwerda, Chapter5).In any numerical method for solving PDE, we must be conce rned withconvergence. The means that as we refine our grid in time and space, thenumerical solution converges to the analytical solution. To this end weexamine the notions of consistency and stability. For precise definitions ofconvergence, consistency, and stability, see Strikwerda, Chapter 1.1. ConsistencyThe errors in approximating the differential operator vanish as 4t,4x → 0. This is typically proved by assuming a sufficiently smoothsolution and using Taylor series expansion.2. StabilityThe solution does not blow up in a finite time. Stability may be uncon-ditional or conditional. In conditional stability, we have restrictionson the value of 4t so that it be sufficiently small.For example, it can be shown that discretizing the scalar, linear, constantcoefficient convection equation above using central differencing in space andforward Euler in time is consistent, but unstable.4According to the Lax-Richtmyer Equivalence Theorem a finite differenceapproximation to a linear partial differential equation is convergent if andonly if it is both consistent and stable.A necessary condition for stability is the Courant-Friedreichs-Lewy con-dition (CFL condition) which asserts that the numerical waves should prop-agate at least as fast as the physical waves. This me ans that the numericalwave speed of 4x/4t must be at least as fast as the physical wave speed|u|, i.e. 4x/4t > |u|. This leads us to the CFL time step restriction of4t <4xmax{|u|}(11)where max{|u|} is chosen to be the largest value of |u| over the entire Carte-sian grid. Equation 11 is usually enforced by choosing a CFL number αwith4tmax{|u|}4x= α (12)and 0 < α < 1. A common near optimal choice is α = .9, and a commonconservative choice is α = .5. A multidimensional CFL condition c an bewritten as4t max|u|4x+|v|4y+|w|4z= α