Exact Inversion of Pentadiagonal Matrix for Semi-Analytic Solution of 2D Poisson Equation

Serigne Bira Gueye

Department of Physics, Faculty of Sciences and Techniques, Cheikh Anta Diop University, Dakar-Fann, Senegal.

DOI: 10.4236/jmp.2022.1312094                

Abstract

This work essentially consists in inverting in an exact, explicit, and original way the pentadiagonal Toeplitz matrix or tridiagonal block matrix resulting from the discretization of the two-dimensional Laplace operator. This method is an algorithm facilitating the resolution of a large number of problems governed by PDEs involving the Laplacian in two dimensions. It guarantees high precision and high efficiency in solving various differential equations.

Keywords

Block Matrix, Pentadiagonal Matrix, 2D Laplacian, Toeplitz Matrix, Inversion

1. Introduction

Partial differential equations govern most phenomena and systems in physics and engineering. Laplace’s operator intervenes there very often. Once discretized, these differential equations are reduced to a matrix form, where the matrix corresponding to the Laplacian of the quantity to be studied takes the form of a pentadiagonal matrix. This matrix can be considered as a tridiagonal block matrix for the two-dimensional case. The inversion of this matrix makes it possible to obtain the solution of several categories of problems. There are different numerical and other approximation methods for inverting this matrix. The originality of this present work is to propose an exact and direct formula for the inverse of this famous pentadiagonal matrix.

Concretely, it will be a matter of determining the exact formula of the following symmetric block matrix M, which is a Toeplitz matrix:

$\underset{:=M}{\underset{︸}{\left(\begin{array}{cccccc}D& A& O_{N_x}& \cdots & \cdots & O_{N_x}\\ A& D& A& \ddots & \cdots & O_{N_x}\\ O_{N_x}& A& D& A& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & O_{N_x}\\ O_{N_x}& \cdots & \ddots & \ddots & \ddots & A\\ O_{N_x}& \cdots & \cdots & O_{N_x}& A& D\end{array}\right)}}$, where $\underset{:=D}{\underset{︸}{\left(\begin{array}{cccccc}d& a& 0& \cdots & \cdots & 0\\ a& d& a& \ddots & \cdots & 0\\ 0& a& d& a& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & 0\\ 0& \cdots & \ddots & \ddots & \ddots & a\\ 0& \cdots & \cdots & 0& a& d\end{array}\right)}}$, and $\underset{:=A}{\underset{︸}{\left(\begin{array}{cccccc}\alpha & 0& 0& \cdots & \cdots & 0\\ 0& \alpha & 0& \ddots & \cdots & 0\\ 0& 0& \alpha & 0& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & 0\\ 0& \cdots & \ddots & \ddots & \ddots & 0\\ 0& \cdots & \cdots & 0& 0& \alpha \end{array}\right)}}$.

D is a tridiagonal, symmetric, and persymmetric Toeplitz matrix of dimension ( $N_x\times N_x$ ). A is a diagonal Toeplitz matrix. We denote by $O_{N_x}$ the null matrix of the same dimension as D. The parameters d, a, and $\alpha$ are reals.

The matrix M, on the other hand, is pentadiagonal. It can be described as a symmetric Toeplitz block matrix. Its dimension is $\left(N_x\times N_y\right)\times \left(N_x\times N_y\right)$.

The inversion method we propose for the block matrix M follows a certain analogy with the approach used in [1] [2] [3] [4] [5] to invert the tridiagonal matrix D. The inverse of M has the following structure:

