AM 205: lecture 16 Last time: hyperbolic PDEs Today: parabolic and - PowerPoint PPT Presentation

AM 205: lecture 16 ◮ Last time: hyperbolic PDEs ◮ Today: parabolic and elliptic PDEs, introduction to optimization

The Wave Equation We now briefly return to the wave equation: u tt − c 2 u xx = 0 In one spatial dimension, this models, say, vibrations in a taut string

The Wave Equation Many schemes have been proposed for the wave equation One good option is to use central difference approximations 1 for both u tt and u xx : U n +1 j + U n − 1 − 2 U n − c 2 U n j +1 − 2 U n j + U n j j j − 1 = 0 ∆ t 2 ∆ x 2 Key points: ◮ Truncation error analysis = ⇒ second-order accurate ◮ Fourier stability analysis = ⇒ stable for 0 ≤ c ∆ t / ∆ x ≤ 1 ◮ Two-step method in time, need a one-step method to “get started” 1 Can arrive at the same result by discretizing the equivalent first order system

Parabolic PDEs

The Heat Equation The canonical parabolic equation is the heat equation u t − α u xx = f ( t , x ) , where α models thermal diffusivity In this section, we shall omit α for convenience Note that this is an Initial-Boundary Value Problem: ◮ We impose an initial condition u (0 , x ) = u 0 ( x ) ◮ We impose boundary conditions on both sides of the domain

The Heat Equation A natural idea would be to discretize u xx with a central difference, and employ the Euler method in time: U n +1 − U n U n j − 1 − 2 U n j + U n j +1 j j − = 0 ∆ x 2 ∆ t Or we could use backward Euler in time: U n +1 U n +1 j − 1 − 2 U n +1 + U n +1 − U n j j +1 j j − = 0 ∆ x 2 ∆ t

The Heat Equation Or we could do something “halfway in between”: U n +1 − U n U n +1 j − 1 − 2 U n +1 + U n +1 U n j − 1 − 2 U n j + U n − 1 − 1 j j +1 j +1 j j = 0 ∆ t 2 ∆ x 2 2 ∆ x 2 This is called the Crank–Nicolson method 2 In fact, it is common to consider a 1-parameter “family” of methods that include all of the above: the θ -method U n +1 − U n − θ U n +1 j − 1 − 2 U n +1 + U n +1 − (1 − θ ) U n j − 1 − 2 U n j + U n j j j j +1 j +1 = 0 ∆ t ∆ x 2 ∆ x 2 where θ ∈ [0 , 1] 2 From a paper by Crank and Nicolson in 1947, note: “Nicolson” is not a typo!

The Heat Equation With the θ -method: ◮ θ = 0 = ⇒ Euler ◮ θ = 1 2 = ⇒ Crank–Nicolson ◮ θ = 1 = ⇒ backward Euler For the θ -method, we can 1. perform Fourier stability analysis 2. calculate the truncation error

The θ -Method: Stability j ( k ) = λ ( k ) n e ik ( j ∆ x ) to get Fourier stability analysis: Set U n λ ( k ) = 1 − 4(1 − θ ) µ sin 2 � 1 � 2 k ∆ x 1 + 4 θµ sin 2 � 1 � 2 k ∆ x where µ ≡ ∆ t / (∆ x ) 2 Here we cannot get λ ( k ) > 1, hence only concern is λ ( k ) < − 1 Let’s find conditions for stability, i.e. we want λ ( k ) ≥ − 1: � 1 � � � 1 �� 1 − 4(1 − θ ) µ sin 2 1 + 4 θµ sin 2 2 k ∆ x ≥ − 2 k ∆ x

The θ -Method: Stability Or equivalently: � 1 � 4 µ (1 − 2 θ ) sin 2 2 k ∆ x ≤ 2 For θ ∈ [0 . 5 , 1] this inequality is always satisfied, hence the θ -method is unconditionally stable (i.e. stable independent of µ ) In the θ ∈ [0 , 0 . 5) case, the “most unstable” Fourier mode is when k = π/ ∆ x , since this maximizes the factor sin 2 � 1 � 2 k ∆ x

The θ -Method: Stability Note that this corresponds to the highest frequency mode that can be represented on our grid, since with k = π/ ∆ x we have e ik ( j ∆ x ) = e π ij = ( e π i ) j = ( − 1) j The k = π/ ∆ x mode: 2 1.5 1 0.5 e πij 0 −0.5 −1 −1.5 −2 0 1 2 3 4 5 6 7 8 9 10 j

The θ -Method: Stability This “sawtooth” mode is stable (and hence all modes are stable) if 1 4 µ (1 − 2 θ ) ≤ 2 ⇐ ⇒ µ ≤ 2(1 − 2 θ ) , Hence for θ ∈ [0 , 0 . 5), the θ -method is conditionally stable

The θ -Method: Stability 45 40 35 30 25 µ 20 15 10 5 0 0 0.1 0.2 0.3 0.4 0.5 θ For θ ∈ [0 , 0 . 5), θ -method is stable if µ is in the “green region,” i.e. approaches unconditional stability as θ → 0 . 5

The θ -Method: Stability Note that if we set θ to a value in [0 , 0 . 5), then stability time-step (∆ x ) 2 restriction is quite severe: ∆ t ≤ 2(1 − 2 θ ) Contrast this to the hyperbolic case where we had ∆ t ≤ ∆ x c This is an indication that the system of ODEs that arise from spatially discretizing the heat equation are stiff

