Lecture 12Wednesday, May 11, 2005Supplementary Reading: CS 205 notes on Poisson’s equation1 Laplace EquationIn 1D the Laplace equation is given bypxx= 0.The solution isp = ax + b.for some constants a and b. In order to find a and b, we need two boundaryconditions. Dirichlet boundary conditions specify the value of p at the boundary,e.g.,p(0) = 0p(1) = 1⇒ p(x) = xNeumann boundary conditions specify the derivatives of the function at theboundary. For example, we might have a Neumann boundary condition atx = 0 and a Dirichlet boundary condition at x = 1,px(0) = 0p(1) = 1⇒ p(x) = 1Recall from the previous lecture that if both boundary conditions are of Neu-mann type, then we either have an inconsistent problem or an underdeterminedproblem. In the case where we have an underdetermined problem, we still donot know the value of the constant b. However, if only pxis actually needed inthe computation, this may be ok.2 Numerical SolutionMore generally, we are interested in numerically solving Poisson’s equationpxx= f(x).1At each grid node, we approximate the equation using the second order centraldifference schemepi+1− 2pi+ pi−14x2= fi.The re sult is a coupled linear system that we need to solve in order to determinep on the entire domain. However, we cannot write this equation as is for thegrid points near the boundary since it will involve points outside of the domain.For example, assume that our domain is the interval [0, 1] and that we have gridpoints 0, 1, . . . , M, M + 1 uniformly spaced on the domain. The equation for p1isp2− 2p1+ p04x2= f1.If we have a Dirichlet boundary condition specified on the left of the domainp(0) = β,then the equation for p1becomesp2− 2p14x2= f1−β4x2.If we have a Neumann boundary condition specified at the half grid point12px(x12) = α,we write the equation for p1asp2−p14x−p1−p04x4x= f1.Sincepx(x12) =p1− p04x+ O(4x2),the equation for p1becomesp2− p14x2= f1+α4x.Let’s look at the matrix equation for the case where we have two Dirichletboundary conditions.−2 11 −2 1.........1 −2 1.........1 −2 11 −2p1...pi−1pipi+1...pM=4x2f1− p04x2f2...4x2fi...4x2fM−14x2fM− pM+12The matrix is symmetric negative definite. This is advantageous bec ausethere are fast linear solvers for such systems, e.g. the conjugate gradients method.In the case with two Neumann boundary conditions, the matrix equation is−1 11 −2 1.........1 −2 1.........1 −2 11 −1p1...pi−1pipi+1...pM=4x2f1− 4x (px)124x2f2...4x2fi...4x2fM−14x2fM− 4x (px)M+12Notice that the matrix has changed. In particular, it is singular since it hasa non-empty null space which is spanned by the vector (1, . . . , 1)T. This isproblematic, but workable. It can be solved for p up to a constant, since for anysolution, ~p, ~p + c(1, . . . , 1)Tis also a solution.In multiple dimension Poisson’s equation is4p = f.In 2D the equation ispxx+ py y= f.We use the second order accurate central difference discretizationpi+1,j− 2pi,j+ pi−1,j4x2+pi,j+1− 2pi,j+ pi,j−14y2= fi,j.In 2D we need boundary conditions specified around the entire domain. Ifat least one boundary condition is Dirichlet, then the resulting matrix will bea banded symmetric positive definite matrix. We can use an iterative solversuch as preconditioned conjugate gradients. If all the boundary conditions areNeumann, then the m atrix will have a null space.3 Computing the PressureThe Laplace and Poisson equations are elliptic partial differential equations.Recall that in solving the Navier-Stokes equations using the projection methodwe derive an elliptic equation for the pressure. To review, the steps in thesolution of the m omentum equation are1. Compute the intermediate velocity field~V?~V?−~Vn4t+~Vn· ∇~Vn= ~g (1)32. Solve an elliptic equation for the pressure4ˆp = ∇ ·~V?(2)3. Compute the divergence free velocity field~Vn+1~Vn+1−~V?+ ∇ˆp = 0 (3)Here we are concerned with ste p 2, the elliptic solve for the pressure. Specif-ically, we address the handling of boundary conditions.Boundary conditions can be applied to either the velocity or the pressure.In order to apply boundary conditions to~Vn+1, we apply them to~V?aftercomputing~V?in equation 1 and before solving equation 2. Then in equation2, we set ∇p ·~N = 0 on the boundary where~N is the local unit normal to theboundary. Note that due to the relation in equation 3, this will result in thecorrect boundary condition for~Vn+1.Recall that the Neumann problem for Poisson’s equation must satisfy thecompatibility condition for a solution to exist. The problem is given by4p = f in Ω∇p · N = g on ∂Ωwhere~N is the unit normal to the boundary. From the equation we have therelationsZΩf dV =ZΩ4p dV =ZΩ∇ · ∇p dV =Z∂Ω∇p · N dS =Z∂Ωg dSwhere the third equality follows from the divergence theorem. The compatibilitycondition isZΩf dV =Z∂Ωg dSIn solving equation 2, f = ∇ ·~V?and g = 0. Therefore, the compatibilitycondition isZΩ∇ ·~V?dV =Z∂Ω~V?·~N dS = 0 (4)where the first equality follows from the divergence theorem. This conditionneeds to be satisfied when specifying the boundary condition on~V?in order toguarantee the existence of a