Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators

Zaman, Mohammad Asif

doi:10.3390/electronics11152365

Open AccessArticle

Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators

by

Mohammad Asif Zaman

Department of Electrical Engineering, Stanford University, Stanford, CA 94305, USA

Electronics 2022, 11(15), 2365; https://doi.org/10.3390/electronics11152365

Submission received: 21 June 2022 / Revised: 24 July 2022 / Accepted: 25 July 2022 / Published: 28 July 2022

Download

Browse Figures

Versions Notes

Abstract

The Poisson equation frequently emerges in many fields of science and engineering. As exact solutions are rarely possible, numerical approaches are of great interest. Despite this, a succinct discussion of a systematic approach to constructing a flexible and general numerical Poisson solver can be difficult to find. In this introductory paper, a comprehensive discussion is presented on how to build a finite difference matrix solver that can solve the Poisson equation for arbitrary geometry and boundary conditions. The boundary conditions are implemented in a systematic way that enables easy modification of the solver for different problems. An image-based geometry-definition approach is also discussed. Python code of the numerical recipe is made publicly available. Numerical examples are presented that show how to set up the solver for different problems.

Keywords:

Poisson equation; matrix operators; finite difference method; partial differential equations; tutorial

1. Introduction

The Poisson equation is an elliptical partial differential equation (PDE) and is ubiquitous in many areas of physics and engineering. Two of its common uses include modeling electrostatic potential and gravitational potential. Some forms of the Poisson equation also appear in fluid flow [1] and heat transfer problems [2,3]. Often, analytic solutions of such problems are not possible for cases of practical interest due to complex geometries. Numerically solving the Poisson equation is the only way to study these systems. As such, numerical solution methods for the Poisson equation are of significant interest [2,4,5].

There are several methods for solving the Poisson equation numerically [6]. The Finite-Difference Method (FDM) is one of the most simple and popular approaches [7,8,9,10]. This method involves replacing the continuous derivative operators with approximate, discrete finite-difference operators [11] that take the form of matrices. The method is intuitive and the solution process involves well-established linear algebraic methods. It suffers from some drawbacks, such as accurately representing complex geometries and boundaries. Methods such as the Finite Element Method (FEM) are better at handling arbitrary geometries. However, due to the ease of implementation and the intuitive nature of FDM, it remains a very popular approach for solving PDEs.

There are several references that discuss solving the Poisson equation using FDM [12]. While these papers cover many details, they usually do not present a systematic approach of constructing the numerical system (i.e., system matrix formation) for arbitrary geometries and arbitrary boundary conditions. Thus, for every variation of the Poisson equation, the linear algebraic system has to be reconstructed from scratch. The goal of this paper is to introduce a systematic approach for solving the Poisson equation for any arbitrary geometry and boundary conditions with minimal modification of the numerical solver/computer code. This would allow one to study a variety of configurations of the physical system without having to reformulate the numerical recipe.

Setting up an FDM solver has some basic steps. The first step is discretization of the solution space. A series of grid points in space where the equation will be solved is defined. At each of these points, the Poisson equation is written in discrete form using difference equations. In the conventional approach, the boundary conditions are incorporated at this step by modifying the difference equations at boundary points. This step depends on the geometry and the type of the boundary conditions. After this adjustment, we are left with a system of linear equations from which a full rank matrix equation can be formed. This matrix is referred to as the system matrix. The system is then solved. Note that when the geometry or the boundary conditions are altered, the system matrix changes and thus has to be reconstructed. Most references do not provide the reader with a systematic approach to do this. In many cases, one is required to study different geometrical variations of the same problem. Thus, one is often left with the cumbersome task of modifying their numerical solver for each configuration.

In this paper, the system matrix is defined as the combination of two matrices. One matrix is independent of the problem geometry and thus remains unchanged for all cases. This is referred to as the derivative operator matrix. For the Poisson equation, this is the second-derivative operator. The other matrix is defined to be dependent on the geometry and the boundary conditions. This is referred to as the boundary operator matrix. Modifying the boundary operator matrix turns out to be relatively simple and intuitive. It is shown that the system matrix can be constructed by compiling rows from these matrices appropriately. Note that when the geometry and/or the boundary conditions change, only the boundary operator matrix and the system matrix construction procedure changes. Thus, different geometry variations (and/or boundary conditions) could be investigated with minor adjustments of the same numerical code. This is the key advantage of the proposed method.

In this paper, a complete numerical recipe of how to implement a flexible Poisson equation solver is presented in detail. In addition, a Python implementation of the algorithm is discussed and is made available to the readers through a GitHub repository [13]. Some example numerical results are also presented that showcase the capabilities of the solver.

The rest of the paper is organized as follows: Section 2 defines the Poisson equation and the types of boundary conditions typically encountered, Section 3 discusses the finite difference formulation for one-dimensional (1D) and two-dimensional (2D) cases; the computer implementation of the algorithm using Python is discussed in Section 4; Section 5 presents numerical results for a few example problems; and concluding remarks are made in Section 6.

2. The Poisson Equation and Boundary Conditions

The Poisson equation consists of three second-order spatial derivatives. The general form of the equation in Cartesian coordinates is:

Δ u (x, y, z) = f (x, y, z), on Ω .

(1)

where

Δ = \nabla^{2} = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}},

(2)

is the Laplacian operator in Cartesian coordinates, x, y, z are the independent Cartesian spatial coordinates, and

Ω

is the solution domain (i.e., the set of x, y, z values for which the equation is to be solved). The dependent variable, u, is the function that is to be solved, and f is known as the source function, which is usually given. In general, it is a function of the space coordinates. In the special case of

f = 0

, the equation reduces to the Laplace equation.

Like any differential equation, boundary conditions must be defined before a solution can be attempted. A complete definition of the problem includes Equation (1) with a defined source function and boundary conditions enforced on u or its first derivatives (i.e.,

\frac{\partial u}{\partial x}

,

\frac{\partial u}{\partial y}

, and

\frac{\partial u}{\partial z}

). Boundary conditions may also be enforced on both u and its derivatives simultaneously. The boundary of the solution domain

Ω

is denoted as

\partial Ω

. The boundary can be subdivided into M number of arbitrary segments or regions,

\partial Ω_{1}, \partial Ω_{2}, \dots, \partial Ω_{M}

, each having its own boundary conditions. Usually, four common types of boundary conditions are encountered:

Dirichlet Boundary Conditions, having the form: ${u |}_{_{\partial Ω_{i}}} = u_{i}$
Neumann Boundary Conditions, having the form: $\frac{\partial u}{\partial η} |_{_{\partial Ω_{i}}} = g_{i}$
Mixed Boundary Conditions, having the form: ${u |}_{_{\partial Ω_{i}}} = u_{i}$ , $\frac{\partial u}{\partial η} |_{_{\partial Ω_{j}}} = g_{j}$
Robin Boundary Conditions, having the form: $(A u + B \frac{\partial u}{\partial η}) |_{_{\partial Ω_{i}}} = g_{i}$

Here

η \in {x, y, z}

represents one of the Cartesian spatial variables,

i, j \in {1, 2, \dots, M}

, A, B are constants, and

u_{i}

,

g_{j}

are given quantities. Note that for Neumann boundary conditions, since the constraint is applied on the derivative, u can only be determined up to a constant. Hence, an additional constraint is required to obtain a unique solution. This constraint can have the following form [14,15]:

\int_{Ω} u d x = 0 .

(3)

The boundary regions can also be specified inside the simulation domain itself and not necessary defined on the outer perimeter (i.e., some

\partial Ω_{i}

among

\partial Ω_{1}, \partial Ω_{2}, \dots, \partial Ω_{M}

may not necessarily reside on the outer edges of

Ω

). This is more common for Dirichlet boundaries [11]. The solution scheme presented here makes no distinction between the inner and outer boundary regions. The method works for any arbitrary