The θ -Method: Accuracy The truncation error analysis is fairly involved, hence we just give the result: u n +1 − u n − θ u n +1 j − 1 − 2 u n +1 + u n +1 − (1 − θ ) u n j − 1 − 2 u n j + u n j j j j +1 j +1 T n ≡ j ∆ x 2 ∆ x 2 ∆ t �� 1 � � ∆ tu xxt − 1 12(∆ x ) 2 u xxxx = [ u t − u xx ] + 2 − θ � 1 24(∆ t ) 2 u ttt − 1 � 8(∆ t ) 2 u xxtt + � 1 � 1 � � ∆ t (∆ x ) 2 u xxxxt − 2 6!(∆ x ) 4 u xxxxxx + 2 − θ + · · · 12 The term u t − u xx in T n j vanishes since u solves the PDE

The θ -Method: Accuracy Key point: This is a first order method, unless θ = 1 / 2, in which case we get a second order method! θ -method gives us consistency (at least first order) and stability (assuming ∆ t is chosen appropriately when θ ∈ [0 , 1 / 2)) Hence, from Lax Equivalence Theorem, the method is convergent

The Heat Equation Note that the heat equation models a diffusive process, hence it tends to smooth out discontinuities Python demo: Heat equation with discontinous initial condition 1 0.8 0.6 0.4 0.2 0 3 10 2 8 6 1 4 2 0 0 t x This is very different to hyperbolic equations, e.g. the advection equation will just transport a discontinuity in u 0

Elliptic PDEs

Elliptic PDEs The canonical elliptic PDE is the Poisson equation In one-dimension, for x ∈ [ a , b ], this is − u ′′ ( x ) = f ( x ) with boundary conditions at x = a and x = b We have seen this problem already: Two-point boundary value problem! (Recall that Elliptic PDEs model steady-state behavior, there is no time-derivative)

Elliptic PDEs In order to make this into a PDE, we need to consider more than one spatial dimension Let Ω ⊂ R 2 denote our domain, then the Poisson equation for ( x , y ) ∈ Ω is u xx + u yy = f ( x , y ) This is generally written more succinctly as ∆ u = f We again need to impose boundary conditions (Dirichlet, Neumann, or Robin) on ∂ Ω (recall ∂ Ω denotes boundary of Ω)

Elliptic PDEs We will consider how to use a finite difference scheme to approximate this 2D Poisson equation First, we introduce a uniform grid to discretize Ω 1 0.9 0.8 0.7 0.6 y 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x

Elliptic PDEs Let h = ∆ x = ∆ y denote the grid spacing Then, ◮ x i = ih , i = 0 , 1 , 2 . . . , n x − 1, ◮ y j = jh , j = 0 , 1 , 2 , . . . , n y − 1, ◮ U i , j ≈ u ( x i , y j ) Then, we need to be able to approximate u xx and u yy on this grid Natural idea: Use central difference approximation!

Elliptic PDEs We have u xx ( x i , y j ) = u ( x i − 1 , y j ) − 2 u ( x i , y j ) + u ( x i +1 , y j ) + O ( h 2 ) , h 2 and u yy ( x i , y j ) = u ( x i , y j − 1 ) − 2 u ( x i , y j ) + u ( x i , y j +1 ) + O ( h 2 ) , h 2 so that u xx ( x i , y j ) + u yy ( x i , y j ) = u ( x i , y j − 1 ) + u ( x i − 1 , y j ) − 4 u ( x i , y j ) + u ( x i +1 , y j ) + u ( x i , y j +1 ) + O ( h 2 ) h 2

Elliptic PDEs Hence we define our approximation to the Laplacian as U i , j − 1 + U i − 1 , j − 4 U i , j + U i +1 , j + U i , j +1 h 2 This corresponds to a “5-point stencil”

Elliptic PDEs As usual, we represent the numerical solution as a vector U ∈ R n x n y We want to construct a differentiation matrix D 2 ∈ R n x n y × n x n y that approximates the Laplacian Question: How many non-zero diagonals will D 2 have? To construct D 2 , we need to be able to relate the entries of the vector U to the “2D grid-based values” U i , j

Elliptic PDEs Hence we need to number the nodes from 1 to n x n y — we number nodes along the “bottom row” first, then second bottom row, etc Let G denote the mapping from the 2D indexing to the 1D indexing. From the above figure we have: G ( i , j ; n x ) = jn x + i , and hence U G ( i , j ; n x ) = U i , j

Elliptic PDEs Let us focus on node ( i , j ) in our F.D. grid, this corresponds to entry G ( i , j ; n x ) of U Due to the 5-point stencil, row G ( i , j ; n x ) of D 2 will only have non-zeros in columns G ( i , j − 1; n x ) = G ( i , j ; n x ) − n x , (1) G ( i − 1 , j ; n x ) = G ( i , j ; n x ) − 1 , (2) G ( i , j ; n x ) = G ( i , j ; n x ) , (3) G ( i + 1 , j ; n x ) = G ( i , j ; n x ) + 1 , (4) G ( i , j + 1; n x ) = G ( i , j ; n x ) + n x (5) ◮ (2), (3), (4), give the same tridiagonal structure that we’re used to from differentiation matrices in 1D domains ◮ (1), (5) give diagonals shifted by ± n x

Elliptic PDEs For example, sparsity pattern of D 2 when n x = n y = 6

Elliptic PDEs Python demo: Solve the Poisson equation − ( x − 0 . 25) 2 − ( y − 0 . 5) 2 � � ∆ u = − exp , for ( x , y ) ∈ Ω = [0 , 1] 2 with u = 0 on ∂ Ω 1 0.06 0.9 0.055 0.8 0.05 0.7 0.045 0.6 0.04 0.5 0.035 y 0.4 0.03 0.3 0.025 0.2 0.02 0.1 0.015 0 0.01 0 0.2 0.4 0.6 0.8 1 x

Nonlinear Equations and Optimization