Discrete Poisson equation

From CFD-Wiki

In mathematics, the discrete Poisson equation is the finite difference analog of the Poisson equation. In it, the discrete Laplace operator takes the place of the Laplace operator. The discrete Poisson equation is frequently used in numerical analysis as a stand-in for the continuous Poisson equation, although it is also studied in its own right as a topic in discrete mathematics.

On a two-dimensional rectangular grid

Using the finite difference numerical method to discretize the 2 dimensional Poisson equation (assuming a uniform spatial discretization) on an m x n grid gives the following formula:

$( {\nabla}^2 u )_{ij} = \frac{1}{dx^2} ( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{ij}) = g_{ij}$

where $2 \le i \le m-1$ and $2 \le j \le n-1$ . The preferred arrangement of the solution vector is to use natural ordering which, prior to removing boundary elements, would look like:

$\begin{bmatrix} U \end{bmatrix} = \begin{bmatrix} u_{11} , u_{21} , \ldots , u_{m1} , u_{12} , u_{22} , \ldots , u_{m3} , \ldots , u_{mn} \end{bmatrix}^{T}$

This will result in an mn x mn linear system:

$\begin{bmatrix} A \end{bmatrix} \begin{bmatrix} U \end{bmatrix} = \begin{bmatrix} b \end{bmatrix}$

where

$A = \begin{bmatrix} ~D & -I & ~0 & ~0 & ~0 & \ldots & ~0 \\ -I & ~D & -I & ~0 & ~0 & \ldots & ~0 \\ ~0 & -I & ~D & -I & ~0 & \ldots & ~0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ ~0 & \ldots & ~0 & -I & ~D & -I & ~0 \\ ~0 & \ldots & \ldots & ~0 & -I & ~D & -I \\ ~0 & \ldots & \ldots & \ldots & ~0 & -I & ~D \end{bmatrix}$

$I$ is the m x m identity matrix, and $D$ , also m x m , is given by:

$D = \begin{bmatrix} ~4 & -1 & ~0 & ~0 & ~0 & \ldots & ~0 \\ -1 & ~4 & -1 & ~0 & ~0 & \ldots & ~0 \\ ~0 & -1 & ~4 & -1 & ~0 & \ldots & ~0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ ~0 & \ldots & ~0 & -1 & ~4 & -1 & ~0 \\ ~0 & \ldots & \ldots & ~0 & -1 & ~4 & -1 \\ ~0 & \ldots & \ldots & \ldots & ~0 & -1 & ~4 \end{bmatrix}$

^[1] For each $u_{ij}$ equation, the columns of $D$ correspond to the $u$ components:

$\begin{bmatrix} u_{1j} , & u_{2j} , & \ldots, & u_{i-1,j} , & u_{ij} , & u_{i+1,j} , & \ldots , & u_{mj} \end{bmatrix}^{T}$

while the columns of $I$ to the left and right of $D$ correspond to the $u$ components:

$\begin{bmatrix} u_{1,j-1} , & u_{2,j-1} , & \ldots, & u_{i-1,j-1} , & u_{i,j-1} , & u_{i+1,j-1} , & \ldots , & u_{m,j-1} \end{bmatrix}^{T}$

and

$\begin{bmatrix} u_{1,j+1} , & u_{2,j+1} , & \ldots, & u_{i-1,j+1} , & u_{i,j+1} , & u_{i+1,j+1} , & \ldots , & u_{m,j+1} \end{bmatrix}^{T}$

respectively.

From the above, it can be inferred that there are $n$ block columns of $m$ in $A$ . It is important to note that prescribed values of $u$ (usually lying on the boundary) would have their corresponding elements removed from $I$ and $D$ . For the common case that all the nodes on the boundary are set, we have $2 \le i \le m - 1$ and $2 \le j \le n - 1$ , and the system would have the dimensions (m - 2) (n - 2) x (m - 2) (n - 2) , where $D$ and $I$ would have dimensions (m-2) x (m-2) .

Example

For a 5 x 5 ( $m = 5$ and $n = 5$ ) grid with all the boundary nodes prescribed, the system would look like:

$\begin{bmatrix} U \end{bmatrix} = \begin{bmatrix} u_{22}, u_{32}, u_{42}, u_{23}, u_{33}, u_{43}, u_{24}, u_{34}, u_{44} \end{bmatrix}^{T}$

with

$A = \begin{bmatrix} ~4 & -1 & ~0 & -1 & ~0 & ~0 & ~0 & ~0 & ~0 \\ -1 & ~4 & -1 & ~0 & -1 & ~0 & ~0 & ~0 & ~0 \\ ~0 & -1 & ~4 & ~0 & ~0 & -1 & ~0 & ~0 & ~0 \\ -1 & ~0 & ~0 & ~4 & -1 & ~0 & -1 & ~0 & ~0 \\ ~0 & -1 & ~0 & -1 & ~4 & -1 & ~0 & -1 & ~0 \\ ~0 & ~0 & -1 & ~0 & -1 & ~4 & ~0 & ~0 & -1 \\ ~0 & ~0 & ~0 & -1 & ~0 & ~0 & ~4 & -1 & ~0 \\ ~0 & ~0 & ~0 & ~0 & -1 & ~0 & -1 & ~4 & -1 \\ ~0 & ~0 & ~0 & ~0 & ~0 & -1 & ~0 & -1 & ~4 \end{bmatrix}$

