A DIMENSION-INDEPENDENT FINITE DIFFERENCE METHOD FOR THE POISSON EQUATION USING AN INDEX MAPPING FUNCTION

ABSTRACT

The numerical solution to the Poisson Problem is widely known and studied in various fields of science for its vast applications. However, most applications consider the case in two dimensions, with fewer studies addressing the problem in dimension three or higher. Most literature texts present a twodimensional case implementation using an approach that makes it difficult to extend to higher dimensions. Our work aims to propose a generalization of the numerical solution to the Poisson problem that can be implemented for any dimension. The strategy used considers an index function that enumerates the mesh points so that, using this index, the implementation is easily extended to any dimension. In addition to the numerical solution extension, we developed the mathematical foundation for the consistency and stability of the solution in arbitrary dimension. The preliminary results consider the implementation in Python and experiments that demonstrate the feasibility of the proposed methodology.

Keywords:
Poisson’s equation; Dirichlet boundary condition; numerical solution

1 INTRODUCTION

It is a known fact that many models exist in the literature for the Poisson Problem in one and two dimensions, both demonstrated with the necessary mathematical rigor and implemented in some programming language. However, there is a smaller amount of work on mathematical modeling for the three-dimensional case of the Poisson Problem. There is also a difficulty when trying to expand the modeling to a new dimension, requiring a new implementation, with non-trivial changes, for each desired dimension. In the present work, we note that the creation of an index function to enumerate the mesh points corresponding to unknown values, not only facilitates the implementation for the two-dimensional and three-dimensional case, but also creates a generalization that facilitates the implementation in any dimension. This approach, using the index function, was not found in any other work during our literature review. This function plays a very important role in the implementation because, by enumerating sampled elements of the domain, it simplifies the construction of the coefficient matrix of the sparse linear system that arises in the modeling of the solution. The index function allows this to be done systematically and efficiently and, moreover, it can be rewritten very easily so that the path to access the elements is not unique and depends on the enumeration one chooses, without this choice interfering with the solution.

In this work, therefore, we generalize the numerical solution to the Poisson problem using the finite difference method in the n dimension, with n ≥ 1. Consider the second-order elliptical problem, in the dimension n, given by

$-\Delta u\left(x\right)=f\left(x\right),x\in \Omega ,$ (1.1)

with the Dirichlet boundary condition

$u\left(x\right)=g\left(x\right),x\in \partial \Omega ,$

where Ω = (r ₁ , s ₁)×... × (r _n , s _n ), r _j < s _j , for j = 1, . . . , n, ∂ Ω is the boundary of Ω and ∆ denotes the Laplacian operator, given by

$\Delta u\left(x\right)=\sum _{j=1}^n\frac{{\partial }^{2}u}{\partial x_j^2}\left(x\right).$

There are different approaches to solving these problems, depending on the geometry of the domain. Some recent studies, including the one, two and three dimensional Poisson equation, present numerical solutions using increasingly accurate and efficient methods ^{19), (18}. There is a large amount of literature about the numerical solution of a Poisson equation, and a good review of results on this important subject can be found in ^{4), (12), (8), (7), (17), (21), (22), (27}.

In this paper, motivated by the two-dimensional case ^{6), (15} and ²⁴, the finite difference method is preferred due to the ease of implementation and computational efficiency for linear problems in rectangular regions of ℝ². In practice, problems may be nonlinear, have nonrectangular domains, and be defined in ℝⁿ with n > 2. Therefore, extending the method to handle these cases is of great importance for solving real-world problems.

In recent decades, numerous applications and extensions of the finite difference method for nonlinear Poisson-type differential equations have emerged. Owing to the nonlinear behavior of the partial differentiable equation under consideration, the theoretical analysis has been proved to be considerably difficult, especially for problems with irregular geometries and non-uniform boundary conditions. To study the nonlinearity in complex solution domain, it is a long history in resorting to numerical solutions. So far, different numerical techniques, including finite difference method ^{1), (5), (16}, exponential finite difference method ²⁰, quasilinear boundary element method ¹⁴, hybrid fundamental solution-based finite element method ²⁶, the method of fundamental solution ^{2), (3}, among others, have been developed to solve nonlinear Poissontype problems. More discussions on this topic, using other numerical techniques, can be found in the literature ^{28), (9}.

There are several studies in the literature in which the finite difference method is applied to solve elliptic equations in the irregular domain ^{11), (13}. In principle, its application is restricted to a domain with a boundary that has a relatively simple geometry ¹⁰, but for a domain with a more general boundary, certain special steps must be taken, with an exploratory approach to overcome this fact, as reported in the literature ²³.

Although the finite difference formulation in ℝⁿ follows classical generalizations, our focus is on an efficient computational abstraction via an index function that is rarely addressed explicitly in the literature. In most cases, the solution goes through the construction of a system of linear equations AU = B, where A is a sparse matrix in which the location of non-zero elements is not trivial, especially as the dimension n of the domain increases. Once the problem is formulated as a linear system, a wide range of existing numerical solvers and optimization techniques can be employed to improve CPU efficiency, particularly given that the resulting system matrix is highly sparse. This opens up opportunities for further performance improvements through the use of specialized sparse matrix libraries and parallel computing techniques.

The main contributions of this paper are the generalization of the numerical solution of the Poisson equation, with Dirichlet boundary conditions, in the dimension n, for all n ≥ 1, and the introduction of an index function I that facilitates the numerical modeling of the solution, and its computational implementation, using the sparse system of linear equations AU = B that results from the numerical problem studied, as seen in Section 2. The index function I, while enumerating the unknown elements of the discrete sample of the domain, facilitates the determination of the numerical solution of the Poisson Problem in any dimension of the Euclidean space and, using it, it is possible to know, for example, in which line and in which column each non-null element of the matrix A of the above system is located.

