Discrete Poisson equation

Using the finite difference numerical method to discretize

the 2-dimensional Poisson equation (assuming a uniform spatial discretization, $\backslash Delta\; x=\backslash Delta\; y$) on an {{math|m × n}} grid gives the following formula:{{citation|title=Numerical Methods for Engineers and Scientists| edition=2nd| first=Joe | last=Hoffman|year=2001|chapter=Chapter 9. Elliptic partial differential equations|publisher=McGraw–Hill| isbn=0-8247-0443-6}}.

$$(\; \{\backslash nabla\}^2\; u\; )\_\{ij\}\; =\; \backslash frac\{1\}\{\backslash Delta\; x^2\}\; (u\_\{i+1,j\}\; +\; u\_\{i-1,j\}\; +\; u\_\{i,j+1\}\; +\; u\_\{i,j-1\}\; -\; 4\; u\_\{ij\})\; =\; g\_\{ij\}$$

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

$$\backslash mathbf\{u\}\; =$$

\begin{bmatrix} u_{11} , u_{21} , \ldots , u_{m1} , u_{12} , u_{22} , \ldots , u_{m2} , \ldots , u_{mn}

\end{bmatrix}^\mathsf{T}

This will result in an {{math|mn × mn}} linear system:

$$A\backslash mathbf\{u\}\; =\; \backslash mathbf\{b\}$$

where

$$A\; =$$

\begin{bmatrix}

~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{bmatrix},

$I$ is the {{math|m × m}} identity matrix, and $D$, also {{math|m × m}}, is given by:Golub, Gene H. and C. F. Van Loan, Matrix Computations, 3rd Ed., The Johns Hopkins University Press, Baltimore, 1996, pages 177–180.

$$D\; =$$

\begin{bmatrix}

~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{bmatrix},

and $\backslash mathbf\{b\}$ is defined by

$$\backslash mathbf\{b\}\; =$$

-\Delta x^2 \begin{bmatrix} g_{11} , g_{21} , \ldots , g_{m1} , g_{12} , g_{22} , \ldots , g_{m2} , \ldots , g_{mn} \end{bmatrix}^\mathsf{T}.

For each $u\_\{ij\}$ equation, the columns of $D$ correspond to a block of $m$ components in $u$:

$$\backslash begin\{bmatrix\}$$

u_{1j} , & u_{2j} , & \ldots, & u_{i-1,j} , & u_{ij} , & u_{i+1,j} , & \ldots , & u_{mj}

\end{bmatrix}^\mathsf{T}

while the columns of $I$ to the left and right of $D$ each correspond to other blocks of $m$ components within $u$:

\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}^\mathsf{T}

and

$$\backslash 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}^\mathsf{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\; \backslash le\; i\; \backslash le\; m\; -\; 1$ and $2\; \backslash le\; j\; \backslash le\; n\; -\; 1$, and the system would have the dimensions {{math|(m − 2)(n − 2) × (m− 2)(n − 2)}}, where $D$ and $I$ would have dimensions {{math|(m − 2) × (m − 2)}}.

For a 3×3 ( $m\; =\; 3$ and $n\; =\; 3$ ) 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}^\mathsf{T}

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{b} =

\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 $u$'s are brought to the right-hand-side of the equation.Cheny, Ward and David Kincaid, Numerical Mathematics and Computing 2nd Ed., Brooks/Cole Publishing Company, Pacific Grove, 1985, pages 443–448. The entire system is {{math|9 × 9}} while $D$ and $I$ are {{math|3 × 3}} and given by:

$$D\; =$$

\begin{bmatrix}

~4 & -1 & ~0 \\

-1 & ~4 & -1 \\

~0 & -1 & ~4 \\

\end{bmatrix}

and

$$-I\; =\; \backslash begin\{bmatrix\}$$

-1 & ~0 & ~0 \\

~0 & -1 & ~0 \\

~0 & ~0 & -1

\end{bmatrix}.

Because $\backslash begin\{bmatrix\}\; A\; \backslash end\{bmatrix\}$ is block tridiagonal and sparse, many methods of solution

have been developed to optimally solve this linear system for $\backslash begin\{bmatrix\}\; U\; \backslash 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(n\; \backslash log(n))$. An optimal $O(n)$ solution can also be computed using multigrid methods.CS267: Notes for Lectures 15 and 16, Mar 5 and 7, 1996, https://people.eecs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html

File:Convergence of Iterative Numerical Methods for Poisson System with 16384 elements.svg

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:

$$\backslash frac\{\; \backslash partial\; v\_x\; \}\{\; \backslash partial\; x\}\; +\; \backslash frac\{\; \backslash partial\; v\_y\; \}\{\; \backslash partial\; y\}\; +\; \backslash frac\{\backslash partial\; v\_z\}\{\backslash 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:

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

where $\backslash nu$ is the kinematic viscosity of the fluid and $V$ is the velocity vector.

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

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. S. P. Meyn and R.L. Tweedie, 2005. [http://probability.ca/MT/ Markov Chains and Stochastic Stability].

Second edition to appear, Cambridge University Press, 2009.S. P. Meyn, 2007. [http://www.meyn.ece.ufl.edu/archive/spm_files/CTCN/CTCN.html Control Techniques for Complex Networks] {{Webarchive|url=https://web.archive.org/web/20141216122500/http://www.meyn.ece.ufl.edu/archive/spm_files/CTCN/CTCN.html |date=December 16, 2014 }}, 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.

{{reflist|30em}}

• 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.
• {{Cite book | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press | publication-place=New York | isbn=978-0-521-88068-8 | chapter=Section 20.4. Fourier and Cyclic Reduction Methods | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=1053 | access-date=August 18, 2011 | archive-date=August 11, 2011 | archive-url=https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=1053 | url-status=dead }}

Category:Finite differences

Category:Numerical differential equations