\partial Ω_{i}

for all i.

3. Finite Difference Discretization and Matrix Operators

All the quantities in Equation (1) are continuous variables. When solving numerically, the equation is solved only at discrete points. The first step is to convert the variables and the equation from continuous domain to discrete domain. The 1D case is considered first for simplicity. Higher-dimensional cases can be built from the one-dimensional case.

3.1. The Poisson Equation in One Dimension

Considering only one space variable, x, Equation (1) can be written as:

\frac{d^{2} u (x)}{d x^{2}} = f (x), x \in Ω .

(4)

where

Ω

is the solution domain defining the possible x values. A set of N discrete x values

x_{1}, x_{2}, \dots, x_{N}

of x is selected where the equation is to be solved. Then, the domain is

Ω = [x_{1}, x_{N}]

. The boundary of the domain is assumed to be the outer edges:

\partial Ω = {x_{1}, x_{N}}

. Some texts use the notation

x_{0}

and

x_{N + 1}

to indicate the boundary points and exclude them from the definition of the solution domain. However, in this paper, the edge points of the solution domain are taken as boundary points unless otherwise stated.

The domain points are referred to as

x_{i}

,

i = 1, 2, \dots, N

. For simplicity, it is assumed that these points are separated by a constant number

Δ x

(i.e.,

x_{i + 1} - x_{i} = Δ x, \forall i

). The corresponding values of f and u are denoted

f_{i}

and

u_{i}

, respectively (i.e.,

f (x_{i}) = f_{i}

,

u (x_{i}) = u_{i}

). It is convenient to express these sets of discrete points as vectors

x = {[x_{1} x_{2} \dots x_{N}]}^{T}

,

f = {[f_{1} f_{2} \dots f_{N}]}^{T}

, and

u = {[u_{1} u_{2} \dots u_{N}]}^{T}

. Note that transpose notations were used to write the column vectors in a compact manner. With vector notation, the differential equation can be imagined as a linear algebra problem where the vector

u

is to be solved for a given

f

. The next step is to convert the second-derivative operator into a matrix operator so that a complete linear algebraic system can be formed.

To convert a continuous derivative operator to a discrete matrix operator, the starting point is the finite difference formulas for calculating the derivatives. The vectors

u^{'} = {[u_{1}^{'} u_{2}^{'} \dots u_{N}^{'}]}^{T}

and

u^{″} = {[u_{1}^{″} u_{2}^{″} \dots u_{N}^{″}]}^{T}

are defined to denote the first and second derivatives, respectively, corresponding to the points

x

