Discrete Poisson equation explained

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 discretizethe 2-dimensional Poisson equation (assuming a uniform spatial discretization,

\Deltax=\Deltay

) on an grid gives the following formula:^[1]

(^2 u)_ = \frac (u_ + u_ + u_ + u_ - 4 u_) = g_

where

2\lei\lem-1

and

2\lej\len-1

. The preferred arrangement of the solution vector is to use natural ordering which, prior to removing boundary elements, would look like:

\mathbf =\begin u_, u_, \ldots, u_, u_, u_, \ldots, u_, \ldots, u_\end^\mathsf

This will result in an linear system: $A\mathbf = \mathbf$ where $A =\begin ~D & -I & ~0 & ~0 & ~0 & \cdots & ~0 \\ -I & ~D & -I & ~0 & ~0 & \cdots & ~0 \\ ~0 & -I & ~D & -I & ~0 & \cdots & ~0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ ~0 & \cdots & ~0 & -I & ~D & -I & ~0 \\ ~0 & \cdots & \cdots & ~0 & -I & ~D & -I \\ ~0 & \cdots & \cdots & \cdots & ~0 & -I & ~D\end,$

is the identity matrix, and

, also, is given by:^[2]

D =\begin ~4 & -1 & ~0 & ~0 & ~0 & \cdots & ~0 \\ -1 & ~4 & -1 & ~0 & ~0 & \cdots & ~0 \\ ~0 & -1 & ~4 & -1 & ~0 & \cdots & ~0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ ~0 & \cdots & ~0 & -1 & ~4 & -1 & ~0 \\ ~0 & \cdots & \cdots & ~0 & -1 & ~4 & -1 \\ ~0 & \cdots & \cdots & \cdots & ~0 & -1 & ~4\end,

and

is defined by

\mathbf =-\Delta x^2 \begin g_, g_, \ldots, g_, g_, g_, \ldots, g_, \ldots, g_ \end^\mathsf.

For each

u_ij

equation, the columns of

correspond to a block of

components in

\begin u_, & u_, & \ldots, & u_, & u_, & u_, & \ldots, & u_\end^\mathsf

while the columns of

to the left and right of

each correspond to other blocks of

components within

\begin u_, & u_, & \ldots, & u_, & u_, & u_, & \ldots, & u_\end^\mathsf

and

\begin u_, & u_, & \ldots, & u_, & u_, & u_, & \ldots, & u_\end^\mathsf

respectively.

From the above, it can be inferred that there are

block columns of

. It is important to note that prescribed values of

(usually lying on the boundary) would have their corresponding elements removed from

and

. For the common case that all the nodes on the boundary are set, we have

2\lei\lem-1

and

2\lej\len-1

, and the system would have the dimensions, where

and

would have dimensions .

Example

For a 3×3 (

m=3

and

n=3

) grid with all the boundary nodes prescribed, the system would look like:

\begin U \end =\begin u_, u_, u_, u_, u_, u_, u_, u_, u_\end^\mathsf

with

A =\left[\begin{array}{ccc|ccc|ccc}
 ~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 \\
 \hline
 -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 \\
 \hline
 ~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{array}\right]

and

$\mathbf =\left[\begin{array}{l}
-\Delta x^2 g_{22} + u_{12} + u_{21} \\
-\Delta x^2 g_{32} + u_{31} ~~~~~~~~ \\
-\Delta x^2 g_{42} + u_{52} + u_{41} \\
-\Delta x^2 g_{23} + u_{13} ~~~~~~~~ \\
-\Delta x^2 g_{33} ~~~~~~~~~~~~~~~~ \\
-\Delta x^2 g_{43} + u_{53} ~~~~~~~~ \\
-\Delta x^2 g_{24} + u_{14} + u_{25} \\
-\Delta x^2 g_{34} + u_{35} ~~~~~~~~ \\
-\Delta x^2 g_{44} + u_{54} + u_{45}
\end{array}\right].$

As can be seen, the boundary

's are brought to the right-hand-side of the equation.^[3] The entire system is while

and

are and given by:

D =\begin ~4 & -1 & ~0 \\ -1 & ~4 & -1 \\ ~0 & -1 & ~4 \\\end

and

-I = \begin -1 & ~0 & ~0 \\ ~0 & -1 & ~0 \\ ~0 & ~0 & -1\end.

Methods of solution

Because

\begin{bmatrix}A\end{bmatrix}

is block tridiagonal and sparse, many methods of solutionhave been developed to optimally solve this linear system for

\begin{bmatrix}U\end{bmatrix}

.Among the methods are a generalized Thomas algorithm with a resulting computational complexity of

O(n^2.5)

, cyclic reduction, successive overrelaxation that has a complexity of

O(n^1.5)

, and Fast Fourier transforms which is

O(nlog(n))

. An optimal

O(n)

solution can also be computed using multigrid methods.^[4]

Applications

In computational fluid dynamics, for the solution of an incompressible flow problem, the incompressibility 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 + \frac + \frac = 0$ where

v_x

is the velocity in the

direction,

v_y

is velocity in

and

v_z

is the velocity in the

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

is the velocity vector.^[5]

The discrete Poisson's equation arises in the theory of Markov chains. It appears as the relative value function for the dynamic programming equation in a Markov decision process, and as the control variate for application in simulation variance reduction.^[6] ^[7] ^[8]

References

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.
Book: Press . WH . Teukolsky . SA . Vetterling . WT . Flannery . BP . 2007 . Numerical Recipes: The Art of Scientific Computing . 3rd . Cambridge University Press . New York . 978-0-521-88068-8 . Section 20.4. Fourier and Cyclic Reduction Methods . http://apps.nrbook.com/empanel/index.html#pg=1053 . August 18, 2011 . August 11, 2011 . https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=1053 . dead .

Notes and References

.
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.
CS267: Notes for Lectures 15 and 16, Mar 5 and 7, 1996, https://people.eecs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html
Fletcher, Clive A. J., Computational Techniques for Fluid Dynamics: Vol I, 2nd Ed., Springer-Verlag, Berlin, 1991, page 334–339.
S. P. Meyn and R.L. Tweedie, 2005. Markov Chains and Stochastic Stability. Second edition to appear, Cambridge University Press, 2009.
S. P. Meyn, 2007. Control Techniques for Complex Networks, Cambridge University Press, 2007.
Asmussen, Søren, Glynn, Peter W., 2007. "Stochastic Simulation: Algorithms and Analysis". Springer. Series: Stochastic Modelling and Applied Probability, Vol. 57, 2007.