The rest of the article is organized as follows: In Section 2, we formulate the finite difference method of 2n + 1 points, with n ≥ 1, to obtain a numerical solution of the Poisson equation, with Dirichlet boundary conditions. In Section 3, we demonstrate the convergence of the method. For this, its consistency and stability is ensured. The computational implementation of the method is discussed in Section 4. In Section 5, we present experimental results to validate the proposed method in the solution of some elliptic problems and, finally, in Section 6, we present the conclusions and future directions of the study.

2 THE FINITE DIFFERENCE METHOD IN DIMENSION N

First, consider u a real function of n real variables sufficiently differentiable, x = (x ₁ , . . . , x _n ) in ℝⁿ and h _j > 0, with j = 1, . . . , n. Using Taylor’s formula, we get the central difference formula

$$\begin{gathered} \frac{{\partial }^{2}u}{\partial x_j^2}\left(x\right)=\frac{u\left(x+h_j{e}_j\right)-2u\left(x\right)+u\left(x-h_j{e}_j\right)}{{\left(h_j\right)}^2} \\ -\frac{{\left(h_j\right)}^2}{12}\frac{{\partial }^{4}u}{\partial x_j^4}\left(x+({\varrho }_j-x_j)e_j\right), \end{gathered}$$ (2.1)

where ρ _j ∈ (x _j −h _j , x _j + h _j ) and e _j = (0,..., 0, 1, 0, ..., 0) is a vector of the canonical basis of ℝⁿ , with the jth coordinate equal to one and the others null, for j = 1, . . . , n.

To approximate second-order partial derivatives $\frac{{\partial }^{2}u}{\partial x_j^2}$, j = 1, . . . , n , by finite differences, we cover the region Ω ∪ ∂Ω with a mesh. The points of this mesh are denoted by $x^{i_1,\dots ,i_n}=(x_1^{i_1},\dots ,x_n^{i_n})$, where $x_j^{i_j}=r_j+i_j.h_j,h_j=\frac{s_j-r_j}{M_j}$ with i _j = 0, 1, . . . , M _j and j = 1, . . . , n.

We denote by Ω_δ the set of mesh points that are interior to Ω and by ∂ Ω_δ the set of mesh points that are on the border of Ω. Then, using the equations (1.1) and (2.1), we get the finite difference method of 2n + 1 points,

${\Delta }_{\delta }{U}_{i_1,\dots ,i_n}=-f\left(x^{i_1,\dots ,i_n}\right)\text{in}{\Omega }_{\delta }$ (2.2)

and

$U_{i_1,\dots ,i_n}=g\left(x^{i_1,\dots ,i_n}\right)\text{on}\partial {\Omega }_{\delta },$ (2.3)

where $U_{i_1,\dots ,i_n}$ is the numerical solution and ∆_δ is the discrete Laplacian operator given by

${\Delta }_{\delta }{U}_{i_1,\dots ,i_n}=\sum _{j=1}^n\frac{U_{i_1,\dots ,i_j+1,\dots ,i_n}-2U_{i_1,\dots ,i_n}+U_{i_1,\dots ,i_j-1,\dots ,i_n}}{{\left(h_j\right)}^2}.$

By substituting into the equation (2.2) each of the (M ₁ − 1) · (M ₂ − 1) · . . . · (M _n − 1) points from Ω_δ and using the equation (2.3), we get the linear system AU = B, with (M ₁ − 1) · (M ₂ − 1) · . . . · (M _n − 1) equations and the same number of unknowns, where U is the column vector given by

$$\begin{gathered} U={(U_{1,1,\dots ,1},\dots ,U_{M_1-1,1,\dots ,1},\dots ,U_{1,M_2-1,\dots ,M_n-1},\dots , \\ {U}_{M_1-1,M_2-1,\dots ,M_n-1})}^T, \end{gathered}$$

$A=\left[\begin{array}{ccccccccc}a& b_1& & b_n& & & & & \\ b_1& a& b_1& & \ddots & & 0& & \\ & \ddots & \ddots & \ddots & 0& \ddots & & & \\ b_n& & \ddots & \ddots & \ddots & & b_n& & \\ & \ddots & 0& \ddots & \ddots & \ddots & & & \\ 0& & \ddots & & b_1& a& b_1& & \\ & & & b_n& & b_1& a& & \\ & & & & & & & & \end{array}\right]$

and

$$\begin{gathered} B={(B_{1,1,\dots ,1},\dots ,B_{M_1-1,1,\dots ,1},\dots ,B_{1,M_2-1,\dots ,M_n-1},\dots , \\ {B}_{M_1-1,M_2-1,\dots ,M_n-1})}^T, \end{gathered}$$

where

$$\begin{gathered} B_{1,1,\dots ,1}=-f\left(x^{1,1,\dots ,1}\right)-\frac{g\left(x^{0,1,\dots ,1}\right)}{(h_1{)}^2}-\frac{g\left(x^{1,0,\dots ,1}\right)}{(h_2{)}^2}-\dots -\frac{g\left(x^{1,1,\dots ,0}\right)}{(h_n{)}^2}, \\ {B}_{M_1-1,1,\dots ,1}=-f\left(x^{M_1-1,1,\dots ,1}\right)-\frac{g\left(x^{M_1,1,\dots ,1}\right)}{(h_1{)}^2}-\frac{g\left(x^{M_1-1,0,\dots ,1}\right)}{(h_2{)}^2}-\dots -\frac{g\left(x^{M_1-1,1,\dots ,0}\right)}{(h_n{)}^2}, \\ {B}_{1,M_2-1,\dots ,M_n-1}=-f\left(x^{1,M_2-1,\dots ,M_n-1}\right)-\frac{g\left(x^{0,M_2-1,\dots ,M_n-1}\right)}{(h_1{)}^2}-\dots -\frac{g\left(x^{1,M_2-1,\dots ,M_n}\right)}{(h_n{)}^2} \end{gathered}$$

and

$$\begin{gathered} B_{M_1-1,M_2-1,\dots ,M_n-1}=-f\left(x^{M_1-1,M_2-1,\dots ,M_n-1}\right)-\frac{g\left(x^{M_1,M_2-1,\dots ,M_n-1}\right)}{(h_1{)}^2}-\dots \\ -\frac{g\left(x^{M_1-1,M_2-1,\dots ,M_n}\right)}{(h_n{)}^2}. \end{gathered}$$