and

$b = \begin{bmatrix} dx^2 g_{22} + u_{12} + u_{21} \\ dx^2 g_{32} + u_{31} ~~~~~~~~ \\ dx^2 g_{42} + u_{52} + u_{41} \\ dx^2 g_{23} + u_{13} ~~~~~~~~ \\ dx^2 g_{33} ~~~~~~~~~~~~~~~~ \\ dx^2 g_{43} + u_{53} ~~~~~~~~ \\ dx^2 g_{24} + u_{14} + u_{25} \\ dx^2 g_{34} + u_{35} ~~~~~~~~ \\ dx^2 g_{44} + u_{54} + u_{45} \\ \end{bmatrix}$

As can be seen, the boundary $u$ 's are brought to the right-hand-side of the equation. ^[2] The entire system is 9 x 9 while $D$ and $I$ are 3 x 3 and given by:

$D = \begin{bmatrix} ~4 & -1 & ~0 \\ -1 & ~4 & -1 \\ ~0 & -1 & ~4 \\ \end{bmatrix}$

and

$-I = \begin{bmatrix} -1 & ~0 & ~0 \\ ~0 & -1 & ~0 \\ ~0 & ~0 & -1 \end{bmatrix}$

Methods of Solution

Because $\begin{bmatrix} A \end{bmatrix}$ is block tridiagonal and sparse, many methods of solution have been developed to optimally solve this linear system for $\begin{bmatrix} U \end{bmatrix}$ . Among the methods are a generalized Thomas algorithm, cyclic reduction, successive overrelaxation, and Fourier transforms. A theoretically optimal $O(n)$ solution can be computed using multigrid methods.

Applications

In Computational fluid dynamics, for the solution of an incompressible flow problem, the incompressiblity condition acts as a constraint for the pressure. There is no explicit form available for pressure in this case due to a strong coupling of the velocity and pressure fields. In this condition, by taking the divergence of all terms in the momentum equation, one obtains the pressure poisson equation. For an incompressible "Flow" this constraint is given by:

$\frac{ \partial v_x }{ \partial x} + \frac{ \partial v_y }{ \partial y} + \frac{\partial v_z}{\partial z} = 0$

where $v_x$ is the velocity in the $x$ direction, $v_y$ is velocity in $y$ and $v_z$ is the velocity in the $z$ direction. Taking divergence of the momentum equation and using the incompressibility constraint, the pressure poisson equation is formed given by:

$\nabla^2 p = f(\nu,V)$

where $\nu$ is the kinematic viscosity of the fluid and $V$ is the velocity vector. ^[3].

Footnotes

↑ Golub, Gene H. and C.F. Van Loan, Matrix Computations, 3rd Ed., The Johns Hopkins University Press, Baltimore, 1996, pages 177-180.
↑ Cheny, Ward and David Kincaid, Numerical Mathematics and Computing 2nd Ed., Brooks/Cole Publishing Company, Pacific Grove, 1985, pages 443-448
↑ Fletcher, Clive A. J., Computational Techniques for Fluid Dynamics: Vol I, 2nd Ed., Springer-Verlag, Berlin, 1991, page 334-339.

References

Cheny, Ward and David Kincaid, Numerical Mathematics and Computing 2nd Ed., Brooks/Cole Publishing Company, Pacific Grove, 1985.

Golub, Gene H. and C.F. Van Loan, Matrix Computations, 3rd Ed., The Johns Hopkins University Press, Baltimore, 1996.

Hoffman, Joe D., Numerical Methods for Engineers and Scientists, 4th Ed., McGraw-Hill Inc., New York, 1992.

Sweet, Roland A., SIAM Journal on Numerical Analysis, Vol. 11, No. 3 , June 1974, 506-520.

[0] Golub, Gene H. and C.F. Van Loan, Matrix Computations, 3rd Ed., The Johns Hopkins University Press, Baltimore, 1996, pages 177-180.

[1] Cheny, Ward and David Kincaid, Numerical Mathematics and Computing 2nd Ed., Brooks/Cole Publishing Company, Pacific Grove, 1985, pages 443-448

[2] Fletcher, Clive A. J., Computational Techniques for Fluid Dynamics: Vol I, 2nd Ed., Springer-Verlag, Berlin, 1991, page 334-339.

[1]

[2]

[3]

Discrete Poisson equation

From CFD-Wiki

Contents

On a two-dimensional rectangular grid

Example

Methods of Solution

Applications

Footnotes

References

Views

My wiki

wiki navigation

wiki search

wiki toolbox