Class 22. PDEs, Part 1• Cf. NRiC §19.Classification of PDEs• A PDE is simply a differential equation of more than one variable (so an ODE is aspecial case of a PDE). PDEs are usually classified into three types:1. Hyperbolic (second or first order in time and space)– Prototype is the wave equation:∂2u∂t2= v2∂2u∂x2(1)(this is the 1-D version), where v = (constant) wave speed and u = amplitude.2. Parabolic(first order in time, second order in space)– Prototype is the diffusion equation:∂u∂t=∂∂xD∂u∂x(2)(1-D), where D = diffusion coefficient, u = amplitude.3. Elliptic (second order in space)– Prototype is the Poisson equation:∂2u∂x2+∂2u∂y2+∂2u∂z2= ρ (3)(3-D), where ρ = density ( if ρ = 0, get Laplace equation).• Note that (1) and (2) define initial value problems. If u(x) (and perhaps ∂u/∂x) definedat t = t0, then equations define how u(x, t) propagates forward in time. ∴ numericalsolutions of (1) and (2) give time evolution of u (e.g., wave amplitude).• On the other hand, (3) defines a boundary value problem. Given static function ρ, findstatic solution u satisfying BCs. ∴ numerical solution of (3) gives space distribution ofu (e.g., gravitational potential).• Distinction between IVPs vs. BVPs more important than distinction between (1) and(2). Often, IVPs are mixture of hyperbolic and para bolic.1Solving Elliptic PDEs (BVP)• Already discussed this at length for PM codes: finite differencing yields lar ge set o fcoupled algebraic equations =⇒ la r ge sparse banded matrix.• Many techniques for solving matrix:1. Relaxation schemes.2. Sparse banded matrix solvers.3. Fourier methods.• Use #3 when you can, #1 or #2 otherwise.Solving Hyperbolic PDEs (IVP)• NRiC §19.1.• Overriding concern is stability of algorithm.Conservative form• Large class of IVP can be put in “ flux-conservative” form:∂u∂t= −∂F(u)∂x, (4)where F = flux of conserved quantity. In multidimensions,∂u∂t= −∇ · F(this is in the f orm of a conservation law).• Fo r example, prototypical hyperbolic PDE∂2u∂t2= v2∂2u∂x2(v constant) can be decomposed into two first-order equations:∂r∂t= v∂s∂x,∂s∂t= v∂r∂x,wherer ≡ v∂u∂x, s ≡∂u∂t.(E.F.T.S.: show that these two equations do indeed combine to give the original second-order equation.) Then letu =rs, F(u) =0 −v−v 0u =−vs−vr.Plugging these into the conservative form (4) gives the decomposed versio n of the PDE.2The scalar advection equation• If we can cast our hyperbolic PDE into conservative form, then all we need to do isdevelop numerical solution strategies for the first-order equations, which can usuallybe written in the form:∂u∂t= −v∂u∂x(5)(v still constant). We happen to a lready know the analytical solution is u = f(x − vt),i.e., f unction f displaced by vt,1x0x0+ vt> 0t= 0tbut we do not necessarily know the exact form of f. Equation (5) is a scalar advectionequation (the quantity u is transported by a “fluid flow” with a speed v) .• Best example of (5) in astrophysics is continuity equation, i.e., conservation law forsome quantity with density ρ. Evolution of ρ (in 1-D) obeys∂ρ∂t+ v∂ρ∂x= 0ifRρdx = constant, i.e., material conserved. Describes how material is mixed in ISM,how mass is tr ansported. One of the equations of fluid dynamics.Forward time centered space (FTCS) scheme• How can we construct a numerical solution to (5)?• Try simple Euler differencing:un+1j− unj∆t= −vunj+1− unj−12∆x. (6)This is first order in time and second o rder in space. Leads to the forward time centeredspace (FTCS) scheme.1To se e this, let w = x − vt and differentiate u = f(w) using the chain rule: ∂f /∂t = (∂f/∂w)(∂w/∂t) =−v(∂f /∂w); −v(∂f /∂x) = −v(∂f /∂w)(∂w/∂x) = −v(∂f /∂w).3• Schematically:new pointknown pointstx−1 +1j jj+1nn• Explicit in time (just solve for un+1j).• What about stability of scheme?von N eumann st ability analysis• To check stability, customary to perform a von Neumann stability analysis.• Treat all coefficients of difference equations as constant in x and t (local analysis).• Then, eigenmodes of difference equations all of formunj= ξneikj∆x, (7)where ξ(k) is the (complex) amplitude.2• The point is that the t dependence of ujis just ξ raised to the nthpower. So if |ξ(k)| > 1for some k, scheme is unstable. ξ is called the amplification factor.• Substitute (7) into (6), divide by ξn, get ( E.F.T.S.):ξ(k) = 1 − iv∆t∆xsin k∆x.Note |ξ(k)| > 1 for all k. ∴ FTCS is unconditionally unstable. Too bad. Simple schemegives garbage.Lax scheme• How do we fix it?• Replace fo r ward Euler time derivative:∂u∂t→un+1j−12(unj−1+ unj+1)∆t,where we have substituted the average value of unj−1and unj+1for unj.2Formally, the eigenmodes can be obtained from Fourier analysis of the finite-difference equations, butthis is beyond our scope.4• Schematically:new pointknown pointstx−1 +1j jj+1nn• FDE becomesun+1j=12(unj−1+ unj+1) −v∆t2∆x(unj+1− unj−1), (8)called the Lax scheme.• von Neumann stability analysis of (8) gives (E.F.T.S.)ξ(k) = cos k∆x − iv∆t∆xsin k∆x,which, for |ξ(k)| ≤ 1, requires|v|∆t∆x≤ 1. (9)• Equation (9) is the Courant condition (or CFL condition, for Courant-Friedrichs-Lewy).• Intuitively, the Courant condition can be thought of as limiting domain over whichinformation can propagate in one timestep to be less than one gridzone, i.e., ∆x ≥|v|∆t:∆x+v −v+v −vUnstable∆Stablet• Simple change in t derivative makes FTCS stable. Why? Write (8) in form of (6) withremainder term:un+1j− unj∆t= −vunj+1− unj−12∆x+12unj+1− 2unj+ unj−1∆t.5But this is just FTCS r epresentatio n of∂u∂t= −v∂u∂x+(∆x)22∆t∂2u∂x2|{z }diffusion term.• Adding diffusion stabilizes scheme: diffusion damps short wavelengths (k∆x ∼ 1),leaves large wavelengths unaffected. This is called numerical dissipation or numericalviscosity.• Damping short scales not as bad a s