The notation X ^T represents the transpose of the X matrix. The matrix A is a sparse matrix, with non-zero elements at 2n + 1 diagonals. The values a, b ₁, . . . , b _n that appear on the diagonals are the coefficients of the discretization of 2n + 1 point, with n ≥ 1, and are given by

$a=-2\sum _{j=1}^n\frac{1}{{\left(h_j\right)}^2}and{b}_j=\frac{1}{{\left(h_j\right)}^2},$

for j = 1, . . . , n. Each number b _j vanishes to zero in some positions of the matrix A. These positions, which depend on the dimension of Euclidean space and the value of each M _j , with j = 1, . . . , n, obey a specific rule and will be seen in Section 4, where an index function enumerates the unknowns of the problem and determines the position of each non-zero element in matrix A.

In the next Section, the convergence of the method will be verified.

3 METHOD CONVERGENCE

To show that the method is convergent, consistency and stability will be analyzed. Thus, for the following definitions and results, given a discrete function V: Ω_δ ∪ ∂Ω_δ → ℝ, consider the operator ∆_δ V: Ω_δ ∪∂Ω_δ → ℝ given by

${\Delta }_{\delta }V\left(x\right)=\sum _{j=1}^n\frac{V(x+h_j{e}_j)-2V\left(x\right)+V(x-h_j{e}_j)}{(h_j{)}^2},$

for all x ∈ Ω_δ ∪ ∂Ω_δ.

Definition 3.1. The local truncation error, denoted by ${\tau }_{i_1,\dots ,i_n}$ , is defined as

${\tau }_{i_1,\dots ,i_n}={\Delta }_{\delta }u\left(x^{i_1,\dots ,i_n}\right)+f\left(x^{i_1,\dots ,i_n}\right).$ (3.1)

The next lemma shows that the local truncation error decreases as we refine the mesh and gives us the convergence rate of this error.

Lemma 3.1.If the solution u of the equation (1.1) is differentiable up to order four in Ω, with limited fourth-order partial derivatives, then the local truncation error satisfies the inequality

$\left|{\tau }_{i_1,\dots ,i_n}\right|\le C\sum _{j=1}^n{\left(h_j\right)}^2,$

where C is a positive constant independent of h_j, with j = 1, . . . , n.

Proof. Applying 2n times the Taylor’s formula, we get the following developments around the point $x^{i_1,\dots ,i_n}$:

$$\begin{gathered} u(x^{i_1,\dots ,i_n}+h_j{e}_j)=u\left(x^{i_1,\dots ,i_n}\right)+h_j\frac{\partial u}{\partial x_j}\left(x^{i_1,\dots ,i_n}\right) \\ +\frac{(h_j{)}^2}{2!}\frac{{\partial }^{2}u}{\partial x_j^2}\left(x^{i_1,\dots ,i_n}\right)+\frac{(h_j{)}^3}{3!}\frac{{\partial }^{3}u}{\partial x_j^3}\left(x^{i_1,\dots ,i_n}\right) \\ +\frac{(h_j{)}^4}{4!}\frac{{\partial }^{4}u}{\partial x_j^4}(x^{i_1,\dots ,i_n}+({\varrho }_j^1-x_j^{i_j}\left)e_j\right) \end{gathered}$$

and

$$\begin{gathered} u(x^{i_1,\dots ,i_n}-h_j{e}_j)=u\left(x^{i_1,\dots ,i_n}\right)-h_j\frac{\partial u}{\partial x_j}\left(x^{i_1,\dots ,i_n}\right) \\ +\frac{(h_j{)}^2}{2!}\frac{{\partial }^{2}u}{\partial x_j^2}\left(x^{i_1,\dots ,i_n}\right)-\frac{(h_j{)}^3}{3!}\frac{{\partial }^{3}u}{\partial x_j^3}\left(x^{i_1,\dots ,i_n}\right) \\ +\frac{(h_j{)}^4}{4!}\frac{{\partial }^{4}u}{\partial x_j^4}(x^{i_1,\dots ,i_n}+({\varrho }_j^2-x_j^{i_j}\left)e_j\right), \end{gathered}$$

where ${\rho }_j^1\in \left(x_j^{i_j},x_j^{i_j}+h_j\right)$ and ${\rho }_j^2\in \left(x_j^{i_j}-h_j,x_j^{i_j}\right)$, with j = 1, . . . , n.

By replacing these expansions into (3.1), simplifying the similar terms and using the fact that the function u satisfies the equation (1.1) at the point $x^{i_1,\dots ,i_n}$, we get

$$\begin{gathered} {\tau }_{i_1,\dots ,i_n}=\frac{1}{4!}\sum _{j=1}^n\left(h_j{)}^2\frac{{\partial }^{4}u}{\partial x_j^4}\right(x^{i_1,\dots ,i_n}+({\varrho }_j^1-x_j^{i_j})e_j) \\ +\frac{1}{4!}\sum _{j=1}^n(h_j{)}^2\frac{{\partial }^{4}u}{\partial x_j^4}(x^{i_1,\dots ,i_n}+({\varrho }_j^2-x_j^{i_j}\left)e_j\right). \end{gathered}$$ (3.2)

Since the fourth-order partial derivatives are limited, it follows from (3.2) that

$\left|{\tau }_{i_1,\dots ,i_n}\right|\le C\sum _{j=1}^n{\left(h_j\right)}^2,$

for some positive constant C, independent of h _j , with j = 1, . . . , n.

Thus, we conclude that the method is of second order. This completes the proof. □

Definition 3.2. The global error, denoted by $e_{i_1,\dots ,i_n}$, is defined as

$e_{i_1,\dots ,i_n}=u\left(x^{i_1,\dots ,i_n}\right)-U_{i_1,\dots ,i_n},$ (3.3)

