# Poisson equation, numerical methods

Methods replacing the original boundary value problem for the Poisson equation

$$\tag{1 } \Delta u ( x) \equiv \ \sum _ { r= } 1 ^ { d } \frac{\partial ^ {2} u ( x) }{\partial x _ {r} ^ {2} } = f ( x) ,\ x \equiv ( x _ {1} \dots x _ {d} ) ,$$

by a system of $N$ linear algebraic equations

$$\tag{2 } L _ {N} ( \mathbf u _ {N} ) = \mathbf f _ {N} ,$$

with a solution $\mathbf u _ {N} \equiv ( u _ {1} \dots u _ {n} )$ which makes it possible to construct a certain approximation $p _ {N} \mathbf u ^ {N}$ to the solution of the original problem as $N \rightarrow \infty$.

One defines such extremely important concepts as the error of the numerical method and bounds for the error (the accuracy) depending on the way one compares the solutions of the original problem (1) and the discrete problem (2). The algebraic properties of the system (2) (discrete analogues of boundary value problems) connected with the stability of solutions (the well-posedness of the discrete problems), the possibility of finding exact or approximate solutions of (2) by some direct or iterative method when carrying out the corresponding computation and corresponding demands on the computer memory (see Minimization of the labour of calculation) give other characteristics of the numerical methods.

Numerical solutions of boundary value problems for the Poisson equation are important not only because these problems often arise in diverse branches of science and technology, but because they frequently are a means for solving more general boundary value problems for both equations and systems of equations of elliptic type as well as for various non-stationary systems. The basic numerical methods for solving the boundary value problems under discussion are projection methods and difference methods (see [1][13]).

## Projection methods.

These encompass a number of methods: variational, least squares, Galerkin, projection-difference, projection-grid, and finite elements methods. In all of these one characteristically reduces the original boundary value problem to an operator equation

$$\tag{3 } L ( u) = f$$

(the operator $L$ acts, for example, from a Hilbert space $H$ into $H$) with subsequent choice of finite-dimensional subspaces $H _ {N}$ and $F _ {N}$( $N \rightarrow \infty$): the problem (3) itself is replaced in these methods by the problem of finding a $\widehat{u} _ {N} \in H _ {N}$ such that for any $v \in F _ {N}$,

$$( L \widehat{u} _ {N} - f , v ) _ {H} = 0 .$$

Then, given bases in $H _ {N}$ and $F _ {N}$, the system (2) is a system for the coefficients in the expansion of $\widehat{u} _ {N}$ with respect to the basis of $H _ {N}$ and $p _ {N} \mathbf u _ {N}$ can be taken as the function $\widehat{u} _ {N}$ itself; it is natural to define the error in the method as $\| u - \widehat{u} _ {N} \| _ {H}$. In the most important cases $H$ is a certain subspace of the Sobolev space $W _ {2} ^ { \prime } ( \Omega )$, $H _ {N} = F _ {N}$ and, if $\psi _ {1} ( x) \dots \psi _ {N} ( x)$ is a basis of $H _ {N}$, then the system (2) takes the form

$$\tag{4 } \sum _ { j= } 1 ^ { N } u _ {j} \int\limits _ \Omega \mathop{\rm grad} \psi _ {j} ( x) \ \mathop{\rm grad} \psi _ {i} ( x) \ d x = \int\limits _ \Omega f ( x) \psi _ {i} ( x) d x ,$$

$$i = 1 \dots N .$$

The error in the method is estimated by the distance in $H$ from the solution of the original problem to the subspace $H _ {N}$( see [1], [5][11]). In modern variants of projection methods the subspaces $H _ {N}$ tend to be chosen so that the functions $\psi _ {i} ( x)$ have local supports and in each equation (4) only a finite number of coefficients are non-zero. Methods of this type are also called projection-grid methods (projection-difference, variational-difference, finite element methods) (see [1], [4], [7][11]). The great merit of these methods is that they can be applied when the domain $\Omega$ of the boundary value problem under consideration is geometrically fairly complicated. Methods which are also close to projection methods but relatively rarely applied are the collocation method, spectral methods and boundary element methods (see [13][15]).

## Difference (finite-difference) methods.

These use a certain grid domain $\Omega _ {N}$ containing $N$ nodes of a grid to approximate the original domain $\Omega$ and usually lead to a system (2) by approximating the Poisson equation and the corresponding boundary conditions by their difference (grid) analogues, using only the values of the function at the chosen nodes (see [1]). The error of the method is usually obtained by comparing the vector $\mathbf u _ {N}$ and the vector obtained by restricting the required solution to the set of nodes under consideration. The well-posedness and approximation can be studied for various choices of norms; it is, in particular, possible to use the maximum principle; convergence is obtained as a corollary of well-posedness and approximation (see [1][4]).

Systems (2) can be derived both by basing them on certain discrete analogues of the corresponding variational problems and by basing them on an approximation to certain integral relations (see [1], [2], [4], [11]); such approaches bring these variants of difference methods somewhat closer to projection-difference methods.