$\underset{:=M^{-1}}{\underset{︸}{\left(\begin{array}{ccccccc}\boxed{B^{11}}& \boxed{B^{12}}& \boxed{B^{13}}& \cdots & \boxed{B^{1l}}& \cdots & \boxed{B^{1N_y}}\\ \boxed{B^{21}}& \boxed{B^{22}}& \boxed{B^{23}}& \cdots & \boxed{B^{2l}}& \cdots & \boxed{B^{2N_y}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \boxed{B^{k1}}& \boxed{B^{k2}}& \boxed{B^{k3}}& \cdots & \boxed{B^{kl}}& \cdots & \boxed{B^{k{N}_y}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ \boxed{B^{N_{y}1}}& \boxed{B^{N_{y}2}}& \boxed{B^{N_{y}3}}& \cdots & \boxed{B^{N_{y}l}}& \cdots & \boxed{B^{N_y{N}_y}}\end{array}\right)}}$

$M^{-1}$ is a symmetric block matrix and the matrices $B^{kl}$ which compose it are of the same dimension as D.

$M^{-1}=\left(B^{kl}\right),B^{kl}=\left(b_{ij}^{kl}\right),k,l=1,\cdots ,N_y;i,j=1,\cdots ,N_x$ (1)

$M^{-1}$ is symmetric and persymmetric. It is composed of symmetric and persymmetric matrices (blocks):

$b_{ij}^{kl}=b_{ij}^{lk}=b_{ij}^{N_y-l+1,N_y-k+1}=b_{ji}^{kl}=b_{N_x-j+1,N_x-i+1}^{kl}$ (2)

Inverting M amounts to determining the blocks $B^{kl}$. To obtain them, we will first recall the formulas for the inversion of D and those for its eigenvalues and its eigenvectors. We will then exploit these parameters of D using them as a key to determine $M^{-1}$. Finally, we will formulate the exact expression for the inverse of M.

2. Inverse of Matrix D

Concerning the matrix D, its inverse is obtained with its cofactors. The inversion of a tridiagonal Toeplitz matrix has been covered in detail in [2] [3] [4] [5].

For the reals d and a, we define the real number ${\Delta }_{da}$, and the complex number ${\theta }_{da}$, as follows:

$\Delta \left(d\mathrm{,}a\right)={\Delta }_{da}=d^2-4a^2\text{and}\theta \left(d\mathrm{,}a\right)={\theta }_{da}=\text{atanh}\left(\frac{\sqrt{{\Delta }_{da}}}{d}\right)\mathrm{,}d\ne 0.$ (3)

Then, the following polynomial $P_n\left(d\mathrm{,}a\right)$ makes it possible to calculate, in an explicit and direct way, the determinants of the sub-matrices $D_{\left(n\times n\right)}$ of D:

$P_n\left(d,a\right)=\{\begin{array}{l}1,n=0\hfill \\ d,n=1\hfill \\ d^n,a=0\hfill \\ a^n\cos \left(n\frac{\pi }{2}\right),d=0\hfill \\ {\left(\mathrm{sgn}\left(d\right)a\right)}^n\left(n+1\right),\left|d\right|=2\left|a\right|\hfill \\ \frac{1}{\sqrt{{\Delta }_{da}}}\left\{{\left[\frac{d+\sqrt{{\Delta }_{da}}}{2}\right]}^{n+1}-{\left[\frac{d-\sqrt{{\Delta }_{da}}}{2}\right]}^{n+1}\right\}={\left(\mathrm{sgn}\left(d\right)a\right)}^n\frac{\sinh \left(\left(n+1\right){\theta }_{da}\right)}{\sinh \left({\theta }_{da}\right)},\text{otherwise}\hfill \end{array}$ (4)

or in another form (for $d\cdot a\ne 0$ ):

$P_n\left(d,a\right)={\displaystyle \underset{i=0}{\overset{\left[n/2\right]}{\sum }}}{\left(-1\right)}^i{C}_i^{n-i}{d}^{n-2i}{a}^{2i}.$ (5)

The determinants of these sub-matrices $D_{\left(n\mathrm{,}n\right)}$ of D are given, for fixed a and b, by $P_n\left(d,a\right)=U_n,n=0,1,\cdots ,N_x$. We denote by B the inverse of D. This matrix B exists when d is different from one of the following $d_i$ values:

$d_i=2a\cos \left(\frac{i\pi }{N_x+1}\right),i=1,2,\cdots ,N_x\cdot$ (6)

Apart from these values $d_i$, the matrix D is regular and its inverse B is symmetric, persymmetric, and is given by the following formula [2] [3] [4] [5]:

$b_{ij}={\left(-1\right)}^{i+j}\cdot a^{j-i}\frac{U_{i-1}\cdot U_{N_x-j}}{U_{N_x}}\mathrm{,}1\le i\le j\le N_x.$ (7)

It yields

$b_{ij}={\left(-1\right)}^{i+j}\cdot a^{j-i}\frac{P_{i-1}\left(d,a\right)\cdot P_{N_x-j}\left(d,a\right)}{P_{N_x}\left(d,a\right)},1\le i\le j\le N_x.$ (8)

Which gives, in case ${\theta }_{da}$ exists and is non-zero:

$b_{ij}={\left(-1\right)}^{i+j}{\left(\mathrm{sgn}\left(d\right)a\right)}^{i-j-1}\frac{\sinh \left(i{\theta }_{da}\right)\sinh \left(\left(N_x-j+1\right){\theta }_{da}\right)}{\sinh \left({\theta }_{da}\right)\sinh \left(\left(N_x+1\right){\theta }_{da}\right)},1\le i\le j\le N_x.$ (9)

This result will be used to determine the inverse of the block matrix M.

3. Block Matrix Inversion

To find the inverse of the block matrix M, one can reason by analogy with the method of inversion of the tridiagonal matrix D. One of the keys to obtain this inverse matrix is to use the polynomial $P_n\left(d\mathrm{,}a\right)$. This time, we are not going to calculate it for real numbers d and a, but rather for matrices D and A. It is then necessary to determine the matrices $P_n\left(D\mathrm{,}A\right)$. The matrix D being Hermitian, its spectral decomposition is: $D=X\Lambda X^t$, where $\Lambda$ is the diagonal matrix containing on its ith row the eigenvalue ${\lambda }_i$ : ${\Lambda }_{ii}={\lambda }_i$.

${\lambda }_i=d+2a\cos \left(\frac{i}{N_x+1}\pi \right),i=1,2,\cdots ,N_x.$ (10)

X is an unitary matrix containing on its ith line the normalized eigenvector ${\stackrel{\to }{x}}_i$ associated with the eigenvalue ${\lambda }_i$. The jth component of vector ${\stackrel{\to }{x}}_i$ is:

$x_{ij}=x_{ji}=\sqrt{\frac{2}{N_x+1}}\cdot \sin \left(\frac{i\cdot j}{N_x+1}\pi \right),i,j=1,2,\cdots ,N_x.$ (11)

Applying the rule $D^n=X{\Lambda }^n{X}^t$, we can calculate $P_n\left(D\mathrm{,}A\right)$ :

$P_n\left(D,A\right)={\displaystyle \underset{i=0}{\overset{\left[n/2\right]}{\sum }}}{\left(-1\right)}^i{C}_i^{n-i}{D}^{n-2i}{A}^{2i}={\displaystyle \underset{i=0}{\overset{\left[n/2\right]}{\sum }}}{\left(-1\right)}^i{C}_i^{n-i}X{\Lambda }^{n-2i}{X}^t{A}^{2i},$ (12)

or

$P_n\left(D,A\right)=X\left({\displaystyle \underset{i=0}{\overset{\left[n/2\right]}{\sum }}}{\left(-1\right)}^i{C}_i^{n-i}{\Lambda }^{n-2i}{A}^{2i}\right)X^t.$ (13)

We can notice the spectral decomposition of the matrix $P_N\left(D\mathrm{,}A\right)$ :

$P_n\left(D\mathrm{,}A\right)=X{P}_n\left(\Lambda \mathrm{,}A\right)X^t\mathrm{,}$ (14)

By analogy with the inversion of the tridiagonal matrix D, the inverse of the block matrix M is given by the matrices $B^{kl}$, with:

$B^{kl}={\left(-1\right)}^{k+l}{\alpha }^{l-k}{P}_{k-1}\left(D,A\right)P_{N_y}{\left(D,A\right)}^{-1}{P}_{N_y-l}\left(D,A\right),1\le k\le l\le N_y$ (15)

We can, therefore, notice that the inverse of M exists when the determinant of $P_{N_y}\left(D,A\right)$ is non-zero.

From Equations (14) and (15), the spectral decomposition of the matrix $B^{kl}$ appears:

$B^{kl}=X\underset{:={\Lambda }^{kl}}{\underset{︸}{{\left(-1\right)}^{k+l}{\alpha }^{l-k}{P}_{k-1}\left(\Lambda ,A\right)P_{N_y}{\left(\Lambda ,A\right)}^{-1}{P}_{N_y-l}\left(\Lambda ,A\right)}}{X}^t,1\le k\le l\le N_y,$ (16)

that is:

$B^{kl}=X{\Lambda }^{kl}{X}^t\mathrm{,}1\le k\le l\le N_y$ (17)

${\Lambda }^{kl}$ is the matrix containing on its ith line the ith eigenvalue of $B^{kl}$. These eigenvalues of $B^{kl}$ are ${\Lambda }_{ii}^{kl}={\lambda }_i^{kl}$, with:

${\lambda }_i^{kl}=\cdot \underset{:={\lambda }_i^{kl}}{\underset{︸}{{\left(-1\right)}^{k+l}{\alpha }^{l-k}\frac{P_{k-1}\left({\lambda }_i,\alpha \right)\cdot P_{N_y-l}\left({\lambda }_i,\alpha \right)}{P_{N_y}\left({\lambda }_i,\alpha \right)}}},1\le k\le l\le N_y,$ (18)

where ${\theta }_{{\lambda }_i\alpha }=\theta \left({\lambda }_i,\alpha \right)$ :

${\lambda }_i^{kl}={\left(-1\right)}^{k+l}{\left(\mathrm{sgn}\left({\lambda }_i\right)\alpha \right)}^{k-l-1}\frac{\sinh \left(k{\theta }_{{\lambda }_i\alpha }\right)\sinh \left(\left(N_y-l+1\right){\theta }_{{\lambda }_i\alpha }\right)}{\sinh \left({\theta }_{{\lambda }_i\alpha }\right)\sinh \left(\left(N_y+1\right){\theta }_{{\lambda }_i\alpha }\right)},1\le k\le l\le N_y.$ (19)

Thus, we obtain all the elements of the matrices $B^{kl}$ composing the block matrix $M^{-1}$ :

$b_{ij}^{kl}={\displaystyle \underset{\nu =1}{\overset{N_x}{\sum }}}{\lambda }_{\nu }^{kl}{x}_{i\nu }{x}_{\nu j},i,j=1,\cdots ,N_x;1\le k\le l\le N_y$ (20)

Equation (20) completely determines the matrix $M^{-1}$. Note that these $B^{kl}$ blocks are symmetric and persymmetric matrices. These different symmetries allow a calculation of the inverse of the pentadiagonal matrix M, in a fast, efficient, precise, and economical way in terms of occupying RAM space. This method allows a very fine implementation of the discretization of the Laplace operator, for 2D problems, even with an inhomogeneous mesh.

4. Conclusion

This work makes it possible to determine the exact formula of the inverse of the symmetric pentadiagonal Toeplitz matrix. This result allows a direct inversion of this matrix of blocks, in an efficient way, precise and economical way, in terms of occupying RAM space. This proposed method will allow simple and accurate resolution of a large category of physics and engineering problems, governed by differential equations involving the two-dimensional Lapalce operator.

Conflicts of Interest

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

References

[1]	Usmani, R.A. (1994) Linear Algebra and its Applications, 212-213, 413-414. https://doi.org/10.1016/0024-3795(94)90414-6
[2]	Gueye, S. (2014) Journal of Electromagnetic Analysis and Applications, 6, 303-308
[3]	Hu, G.Y. and O’Connell, R.F. (1996) Journal of Physics A: Mathematical and General, 29, 1511-1513. https://doi.org/10.1088/0305-4470/29/7/020
[4]	Da Fonseca, C.M. and Petronilho, J. (2001) Linear Algebra and Its Applications, 325, 7-21. https://doi.org/10.1016/S0024-3795(00)00289-5
[5]	Gueye, S.B., Mbow, C. and Diagana, M.F. (2021) Journal of Electromagnetic Analysis and Applications, 13, 135-143.