where $U_{i_1,\dots ,i_n}$ denotes the numerical solution of the Poisson problem calculated at the point $x^{i_1,\dots ,i_n}$.

In the following theorem, we demonstrate that the numerical method is stable. For this, we use the Discrete Maximum Principle ²⁵.

Theorem 3.1. (Discrete Maximum Principle) Consider V: Ω_δ ∪ ∂Ω_δ → ℝ a discrete function.

(i) If ∆_δ V (x) ≥ 0, for all x ∈ Ω_δ , then

$\underset{x\in {\Omega }_{\delta }}{\max }V\left(x\right)\le \underset{x\in \partial {\Omega }_{\delta }}{\max }V\left(x\right).$

(ii) If ∆_δ V (x) ≤ 0, for all x ∈ Ω_δ , then

$\underset{x\in {\Omega }_{\delta }}{\min }V\left(x\right)\ge \underset{x\in \partial {\Omega }_{\delta }}{\min }V\left(x\right).$

Proof.

(i) Suppose the maximum of V does not occur on the boundary of Ω_δ . This implies that there is P ₀ ∈ Ω_δ so that V (P ₀) = M ₀, with V (P) ≤ M ₀, for all P ∈ Ω_δ , and V (P) < M ₀, for any P ∈ ∂ Ω_δ . Now consider the points given by

$P_j^1=P_0+h_j{e}_j\text{and}{P}_j^2=P_0-h_j{e}_j,$

where j = 1, . . . , n. Hence,

${\Delta }_{\delta }V\left(P_0\right)=\sum _{j=1}^n\left(h_j{)}^{-2}\right[V\left(P_j^1\right)+V\left(P_j^2\right)]-2V(P_0\left)\sum _{j=1}^n\right(h_j{)}^{-2}.$

As ∆_δ V (P ₀) ≥ 0, we conclude that

$2V\left(P_0\right)\sum _{j=1}^n(h_j{)}^{-2}\le \sum _{j=1}^n(h_j{)}^{-2}\left[V\right(P_j^1)+V(P_j^2\left)\right].$

Consequently,

$M_0\le {\left(\sum _{j=1}^n(h_j{)}^{-2}\right)}^{-1}\sum _{j=1}^n(h_j{)}^{-2}\frac{V\left(P_j^1\right)+V\left(P_j^2\right)}{2}.$

Since V(Q) ≤ M ₀, for all Q ∈ Ω_δ ∪ ∂Ω_δ , we conclude that $V\left(P_j^1\right)=M_0$ and $V\left(P_j^2\right)=M_0$, for all j = 1, . . . , n. In fact, suppose $V\left(P_j^1\right)<M_0$ or $V\left(P_j^2\right)<M_0$, for some j = 1, . . . , n. It follows from (3.4) that

$M_0<{\left(\sum _{j=1}^n(h_j{)}^{-2}\right)}^{-1}\sum _{j=1}^n(h_j{)}^{-2}\frac{M_0+M_0}{2}=M_0.$

Which is absurd. Therefore, $V\left(P_j^1\right)=M_0$ and $V\left(P_j^2\right)=M_0$, for all j = 1, . . . , n.

This argument is repeated for each of the interior points $P_j^1$ and $P_j^2$ instead of P ₀. By repetition, each point of Ω_δ ∪ ∂Ω_δ appears as one of $P_j^1$ and $P_j^2$, for some corresponding P ₀.

Thus, we conclude that

$V\left(P\right)=M_0\text{, for all}P\in {\Omega }_{\delta }\cup \partial {\Omega }_{\delta }.$

But this contradicts the fact that V(P) < M ₀, for all P ∈ ∂ Ω_δ . Thus, the item (i) is established.

(ii) Firstly, note that

$max\left[-V\left(x\right)\right]=-minV\left(x\right)\text{and}{\Delta }_{\delta }\left(-V\right)=-{\Delta }_{\delta }V.$

Suppose ∆_δ V(x) ≤ 0, for all x ∈ Ω_δ . So,

${\Delta }_{\delta }(-V(x\left)\right)\ge 0,$

for all x ∈ Ω_δ . By the item (i),

$\underset{x\in {\Omega }_{\delta }}{\max }\left[-V\left(x\right)\right]\le \underset{x\in \partial {\Omega }_{\delta }}{\max }\left[-V\left(x\right)\right]$

and consequently,

$\underset{x\in {\Omega }_{\delta }}{\min }V\left(x\right)\ge \underset{x\in \partial {\Omega }_{\delta }}{\min }V\left(x\right).$

Which completes the proof. □

The next theorem gives us a bound for the solution of the equation (2.2).

Theorem 3.2. (A Priori estimate)Consider V: Ω_δ ∪ ∂Ω_δ → ℝ a discrete function. Then,

$\underset{x\in {\Omega }_{\delta }}{\max }\left|V\right(x\left)\right|\le \underset{x\in \partial {\Omega }_{\delta }}{\max }\left|V\right(x\left)\right|+\frac{r_1^2+s_1^2}{2}\underset{x\in {\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }V\right(x\left)\right|.$ (3.5)

Proof. Consider ψ : Ω_δ ∪ ∂Ω_δ → ℝ a discrete function defined by $\psi \left(x\right)=\frac{1}{2}{x}_1^2$, for all x = (x ₁, . . . , x _n ) ∈ Ω_δ ∪ ∂Ω_δ. Note that, for all x ∈ Ω_δ ∪ ∂Ω_δ,

$0\le \psi \left(x\right)\le \frac{r_1^2+s_1^2}{2}and{\Delta }_{\delta }\psi \left(x\right)=1.$

Consider V₊: Ω_δ ∪ ∂Ω_δ → ℝ and V₋: Ω_δ ∪ ∂Ω_δ → ℝ two discrete functions defined by

$V_{+}\left(x\right)=V\left(x\right)+N_0\psi \left(x\right)and{V}_{-}\left(x\right)=-V\left(x\right)+N_0\psi \left(x\right),$

where

$N_0=\underset{x\in {\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }V\right(x\left)\right|.$

So, for every x ∈ Ω_δ, we have that ∆_δ V ₊(x) = ∆_δ V(x) + N ₀ ≥ 0 and ∆_δ V ₋(x) = −∆_δ V(x) + N ₀ ≥ 0.

By applying the item (i) of the Theorem 3.1 to V ₊, we get

$$\begin{gathered} V\left(x\right)\le \underset{x\in {\Omega }_{\delta }}{\max }{V}_{+}\left(x\right)\le \underset{x\in \partial {\Omega }_{\delta }}{\max }\left[V\right(x)+N_0\psi (x\left)\right] \\ \le \underset{x\in \partial {\Omega }_{\delta }}{\max }V\left(x\right)+N_0\frac{r_1^2+s_1^2}{2} \\ \le \underset{x\in \partial {\Omega }_{\delta }}{\max }\left|V\right(x\left)\right|+\frac{r_1^2+s_1^2}{2}\underset{x\in {\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }V\right(x\left)\right|, \end{gathered}$$ (3.6)

for all x ∈ Ω_δ.

Using the Theorem 3.1, item (i), again, we deduce that

$-V\left(x\right)\le \underset{x\in \partial {\Omega }_{\delta }}{\max }\left|V\right(x\left)\right|+\frac{r_1^2+s_1^2}{2}\underset{x\in {\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }V\right(x\left)\right|,$ (3.7)

for all x ∈ Ω_δ. As a consequence of (3.6) and (3.7), we get (3.5). □

Note 3.1. In the proof of Theorem 3.2, we can replace $\frac{r_1^2+s_1^2}{2}$ in (3.5) for $\frac{r_j^2+s_j^2}{2}$ with j = 2, . . . , n, as long as we use one of the discrete functions ψ_j , defined by ${\psi }_j\left(x\right)=\frac{x_j^2}{2}$, for x = (x ₁, . . . , x _n ), in place of ψ.

Now, let’s get an estimate for the global error. Then, consider the discrete function $e_{i_1,\dots ,i_n}$ defined in (3.3). By the Theorem 3.2 , we deduce that

$\left|e_{i_1,\dots ,i_n}\right|\le \underset{\partial {\Omega }_{\delta }}{\max }\left|e_{i_1,\dots ,i_n}\right|+\frac{r_1^2+s_1^2}{2}\underset{{\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }{e}_{i_1,\dots ,i_n}\right|.$

Once we have, by the equation (2.3),

$e_{i_1,\dots ,i_n}=u\left(x^{i_1,\dots ,i_n}\right)-g\left(x^{i_1,\dots ,i_n}\right)=0$

over the boundary of Ω_δ, we get

$\left|e_{i_1,\dots ,i_n}\right|\le \frac{r_1^2+s_1^2}{2}\underset{{\Omega }_{\delta }}{\max }\left|{\Delta }_{\delta }{e}_{i_1,\dots ,i_n}\right|.$ (3.8)

From (2.2) and (3.1), we conclude that

$$\begin{gathered} {\Delta }_{\delta }{e}_{i_1,\dots ,i_n}={\Delta }_{\delta }u\left(x^{i_1,\dots ,i_n}\right)-{\Delta }_{\delta }{U}_{i_1,\dots ,i_n} \\ ={\Delta }_{\delta }u\left(x^{i_1,\dots ,i_n}\right)+f\left(x^{i_1,\dots ,i_n}\right) \\ ={\tau }_{i_1,\dots ,i_n}. \end{gathered}$$

Using this, (3.8) and the Lemma 3.1, we have that

$\left|e_{i_1,\dots ,i_n}\right|\le C\sum _{j=1}^n{\left(h_j\right)}^2,$ (3.9)

for some positive constant C, independent of h _j , with j = 1, . . . , n.

Note 3.2. We conclude, from Lemma 3.1 and from (3.9), that the numerical method is convergent.

The next result establishes the uniqueness of the solution of the system AU = B and is a consequence of the Discrete Maximum Principle.

Corollary 3.1. The resulting system of linear equations AU = B has a unique solution.

Proof. It is sufficient to show that the only solution of the homogeneous linear system AU = 0 is the trivial solution. For this, consider the homogeneous problem

$\left\{\begin{array}{l}-\Delta u=0,\text{in}\Omega ,\\ u=0,\text{on}\partial \Omega .\end{array}\right.$ (3.10)

The function u = 0 is the only solution to the above problem. Discretizing the problem (3.10), we obtain a homogeneous linear system for the unknowns $U_{i_1,\dots ,i_n}$. Since $U_{i_1,\dots ,i_n}$ is a solution of the difference equation (2.2), with f = 0, we conclude that

${\Delta }_{\delta }{U}_{i_1,\dots ,i_n}=0\text{in}{\Omega }_{\delta }\text{and}{U}_{i_1,\dots ,i_n}=0\text{on}\partial {\Omega }_{\delta }.$

By theorem 3.1, item (i), we have

$\underset{{\Omega }_{\delta }}{\max }{U}_{i_1,\dots ,i_n}\le \underset{\partial {\Omega }_{\delta }}{\max }{U}_{i_1,\dots ,i_n}=0.$

Applying again the Theorem 3.1, item (ii), we obtain

$\underset{{\Omega }_{\delta }}{\min }{U}_{i_1,\dots ,i_n}\ge \underset{\partial {\Omega }_{\delta }}{\min }{U}_{i_1,\dots ,i_n}=0.$

Thus, $U_{i_1,\dots ,i_n}=0$ in Ω_δ. Therefore, U = 0 is the only solution for the linear system AU = 0 and the proof of the corollary is complete. □

In the next section the numerical method, given by the equations (2.2) and (2.3), will be implemented computationally.

4 COMPUTER IMPLEMENTATION

We consider the function u with domain [r ₁ , s ₁] × · · · × [r _n , s _n ] ⊂ ℝⁿ and the discrete sample of (M ₁ + 1) · (M ₂ + 1) · . . . · (M _n + 1) domain points. Along each axis x _j , for j = 1, . . . , n, uniform sampling of points defines segments of length

$d{x}_j=\frac{s_j-r_j}{M_j}.$

Initially, we build a multidimensional matrix G with the dimension (M ₁ + 1)· (M ₂ + 1)·...· (M _n + 1), whose elements must be images from the function u, known on the domain boundary, that is,

$G[i_1,i_2,\dots ,i_n]=u(r_1+i_{1}d{x}_1,r_2+i_{2}d{x}_2,\dots ,r_n+i_{n}d{x}_n).$

Note that only the values of u on the domain boundary are known, that is, we know G[i ₁ , i ₂ , . . . , i _n ] only when i _j = 0 or i _j = M _j , for some j ∈ {1, 2, . . . , n}. Within the domain, when 0 < i _j < M _j for all j, the values G[i ₁ , i ₂ , . . . , i _n ] are unknown.

The next step is the construction of a multidimensional matrix F with dimension (M ₁ + 1)· (M ₂ + 1) · . . . · (M _n + 1), whose elements must be values of the Laplacian ∆u, that is,

$F[i_1,i_2,\dots ,i_n]=\Delta u(r_1+i_{1}d{x}_1,r_2+i_{2}d{x}_2,\dots ,r_n+i_{n}d{x}_n).$

In this case, only the values of ∆u in the domain interior are known, that is, we know F[i ₁ , i ₂ , . . . , i _n ] only when 0 < i _j < M _j , for all j ∈ {1, 2, . . . , n}.

Finally, each known element of matrix F is related to unknown elements of matrix G by the equation,

$$\begin{gathered} F[i_1,...,i_n]=\sum _{j=1}^n\left\{b_{j}G[i_1,...,i_j-1,...,i_n]+b_{j}G[i_1,...,i_j+1,...,i_n]\right\} \\ +aG[i_1,i_2,...,i_n], \end{gathered}$$ (4.1)

where $a=-2\sum _{j=1}^n\frac{1}{d{x}_j^2}$ and $b_j=\frac{1}{d{x}_j^2}$, for j = 1, ..., n.

The same equation (4.1) can be rewritten considering all the elements of the matrix G, in order to obtain an equation with (M ₁ − 1)...(M _n − 1) unknowns, given by

$F[i_1,...,i_n]=\sum _{l_1=1}^{M_1-1}...\sum _{l_n=1}^{M_n-1}\alpha (l_1,...,l_n)G[l_1,...,l_n],$ (4.2)

with i _j = 1, . . . , M _j − 1 and j = 1, . . . , n, where the coefficients α(l ₁ , ..., l _n ) are null, except for the following cases:

$\left\{\begin{array}{l}\alpha (i_1,i_2,...,i_n)=a;\\ \alpha (i_1,...,i_j-1,...,i_n)=b_j;\\ \alpha (i_1,...,i_j+1,...,i_n)=b_j.\end{array}\right.$ (4.3)

To construct a system of linear equations, the unknowns G[i ₁ , i ₂ , . . . , i _n ] are enumerated by an index function

$I:\{1,\dots ,M_1-1\}\times \{1,\dots ,M_2-1\}\times ...\times \{1,\dots ,M_n-1\}\to \mathrm{\mathbb{N}}$

defined by

$$\begin{gathered} I(i_1,\dots ,i_n)=(i_1-1)+(M_1-1)(i_2-1)+(M_1-1)(M_2-1)(i_3-1)+\dots \\ +(M_1-1)(M_2-1)\dots (M_{n-1}-1)(i_n-1). \end{gathered}$$ (4.4)

Thus, we have (M ₁ − 1) . . . (M _n − 1) equations, with the same number of unknowns, determining a linear system AU = B, whose solution vector U provides the searched unknowns, that is,

$G[i_1,i_2,...,i_n]=U\left[I\right(i_1,i_2,...,i_n\left)\right].$ (4.5)

The matrix A is a sparse matrix with null elements, except for:

$\left\{\begin{array}{l}A\left[I\right(i_1,i_2,...,i_n),I(i_1,i_2,...,i_n\left)\right]=a;\\ A\left[I\right(i_1,i_2,...,i_n),I(i_1,...,i_j-1,...,i_n\left)\right]=b_j\text{if}{i}_j>1\text{and}\\ A\left[I\right(i_1,i_2,...,i_n),I(i_1,...,i_j+1,...,i_n\left)\right]=b_j\text{if}{i}_j<M_j-1,\\ \text{for}j=1,...,n.\end{array}\right.$ (4.6)

Note that the index I(i ₁ , . . . , i _n ) indicates that the coefficient α(i ₁ , i ₂ , ..., i _n ) = a of the equation (4.2) associated with the unknown G[i ₁ , . . . , i _n ] will appear in the line I(i ₁ , . . . , i _n ) and column I(i ₁ , . . . , i _n ) of the matrix A. The indices I(i ₁ , . . . , i _j − 1, . . . , i _n ) and I(i ₁ , . . . , i _j + 1, . . . , i _n ) indicate that the coefficients α(i ₁ , . . . , i _j −1,..., i _n ) = b _j and α(i ₁ , . . . , i _j +1,..., i _n ) = b _j , j = 1, . . . , n, will appear in the line I(i ₁ , . . . , i _n ) and columns I(i ₁ , . . . , i _j − 1, . . . , i _n ) and I(i ₁ , . . . , i _j + 1, . . . , i _n ) of the matrix A, respectively. In the other positions of the line I(i ₁ , . . . , i _n ) the coefficients α(l ₁ , ..., l _n ) of the equation (4.2) are all null. Figure 1 illustrates the final position of each element in matrix A, for the case n = 3.

Figure 1:
Illustration of the coefficient matrix of the system AU = B with indication of the diagonals corresponding to the coefficients a, b ₁ , b ₂ , b ₃.

Now, notice that the vector B, which has (M ₁ − 1) . . . (M _n − 1) independent terms, is not given just by values of the Laplacian F[i ₁ , i ₂ , . . . , i _n ], but known boundary values G[i ₁ , i ₂ , . . . , i _n ], for i _j = 0 or i _j = M _j , are incorporated into the independent terms as follows:

$$\begin{gathered} B\left[I\right(i_1,i_2,\dots ,i_n\left)\right]=F[i_1,i_2,\dots ,i_n] \\ -\sum _{j=1}^n{s}_{j}G[i_1,\dots ,i_j-1,\dots ,i_n] \\ -\sum _{j=1}^n{t}_{j}G[i_1,\dots ,i_j+1,\dots ,i_n], \end{gathered}$$

where the coefficients s _j and t _j , for j = 1, . . . , n, are all null, except in the cases:

$\left\{\begin{array}{l}{s}_j=b_j\text{if}{i}_j=1;\\ t_j=b_j\text{if}{i}_j=M_j-1.\end{array}\right.$

Using the implementation details above, below is a code example, in Python, for solving the following Problem in ℝ³:

$\left\{\begin{array}{l}\Delta u=12x^2+12y^2+12z^2,\text{in}\Omega ,\\ u=x^4+y^4+z^4,\text{on}\partial \Omega ,\end{array}\right.$

where Ω = (−1,1)×(−1,1)×(−1,1).

5 EXPERIMENTS

In this section, we present numerical experiments to evaluate the effectiveness of the proposed method. To illustrate the implementation and computational performance, we consider a set of test cases. In each case, the integration domain is a rectangular block discretized using a uniform mesh. The tables below display the average error, (e), calculated according to the following formula:

$e=\frac{{\displaystyle \sum _{i_1=0}^{M_1}}{\displaystyle \sum _{i_2=0}^{M_2}}\dots {\displaystyle \sum _{i_n=0}^{M_n}}\left|G(i_1,\dots ,i_n)-u(r_1+i_{1}d{x}_1,\dots ,r_n+i_{n}d{x}_n)\right|}{\left(M_1+1\right).\left(M_2+1\right)\dots \left(M_n+1\right)},$ (5.1)

where G is the approximate solution by the method and u is the exact solution.

We also show the evolution of the error in the numerical solution, for each problem, as the resolution of the considered mesh increases. The solution for the linear system AU = B was obtained using spsolve package for sparse matrix in Python. The numerical solution for the examples 1 and 2, considered next, is obtained using the procedure in the section 4.

Example 1. Consider the following problem in ℝ²

$\left\{\begin{array}{l}\Delta u=-sin\left(x\right)-sin\left(y\right),\text{in}\Omega ,\\ \\ u=sin\left(x\right)+sin\left(y\right),\text{on}\partial \Omega ,\end{array}\right.$

where Ω = (−π, π) × (−π, π). The exact solution to this problem is given by u(x, y) = sin(x) + sin(y).

In Table 1 we present values of the average error, given by the formula (5.1), for some resolutions M × M of the mesh, with M = 4, 8, 16, 32, 64, 128, 256, 512, 1024. In this table one can observe the convergence of the error. The curve fitting for this data indicates that the convergence is approximately of order O(1.978), which is close to quadratic, as expected.

Table 1:
Error in the solution of the problem Example 1, for some resolutions.

Resolution	Error (e)
4x4	9.1810e-02
8x8	2.6437e-02
16x16	6.8522e-03
32x32	1.7279e-03
64x64	4.3293e-04
128x128	1.0829e-04
256x256	2.7080e-05
512x512	6.7730e-06
1024x1024	1.6960e-06

Example 2. Consider the following problem in ℝ³

$\left\{\begin{array}{l}\Delta u=12x^2+12y^2+12z^2,\text{in}\Omega ,\\ \\ u=x^4+y^4+z^4,\text{on}\partial \Omega ,\end{array}\right.$

where Ω = (−1, 1) × (−1, 1) × (−1, 1).

In this case, with the domain in ℝ³, the resolution M ×M ×M leads to the sampling of M ³ points from the domain. Due to the limitation of memory allocation capacity, we only considered up to resolution 50 × 50 × 50, but it is already possible to notice the convergence (Table 2). For this data, the convergence order is approximately O(1.915), which is also close to quadratic, as expected. Note that in this case the matrix A, from the system AU = B, has 49³ = 117, 649 rows and columns.

Table 2:
Error in the solution of the problem Example 2, for some resolutions.

Resolution	Error (e)
5 x 5 x 5	6.1741e-02
10 x 10 x 10	1.8259e-02
15 x 15 x 15	8.3816e-03
20 x 20 x 20	4.7689e-03
30 x 30 x 30	2.1370e-03
40 x 40 x 40	1.2055e-03
50 x 50 x 50	7.7261e-04

Figure 1 illustrates the matrix A of the system AU = B for the mesh with resolution 5 × 5 × 5. In this case A has 4³ = 64 rows and columns. The diagonals corresponding to the coefficients a, b ₁ , b ₂ , b ₃ are indicated in the figure, as automatically defined in Equation 4.6.

Numerical experiments confirm the consistency and stability of the method, ie, the convergence of the numerical solution to the exact solution of the problem while the mesh is refined as demonstrated in the section 3. The code in Python presented in the section 4 uses the index function to solve the problem of Example 2, in the dimension n = 3, and can be modified to consider larger dimensions, as detailed in the same Section.

6 CONCLUSIONS

This article presents the numerical solution for the n-dimensional Poisson problem, for all n ≥ 1. We propose a new approach based on an index function that plays a central role in obtaining the numerical solution, as it facilitates generalization and supports computational implementation in any dimension. The index function simplifies the construction of both the coefficient matrix and the right-hand side vector of the sparse linear system that arises in the modeling process. Moreover, this function can be easily adapted, since the access path to array elements is not unique and depends on the chosen enumeration, without affecting the final solution. This concept can be extended to obtain numerical solutions for other elliptic problems.

As future work, we intend to investigate how the use of the index function can be adapted and extended to more complex scenarios, including problems with irregular boundaries and nonlinear Poisson-type equations.

Acknowledgments

The authors are grateful to the support of Unimontes and IFNMG.

REFERENCES

• ¹ B.M. Averick & J.M. Ortega. Fast solution of nonlinear Poisson-type equations. SIAM Journal on Scientific Computing, 14(1) (1993), 44-48.
• ² K. Balakrishnan & P.A. Ramachandran. A particular solution Trefftz method for non-linear Poisson problems in heat and mass transfer. Journal of Computational Physics, 150(1) (1999), 239-267.
• ³ K. Balakrishnan & P.A. Ramachandran. Osculatory interpolation in the method of fundamental solution for nonlinear Poisson problems. Journal of Computational Physics , 172(1) (2001), 1-18.
• ⁴ R.F. Boisvert. Families of high order accurate discretizations of some elliptic problems. SIAM Journal on Scientific and Statistical Computing, 2(3) (1981), 268-284.
• ⁵ W.R. Bowen & P.M. Williams. Finite difference solution of the 2-dimensional Poisson-Boltzmann equation for spheres in confined geometries. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 204(1-3) (2002), 103-115.
• ⁶ J.A. Cuminato & M.M. Junior. “Discretização de equações diferenciais parciais: técnicas de diferenças finitas”. SBM, Rio de Janeiro, 1 ed. (2013).
• ⁷ L.P. da Silva, M.A.V. Pinto & L.K. Araki. Higher-order methods for the Poisson equation obtained with geometric multigrid and completed Richardson extrapolation. Computational and Applied Mathematics, 43(7) (2024), 395.
• ⁸ L.P. da Silva, B.B. Rutyna, A.R. Santos Righi & M.A. Villela Pinto. High order of accuracy for Poisson equation obtained by grouping of repeated Richardson extrapolation with fourth order schemes. Computer Modeling in Engineering & Sciences, 128(2) (2021), 699-715.
• ⁹ G.E. Fasshauer. Newton iteration with multiquadrics for the solution of nonlinear PDEs. Computers & Mathematics with Applications, 43(3-5) (2002), 423-438.
• ¹⁰ J.H. Ferziger & M. Perić. “Computational methods for fluid dynamics”. Springer (2002).
• ¹¹ M.M. Gupta, R.P. Manohar & J.W. Stephenson. A single cell high order scheme for the convectiondiffusion equation with variable coefficients. International Journal for Numerical Methods in Fluids, 4(7) (1984), 641-651.
• ¹² G.M.A.M.L. Guta. Solution of Two Dimensional Poisson Equation Using Finite Difference Method with Uniform and Non-uniform Mesh Size. (2019).
• ¹³ H. Johansen & P. Colella. A Cartesian grid embedded boundary method for Poisson’s equation on irregular domains. Journal of Computational Physics , 147(1) (1998), 60-85.
• ¹⁴ J.J. Kasab, S.R. Karur & P. Ramachandran. Quasilinear boundary element method for nonlinear Poisson type problems. Engineering Analysis with Boundary Elements, 15(3) (1995), 277-282.
• ¹⁵ L. Lapidus & G.F. Pinder. “Numerical solution of partial differential equations in science and engineering”. John Wiley & Sons (2011).
• ¹⁶ Z. Li, C. Pao & Z. Qiao. A finite difference method and analysis for 2D nonlinear Poisson-Boltzmann equations. Journal of Scientific Computing, 30 (2007), 61-81.
• ¹⁷ B. Mebrate, P.R. Koya et al. Numerical solution of a two dimensional poisson equation with dirichlet boundary conditions. American Journal of Applied Mathematics, 3(6) (2015), 297-304.
• ¹⁸ H. Moghaderi, M. Dehghan & M. Hajarian. A fast and efficient two-grid method for solving d-dimensional poisson equations. Numerical Algorithms, 72 (2016), 483-537.
• ¹⁹ F.M. Okoro & A.E. Owoloko. Compact finite difference schemes for Poisson equation using direct solver. Journal of Mathematics and Technology, ISSN, (2010), 2078-0257.
• ²⁰ P. Pandey. A higher accuracy exponential finite difference method for the numerical solution of second order elliptic partial differential equations. J. Math. Comput. Sci., 3(5) (2013), 1325-1334.
• ²¹ J.B. Rosser. Nine-point difference solutions for Poisson’s equation. Computers & Mathematics with Applications , 1(3-4) (1975), 351-360.
• ²² J.B. Rosser. Finite-difference solution of Poisson’s equation in rectangles of arbitrary proportions. Zeitschrift für angewandte Mathematik und Physik ZAMP, 28 (1977), 185-196.
• ²³ G.H. Shortley & R. Weller. The numerical solution of Laplace’s equation. Journal of Applied Physics, 9(5) (1938), 334-348.
• ²⁴ G.D. Smith. “Numerical solution of partial differential equations: finite difference methods”. Oxford university press (1985).
• ²⁵ J.C. Strikwerda. “Finite difference schemes and partial differential equations”. SIAM (2004).
• ²⁶ H. Wang, Q.H. Qin & X.P. Liang. Solving the nonlinear Poisson-type problems with F-Trefftz hybrid finite element model. Engineering Analysis with Boundary Elements , 36(1) (2012), 39-46.
• ²⁷ Y. Wang & J. Zhang. Sixth order compact scheme combined with multigrid method and extrapolation technique for 2D Poisson equation. Journal of Computational Physics , 228(1) (2009), 137-146.
• ²⁸ T. Zhu, J. Zhang & S.N. Atluri. A meshless local boundary integral equation (LBIE) method for solving nonlinear problems. Computational mechanics, 22(2) (1998), 174-186.

Data availability

Do not apply.

*

Corresponding author: Antonio Wilson Vieira - E-mail: [e-mail]

Associate editor:

Pedro Lima

Publication Dates

• Publication in this collection
  28 July 2025
• Date of issue
  2025

History

• Received
  12 Apr 2022
• Accepted
  05 June 2025