The most intensively studied methods for solving systems of grid equations (2) are the simplest difference analogues on a parallelepiped grid (see [1], [10]). In the case of two variables and when the domain $\Omega$ in the plane is rectangular, direct methods are often applied for a number of boundary conditions; these direct methods enable one to find the solution of (2) at the cost of $O ( N \mathop{\rm ln} N )$ arithmetic operations. Such a method is that of separation of variables, which uses the discrete Fourier transform and the reduction method (cf. Separation of variables, method of); methods are also known with a bound $O ( N)$ on the numbers of operations (see [11][13]). When $d \geq 2$ and one may use separation of variables (in this case $\Omega$ is a parallelepiped), the solution of (2) can be found with accuracy $\epsilon > 0$ at the cost of $O ( N \mathop{\rm ln} N | \mathop{\rm ln} \epsilon | )$ operations by an iterative alternating-directions method (see [1], [11]); iterative methods with factorizable operators (methods of successive overrelaxation with symmetrization, incomplete matrix factorization, alternating-triangular methods) enable one to find the solution of (2) with accuracy $\epsilon$ at a cost of $O ( N ^ {1+} 1/2d | \mathop{\rm ln} \epsilon | )$ operations (see [11], [12]) in fairly general situations.

In the case of simply-connected and multiply-connected domains $\Omega$ in the plane which are composed of a finite number of rectangles, the solution of the system (2) with accuracy $\epsilon$ can be found at the cost of $O ( N \mathop{\rm ln} N + N ^ {3/4} | \mathop{\rm ln} \epsilon | \mathop{\rm ln} N )$ operations by cutting $\Omega$ up into rectangles (see [11], [13]). Similar asymptotic solutions for discrete analogues of the Dirichlet problem can be obtained for certain domains using the methods of capacity and of fictitious unknowns (see [11], [13]). For a number of systems (2) which are projection-difference analogues of the original problem, computation lengths of the order $O ( N \mathop{\rm ln} N | \mathop{\rm ln} \epsilon | )$ and sometimes also $O ( N \mathop{\rm ln} N )$( when $\epsilon$ is of order $N ^ {- \alpha }$, $\alpha > 0$) can be achieved by iterative methods using the spectral equivalence of operators (see [11], [13]). In a number of cases the use of a sequence of grids enables one to obtain methods which give the solution to (2) with an accuracy of order $N ^ {- \alpha }$, $\alpha > 0$, and with asymptotically-minimal computational work (the number of operations is $O ( N)$) (see for example [1], [10], [11]).

#### References

 [1] S.K. Godunov, V.S. Ryaben'kii, "The theory of difference schemes" , North-Holland (1964) (Translated from Russian) [2] O.A. Ladyzhenskaya, "The boundary value problems of mathematical physics" , Springer (1986) (Translated from Russian) [3] G.I. Marchuk, V.V. Shaidurov, "Difference methods and their extrapolations" , Springer (1983) (Translated from Russian) [4] A.A. Samarskii, V.B. Andreev, "Méthodes aux différences pour équations elliptiques" , MIR (1978) (Translated from Russian) [5] S.G. Mikhlin, "The numerical performance of variational methods" , Wolters-Noordhoff (1971) (Translated from Russian) [6] M.A. Krasnosel'skii, G.M. Vainikko, P.P. Zabreiko, et al., "Approximate solution of operator equations" , Wolters-Noordhoff (1972) (Translated from Russian) [7] J.P. Aubin, "Approximation of elliptic boundary value problems" , Wiley (1972) [8] G. Strang, J. Fix, "An analyse of the finite element method" , Prentice-Hall (1973) [9] P.G. Ciarlet, "The finite element method for elliptic problems" , North-Holland (1978) [10] W. Hackbush, "Multigrid methods and applications" , Springer (1985) [11] E.G. D'yakonov, "Minimization of computational work. Asymptotically-optimal algorithms" , Moscow (1989) (In Russian) [12] A.A. Samarskii, E.S. Nikolaev, "Numerical methods for grid equations" , 1–2 , Birkhäuser (1989) (Translated from Russian) [13] R. Glowinski (ed.) G.H. Golub (ed.) G.A. Meurant (ed.) J. Periaux (ed.) , First Internat. Symp. Domain Decomposition Methods For Partial Diff. Eqs. , SIAM (1988) [14] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, "Spectral methods in fluid dynamics" , Springer (1988) [15] P.K. Banerjee, R. Butterfield, "Boundary element methods in engineering science" , McGraw-Hill (1981)

During the last few years a new and very powerful class of methods has come up which have essentially an $O( N)$ complexity. These so-called multigrid methods are based on employing smoothing effects of higher frequencies on courser grids and damping of lower frequencies on finer grids. They must be considered to be among the most powerful tools available today (1990) and can be used for finite-difference and finite-element discretization alike.