(i.e.,

{u_{i}}^{'}

=

u^{'}

(x_{i})

, and

{u_{i}}^{″}

=

u^{″}

(x_{i})

). Using the second-order center-difference formulas, the first and second derivatives can be expressed as [16]:

\begin{matrix} {\frac{d u}{d x}|}_{x = x_{i}} & \approx & u_{i}^{'} = \frac{- u_{i - 1} + u_{i + 1}}{2 Δ x}, \end{matrix}

(5)

\begin{matrix} {\frac{d^{2} u}{d x^{2}}|}_{x = x_{i}} & \approx & u_{i}^{″} = \frac{u_{i - 1} - 2 u_{i} + u_{i + 1}}{{(Δ x)}^{2}} . \end{matrix}

(6)

Note that there are other higher-order finite difference formulas for the above derivatives. The widely used second-order center-difference formulas are selected for the current analysis. Now, from Equations (5) and (6), focusing only on three consecutive quantities

u_{i - 1}

,

u_{i}

,

u_{i + 1}

, the following matrix equation can be written:

\begin{matrix} u_{i}^{'} & = & \frac{1}{2 Δ x} [\begin{matrix} - 1 & 0 & 1 \end{matrix}] \cdot [\begin{matrix} u_{i - 1} \\ u_{i} \\ u_{i + 1} \end{matrix}], \end{matrix}

(7)

\begin{matrix} u_{i}^{″} & = & \frac{1}{{(Δ x)}^{2}} [\begin{matrix} 1 & - 2 & 1 \end{matrix}] \cdot [\begin{matrix} u_{i - 1} \\ u_{i} \\ u_{i + 1} \end{matrix}] . \end{matrix}

(8)

Equations (7) and (8) work fine for

i = 2

to

N - 1

. However, for the boundary points

i = 1

and

i = N

, the subscripts overflow the bounds, requiring values of

u_{0}

and

u_{N + 1}

, respectively. Thus, the equation is not valid for those two cases (note that the outer boundary points are

x_{1}

and

x_{N}

). We resort to using forward (for

i = 1

) and backward (for

i = N

) finite difference formulas for these two cases [16]:

\begin{matrix} u_{1}^{'} & = & \frac{- 3 u_{1} + 4 u_{2} - u_{3}}{2 Δ x} = \frac{1}{2 Δ x} [\begin{matrix} - 3 & 4 & - 1 \end{matrix}] \cdot {[\begin{matrix} u_{1} & u_{2} & u_{3} \end{matrix}]}^{T}, \end{matrix}

(9)

\begin{matrix} u_{1}^{″} & = & \frac{2 u_{1} - 5 u_{2} + 4 u_{3} - u_{4}}{{(Δ x)}^{2}} = \frac{1}{{(Δ x)}^{2}} [\begin{matrix} 2 & - 5 & 4 & - 1 \end{matrix}] \cdot {[\begin{matrix} u_{1} & u_{2} & u_{3} & u_{4} \end{matrix}]}^{T}, \end{matrix}

(10)

u_{N}^{'} = \frac{u_{N - 2} - 4 u_{N - 1} + 3 u_{N}}{2 Δ x} = \frac{1}{2 Δ x} [\begin{matrix} 1 & - 4 & 3 \end{matrix}] \cdot {[\begin{matrix} u_{N - 2} & u_{N - 1} & u_{N} \end{matrix}]}^{T},

(11)

\begin{matrix} u_{N}^{″} & = & \frac{- u_{N - 3} + 4 u_{N - 2} - 5 u_{N - 1} + 2 u_{N}}{{(Δ x)}^{2}} \\ = & \frac{1}{{(Δ x)}^{2}} [\begin{matrix} - 1 & 4 & - 5 & 2 \end{matrix}] \cdot {[\begin{matrix} u_{N - 3} & u_{N - 2} & u_{N - 1} & u_{N} \end{matrix}]}^{T} . \end{matrix}

(12)

Now, using Equations (7) to (12), the first and second derivatives at all points in the solution space can be calculated. These equations for all i values can be compiled to form the following matrix equations:

u^{'} = [\begin{matrix} u_{1}^{'} \\ u_{2}^{'} \\ u_{3}^{'} \\ ⋮ \\ u_{N - 1}^{'} \\ u_{N}^{'} \end{matrix}] = \underset{D_{x}}{\underset{︸}{\frac{1}{2 Δ x} [\begin{matrix} - 3 & 4 & - 1 & 0 & \dots & 0 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 1 & \dots & 0 & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & \dots & 0 & - 1 & 0 & 1 \\ 0 & 0 & 0 & 0 & \dots & 0 & 1 & - 4 & 3 \end{matrix}]}} \cdot [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \\ u_{N - 1} \\ u_{N} \end{matrix}],

(13)

u^{''} = [\begin{matrix} u_{1}^{″} \\ u_{2}^{″} \\ u_{3}^{″} \\ ⋮ \\ u_{N - 1}^{″} \\ u_{N}^{″} \end{matrix}] = \underset{D_{2 x}}{\underset{︸}{\frac{1}{{(Δ x)}^{2}} [\begin{matrix} 2 & - 5 & 4 & - 1 & \dots & 0 & 0 & 0 & 0 \\ 1 & - 2 & 1 & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & 1 & - 2 & 1 & \dots & 0 & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & \dots & 0 & 1 & - 2 & 1 \\ 0 & 0 & 0 & 0 & \dots & - 1 & 4 & - 5 & 2 \end{matrix}]}} \cdot [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \\ u_{N - 1} \\ u_{N} \end{matrix}] .

(14)

It should be noted that for nonuniform grids, the term

Δ x

would no longer be constant. As a result, in the place of the common

Δ x

factor in the matrices, every matrix element would have a denominator corresponding to its own

Δ x

value. Every other step discussed in this manuscript would remain unchanged. With the help of these matrices, the derivative operations can be expressed in linear algebra form as:

\begin{matrix} u^{'} & = & D_{x} u, \end{matrix}

(15)

\begin{matrix} u^{″} & = & D_{2 x} u . \end{matrix}

(16)

Here

D_{x}

and

D_{2 x}

are

N \times N

matrices, as shown in Equations (13) and (14). These are referred to as matrix operators for the first and second derivatives, respectively, as they operate on a vector to calculate its derivatives. Equation (4) can be converted to a discrete linear algebraic system using these matrix operators:

D_{2 x} u = f, on Ω .

(17)

Now the boundary conditions must be implemented in such a way that they can be easily integrated with Equation (17).

Let us consider a general implementation of boundary conditions in matrix form:

B u = g, on \partial Ω .

(18)

Here B is referred to as the boundary operator matrix, and

g

is the known/given boundary function vector. For a given problem, the boundary operator matrix is determined by the type of the boundary conditions. For Dirichlet boundary conditions, the boundary function values are directly assigned to the boundary points. In such a case, B would be an identity matrix, and

g

would list the values of

u

at the boundary. For Neumann boundary conditions, the value of the first derivative is assigned at the boundary points. For this case, B would be the first-derivative matrix operator

D_{x}

, and

g

would list the values of

u^{'}

at the boundary. Thus,

B = \{\begin{matrix} I_{N}, & for Dirichlet boundary conditions \\ D_{x}, & for Neumann boundary conditions \end{matrix}

(19)

Here

I_{N}

is the

N \times N

identity matrix, and

D_{x}

is expressed in Equation (13). Note that although Equation (19) relates the quantities in such a way that B operates on the entire

u

, the equation is valid only at the boundary

\partial Ω

. In fact,

g

elements at non-boundary points are unknown. In other words, the locations of the boundaries are not defined explicitly within B. However, to understand which entries of B come from

I_{N}

and which entries come from

D_{x}

, information about the geometry of the system is required. In the following paragraphs, we will discuss how to compile the linear algebraic system so that only the relevant rows of B are used. This is where the remaining geometry information will be encoded. Before that, it is important to discuss how a multi-boundary system would be treated. As previously mentioned,

\partial Ω

can be subdivided into smaller regions

\partial Ω_{1}, \partial Ω_{2} \dots \partial Ω_{M}

each with its own boundary conditions. Each

\partial Ω_{i}

corresponds to a set of specific rows in the boundary operator matrix. For such cases,

g

would be defined appropriately in a piecewise manner. For mixed boundary conditions, the rows of B would be a compilation of rows taken from either

I_{N}

or

D_{x}

depending on which boundary regime the corresponding equation falls under. It should also be noted that for Neumann boundary conditions, an additional constraint such as Equation (3) would be needed to obtain a unique solution.

After defining the matrix forms for the PDE and the boundary conditions, the process of integrating those into a single linear algebraic equation can be started. It is possible to simply solve Equation (17) under the constraint defined in Equation (19). However, an alternative and somewhat more elegant approach would be to form a single system matrix that incorporates both the PDE and the boundary conditions in a compact form. By observing the full form of the matrix and vectors in Equation (17), it can be noted that each row of the matrix defines an equation containing a few elements of

u

. Since some of the elements of

u

are on the boundary and follow Equation (19), replacing the corresponding rows of

D_{2 x}

with those of B can lead to the formation of a system matrix,

L

. This can be mathematically described as:

L = \{\begin{matrix} B, & on \partial Ω \\ D_{2 x}, & otherwise \end{matrix} .

(20)

Of course, similar changes to the right-hand vector would need to be made as well. The system right-hand vector,

b

, can be defined as:

b = \{\begin{matrix} g, & on \partial Ω \\ f, & otherwise \end{matrix} .

(21)

According to this definition, if the k-th element of

u

,

u_{k}

, corresponds to a Dirichlet boundary point, then the k-th row of

L

would be the k-th row of

I_{N}

. Further, the k-th element of

b

would be the k-th row of

g

. All rows corresponding to non-boundary points will be identical to the rows of

D_{2 x}

. Now, the system equation can finally be written as:

L u = b .

(22)

Dimension-wise,

L

is an

N \times N

matrix, and

b

is an

N \times 1

vector. This is a fully defined linear system from which

u

can be easily solved. Note that the derivative operator, the boundary conditions, and the problem geometry are all encoded within

L

and b. Let us consider the 1D case with

x_{1}

and

x_{N}

being the boundary points. For Dirichlet boundary conditions, this correspond to

{[u_{1} u_{N}]}^{T} = {[g_{1} g_{N}]}^{T}

, where

g_{1}

and

g_{N}

are given. Using Equations (20) and (21), the matrix system is constructed and solved. Note that

g_{i}, \forall i \neq 1, N

, are not defined and are not needed to form

b

.

One notable feature of the proposed system matrix construction procedure is that part of it is problem independent. Note that only the B matrix and the

g

vector depend on the boundary condition types at different boundaries. The

D_{2 x}

matrix is the same for all problems. The construction of

L

involves swapping out certain rows of

D_{2 x}

with B. This part depends on the problem geometry (as points on

\partial Ω

need to be identified). Thus, modifying the solver for different geometries and boundary conditions involves selecting the appropriate form of B and modifying the row-swapping operation of

L

only (the same operation must be done for the right-hand vector, also). This is a major benefit over conventional procedures for which the system matrix is constructed from scratch for every geometry/boundary condition.

It is obvious that the size of the matrix operators, having size

N \times N

, can be very large for moderately large values of N. As will be discussed later, the sizes of the matrices are even larger for cases with higher dimensions. Working with such large matrices in full form can be computationally demanding. Fortunately, these matrices are sparse in nature (i.e., the number of non-zero elements is small compared to the full size of the matrix), as can be seen from the forms of

D_{x}

and

D_{2 x}

. As such, when solving numerically, they are always defined as sparse matrices, which significantly reduces memory requirements and speeds up the operation. Almost all programming languages commonly used for numerical analysis have sparse-matrix libraries that can be used for this (e.g., Python has scipy.sparse package).

3.2. The Poisson Equation in Two Dimensions

The 2D Poisson equation in a Cartesian

x y

plane is given by:

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = f (x, y), x, y \in Ω .

(23)

Here the solution

u = u (x, y)

spans a 2D space. The matrix analysis for the 1D Poisson equation discussed in Section 3.1 can be extended for this 2D case. The matrix operators developed for the 1D case can be used as the building blocks. The key idea is to convert the 2D problem into an equivalent 1D problem. This is done by mapping the 2D solution grid into a 1D vector, solving the 1D vector, and then remapping it to the original 2D space. The main tasks are developing the mapping algorithm, figuring out how the boundary conditions are mapped from 2D to 1D space, and assembling the corresponding matrix operators.

The 2D solution space in the

x y

-Cartesian plane is defined as

Ω = [x_{1}, x_{N_{x}}] \times [y_{1}, y_{N_{y}}]

. This denotes a continuous rectangular region with x values spanning from

x_{1}

to

x_{N_{x}}

and y values spanning from

y_{1}

to

y_{N_{y}}

. Note that the analysis would hold for any shape of solution space. In discrete space, a grid of equally spaced

N_{y} \times N_{x}

points

(y_{k}, x_{j})

are considered, where

j = 1, 2, \dots, N_{x}

and

k = 1, 2, \dots, N_{y}

. The x and y spacing are denoted

Δ x

and

Δ y

, respectively. These points can be imagined as elements of a matrix with

N_{y}

rows and

N_{x}

columns, as shown in Figure 1. Any point can be addressed by the matrix coordinate

(k, j)

, where k refers to the row number and j refers to the column number. The solution u at point

(y_{k}, x_{j})

is denoted as

u (k, j) = u_{k j}

. The order of the indices j and k is selected in this way to be consistent with syntax of common numerical coding languages (e.g., Python and MATLAB).

Figure 1. Process of mapping a 2D grid to a 1D vector and vice-versa.

Now,

u_{k j}

, being 2D, is not an ideal input for a matrix operator. The elements of the 2D grid are rearranged by stacking the columns on top of each other to form one large

N_{x} N_{y} \times 1

vector, referred to as

u_{1 D}

. This is shown in Figure 1 [17]. Each element of this

u_{1 D}

corresponds to a specific grid point

(y_{k}, x_{j})

and the corresponding solution

u_{k j}

. The elements of

u_{1 D}

can be addressed by a single index

i = 1, 2, \dots, N_{x} N_{y}

. Note that there are many other valid approaches to mapping a 2D matrix into an 1D vector. This discussion is limited to the column-stacking method, which is relatively simple and intuitive. For this mapping, the 2D and 1D indices are related as:

\begin{matrix} i = N_{y} (j - 1) + k \end{matrix}

(24)

The opposite transformation is:

\begin{matrix} k & = & i \mod N_{y}, \end{matrix}

(25)

\begin{matrix} j & = & \frac{i - k}{N_{y}} + 1 . \end{matrix}

(26)

Here mod refers to the modulo operation. Using these relationships, it is possible to transfer back and forth between the 2D representation and the equivalent 1D representation. This process is referred to as mapping. The mapping process is shown in Figure 1.

After obtaining

u_{1 D}

from the mapping process, the corresponding derivative matrix operators need to formulated. The matrix operators for the first and second partial derivatives with respect to x are expressed as

D_{x}^{2 D}

and

D_{2 x}^{2 D}

, respectively. Similar quantities with respect to y are expressed as

D_{y}^{2 D}

and

D_{2 y}^{2 D}

, respectively. The superscripts denote that these are 2D matrix operators. For the 1D case discussed in Section 3.1, the consecutive elements of the vector

u

corresponded to consecutive points along a 1D spatial axis. For the 2D case, however, the relation between the elements of

u_{1 D}

and the corresponding spatial points is more complex. Every

N_{y}

block of elements in

u_{1 D}

(shown in different colors in Figure 1) represents a column in the 2D representation. They have consecutive y values for a given x value. Within this block, matrix operators of similar form as Equation (13) and (14) would work for the y derivative. As one moves from one of these blocks to the next (moving from one color to another in Figure 1), the y values restart from the beginning and the x values increment by one step. Our x derivative operator must operate on this jumbled space and accurately pick out the correct elements for a differentiation operation.

First, the y derivatives are considered, as they are simpler to understand than the x derivatives for this mapping process. The y partial derivatives in vector form are denoted

u_{y, 1 D}

and

u_{2 y, 1 D}

. The 2D matrix operator for the y partial derivative must be such that it sifts each

N_{y}

block of elements and multiplies it by the 1D-derivative matrix operator. As it goes through each of the columns, the results should be stacked on a column to maintain the same spatial mapping as

u_{1 D}

. This can be done by forming a block diagonal matrix whose diagonal blocks are copies of the 1D

D_{y}

matrix operator [8,17] (

D_{y}

has the same form as

D_{x}

from Equation (13), with

Δ x

replaced by

Δ y

). The process is shown graphically in Figure 2. Each rectangular block in

D_{y}^{2 D}

is an

N_{y} \times N_{y}

matrix. The zero blocks are zero matrices of the same size. Each block operates on an

N_{y} \times 1

sized block of

u_{1 D}

(represented as different-colored columns and labeled

C o l_{1}, C o l_{2}, \dots, C o l_{N_{x}}

). This is shown in the figure by connecting a

D_{y}^{2 D}

matrix block with the corresponding

u_{1 D}

column block by curved lines. For a given block-row (the rows drawn in the figure with black lines) of the matrix, all elements of

u_{1 D}

beside a specific block of

N_{y}

consecutive elements are ignored (i.e., multiplied by zeros). Thus, the derivative operation is performed for each of these column blocks as one moves from one block-row of the matrix to another. This matrix multiplication produces a column vector representing the y derivatives mapped in the same format as

u_{1 D}

. The second-derivative matrix,

D_{2 y}^{2 D}

, can be similarly constructed by replacing

D_{y}

with

D_{2 y}

in each block.

Figure 2. Matrix first-derivative operator (with respect to y) for the two-dimensional case.

D_{y}

has the same form as

D_{x}

which is defined in Equation (13).

Let us now consider the matrix operators for the x partial derivatives,

D_{x}^{2 D}

and

D_{2 x}^{2 D}

. When multiplied with

u_{1 D}

, these matrices produce the partial derivatives in vector form denoted

u_{x, 1 D}

and

u_{2 x, 1 D}

, respectively. First, the locations of the elements representing consecutive x values are identified and listed in the 1D-mapped vectors. The x values increase by one when moving from left to right in a given row (move by one column) in the 2D representation (Figure 1). In the 1D map, movement of one space along a row translates into movement by

N_{y}

elements in the corresponding 1D-mapped vector. A block matrix construction is again used to form

D_{x}^{2 D}

from

D_{x}

. This time, however, each row of

D_{x}

must operate on elements of

u_{1 D}

that are

N_{y}

spaces apart. This can be done by adding

N_{y}

zeros after each element in a row of

D_{x}

and then copying and arranging the rows appropriately [17]. The required configuration is shown in Figure 3. The elements of

D_{x}

are repeated in the diagonal of each block in

D_{x}^{2 D}

. Thus, each block is a diagonal matrix. Each block sifts out certain elements of

u_{1 D}

, as shown by the curved lines in the figure. Careful observation shows that these elements are indeed the ones required to calculate the derivatives at that point. A similar block matrix can be constructed for

D_{2 x}^{2 D}

by replacing

D_{x}

elements with

D_{2 x}

elements.

Figure 3. Matrix first-derivative operator (with respect to x) for the two-dimensional case.

D_{x} (j, k)

refers to the element from row j and column k of the

D_{x}

matrix defined in Equation (13).

Now that the matrix operators have been defined, their construction using simple mathematical operations is discussed. A powerful tool in constructing block matrices is the Kronecker product [4,18]. The Kronecker product between two matrices A (size

p \times q

) and B (size

m \times n

) is defined as:

A \otimes B = [\begin{matrix} A_{11} B & A_{12} B & \dots & A_{1 q} B \\ A_{21} B & A_{22} B & \dots & A_{2 q} B \\ ⋮ & ⋮ & ⋱ & ⋮ \\ A_{p 1} B & A_{p 2} B & \dots & A_{p q} B \end{matrix}]

(27)

where

A_{i j}

represents the

(i, j)

element of matrix A, and ⊗ represents Kronecker product operation. Note that each term in Equation (27) is a matrix itself. The

A \otimes B

block matrix is formed by tiling B matrices after multiplying them with an element of the A matrix. It is not difficult to see that

D_{y}^{2 D}

can be constructed from the Kronecker product of

D_{y}

and an identity matrix of size

N_{x}

(denoted as

I_{N_{x}}

):

D_{y}^{2 D} = I_{N_{x}} \otimes D_{y} .

(28)

Careful observation of Figure 3 and Equation (27) suggests that

D_{x}^{2 D}

can also be computed using the Kronecker product of

D_{x}

and an identity matrix of size

N_{y} \times N_{y}

(denoted

I_{N_{y}}

):

D_{x}^{2 D} = D_{x} \otimes I_{N_{y}} .

(29)

The second-derivative matrices are calculated identically:

\begin{matrix} D_{2 y}^{2 D} & = & I_{N_{x}} \otimes D_{2 y}, \end{matrix}

(30)

\begin{matrix} D_{2 x}^{2 D} & = & D_{2 x} \otimes I_{N_{y}} . \end{matrix}

(31)

The overall sizes of these matrices are

N_{x} N_{y} \times N_{x} N_{y}

. This can be very large even for moderate values of

N_{x}

and

N_{y}

. However, it is noted that the sparsity of the matrices from the 1D case is maintained in the 2D case as well. Even more so than with the 1D case, any numerical implementation of the 2D Poisson equation should employ sparse-matrix algorithms.

The step after constructing the differentiation operators is formulation of the boundary operator matrix. The 2D boundary operator matrix,

B^{2 D}

, would be identical to the one discussed in the previous section for the one-dimensional case. Equation (19) would be slightly modified to accommodate for the larger two-dimensional matrices:

B^{2 D} = \{\begin{matrix} I_{N_{x} N_{y}}, & for Dirichlet boundary conditions \\ D_{x}^{2 D} or D_{y}^{2 D} & for Neumann boundary conditions \end{matrix}

(32)

where

I_{N_{x} N_{y}}

is an identity matrix of size

N_{x} N_{y} \times N_{x} N_{y}

. As there are two independent variables, there are two possible Neumann boundary conditions that can be applied. Which derivative to use for the Neumann boundary is determined by the problem statement. Having defined the boundary operator, construction of the system matrix operator is straightforward:

L = \{\begin{matrix} B^{2 D}, & on \partial Ω \\ D_{2 x}^{2 D}, & otherwise \end{matrix} .

(33)

The right-hand vector also follows a similar definition as the 1D case:

b = \{\begin{matrix} g_{1 D}, & on \partial Ω \\ f_{1 D}, & otherwise \end{matrix},

(34)

where

f_{1 D}

and

g_{1 D}

are 1D-mapped versions of the corresponding two-dimensional quantities,

f (y_{k}, x_{j})

and

g (y_{k}, x_{j})

, respectively. Obviously, the algorithm used for the 1D mapping here should be identical to the one used for calculating

u_{1 D}

. Finally, the system equation is constructed as:

L u_{1 D} = b .

(35)

Just like the one-dimensional case, this is a simple linear algebraic equation that can be solved easily. It should be noted that the solution is in a 1D-mapped vector form. Thus, an additional step of converting it back to the 2D form is necessary. This can be done using Equations (25) and (26).

It should be noted that only the

B^{2 D}

matrix, the

g_{1 D}

vector, and the construction of

L

and

b

are geometry and boundary condition dependent. The rest of the matrices and vectors are problem independent. Thus, configuring the solver for different geometries and boundary conditions is straightforward.

3.3. Extending to Three Dimensions

In the previous subsection, it was discussed how the two-dimensional system can be constructed using the one-dimensional operators as building blocks. A similar approach can be used to construct a three-dimensional system from the two-dimensional matrix operators. The process will involve Kronecker products of the identity matrices and the two-dimensional matrix operators. As this is intended to be an introductory paper, the details of the three-dimensional case are omitted for the sake of brevity.

4. Python Implementation

Python code for solving the Poisson equation in two-dimensions has been developed. It is publicly available on GitHub [13]. The code uses the previously discussed approach of constructing the matrix operators. The code is general and can solve any system when the boundary conditions are appropriately defined. In addition to the basic numerical recipe, a few additional things had to be taken into account for computational efficiency, with the most-necessary thing being the use of sparse matrices. The sparse matrix sub-package of SciPy (Scientific Python) [19] is used for constructing matrix operators. In this way, only the non-zero elements of the matrices are indexed and stored. This makes it possible to store large sparse matrices without wasting computer memory on storing the zeros. Further, mathematical operations can be performed more efficiently. For this work, the sparse linear algebra solver of SciPy, spsolve, was used.

In the following paragraphs, how the Python code should be used for specific problems is discussed. The focus is mainly on how geometric features and boundary conditions are defined inside the code. A few key features of the code are also highlighted.

The code requires the user to define the source function and the geometry of the problem. From those, the code automatically creates the boundary operator matrices, the system matrix, and the right-hand vectors. How the geometry is defined inside the code is discussed first. Two arrays that denote the x and y independent variables are declared. The range and the resolution of the variables are also defined here. An example code snippet is given below:

import numpy as np

# Define independent variables.

Nx = 300                    # No. of grid points along x direction

Ny = 200                    # No. of grid points along y direction

x = np.linspace(-6,6,Nx)     # x variable in 1D

y = np.linspace(-3,3,Ny)     # y variable in 1D

The size of the arrays determines the size of the matrix operators. Note that

Δ x

and

Δ y

are implicitly defined here. The 1D arrays are converted into 2D meshgrid data format, which is then unraveled into an equivalent 1D structure for ease of indexing later. The code for these operations is:

X,Y = np.meshgrid(x,y)       # 2D meshgrid

# 1D indexing

Xu = X.ravel()               # Unravel 2D meshgrid to 1D array

Yu = Y.ravel()

After the space variables, the source function needs to be defined. For a source term having a closed-form expression, a simple for loop can be used to do this operation, as shown below:

# Source term for the diode example

f = np.zeros(Nx∗Ny)

for m in range(Nx∗Ny):

    if np.abs(Xu[m]) < 3.8 and np.abs(Yu[m])<1:

        f[m] = -np.tanh(Xu[m])/np.cosh(Xu[m])

This test source function is used in the semiconductor p–n junction example problem discussed later in Section 5.2. In this example, the source function is set equal an analytical functional at specific regions of space and set to zero elsewhere. Instead of an analytical function, tabulated values can also be used by utilizing an interpolation function. Using such piecewise definitions, source functions for most practical problems can be implemented.

The next step is to define outer boundary regions. The definition consists of the boundary values at the four outer walls (for 2D problems) as well as the types of boundary conditions. An example code snippet is given below:

# Dirichlet/Neumann boundary conditions at outer walls

uL, uR, uT, uB = 0, 0, 0, 0

ub_o = [uL, uR, uT, uB]

# Type of outer boundary conditions. 0 = Dirichlet,

# 1 = Neumann x derivative, 2 = Neumann y derivative

B_type_o = [1,1,2,2]

# Order: element 1 = left, element 2 = right, element 3 = top,

# element 4 = bottom

Now we discuss the non-trivial part of the geometry: the internal regions where boundary conditions are enforced. For now, the discussion is limited to rectangular regions. Two Python lists containing the high and low limits of

x, y

values are used to define these inner boundary regions. Then, the boundary type and boundary values are set as previously seen. The following code snippet defines two rectangular regions, each with Dirichlet boundary conditions (one with

- 2

and one with 2). This code snippet is from the semiconductor p–n junction diode example discussed later in Section 5.2.

# lower and upper limits of x defining the inner boundary region

# Format:[ [xlow1,xhigh1], [xlow2,xhigh2],.... ]

xb_i = [[-4,-3.8],[3.8,4]]

yb_i = [[-1,1],[-1,1]]

# Type of inner boundary conditions. 0 = Dirichlet,

# 1 = Neumann x derivative, 2 = Neumann y derivative

B_type_i = [0,0]

ub_i = [-2,2]            # boundary values at inner region

The points that are within the rectangular region defined by the limits of x and y are automatically calculated by the code. Using this method, any number of arbitrary inner boundary regions can be defined as long as they are rectangular in shape.

For non-rectangular geometric features, the regions can be difficult to define by code algebraically. To tackle this, an image-based geometry-definition feature is implemented. The idea is that the geometric features can be extracted from a bitmap image file by code. Different colors in an image represent different regions (some of these will be boundary regions). A user can draw a bitmap image using any graphical editor. The graphic editor GraphicsGale [20], which is a freeware, is recommended. The pixels of the image will be treated as a grid point of the solution domain. With this approach, a user can define any arbitrary geometry by drawing it in an image editor and then importing it in the solver. The code snippet of the image import feature is given below:

# Read image file

im = Image.open(im_file_name)

boundary_region_color_threshold = 10    # color value that identifies boundary

strct_map = np.flipud(np.array(im))     # flipud is done to match orientation

strct_map_u = strct_map.ravel()

regions = np.unique(strctmapu)           # Find distinct color regions

n_regions = np.size(regions)             # Number of distinct color~regions

b_regions_indices_u = []

o_regions_indices_u = []

for m in range(nregions):

    ind_u_temp = np.squeeze(np.where(strctmapu== regions[m])

    if regions[m] <= boundary_region_color_threshold: # boundary region

        b_regions_indices_u.append(ind_u_temp)

    else:

        o_regions_indices_u.append(ind_u_temp)

The code reads an image file and identifies regions with distinct colors. An arbitrary color threshold value is defined. Any color value less than the threshold is considered a boundary region. The other colors are treated as normal regions where the solution is to be calculated. The indices of all the boundary (and non-boundary) regions are stored in a list. These indices correspond to the x and y variable indices.

Section 5.3 discusses an example problem using the image import feature. More details about the code can be found within the comments of the code file on the GitHub repository [13].

5. Numerical Results

Using the developed code, a few simple example problems are solved. The examples are limited to electrostatics. The electrostatic potential distribution in a two-dimensional space is determined by the Poisson equation in the following form:

\nabla^{2} u (x, y) = f (x, y) = \frac{ρ}{ϵ}, on Ω .

(36)

Here

ρ

is the charge density,

ϵ

is the permittivity, and u represents the electric potential or voltage. The electric field,

E

, is given by the negative gradient of u:

E = - \nabla u (x, y) .

(37)

Many interesting electrostatics problem involve solving the electric potential and electric field distribution using the Poisson equation. A few example cases are discussed here. As the focus is on the numerical approach rather than the physical system, simple source functions and boundary conditions are used that are not necessarily practical. Units of the variables are also ignored.

It should be noted that all three examples are solved using the same numerical solver with very little modification of the Python code. Only the definition of the geometry has to be redone. The rest of the code works without any further modifications. This is the distinct advantage of the proposed approach.

5.1. Example 1: Parallel Plate Capacitor

A parallel plate capacitor consists of two parallel conducting plates separated by an insulating region. When voltage is applied between the plates, the capacitor stores energy in the gap between the plates in the form of electric fields. Usually, the charge density is taken to be zero inside a parallel plate capacitor (i.e.,

ρ = 0

), which implies

f (x, y) = 0

. Thus, with a zero source term, the equation reduces to a Laplace equation:

\nabla^{2} u (x, y) = 0, on Ω .

(38)

A parallel plate capacitor bounded within a larger simulation space with Neumann boundary conditions is simulated. This example is selected to show the versatility of the numerical solver in implementation of both inner and outer boundary conditions. The solution domain is selected to be be a rectangle bounded by

- 6 \leq x \leq 6

,

- 3 \leq y \leq 3

(i.e.,

Ω ≜ [- 6, 6] \times [- 3, 3]

). The outer edge boundaries are defined as:

\begin{matrix} \partial Ω_{1} & ≜ & \{(x, y) | x = - 6\} (left outer edge), \end{matrix}

(39)

\begin{matrix} \partial Ω_{2} & ≜ & \{(x, y) | x = 6\} (right outer edge), \end{matrix}

(40)

\begin{matrix} \partial Ω_{3} & ≜ & \{(x, y) | y = - 3\} (bottom outer edge), \end{matrix}

(41)

\begin{matrix} \partial Ω_{4} & ≜ & \{(x, y) | y = 3\} (top outer edge) . \end{matrix}

(42)

The boundary conditions at the outer edges are taken to be:

\begin{matrix} \frac{\partial u}{\partial x} & = & 0, at \partial Ω_{1} and \partial Ω_{2}, \end{matrix}

(43)

\begin{matrix} \frac{\partial u}{\partial y} & = & 0, at \partial Ω_{3} and \partial Ω_{4} . \end{matrix}

(44)

These are Neumann boundaries that force the electric field to be zero at the outer edges. Again, these boundaries are not usually important, as only the solution within the region bounded by the parallel plates is of interest. These extra boundary conditions are added to create a more complex demonstration problem.

Next, two inner boundary regions are defined by two rectangles representing the parallel plates:

\begin{matrix} \partial Ω_{5} & ≜ & \{(x, y) | - 2 \leq x \leq 2 and 0.5 \leq y \leq 0.6\} (top plate) \end{matrix}

(45)

\begin{matrix} \partial Ω_{6} & ≜ & \{(x, y) | - 2 \leq x \leq 2 and - 0.6 \leq y \leq - 0.5\} (bottom plate) \end{matrix}

(46)

Dirichlet boundary conditions are set at these boundaries:

\begin{matrix} u & = & 2, at \partial Ω_{5}, \end{matrix}

(47)

\begin{matrix} u & = & - 2, at \partial Ω_{6} . \end{matrix}

(48)

These values represent the voltage applied to these plates. The solution domain with all the boundary regions is shown in Figure 4.

Figure 4. Solution domain for the two-dimensional parallel plate capacitor problem. The boundaries defined by the red points are Neumann boundaries. The gray regions indicate the parallel plates, where Dirichlet boundary conditions are set. The light green regions have no boundary conditions enforced on them.

Now that the problem, the solution domain, and the boundary conditions are defined, the process of constructing the boundary operator matrix, the system matrix, and the right-hand vector can be started. The boundary operator follows Equation (19):

B = \{\begin{matrix} D_{x}^{2 D}, & at \partial Ω_{1}, \partial Ω_{2} \\ D_{y}^{2 D}, & at \partial Ω_{3}, \partial Ω_{4} \\ I_{N}, & at \partial Ω_{5}, \partial Ω_{6} \end{matrix}

(49)

The system matrix has the same form as defined in Equation (35). The right hand vector follows:

b = \{\begin{matrix} 2, & on \partial Ω_{5} \\ - 2, & on \partial Ω_{6} \\ 0, & otherwise \end{matrix}

(50)

With all the mathematical quantities defined, the linear system of the equation is solved using a sparse-matrix solver. After mapping the solution back into 2D space,

u (x, y)

and

- \nabla u

are plotted in Figure 5. From the

u (x, y)

contour plot, it can be seen that the top and bottom electrodes have values of

+ 2

and

- 2

, respectively, and the potential gets distributed to other points in space. The

- \nabla u

distribution shows how the electric field lines point from the positive voltage side to the negative voltage side. The lines inside the capacitor are oriented vertically, as one would expect from an ideal parallel plate capacitor. It is also noted that the fringe fields at the edge of the plates are curved. The field magnitude is large at the edges due to the sharp potential transition. The electric field lines outside the capacitor region are mostly determined by the Neumann boundary conditions enforced at the outer edges of the solution domain. As previously stated, the solution at these locations is not usually relevant when studying a parallel plate capacitor unless one wishes to study the fringe fields. It is possible to limit the solution domain to the inside of the capacitor only, which would require setting only outer-edge boundaries. The slightly more-complicated boundaries were enforced for demonstration purposes. The proposed scheme can model physical systems that require boundary conditions to be set both at the outer edges as well as at inner regions.

Figure 5. Solution of the Poisson equation for a parallel plate capacitor: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

5.2. Example 2: Semiconductor p–n Junction

For the second example, a highly simplified semiconductor p–n junction diode structure is considered. The diode consists of two different types of semiconductors that are sandwiched together. Voltage is applied at the two ends of the semiconductors. Such a device can act as a rectifier (allowing electric current to flow in only one direction). More details about device operation are avoided as they are not in the scope of this text. The potential profile along the junction of the device follows the Poisson equation as expressed in Equation (36). The charge profile along the junction can be modeled in various ways. A common model of charge density is [21]:

ρ = 2 ρ_{0} sech (\frac{x}{a}) tanh (\frac{x}{a}) .

(51)

Here

ρ_{0}

and a are parameters related to the geometrical and material properties of the device. It is assumed that the two semiconductor materials are placed along the x axis to form a junction along the line

x = 0

. Being consistent with this form of charge density profile, the source function,

f (x, y)

, is defined as:

f (x, y) = \{\begin{matrix} sech (x) tanh (x), & for | x | < 3.8, | y | < 1 \\ 0, & otherwise \end{matrix} .

(52)

Here convenient values of parameters that produce a simple expression are used. Note that for this demonstration problem, only the same functional variation is maintained. For simplicity, realistic values of the functions are not used. In addition to the source function, two Dirichlet boundary regions,

\partial Ω_{5}

and

\partial Ω_{6}

, are defined. These represent metal contacts on either side of the semiconductors where voltage can be applied:

\begin{matrix} \partial Ω_{5} & ≜ & \{(x, y) | 3.8 \leq x \leq 4 and - 1 \leq y \leq 1\} (right contact) \end{matrix}

(53)

\begin{matrix} \partial Ω_{6} & ≜ & \{(x, y) | - 4 \leq x \leq 3.8 and - 1 \leq y \leq 1\} (left contact) \end{matrix}

(54)

The four outer boundary edges (i.e.,

\partial Ω_{1}

to

\partial Ω_{4}

) have the same definitions as in the previous example (Equations (39)–(42)). The simulation domain and the source function are shown in Figure 6. As with the previous example, the region outside the device (where

f (x, y) = 0

) is not meaningful. It is simulated only as a test for the numerical solver.

Figure 6. (a) Solution domain for the two-dimensional p–n junction problem. The red points are Neumann boundaries. The gray regions indicate the inner Dirichlet boundaries. (b) The source function

f (x, y)

.

The system matrix and the right-hand vector can be constructed using the same process as discussed in Section 5.1. The exact same boundary values are used, and a voltage of

u = - 2

and

u = 2

is applied to the left and right contact regions, respectively. Please note that the voltages applied to the metal contacts affect the charge density profile (and hence the source function

f (x, y)

). It is assumed that the

f (x, y)

function already contains the charge profile that would exist when the aforementioned voltage is applied. More-complex modeling of the charge distribution would be needed to study the p–n junction under a variety of voltage conditions. The results for our simplified case are shown in Figure 7. It can be noted that the electric field (i.e.,

- \nabla u

) is maximum near the junction along

x = 0

, which is to be expected for this configuration [21].

Figure 7. Solution of the Poisson equation for a simplified p–n junction: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

5.3. Example 3: Complex Arbitrary Geometries with Image Import

There are cases where the geometry of a problem can be more complex than that of the two examples discussed above. For those cases, it can be difficult to define the geometry through text data or code. For example, any region with non-rectangular shapes can be cumbersome to define by code alone. For such cases, the image import feature discussed in Section 4 is used. This feature is highlighted in this example. First, an arbitrary geometry is drawn in a 300 pixel × 200 pixel grid using the free image editor GraphicsGale [20]. Note that this can be done using any image editor. Each separate region of the geometry is drawn using a separate color. For this example, the geometry shown in Figure 8a is selected. The image is imported into the solver using the developed code, and the boundary regions are highlighted in Figure 8b. In the following, all colors are referred to with respect to Figure 8b. The orange and blue pixels represent zero Neumann boundary conditions (x- and y-derivative Neumann boundary condition, respectively) at the outer walls. The yellow rectangular regions are a Dirichlet boundary of value 2. The green circular regions are a Dirichlet boundary of value

- 2

, and the red circular region is a Dirichlet boundary of value 1. Note that some colored regions of Figure 8a are not reflected in Figure 8b because they are not defined as boundary regions. Those regions are considered simple solution regions instead. Note that there can be more than one of these regions (as in this example). Having multiple such regions can help us define complex source functions that have different values in different regions. For simplicity, a source function that has zero values everywhere is assumed.

Figure 8. (a) Solution domain for the two-dimensional arbitrary problem with external-image-defined geometry. (b) The boundary regions identified by the code. Different colors represent different boundary conditions to be enforced.

The Poisson equation is solved for the aforementioned geometry, and the results are shown in Figure 9. Like the previous examples, these results are interpreted in context of electric potential and electric field. Both the potential and the electric field have the expected distribution.

Figure 9. Solution of the Poisson equation for an arbitrary image-defined geometry: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

6. Conclusions

A method for solving the Poisson equation in 1D and 2D using a finite difference approach is presented. The main focus was on the systematic construction of the solver so that it can be applied to a variety of problems involving the Poisson equation. First, the matrix operators for the derivative operation and the matrix operators for boundary condition enforcement were constructed. The two matrices were combined to form a system matrix operator. A similar operation was performed for the right-hand side vector. Thus, a systematic approach of converting the Poisson equation into a linear algebraic/matrix equation was developed. The boundary conditions with geometry information were coded into a matrix separately from the derivative operator, making the approach easy to adjust for different geometries and boundary conditions. Three example problems were solved using the developed approach. The key benefit of the proposed method is that very little modifications of the numerical solver are needed to handle three problems that are significantly different from each other. The procedure for setting up the solver for these problems is discussed thoroughly. Further, the Python code for the algorithm is available to the readers though a GitHub repository. It should be noted that the proposed method has some of the drawbacks inherent with all finite difference solvers, for example, meshing constraints when handling complex geometries and boundaries.

Funding

This research received no external funding.

Data Availability Statement

The data that support the findings of this study are openly available in the GitHub repository at https://github.com/zaman13/Poisson-solver-2D, Accessed date: 1 July 2022.

Acknowledgments

The author thanks Paul C. Hansen.

Conflicts of Interest

The author declares no conflict of interest.

List of Symbols

$\nabla^{2}$	Laplacian operator in Cartesian coordinates
u	Dependent variable of the Poisson equation (to be solved)
f	Source function of the Poisson equation
g	Boundary function of the Poisson equation
$Ω$	The solution domain
$\partial Ω$	Boundary of the solution domain
$u$	Vectorized form of u on the solution grid
$f$	Vectorized form of f on the solution grid
$g$	Vectorized form of g on the solution boundary
$b$	The right hand vector
$u_{1 D}$	1D-mapped vectorized form of u for the 2D problem
$f_{1 D}$	1D-mapped vectorized form of f for the 2D problem
$g_{1 D}$	1D-mapped vectorized form of g for the 2D problem
$D_{x}$	First-derivative matrix operator with respect to x in 1D
$D_{y}$	First-derivative matrix operator with respect to y in 1D
$D_{2 x}$	Second-derivative matrix operator with respect to x in 1D
$D_{2 y}$	Second-derivative matrix operator with respect to y in 1D
$D_{x}^{2 D}$	First-derivative matrix operator with respect to x in 2D
$D_{y}^{2 D}$	First-derivative matrix operator with respect to y in 2D
$D_{2 x}^{2 D}$	Second-derivative matrix operator with respect to x in 2D
$D_{2 y}^{2 D}$	Second-derivative matrix operator with respect to y in 2D
B	Boundary operator matrix
$I_{N}$	Identity matrix of dimension $N \times N$
$L$	System matrix
⊗	Kronecker product
$E$	The electric field
$ρ$	Charge density
$ϵ$	Permittivity of the medium

References

Mori, Y.; Takabatake, K.; Tsugeno, Y.; Sakai, M. On artificial density treatment for the pressure Poisson equation in the DEM-CFD simulations. Powder Technol. 2020, 372, 48–58. [Google Scholar] [CrossRef]
Radhakrishnan, A.; Xu, M.; Shahane, S.; Vanka, S.P. A Non-Nested Multilevel Method for Meshless Solution of the Poisson Equation in Heat Transfer and Fluid Flow. arXiv 2021, arXiv:2104.13758. [Google Scholar]
Baranyi, L. Computation of unsteady momentum and heat transfer from a fixed circular cylinder in laminar flow. J. Comput. Appl. Mech. 2003, 4, 13–25. [Google Scholar]
Brewer, J. Kronecker products and matrix calculus in system theory. IEEE Trans. Circuits Syst. 1978, 25, 772–781. [Google Scholar] [CrossRef]
Schroeder, J.; Muller, R. IGFET Analysis through numerical solution of Poisson’s equation. IEEE Trans. Electron Devices 1968, 15, 954–961. [Google Scholar] [CrossRef]
Matsuura, T.; Saitoh, S.; Trong, D. Numerical solutions of the poisson equation. Appl. Anal. 2004, 83, 1037–1051. [Google Scholar] [CrossRef]
Causon, D.; Mingham, C. Introductory Finite Difference Methods for PDEs; Bookboon: London, UK, 2010. [Google Scholar]
LeVeque, R.J. Finite Difference Methods for Ordinary and Partial Differential Equations; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2007. [Google Scholar] [CrossRef]
Press, W.H.; Flannery, B.P. Numerical Recipes: The Art Of Scientific Computing, 3rd ed.; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
Burden, A.M.; Burden, R.L.; Faires, J.D. Numerical Analysis, 10th ed.; Cengage: Boston, MA, USA, 2016. [Google Scholar] [CrossRef]
Nagel, J.R. Numerical Solutions to Poisson Equations Using the Finite-Difference Method [Education Column]. IEEE Antennas Propag. Mag. 2014, 56, 209–224. [Google Scholar] [CrossRef]
Abdallah, S. Numerical solutions for the pressure Poisson equation with Neumann boundary conditions using a non-staggered grid, I. J. Comput. Phys. 1987, 70, 182–192. [Google Scholar] [CrossRef]
Zaman, M. Poisson-sover-2D, 2022. Available online: https://github.com/zaman13/Poisson-solver-2D (accessed on 1 July 2022).
Sapa, L.; Bożek, B.; Danielewski, M. Difference Methods To One And Multidimensional Interdiffusion Models With Vegard Rule. Math. Model. Anal. 2019, 24, 276–296. [Google Scholar] [CrossRef]
Sapa, L.; Bożek, B.; Tkacz–Śmiech, K.; Zajusz, M.; Danielewski, M. Interdiffusion in many dimensions: Mathematical models, numerical simulations and experiment. Math. Mech. Solids 2020, 25, 2178–2198. [Google Scholar] [CrossRef]
Dunn, S.M.; Constantinides, A.; Moghe, P.V. Finite Difference Methods, Interpolation and Integration. In Numerical Methods in Biomedical Engineering; Elsevier: Amsterdam, The Netherlands, 2006; pp. 163–208. [Google Scholar] [CrossRef]
Basilisk Software. Available online: http://basilisk.fr/sandbox/easystab/diffmat_2D.m (accessed on 12 April 2022).
Regalia, P.A.; Sanjit, M.K. Kronecker Products, Unitary Matrices and Signal Processing Applications. SIAM Rev. 1989, 31, 586–613. [Google Scholar] [CrossRef]
Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [Green Version]
GraphicsGale Software. Version 2.08.20. 2020. Available online: https://graphicsgale.com/us/ (accessed on 12 April 2022).
Hayt, W. Engineering Electromagnetics; McGraw-Hill: New York, NY, USA, 2012. [Google Scholar]

Figure 1. Process of mapping a 2D grid to a 1D vector and vice-versa.

Figure 2. Matrix first-derivative operator (with respect to y) for the two-dimensional case.

D_{y}

has the same form as

D_{x}

which is defined in Equation (13).

Figure 3. Matrix first-derivative operator (with respect to x) for the two-dimensional case.

D_{x} (j, k)

refers to the element from row j and column k of the

D_{x}

matrix defined in Equation (13).

Figure 4. Solution domain for the two-dimensional parallel plate capacitor problem. The boundaries defined by the red points are Neumann boundaries. The gray regions indicate the parallel plates, where Dirichlet boundary conditions are set. The light green regions have no boundary conditions enforced on them.

Figure 5. Solution of the Poisson equation for a parallel plate capacitor: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

Figure 6. (a) Solution domain for the two-dimensional p–n junction problem. The red points are Neumann boundaries. The gray regions indicate the inner Dirichlet boundaries. (b) The source function

f (x, y)

.

Figure 7. Solution of the Poisson equation for a simplified p–n junction: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

Figure 8. (a) Solution domain for the two-dimensional arbitrary problem with external-image-defined geometry. (b) The boundary regions identified by the code. Different colors represent different boundary conditions to be enforced.

Figure 9. Solution of the Poisson equation for an arbitrary image-defined geometry: (a) potential distribution, u, and (b) electric field distribution. The colors in (b) represent the magnitude of the electric field,

| - \nabla u |

, and the arrow lines are electric field lines showing the direction of

- \nabla u

.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Zaman, M.A. Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators. Electronics 2022, 11, 2365. https://doi.org/10.3390/electronics11152365

AMA Style

Zaman MA. Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators. Electronics. 2022; 11(15):2365. https://doi.org/10.3390/electronics11152365

Chicago/Turabian Style

Zaman, Mohammad Asif. 2022. "Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators" Electronics 11, no. 15: 2365. https://doi.org/10.3390/electronics11152365

APA Style

Zaman, M. A. (2022). Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators. Electronics, 11(15), 2365. https://doi.org/10.3390/electronics11152365

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Citations

Crossref

8

Scopus

10

Web of Science

8

Google Scholar

[click to view]

Article Access Statistics

For more information on the journal statistics, click here.

Multiple requests from the same IP address are counted as one view.

Article Menu

Numerical Solution of the Poisson Equation Using Finite Difference Matrix Operators

Abstract

1. Introduction

2. The Poisson Equation and Boundary Conditions

3. Finite Difference Discretization and Matrix Operators

3.1. The Poisson Equation in One Dimension

3.2. The Poisson Equation in Two Dimensions

3.3. Extending to Three Dimensions

4. Python Implementation

5. Numerical Results

5.1. Example 1: Parallel Plate Capacitor

5.2. Example 2: Semiconductor p–n Junction

5.3. Example 3: Complex Arbitrary Geometries with Image Import

6. Conclusions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

List of Symbols

References

Share and Cite

Article Metrics

Citations

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI