The Numerical Solution of Poisson Equation with Dirichlet Boundary Conditions

Peng Guo

The Business School, Shanghai Dianji University, Shanghai, China.

DOI: 10.4236/jamp.2021.912194                

Abstract

This work mainly focuses on the numerical solution of the Poisson equation with the Dirichlet boundary conditions. Compared to the traditional 5-point finite difference method, the Chebyshev spectral method is applied. The numerical results show the Chebyshev spectral method has high accuracy and fast convergence; the more Chebyshev points are selected, the better the accuracy is. Finally, the error of two numerical results also verifies that the algorithm has high precision.

Keywords

Poisson Equation, Finite Difference Method, Chebyshev, Spectral Method

1. Introduction

The Laplace equation $\nabla u=0$ can be dated back to 1782 when the French mathematician P.S. Laplace discussed the problem of the gravitational field [1] [2]. About thirty years later, S.D. Poisson pointed out that if the density of gravitational field is considered, the Laplace equation should have a new form i.e. $\nabla u=f$, which we called Poisson equation. Now Poisson equation has been applied to modelling the temperature distribution of stable temperature field with a stable heat source or without internal heat source, the stable non-rotating flow of incompressible fluid in hydrodynamics, etc. In recent several years, some researchers find that the 2-dimensional Poisson equation with Dirichlet boundary condition is a good tool to cope with seamless image composite problems [3] [4] [5] [6] [7].

In the following part, we will focus on Poisson equation with the Dirichlet boundary condition [8].

$\{\begin{array}{ll}-\nabla u=f\left(x,y\right)\hfill & \left(x,y\right)\in \Omega \hfill \\ u\left(x,y\right)=g\left(x,y\right)\hfill & \left(x,y\right)\in \partial \Omega \hfill \end{array}$ (1)

The above Equation (1) is the 2-dimensional Poisson equation, where $\nabla u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$, $g\left(x\mathrm{,}y\right)$ is the boundary condition.

In general, the Poisson equation is hard to get the analytical solution, only a few can find the exact solution. Therefore, the numerical algorithm is a good way to deal with this problem. For the second-order parabolic type differential equation, the finite difference method is a popular discretization method. The precision of the traditional 5-point difference method is not good enough. In order to improve the accuracy, some scholars proposed a new method based on 5-point difference scheme, such as 9-point difference scheme. The error of the 9-point difference scheme is fourth-order. But if we want to continue to improve the accuracy of the algorithm based on the 9-point difference method, it will be very difficult. Compared to the finite difference method the spectral method has high accuracy and less computation. Especially for multidimensional problems, such as 2-dimensional or 3-dimensional Poisson equation if we use the finite difference method, we need to calculate so many nodes the amount of calculation is very large. In addition to the finite difference method and spectral method, the finite element method is also an effective method to deal with differential equations. For the in-depth discussion of these methods, we can refer to [8] - [23].

The spectral method is a well-developed algorithm, which has infinite order accuracy and exponential order convergence speed in theory [17] [18] [19]. In this work, we will use Chebyshev spectral method to find the numerical solution of 2-dimensional Poisson equation.

2. Preliminaries

Some basic contents including Chebyshev polynomials and Chebyshev points will be introduced in this part. The Chebyshev points will be used to construct the differentiation matrices, which is the key point to obtain high-precision solutions.

Definition 1. In the interval $\left[-\mathrm{1,1}\right]$, with the weight $\omega \left(x\right)=\frac{1}{\sqrt{1-x^2}}$, the polynomials with the following relations are called Chebyshev polynomials [20].

$\{\begin{array}{l}{T}_0\left(x\right)=1,\\ T_1\left(x\right)=x,\\ T_{n+1}\left(x\right)=2x{T}_n\left(x\right)-T_{n-1}\left(x\right),n\ge 1.\end{array}$ (2)

where $T_n\left(x\right)\mathrm{,}n\ge 0$ are the Chebyshev polynomials.

A big difference between Chebyshev polynomials and other polynomials is that Chebyshev polynomials have an explicit relation, i.e. $T_n\left(x\right)=\cos \left(n\cdot \arccos \left(x\right)\right)$. The zeros of Chebyshev polynomials $x_k=\cos \frac{k\pi }{n},n=0,1,\cdots ,n$ are called Chebyshev points. These Chebyshev points are highly clustered on both sides of the interval. This feature of the Chebyshev points can help one overcome Runge phenomenon in interpolation. The Chebyshev polynomials satisfy the following orthogonal properties with weight $\frac{1}{\sqrt{1-x^2}}$

Property 1

${\displaystyle {\int }_{-1}^1}\frac{1}{\sqrt{1-x^2}}{T}_m\left(x\right)T_n\left(x\right)\text{d}x=\{\begin{array}{ll}0,\hfill & m\ne n,\hfill \\ \frac{\pi }{2},\hfill & m=n=0,\hfill \\ \pi ,\hfill & m=n\ne 0.\hfill \end{array}$ (3)

First-order derivative of Chebyshev polynomials ${T^{\prime }}_n\left(x\right)$ and the original Chebyshev polynomials $T_n\left(x\right)$ satisfy the following relations.

Property 2

$\{\begin{array}{l}{T}_0\left(x\right)={T^{\prime }}_1\left(x\right),\\ T_1\left(x\right)=\frac{1}{4}{T^{\prime }}_2\left(x\right),\\ T_n\left(x\right)=\frac{1}{2\left(n+1\right)}{T^{\prime }}_{n+1}\left(x\right)-\frac{1}{2\left(n-1\right)}{T^{\prime }}_{n-1}\left(x\right),n\ge 2.\end{array}$ (4)

There are some other polynomials such as Legendre polynomials, Jacobi polynomials, etc. For a more detailed discussion about the properties of Chebyshev and other polynomials, one can refer to [17] [19].

3. Numerical Algorithms for 2-Dimensional Poisson Equation

In this part, we will show two different schemes to solve the 2-dimensional Poisson equation. The finite difference method with five points will be applied first, then the Chebyshev spectral method will be considered. The 2-dimensional Poisson equation we will discuss is as follows,

$\{\begin{array}{l}-\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=f\left(x,y\right),\left(x,y\right)\in \left[a,b\right]\times \left[c,d\right]\\ u\left(a,y\right)=g_1\left(y\right),\\ u\left(b,y\right)=g_2\left(y\right),\\ u\left(x,c\right)=g_3\left(x\right),\\ u\left(x,d\right)=g_4\left(x\right).\end{array}$ (5)

1) 5-point finite difference method

First, we will give a discretization in x and y directions. In x direction we choose $N_1+1$ points and give an equidistant discretization, that is $h_1=\frac{b-a}{N_1}$. Then the points in the interval $\left[a\mathrm{,}b\right]$ can be expressed as $x_i=a+i{h}_1$, where $i=0,1,2,\cdots ,N_1$. Similarly, in y direction, the points in $\left[c\mathrm{,}d\right]$ can be expressed as $y_j=a+j{h}_2,j=0,1,2,\cdots ,N_2$, where $h_2=\frac{d-c}{N_2}$.

Both in x and y directions we use the central difference method respectively, as shown in Figure 1, we have

Figure 1. 5-point finite difference method.

${\left(\frac{{\partial }^{2}u}{\partial x^2}\right)}_{\left(x_i,y_j\right)}=\frac{u_{i+1,j}+u_{i-1,j}-2u_{i,j}}{h_1^2}+O\left(h_1^2\right),$ (6)

${\left(\frac{{\partial }^{2}u}{\partial y^2}\right)}_{\left(x_i,y_j\right)}=\frac{u_{i,j+1}+u_{i,j-1}-2u_{i,j}}{h_2^2}+O\left(h_2^2\right).$ (7)

If we choose the same step size, i.e. $h_1=h_2=h$, then for the problem we have the following algorithm,

$-\frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}}{h^2}=f_{i,j},$ (8)

where the corresponding truncation error is $O\left(h_1^2+h_2^2\right)$, if $h_1=h_2=h$, the truncation error can be written as $O\left(h^2\right)$.

The size of the nodes on grids is $\left(N_1+1\right)\times \left(N_2+1\right)$. Considering that the boundary values are known, we arrange the remaining nodes in columns to form a new column vector. Let’s $U={\left(u_{1,1},u_{2,1},\cdots ,u_{N_1-1,1},\cdots ,u_{N_1-1,N_2-1}\right)}^{\text{T}}$, the 2-dimensional Poisson equation can be written in a matrix form $AU=F$, where $A_{N_1-\mathrm{1,}{N}_2-1}$ is the coefficient matrix, F is a matrix on the right-hand side of the relation (8). The detailed expression of the coefficient matrix A is as follows,

$A=h^{-2}\left(\begin{array}{ccccc}B& -I& & & \\ -I& B& -I& & \\ & \ddots & \ddots & \ddots & \\ & & -I& B& -I\\ & & & -I& B\end{array}\right),$ (9)

where $B_{N_1-\mathrm{1,}{N}_1-1}$ has the following expression,

$B=\left(\begin{array}{ccccc}4& -1& & & \\ -1& 4& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 4& -1\\ & & & -1& 4\end{array}\right)\mathrm{.}$ (10)

2) Chebyshev spectral method

The Chebyshev points are $x_k=-\cos \frac{k\pi }{n},n=0,1,\cdots ,n$. Obviously, all the Chebyshev points are in the interval $\left[-\mathrm{1,1}\right]$. Before we use the Chebyshev spectral method, we should map the interval $\left[-\mathrm{1,1}\right]$ into the interval $\left[a\mathrm{,}b\right]$ with the following transformation,

$x=\frac{b+a}{2}+\frac{b-a}{2}{x}_k.$ (11)

Similarly, for $y\in \left[c\mathrm{,}d\right]$, if we use the transformation

$y=\frac{d+c}{2}+\frac{d-c}{2}{x}_k,$ (12)

the interval $\left[-\mathrm{1,1}\right]$ can also be mapped into $\left[c\mathrm{,}d\right]$. The differentiation matrices $D^1$ for the interval $\left[-\mathrm{1,1}\right]$ can be derived directly from the explicit expression of Chebyshev polynomials. But the expression of the first-order differentiation matrices is slightly different, i.e. $D_{\left[a\mathrm{,}b\right]}^1=\frac{b-a}{2}{D}^1$, $D_{\left[c\mathrm{,}d\right]}^1=\frac{d-c}{2}{D}^1$. For more detailed proofs or properties, one can refer to [17].

$\{\begin{array}{l}{D}_{0,0}^1=-D_{N,N}^1=\frac{2N^2+1}{6},\\ D_{i,i}^1=\frac{-x_i}{2\left(1-x_i^2\right)},i=1,2,\cdots ,N-1,\\ D_{i,j}^1=\frac{c_i{\left(-1\right)}^{i+j}}{c_j\left(x_i-x_j\right)},i\ne j,i,j=0,1,\cdots ,N,\end{array}$ (13)

where

$c_i=\{\begin{array}{l}2,i=0,i=N\\ 1,\text{otherwise}.\end{array}$ (14)

The second-order derivative can be computed easily by $D^2={\left(D^1\right)}^2$. If we arrange the matrixu in a column vector, $\frac{{\partial }^{2}u}{\partial x^2}$ can be writen in the matrix form $\left(D^2\otimes I\right)u$, and $\frac{{\partial }^{2}u}{\partial y^2}$ can be written in the matrix form $\left(I\otimes D^2\right)u$. Let’s L denote $\left(D^2\otimes I+I\otimes D^2\right)$, where $\otimes$ is the Kronecker product. Then for the 2-dimensional Poisson equation, the Chebyshev spectral algorithm can be given as follows,

$-LU=F.$ (15)

where U and F are similar as given in 5-point difference method. The numerical results can be computed directly. Jie et al. have proved the error estimate for the Chebyshev spectral method with Dirichlet boundary conditions [17].

4. Numerical Examples

In this part, the proposed algorithm will be employed to solve 2-dimensional Poisson equation with Dirichlet boundary conditions.

Example 1. Consider the following 2-dimensional Poisson equation

$\{\begin{array}{ll}-\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)=-{\text{e}}^{xy}\left(x^2+y^2\right),\hfill & \left(x,y\right)\in \left[0,1\right]\times \left[0,1\right],\hfill \\ u\left(x,0\right)=1,\hfill & 0\le x\le 1\hfill \\ u\left(x,1\right)={\text{e}}^{-x},\hfill & 0\le x\le 1\hfill \\ u\left(0,y\right)=1,\hfill & 0\le y\le 1\hfill \\ u\left(1,y\right)={\text{e}}^{-y},\hfill & 0\le y\le 1\hfill \end{array}$ (16)

where the exact solution is $u\left(x\mathrm{,}y\right)={\text{e}}^{-xy}$.

Figure 2 shows the numerical solution of Example 1 with five points finite difference method. Figure 3 shows the absolute error at different points in the interval $\left[\mathrm{0,1}\right]\times \left[\mathrm{0,1}\right]$. The biggest error is about 8.5 × 10⁻⁶. Table 1 is the detailed error when x and y take different values.

Table 2 shows the absolute error when x and y take different values. The biggest absolute error is 3.7 × 10⁻⁸. From the absolute error tables obtained by the two methods, we find that the accuracy of the spectral method is better than that of the finite difference method. The more the Chebyshev points we use, the higher accuracy of the numerical results we get. Theoretically, if we use enough Chebyshev points, we can obtain the numerical solutions with arbitrary accuracy.

Figure 4 is the numerical solution that we get with the Chebyshev spectral method. Figure 4 is very similar to Figure 2, but from the above two error tables, we find that the solution obtained by the spectral method is more consistent with the exact solution.

Figure 2. The numerical solution with finite difference method.

Figure 3. The absolute error with finite difference method.

Table 1. The errors of Example 1 with finite difference method.

Table 2. The errors of Example 1 with spectral method.

Figure 4. The numerical solution with Chebyshev spectral method.

Figure 5 is the error of the solution obtained by the Chebyshev spectral method. From this figure, we can see that the overall error of the solution obtained by the Chebyshev spectrum method is relatively small compared with the difference method.

Example 2. Consider the following 2-dimensional Laplace equation,

$\{\begin{array}{ll}\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=0,\hfill & \left(x,y\right)\in \left[0,1\right]\times \left[0,1\right],\hfill \\ u\left(x,0\right)=\sin \pi x,\hfill & 0\le x\le 1\hfill \\ u\left(x,1\right)={\text{e}}^{-x}\sin \pi x,\hfill & 0\le x\le 1\hfill \\ u\left(0,y\right)=0,\hfill & 0\le y\le 1\hfill \\ u\left(1,y\right)=0,\hfill & 0\le y\le 1\hfill \end{array}$ (17)

where the exact solution is $u\left(x\mathrm{,}y\right)={\text{e}}^{-\pi y}\sin \pi x$.

Figure 6 and Figure 7 show the absolute error of example 2 with finite difference method and Chebyshev spectral method respectively. When we choose the nodes $n_x=n_y=51$ in x and y directions, for the finite difference method the maximum error is 1.16 × 10⁻⁴ as shown in Figure 6. For the same example, if we use spectral method with 51 × 51 Chebyshev points in the interval $\left[\mathrm{0,1}\right]\times \left[\mathrm{0,1}\right]$, the maximum error is 2.33 × 10⁻⁷ as shown in Figure 7. From the absolute error of example 2, we find that the accuracy of spectral method is much higher than that of finite difference method.

Figure 5. The absolute error of Chebyshev spectral method.

Figure 6. The absolute error with finite difference method.

The second example is the Laplace equation which is a special case of Poisson equation. Table 3 shows the absolute error of example with finite difference method and Chebyshev spectral method. $Erro{r}_f$ represents the maximum absolute error with finite difference method and $Erro{r}_s$ represents the maximum absolute error with Chebyshev spectral method. We select the same number of nodes, i.e. $n_x=n_y=10,30,50,70,90$, and then compare the absolute errors

Figure 7. The absolute error with Chebyshev spectral method.

Table 3. The errors of Example 2.

of the two methods. From Table 3, we can clearly find that the accuracy of spectral method is better than that of finite difference method.

5. Conclusion

In this work, we study the 2-dimensional Poisson equation with Dirichlet boundary conditions. We use the five-point difference method and Chebyshev spectral method to solve the corresponding two-dimensional Poisson equation. For each method, we give the corresponding numerical algorithm. Finally, we give the numerical solution of the corresponding algorithm through a numerical case. We also obtain the absolute error respectively. The absolute error of two methods reveals that the accuracy of the spectral method is better than that of the difference method.

Acknowledgements

This research was funded by the Humanity and Social Science Youth Foundation of Ministry of Education (No. 18YJC630120), the Applied Mathematics of Shanghai DianJi University (No. 16JCXK02) and the Key Projects of Shanghai Soft Science Research Program (No. 18692106600).

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Abramowitz, M. and Stegun, I.A. (1972) Handbook of Mathematical Functions. Dover, New York.
[2]	Wang, Z.X. and Guo, D.R. (1982) Special Functions. World Scientific Press, Singapore.
[3]	Hao, W. and Dan, X. (2009) Improved Poisson Image Editing and Its Implementation on GPU. 2009 IEEE 10th International Conference on Computer-Aided Industrial Design & Conceptual Design, Wenzhou, 26-29 November 2009, 1044-1048. https://doi.org/10.1109/CAIDCD.2009.5375343
[4]	Morel, J.M., Petro, A.B. and Sbert, C. (2012) Fourier Implementation of Poisson Image Editing. Pattern Recognition Letters, 33, 342-348. https://doi.org/10.1016/j.patrec.2011.10.010
[5]	Martino, J.M.D., Facciolo, G. and Meinhardt-Llopis, E. (2016) Poisson Image Editing. Image Processing On Line, 5, 300-325. https://doi.org/10.5201/ipol.2016.163
[6]	Tian, C., Li, H. and Gao, X. (2018) Photo-Realistic 2D Expression Transfer Based on FFT and Modified Poisson Image Editing. Neurocomputing, 309, 1-10. https://doi.org/10.1016/j.neucom.2018.03.045
[7]	Kruse, A.W., Alenin, A.S. and Tyo, S. (2020) Contrast-Based Phase Angle of Polarization Image Construction Using Poisson Image Editing. Conference: Computational Optical Sensing and Imaging, Washington DC, 22-26 June 2020, CTh4C.5. https://doi.org/10.1364/COSI.2020.CTh4C.5
[8]	Sauer, T. (2012) Numerical Analysis. 2nd Edition, Addison Wesley Longman, Boston.
[9]	Dou, M.Y., Li, H.Y. and Liu, X.W. (2017) An Inertial Proximal Peaceman-Rachford Splitting Method. Scientia Sinica: Mathematica, 47, 333-348. https://doi.org/10.1360/N012016-00134
[10]	Sweilam, N.H., Moharram, H.M. and Ahmed, S. (2012) On the Parallel Iterative Finite Difference Algorithm for 2-D Poisson equation with MPI Cluster. International Conference on Informatics & Systems IEEE, Giza, 14-16 May 2012, MM-78-MM-85.
[11]	Yang, F. and Fu, C.L. (2010) The Inverse Problem of the Identification of the Heat Source for the Poisson equation. Acta Mathematica Scientia, 30, 1080-1087.
[12]	Barnett, A., Wu, B. and Veerapaneni, S. (2015) Spectrally Accurate Quadratures for Evaluation of Layer Potentials Close to the Boundary for the 2D Stokes and Laplace Equations. SIAM Journal on Scientific Computing, 37, 519-542. https://doi.org/10.1137/140990826
[13]	Heinrich, B. (1996) The Fourier-Finite-Element Method for Poisson equation in Axisymmetric Domains with Edges. SIAM Journal on Numerical Analysis, 33, 1885-1911. https://doi.org/10.1137/S0036142994266108
[14]	Nkemzi, B. and Jung, M. (2020) The Fourier-Finite-Element Method for Poisson equation in Three-Dimensional Axisymmetric Domains with Edges: Computing the Edge Flux Intensity Functions. Journal of Numerical Mathematics, 28, 75-98. https://doi.org/10.1515/jnma-2019-0002
[15]	Febrianto, E., Ortiz, M. and Cirak, F. (2021) Mollified Finite Element Approximants of Arbitrary Order and Smoothness. Computer Methods in Applied Mechanics and Engineering, 373, Article ID: 113513. https://doi.org/10.1016/j.cma.2020.113513
[16]	Et., A. (2021) A New Finite Element Method Using Matlab for Poisson Equation in Axisymmetric Domain. Turkish Journal of Computer and Mathematics Education, 12, 671-678. https://doi.org/10.17762/turcomat.v12i1S.1963
[17]	Jie, S. and Tao, T. (2006) Spectral and High-Order Methods with Applications. Science Press, Beijing.
[18]	Jie, S., Tao, T. and Wang, L.L. (2011) Spectral Methods: Algorithms, Analysis and Applications. Springer-Verlag, Berlin.
[19]	Boyd, J.P. (2001) Chebyshev and Fourier Spectral Methods. Dover, New York.
[20]	Trefethen, L.N. (2000) Spectral Methods in MATLAB. SIAM, Philadelphia. https://doi.org/10.1137/1.9780898719598
[21]	Kosloff, D. (2012) Solution of the Equations of Dynamic Elasticity by a Chebychev Spectral Method. Geophysics, 55, Article No. 734. https://doi.org/10.1190/1.1442885
[22]	Long, Y. and Zeng, B. (2021) Sign-Changing Solutions for Discrete Dirichlet Boundary Value Problem. Journal of Applied Mathematics and Physics, 5, 2228-2243. https://doi.org/10.4236/jamp.2017.511182
[23]	Mangoubi, J., Moukoko, D., Moukamba, F., et al. (2016) Existence and Uniqueness of Solution for Cahn-Hilliard Hyperbolic Phase-Field System with Dirichlet Boundary Condition and Regular Potentials. Applied Mathematics, 7, 1919-1926. https://doi.org/10.4236/am.2016.716157