- Email: [email protected]
A fast numerical algorithm for the inverse of a tridiagonal and pentadiagonal matrix
Available online at www.sciencedirect.com
Applied Mathematics and Computation 202 (2008) 441–445 www.elsevier.com/locate/amc
A fast numerical algorithm for the inverse of a tridiagonal and pentadiagonal matrix Ahmed Driss Aiat Hadj *, Mohamed Elouafi Preparatory Classes to High School Engineering, BP 3117, Souanni, Tangier, Morocco
Abstract In the current paper, we present a computationally efficient algorithm for obtaining the inverse of a tridiagonal matrix and a pentadiogonal matrix. Few conditions are required, and the algorithm is suited for implementation using computer algebra systems. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Determinants; Inverse matrix; Tridiagonal and pentadiagonal matrices
1. Introduction Tridiagonal and pentadiagonal matrices frequently arise from boundary value problems involving fourth order derivatives, and a fast computational formula for the inverse is required to test efficiently the existence of unique solutions of the PDEs. For such systems, the solution can be obtained in OðnÞ operations instead of Oðn3 Þ required by Gaussian elimination. The methods usually used such as the Thomas algorithm and the well based LU decomposition methods require some restrictive condition on the matrix to be inversed [2,6]. Recently explicit formula for the elements of the inverse of a general tridiagonal matrix is given by Mallik [4] and Kilic [5]. In this paper, we present a new algorithm with totally different approach for numerically inversing a tridiagonal and a pentadiagonal matrix. The conditions required are less few and the algorithm always works. Moreover, the algorithm is suited for implementation in a parallel machine. The rest of this paper is organized as follows: in the next section, we give recurrence formulas for the columns of the inverse of a tridiagonal matrix. In Section 3, a similar approach is applied to a pentadiagonal matrix.
*
Corresponding author. E-mail addresses: [email protected] (A.D.A. Hadj), med_elouafi@hotmail.com (M. Elouafi).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.02.026
442
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
2. Inverse of a tridiagonal matrix In this section, we give recurrence formulas for the columns of the inverse of a tridiagonal matrix. We consider the tridiagonal matrix T ¼ ðaij Þ16i;j6n 2 Mn;n ðKÞ over an arbitrary field K. Thus, aii ¼ ai ;
aiiþ1 ¼ ci ;
aii1 ¼ bi ;
and
aij ¼ 0
where fai g16i6n ; fbi g16i6n1 ; and fci g16i6n1 are i ¼ 1; 2; . . . ; n 1; cn ¼ 1. Assume that T is non-singular and denote
for ji jj > 1;
sequences
of numbers
such
that
ci 6¼ 0
for
all
T 1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ; where C j is the jth column of the inverse T 1 . From the relation T 1 T ¼ I n , where I n denotes the identity matrix of order n, we deduce the relations C n1 ¼ C j1
1
ðEn an C n Þ; cn1 1 ¼ ðEj aj C j bjþ1 C jþ1 Þ cj1
ð1Þ for j ¼ n 1; . . . ; 2;
t
where Ej ¼ ½ðdi;j Þ16i6n 2 Kn is the vector of order j of the canonical basis of Kn ðdi;j is the Kronecker symbolÞ. It follow from (1) that knowing the last column C n we can recursively determine all other columns C n1 ; C n2 ; . . . ; C 1 . We give straightforward recurrence formulas for computing C n . We consider the sequence of numbers ðAi Þ06i6n defined by A0 ¼ 1; ð2Þ a1 A0 þ c1 A1 ¼ 0; and biþ1 Ai1 þ aiþ1 Ai þ ciþ1 Aiþ1 ¼ 0
for 1 6 i 6 n 1;
where fai g16i6n ; fbi g16i6n1 and fci g16i6n1 are sequences j ¼ 1; 2; . . . ; n 1; cn ¼ 1. We can give a matrix form to this term recurrence (2) T A ¼ An En ;
of
numbers
such
that
ci 6¼ 0
for
all
ð3Þ
t
where A ¼ ½A0 ; A1 ; . . . ; An1 . Lemma 2.1. If An ¼ 0, the matrix T is singular. Proof. If An ¼ 0, then A is non-null and by (3) T A ¼ 0. We conclude that T is singular. Let’s define i
Di ¼ ð1Þ
i Y
h
! c k Ai ;
k¼1
for i ¼ 1; . . . ; n and D0 ¼ 1. Then fDi g verifies the three-term recurrence relation Di ¼ ai Di1 ci1 bi Di2 ;
i ¼ 2; 3; . . . ; n;
with initial conditions D0 ¼ 1 and D1 ¼ a1 . It is not difficult to show that Di is the ith principal minor of T, in particular, Dn ¼ detðT Þ. Theorem 2.1. Suppose that An 6¼ 0, then T is invertible and if T 1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ is the inverse of T then t A0 A1 An1 ; ;...; ; Cn ¼ An An An and C j can be calculated as in (1) for j ¼ n 1; . . . ; 1.
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
443
n Qn Proof. Since detðT Þ ¼ Dn ¼ ð1Þ k¼1 ck An 6¼ 0, the matrix T is invertible. From relation (3) we have TC n ¼ En , the proof is completed. h Remark 2.1. Suppose that we have ci ¼ 0 for some i 2 f1; 2; . . . ; n 1g. Let T ðeÞ be the matrix obtained from T by replacing ci such that ci ¼ 0 with e > 0. We have lim!0 T ðeÞ ¼ T . If T is invertible then for some g > 0; T ðeÞ is invertible for e < g too and we can calculate the inverse ðT ðeÞÞ1 with the algorithm above. Then, we can compute T 1 by tending e to 0 in ðT ðeÞÞ1 . 3. Inverse of a pentadiagonal matrix In this section, we give recurrence formulas for the columns of the inverse of a pentadiagonal matrix. We consider the pentadiagonal matrix 0
a1
B Ba B 2 B B B b3 P¼B B B0 B B. B. @. 0
b1 a2 a3 .. . ..
.
.. . b2 c 2 .. .. .. . . . .. .. .. . . . .. .. .. . . . c1
0
0
bn
an
0 .. .
1
C C C C C 0 C C 2 Mn;n ðKÞ; C cn2 C C C C bn1 A an
where fai g16i6n ; fbi g16i6n ; fci g16i6n ; fai g26i6n and fbi g36i6n are sequences of numbers such that ci 6¼ 0; cn ¼ cn1 ¼ 1 and bn ¼ 0. Assume that P is non-singular and denote P1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ; where C j is the jth column of the inverse P1 . From the relation P1 P ¼ I n , where I n denotes the identity matrix of order n. we deduce the relation 1 ðEn an C n bn1 C n1 Þ; cn2 1 C n3 ¼ ðEn1 bn2 C n2 an1 C n1 an C n Þ; cn3 1 ðEj bj1 C j1 aj C j ajþ1 C jþ1 bjþ2 C jþ2 Þ C j2 ¼ cj2
C n2 ¼
ð4Þ for 3 6 j 6 n 2;
t
where Ej ¼ ½ðdi;j Þ16i6n 2 Kn is the vector of order j of the canonical basis of Kn . It follows from (4) that knowing the last columns C n and C n1 we can recursively determine all other columns C n2 ; . . . ; C 1 . As in the case of a tridiagonal matrix, we will give recurrence formulas for computing C n and C n1 . Consider the sequence of numbers fAi giP0 and fBi giP0 characterized by a term recurrence relation A0 ¼ 0; A1 ¼ 1; a1 A0 þ b1 A1 þ c1 A2 ¼ 0; a2 A0 þ a2 A1 þ b2 A2 þ c2 A3 ¼ 0; bi Ai3 þ ai Ai2 þ ai Ai1 þ bi Ai þ ci Aiþ1 ¼ 0
ð5Þ for i P 3;
444
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
and B0 ¼ 1; B1 ¼ 0; ð6Þ
a1 B0 þ b1 B1 þ c1 B2 ¼ 0; a2 B0 þ a2 B1 þ b2 B2 þ c2 B3 ¼ 0; bi Bi3 þ ai Bi2 þ ai Bi1 þ bi Bi þ ci Biþ1 ¼ 0
for i P 3:
We can give a matrix form to this term recurrence (2) PA ¼ An En1 Anþ1 En ; PB ¼ Bn En1 Bnþ1 En ;
ð7Þ ð8Þ
where A ¼ ½A0 ; A1 ; . . . ; An1 t , B ¼ ½B0 ; B1 ; . . . ; Bn1 t . Let’s define for 0 6 i 6 n þ 1 Anþ1 Ai An Ai X i ¼ det and Y i ¼ det : Bnþ1 Bi Bn Bi
ð9Þ
We have immediately PX ¼ X n En1 ;
ð10Þ
PY ¼ Y nþ1 En ; t
t
where X ¼ ½X 0 ; X 1 ; . . . ; X n1 , Y ¼ ½Y 0 ; Y 1 ; . . . ; Y n1 . We remark that X n ¼ Y nþ1 . Lemma 3.1. If X n ¼ 0, the matrix P is singular. Proof. We have X 0 ¼ Anþ1 ; Y 0 ¼ An . If Anþ1 6¼ 0 then X is non-null and PX ¼ 0. If An 6¼ 0 then Y is non-null and PY ¼ 0. If Anþ1 ¼ An ¼ 0, then A is non-null and PA ¼ 0. The proof is completed.
h
Theorem 3.1. Suppose that X n 6¼ 0, then P is invertible and
t Y 0 Y 1 Y n1 ; ;...; ; Y nþ1 Y nþ1 Y nþ1 ð11Þ t X 0 X 1 X n1 C n1 ¼ ; ;...; : Xn Xn Xn Q n2 c Proof. We have detðPÞ ¼ X n 6¼ 0 (see [1]). So P is invertible. From the relations Eqs. (10) and (11) k k¼1 we have PC n ¼ En ; PC n1 ¼ En1 . The proof is completed. h Cn ¼
4. Conclusion One can apply the algorithms above for solving tridiagonal and pentadiagonal systems and show that they require only OðnÞ number operations for a matrix of size n. In Table 1 we give the running time of the algorithm for a pentadiagonal system in relation with the order n of the system. The algorithm was tested in a Pentium 4, 2.4 GH.
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
445
Table 1 Running time of the algorithm n
The running time
103 104 105
4 ms 6s 8 min
References [1] D. Aiat Hadj, M. Elouafi, On the characteristic polynomial, eigenvectors and determinant of a pentadiagonal matrix, Appl. Math. Comput. 198 (2) (2008) 634–642. [2] M. El-Mikkawy, A. Karawia, Inversion of general tridiagonal matrices, Appl. Math. Lett. 19 (2006) 712–720. [4] Ranjan K. Mallik, The inverse of a tridiagonal matrix, Linear Algebra Appl. 325 (1–3) (2001) 109–139. [5] Emrah Kilic, Explicit formula for the inverse of a tridiagonal matrix by backward continued fractions, Appl. Math. Comput. 197 (1) (2008) 345–357. [6] Janina Losiak, E. Neuman, Jolanta Nowak, The inversion of cyclic tridiagonal matrices, Zastos. Mat. 20 (1) (1988) 93–102.
Further reading [3] Samuel Daniel Conte, Carl W. De Boor, Elementary Numerical Analysis: An Algorithmic Approach, McGraw-Hill Higher Education, 1980.
Applied Mathematics and Computation 202 (2008) 441–445 www.elsevier.com/locate/amc
A fast numerical algorithm for the inverse of a tridiagonal and pentadiagonal matrix Ahmed Driss Aiat Hadj *, Mohamed Elouafi Preparatory Classes to High School Engineering, BP 3117, Souanni, Tangier, Morocco
Abstract In the current paper, we present a computationally efficient algorithm for obtaining the inverse of a tridiagonal matrix and a pentadiogonal matrix. Few conditions are required, and the algorithm is suited for implementation using computer algebra systems. Ó 2008 Elsevier Inc. All rights reserved. Keywords: Determinants; Inverse matrix; Tridiagonal and pentadiagonal matrices
1. Introduction Tridiagonal and pentadiagonal matrices frequently arise from boundary value problems involving fourth order derivatives, and a fast computational formula for the inverse is required to test efficiently the existence of unique solutions of the PDEs. For such systems, the solution can be obtained in OðnÞ operations instead of Oðn3 Þ required by Gaussian elimination. The methods usually used such as the Thomas algorithm and the well based LU decomposition methods require some restrictive condition on the matrix to be inversed [2,6]. Recently explicit formula for the elements of the inverse of a general tridiagonal matrix is given by Mallik [4] and Kilic [5]. In this paper, we present a new algorithm with totally different approach for numerically inversing a tridiagonal and a pentadiagonal matrix. The conditions required are less few and the algorithm always works. Moreover, the algorithm is suited for implementation in a parallel machine. The rest of this paper is organized as follows: in the next section, we give recurrence formulas for the columns of the inverse of a tridiagonal matrix. In Section 3, a similar approach is applied to a pentadiagonal matrix.
*
Corresponding author. E-mail addresses: [email protected] (A.D.A. Hadj), med_elouafi@hotmail.com (M. Elouafi).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.02.026
442
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
2. Inverse of a tridiagonal matrix In this section, we give recurrence formulas for the columns of the inverse of a tridiagonal matrix. We consider the tridiagonal matrix T ¼ ðaij Þ16i;j6n 2 Mn;n ðKÞ over an arbitrary field K. Thus, aii ¼ ai ;
aiiþ1 ¼ ci ;
aii1 ¼ bi ;
and
aij ¼ 0
where fai g16i6n ; fbi g16i6n1 ; and fci g16i6n1 are i ¼ 1; 2; . . . ; n 1; cn ¼ 1. Assume that T is non-singular and denote
for ji jj > 1;
sequences
of numbers
such
that
ci 6¼ 0
for
all
T 1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ; where C j is the jth column of the inverse T 1 . From the relation T 1 T ¼ I n , where I n denotes the identity matrix of order n, we deduce the relations C n1 ¼ C j1
1
ðEn an C n Þ; cn1 1 ¼ ðEj aj C j bjþ1 C jþ1 Þ cj1
ð1Þ for j ¼ n 1; . . . ; 2;
t
where Ej ¼ ½ðdi;j Þ16i6n 2 Kn is the vector of order j of the canonical basis of Kn ðdi;j is the Kronecker symbolÞ. It follow from (1) that knowing the last column C n we can recursively determine all other columns C n1 ; C n2 ; . . . ; C 1 . We give straightforward recurrence formulas for computing C n . We consider the sequence of numbers ðAi Þ06i6n defined by A0 ¼ 1; ð2Þ a1 A0 þ c1 A1 ¼ 0; and biþ1 Ai1 þ aiþ1 Ai þ ciþ1 Aiþ1 ¼ 0
for 1 6 i 6 n 1;
where fai g16i6n ; fbi g16i6n1 and fci g16i6n1 are sequences j ¼ 1; 2; . . . ; n 1; cn ¼ 1. We can give a matrix form to this term recurrence (2) T A ¼ An En ;
of
numbers
such
that
ci 6¼ 0
for
all
ð3Þ
t
where A ¼ ½A0 ; A1 ; . . . ; An1 . Lemma 2.1. If An ¼ 0, the matrix T is singular. Proof. If An ¼ 0, then A is non-null and by (3) T A ¼ 0. We conclude that T is singular. Let’s define i
Di ¼ ð1Þ
i Y
h
! c k Ai ;
k¼1
for i ¼ 1; . . . ; n and D0 ¼ 1. Then fDi g verifies the three-term recurrence relation Di ¼ ai Di1 ci1 bi Di2 ;
i ¼ 2; 3; . . . ; n;
with initial conditions D0 ¼ 1 and D1 ¼ a1 . It is not difficult to show that Di is the ith principal minor of T, in particular, Dn ¼ detðT Þ. Theorem 2.1. Suppose that An 6¼ 0, then T is invertible and if T 1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ is the inverse of T then t A0 A1 An1 ; ;...; ; Cn ¼ An An An and C j can be calculated as in (1) for j ¼ n 1; . . . ; 1.
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
443
n Qn Proof. Since detðT Þ ¼ Dn ¼ ð1Þ k¼1 ck An 6¼ 0, the matrix T is invertible. From relation (3) we have TC n ¼ En , the proof is completed. h Remark 2.1. Suppose that we have ci ¼ 0 for some i 2 f1; 2; . . . ; n 1g. Let T ðeÞ be the matrix obtained from T by replacing ci such that ci ¼ 0 with e > 0. We have lim!0 T ðeÞ ¼ T . If T is invertible then for some g > 0; T ðeÞ is invertible for e < g too and we can calculate the inverse ðT ðeÞÞ1 with the algorithm above. Then, we can compute T 1 by tending e to 0 in ðT ðeÞÞ1 . 3. Inverse of a pentadiagonal matrix In this section, we give recurrence formulas for the columns of the inverse of a pentadiagonal matrix. We consider the pentadiagonal matrix 0
a1
B Ba B 2 B B B b3 P¼B B B0 B B. B. @. 0
b1 a2 a3 .. . ..
.
.. . b2 c 2 .. .. .. . . . .. .. .. . . . .. .. .. . . . c1
0
0
bn
an
0 .. .
1
C C C C C 0 C C 2 Mn;n ðKÞ; C cn2 C C C C bn1 A an
where fai g16i6n ; fbi g16i6n ; fci g16i6n ; fai g26i6n and fbi g36i6n are sequences of numbers such that ci 6¼ 0; cn ¼ cn1 ¼ 1 and bn ¼ 0. Assume that P is non-singular and denote P1 ¼ ðC 1 ; C 2 ; . . . ; C n Þ; where C j is the jth column of the inverse P1 . From the relation P1 P ¼ I n , where I n denotes the identity matrix of order n. we deduce the relation 1 ðEn an C n bn1 C n1 Þ; cn2 1 C n3 ¼ ðEn1 bn2 C n2 an1 C n1 an C n Þ; cn3 1 ðEj bj1 C j1 aj C j ajþ1 C jþ1 bjþ2 C jþ2 Þ C j2 ¼ cj2
C n2 ¼
ð4Þ for 3 6 j 6 n 2;
t
where Ej ¼ ½ðdi;j Þ16i6n 2 Kn is the vector of order j of the canonical basis of Kn . It follows from (4) that knowing the last columns C n and C n1 we can recursively determine all other columns C n2 ; . . . ; C 1 . As in the case of a tridiagonal matrix, we will give recurrence formulas for computing C n and C n1 . Consider the sequence of numbers fAi giP0 and fBi giP0 characterized by a term recurrence relation A0 ¼ 0; A1 ¼ 1; a1 A0 þ b1 A1 þ c1 A2 ¼ 0; a2 A0 þ a2 A1 þ b2 A2 þ c2 A3 ¼ 0; bi Ai3 þ ai Ai2 þ ai Ai1 þ bi Ai þ ci Aiþ1 ¼ 0
ð5Þ for i P 3;
444
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
and B0 ¼ 1; B1 ¼ 0; ð6Þ
a1 B0 þ b1 B1 þ c1 B2 ¼ 0; a2 B0 þ a2 B1 þ b2 B2 þ c2 B3 ¼ 0; bi Bi3 þ ai Bi2 þ ai Bi1 þ bi Bi þ ci Biþ1 ¼ 0
for i P 3:
We can give a matrix form to this term recurrence (2) PA ¼ An En1 Anþ1 En ; PB ¼ Bn En1 Bnþ1 En ;
ð7Þ ð8Þ
where A ¼ ½A0 ; A1 ; . . . ; An1 t , B ¼ ½B0 ; B1 ; . . . ; Bn1 t . Let’s define for 0 6 i 6 n þ 1 Anþ1 Ai An Ai X i ¼ det and Y i ¼ det : Bnþ1 Bi Bn Bi
ð9Þ
We have immediately PX ¼ X n En1 ;
ð10Þ
PY ¼ Y nþ1 En ; t
t
where X ¼ ½X 0 ; X 1 ; . . . ; X n1 , Y ¼ ½Y 0 ; Y 1 ; . . . ; Y n1 . We remark that X n ¼ Y nþ1 . Lemma 3.1. If X n ¼ 0, the matrix P is singular. Proof. We have X 0 ¼ Anþ1 ; Y 0 ¼ An . If Anþ1 6¼ 0 then X is non-null and PX ¼ 0. If An 6¼ 0 then Y is non-null and PY ¼ 0. If Anþ1 ¼ An ¼ 0, then A is non-null and PA ¼ 0. The proof is completed.
h
Theorem 3.1. Suppose that X n 6¼ 0, then P is invertible and
t Y 0 Y 1 Y n1 ; ;...; ; Y nþ1 Y nþ1 Y nþ1 ð11Þ t X 0 X 1 X n1 C n1 ¼ ; ;...; : Xn Xn Xn Q n2 c Proof. We have detðPÞ ¼ X n 6¼ 0 (see [1]). So P is invertible. From the relations Eqs. (10) and (11) k k¼1 we have PC n ¼ En ; PC n1 ¼ En1 . The proof is completed. h Cn ¼
4. Conclusion One can apply the algorithms above for solving tridiagonal and pentadiagonal systems and show that they require only OðnÞ number operations for a matrix of size n. In Table 1 we give the running time of the algorithm for a pentadiagonal system in relation with the order n of the system. The algorithm was tested in a Pentium 4, 2.4 GH.
A.D.A. Hadj, M. Elouafi / Applied Mathematics and Computation 202 (2008) 441–445
445
Table 1 Running time of the algorithm n
The running time
103 104 105
4 ms 6s 8 min
References [1] D. Aiat Hadj, M. Elouafi, On the characteristic polynomial, eigenvectors and determinant of a pentadiagonal matrix, Appl. Math. Comput. 198 (2) (2008) 634–642. [2] M. El-Mikkawy, A. Karawia, Inversion of general tridiagonal matrices, Appl. Math. Lett. 19 (2006) 712–720. [4] Ranjan K. Mallik, The inverse of a tridiagonal matrix, Linear Algebra Appl. 325 (1–3) (2001) 109–139. [5] Emrah Kilic, Explicit formula for the inverse of a tridiagonal matrix by backward continued fractions, Appl. Math. Comput. 197 (1) (2008) 345–357. [6] Janina Losiak, E. Neuman, Jolanta Nowak, The inversion of cyclic tridiagonal matrices, Zastos. Mat. 20 (1) (1988) 93–102.
Further reading [3] Samuel Daniel Conte, Carl W. De Boor, Elementary Numerical Analysis: An Algorithmic Approach, McGraw-Hill Higher Education, 1980.