- Email: [email protected]
A new algorithm for solving nearly penta-diagonal Toeplitz linear systems
Computers and Mathematics with Applications 63 (2012) 1238–1243
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A new algorithm for solving nearly penta-diagonal Toeplitz linear systems✩ Jiteng Jia a,∗ , Qiongxiang Kong b , Tomohiro Sogabe c a
Department of Mathematical Sciences, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
b
Department of Building Environment and Services Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
c
Graduate School of Information Science and Technology, Aichi Prefectural University, 1522-3 Kumabari, Nagakute, Aichi 480-1198, Japan
article
abstract
info
Article history: Received 20 September 2011 Received in revised form 13 December 2011 Accepted 16 December 2011
A new numerical algorithm for solving nearly penta-diagonal Toeplitz linear systems is presented. The algorithm is suited for implementation using Computer Algebra Systems (CASs) such as MATLAB, MACSYMA and MAPLE. Numerical examples are given in order to illustrate the efficiency of our algorithm. © 2011 Elsevier Ltd. All rights reserved.
Keywords: Nearly penta-diagonal matrices Sherman–Morrison–Woodbury formula Toeplitz matrices Linear systems
1. Introduction Linear systems of equations with nearly penta-diagonal Toeplitz coefficient matrices appear in many applications such as differential equations, parallel computing, matrix algebra, data interpolation, boundary value problems, etc [1–12]. In this paper, we consider an n × n nearly penta-diagonal Toeplitz system of linear equations given by Ax = f,
(1.1)
where
a 0 a1 a2 0 A= . .. 0
an − 2 an − 1
a− 1 a0 a1
a− 2 a− 1 a0
0 a− 2 a− 1
···
0
0 a− 2
···
.
.
.
.
.
..
..
. ···
..
..
.
0
..
..
.
0
···
a2 0
an − 2
0
···
a2 − n 0
.
..
.
··· .. . .. .
a1 a2 0
a0 a1 a2
a− 1 a0 a1
..
..
0
..
a1−n a2−n 0
.. . , 0 a− 2
(1.2)
a− 1 a0
✩ This work was supported by the Natural Science Foundation of China (NSFC) under grant 11071192 and the International Science and Technology Cooperation Program of China under grant 2010DFA14700. ∗ Corresponding author. E-mail addresses: [email protected] (J. Jia), [email protected] (T. Sogabe).
0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2011.12.044
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1239
x = (x1 , x2 , . . . , xn )T , f = ( f1 , f2 , . . . , fn )T . Xiao-Guang Lv and Jiang Le [13] presented a fast computational algorithm for solving nearly penta-diagonal linear systems based on the use of any penta-diagonal linear solver. By using LU decomposition, Neossi Nguetchue and Abelman [14] gave a numerical algorithm to obtain the solution of nearly pentadiagonal linear systems. Recently, Xiangjian Xu [15] has proposed another efficient algorithm to solve symmetric Toeplitz penta-diagonal linear systems based on the Sherman–Morrison–Woodbury formula. In this study, we derive a numerical algorithm for solving the nearly penta-diagonal Toeplitz linear systems (1.1) and show that the total number of operations is less than those of the two algorithms in [13,14]. The rest of this paper is organized as follows: In Section 2, we describe a more efficient algorithm for solving nearly pentadiagonal Toeplitz linear systems. In Section 3, numerical examples are provided to show the performance and efficiency of our algorithm. Finally, we make some concluding remarks in Section 4. 2. Main results We present two algorithms for solving penta-diagonal linear systems and nearly penta-diagonal Toeplitz linear systems respectively. 2.1. An algorithm for solving penta-diagonal linear systems Consider the penta-diagonal system of linear equations A′ x = f given by
b c a2
a−1 b a1
a− 2 a−1 a0
a− 2 a− 1
a− 2
.
.
.
.
..
..
..
a2
a1 a2
..
..
a0 a1 a2
a−2 e
.
a− 1 d a1
x1 x2 x3
.. . .. .
x n −1 xn
d
f1 f2 f3 .. = . .. .
fn−1 fn
.
(2.1)
Without loss of generality, suppose that the penta-diagonal matrix A′ is nonsingular and that the element a−2 ̸= 0. We now introduce some notations for later use.
V =
b a− 1
U =
c b
··· ···
0 0
a2 a1 0 0
0 a2 a2 0
ˆf = ( f1 , f2 , . . . , fn−2 )T , xˆ = (x1 , x2 )T ,
T
···
0
0
···
a1 a2
d a1
e d
0 0
∈ R(n−2)×2 ,
∈ R2×(n−2) ,
˜f = ( fn−1 , fn )T ,
x˜ = (x3 , x4 , . . . , xn )T ,
ˆf f= ˜f , xˆ x= ˜ . x
According to the above notation, linear system (2.1) can be written in the form
V O
L U
ˆf xˆ = ˜ , x˜ f
(2.2)
where L is a lower triangular Toeplitz matrix of size (n − 2) × (n − 2). Thus (2.2) is equivalent to
V xˆ + Lx˜ = ˆf, U x˜ = ˜f.
(2.3)
Since det L = (a−2 )n−2 ̸= 0, L is invertible. Assume that UL−1 V is nonsingular, it is easy to deduce that xˆ = (UL−1 V )−1 (UL−1 ˆf − ˜f), x˜ = L−1 (ˆf − V xˆ ).
Proposition 2.1. The inverse matrix of a lower triangular Toeplitz matrix is a lower triangular Toeplitz matrix. Proof. See [16].
(2.4)
1240
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
Proposition 2.2. Let (t1 , t2 , . . . , tn )T be the first column of an n × n lower triangular Toeplitz matrix T , (ξ1 , ξ2 , . . . , ξn )T be the first column of T −1 , where T −1 is the inverse matrix of T . Then
ξ1 ξ2 ti , . . . , t2 , t1 .. = 0, . ξi
t1 ξ1 = 1,
i = 2, . . . , n.
Proof. Since (ξ1 , ξ2 , . . . , ξn )T is the first column of T −1 , we have t1 t2
t1
. . .
..
..
. ···
tn
. ···
1 ξ1 ξ2 0 . = . . . . . . ξn 0 t1
By simple algebraic manipulations, the conclusion is true. Obviously, L
−1
is a lower triangular Toeplitz matrix. Let the first column of L−1 be (η1 , η2 , . . . , ηn−2 )T . Then we have
xˆ =
br1 + cr2 + a2 r3 bs1 + cs2 + a2 s3
n −2 ri fi − fn−1 −1 a−1 r1 + br2 + a1 r3 + a2 r4 i =1 , a−1 s1 + bs2 + a1 s3 + a2 s4 n −2 si fi − fn
(2.5)
i =1
where ri = a2 ηn−i−4 + a1 ηn−i−3 + dηn−i−2 + eηn−i−1 for i = 0, 1, . . . , n − 2, sj = rj−1 − eηn−j for j = 1, 2, . . . , n − 2, and ηk = 0 for k ≤ 0 or k ≥ n − 2. After determining xˆ , set f1′ = f1 − bx1 − a−1 x2 , f2′ = f2 − cx1 − bx2 , f3′ = f3 − a2 x1 − a1 x2 , f4′ = f4 − a2 x2 and fi′ = fi for i = 5, 6, . . . , n − 2, then by using backward substitution, we can obtain x˜ as follows: x˜ = (x3 , . . . , xn−1 , xn )T , where x3 =
f1′ a−2
, xi =
1 a−2
′ ( fi− 2 −
(2.6)
i−3
j =1
ai−j−4 xj+2 ) for i = 4, 5, . . . , n, and ak = 0 for k > 2.
Remark 2.1. If ˆf is a standard unit vector, we can obtain x˜ from the second equation of (2.4) directly, and the computational costs are much less than those of (2.6). Thus, it is important to compute the first column elements ηi (i = 1, 2, . . . , n − 2) of L−1 . From Proposition 2.2, we can obtain the general recurrence formula
ηi = − where η1 =
i−1 1
a− 2 j = 1
1 a−2
ai−j−2 ηj ,
i = 2, 3, . . . , n − 2,
(2.7)
and ak = 0 for 2 < k < n − 2.
In the following, we state the algorithm for solving the linear system (2.1): Algorithm 2.1 Step 1. Input a−2 , a−1 , a0 , a1 , a2 , b, c , d, e, n, f. Step 2. Set η1 = a 1 , ak = 0 for 2 < k < n − 2, then compute η2 , η3 , . . . , ηn−2 by using (2.7). −2
Step 3. Step 4.
Step 5.
Step 6. Step 7. Step 8.
Set ηk = 0 for k ≤ 0 or k > n − 2. For i = 0, 1, . . . , n − 2, compute ri = a2 ηn−i−4 + a1 ηn−i−3 + dηn−i−2 + eηn−i−1 . End. For j = 1, 2, . . . , n − 2, compute sj = rj−1 − eηn−j . End. Compute xˆ = (x1 , x2 )T by using (2.5). Set f1′ = f1 − bx1 − a−1 x2 , f2′ = f2 − cx1 − bx2 , f3′ = f3 − a2 x1 − a1 x2 , f4′ = f4 − a2 x2 and fi′ = fi for i = 5, 6, . . . , n − 2, then compute x˜ = (x3 , x4 , . . . , xn−2 )T from the second equation of (2.4) or by using (2.6). xˆ Output the solution of system: x = . x˜
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1241
2.2. Solving nearly penta-diagonal Toeplitz linear systems Now, we consider the nearly pentadiagonal Toeplitz linear system (1.1). The matrix A of the form (1.2) can be decomposed as follows: A = A1 + U ′ V ′ ,
(2.8)
where
a − a 0 n −2 a1 − an − 1 a2 A1 =
U = (u1 , u2 ) = ′
V′ =
an − 2 an − 1
a− 1 a0 − an − 2 a1
..
a− 2 a− 1 a0
a− 2 a− 1
a− 2
.
.
.
..
.
..
a2
1 0
0 an − 2
0 1
···
0
0
···
···
0
0
···
..
a1 a2 1 0
a2−n 0
a0 a1 a2
T
0 1
a1−n a2−n
..
.
a− 1 a0 − a2 − n a1
a− 2 a− 1 − a1 − n a0 − a2 − n
,
, .
Proposition 2.3 (Sherman–Morrison–Woodbury Formula). Let A be an n × n nonsingular matrix, X and Y be n × m (n ≥ m) matrices. If Im + Y T A−1 X is nonsingular, then A + XY T is nonsingular. In addition, the inverse of the matrix A + XY T can be explicitly given by
(A + XY T )−1 = A−1 − A−1 X (Im + Y T A−1 X )−1 Y T A−1 . Proof. See [17].
By the Sherman–Morrison–Woodbury formula, the matrix A of (1.2) is invertible if and only if 1 ′ M = I2 + V ′ A − 1 U ,
is invertible, and 1 −1 ′ ′ −1 ′ −1 ′ −1 A−1 = A− V A1 . 1 − A1 U (I2 + V A1 U )
(2.9)
Let y = (y1 , y2 , . . . , yn )T , z = (z1 , z2 , . . . , zn )T , and w = (w1 , w2 , . . . , wn )T be solutions of the following equations A1 y = f,
A1 z = u1 ,
A1 w = u2 ,
respectively. Then, we have M = I2 + V (z, w) = ′
m1 m3
m2 m4
,
(2.10)
where m1 = an−2 z1 + a2−n zn−1 + a1−n zn + 1, m2 = an−2 w1 + a2−n wn−1 + a1−n wn , m3 = an−1 z1 + an−2 z2 + a2−n zn , m4 = an−1 w1 + an−2 w2 + a2−n wn + 1. On both sides of Eq. (2.9), multiply by f on the right to obtain x = y − (z, w)
1 d
m4 −m 3
−m2 m1
an−2 y1 + a2−n yn−1 + a1−n yn an−1 y1 + an−2 y2 + a2−n yn
= y − d1 z − d2 w, where d1 = d2 =
1 d 1 d
[(m4 an−2 − m2 an−1 )y1 + (m4 a1−n − m2 a2−n )yn + m4 a2−n yn−1 − m2 an−2 y2 ], [(m1 an−1 − m3 an−2 )y1 + (m1 a2−n − m3 a1−n )yn − m3 a2−n yn−1 + m1 an−2 y2 ].
1242
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243 Table 1 Total operations for solving the nearly penta-diagonal Toeplitz linear systems.
Operations
Algorithm 1 [13]
NPENTA algorithm [14]
Our algorithm
58n − 151
53n − 142
48n − 24
The resulting algorithm is summarized as follows: Algorithm 2.2 Step 1. Input A1 , U ′ , V ′ , n, f. Step 2. Using Algorithm 2.1, solve equations: A1 y = f, A1 z = u1 , A1 w = u2 , (These systems can be solved in parallel). Step 3. Compute m1 , m2 , m3 and m4 by using (2.10). Step 4. Compute d = m1 m4 − m2 m3 , d1 = 1d [(m4 an−2 − m2 an−1 )y1 + (m4 a1−n − m2 a2−n )yn + m4 a2−n yn−1 − m2 an−2 y2 ], Step 5.
d2 = 1d [(m1 an−1 − m3 an−2 )y1 + (m1 a2−n − m3 a1−n )yn − m3 a2−n yn−1 + m1 an−2 y2 ]. Output the solution of system: x = y − d1 z − d2 w.
Let us look at the Algorithm 2.2 and consequently determine its computational cost. Step1 costs four operations. In step2, obtaining y, z and w, takes 30n − 55, 7n − 15 and 7n − 15 operations respectively. Step3 and Step4 amount to 57 operations. Step5 costs 4n operations. On the whole, we need 48n − 24 operations to compute the solution of the system (1.1). The cost of our algorithm is less than those of the two algorithms, when n ≥ 24. Now, we compare computational costs for Lv and Le’s algorithm [13], NPENTA algorithm [14], and our algorithm in Table 1. 3. Examples We provide the results of numerical examples to illustrate the effectiveness of our algorithm. All tests are performed in MATLAB 7.0 using double precision arithmetic. Example 3.1. First, we consider the following 8 × 8 nearly penta-diagonal Toeplitz linear system
−1
1 1 3 0 0 0 −1 1
2 −1 1 1 3 0 0 0
1 1 3 0 0 0 −1
0 2 −1 1 1 3 0 0
0 0 2 −1 1 1 3 0
0 0 0 2 −1 1 1 3
1 0 0 0 2 −1 1 1
1 x1 −2 1 x2 4 0 x3 6 0 x4 = 6 . 0 x5 6 6 2 x6 3 x7 −1 5
x8
1
Step 1. 2 0 3 0 A1 = 0 0 0 0
U′ =
1 0
−1
2 −1 1 1 3 0 0 0
2 1 3 0 0 0 0 0 1
0 0
0 0
0 2 −1 1 1 3 0 0 0 0
ˆf = (1, 4, 6, 6, 6, 6)T ,
0 0 2 −1 1 1 3 0 0 0
1 0
0 0 0 2 −1 1 1 3
T
0 1
0 0 0 0 2 −1 0 1
,
0 0 0 0 , 0 2 1 0
V′ =
−1 1
0 −1
0 0
0 0
0 0
0 0
1 0
˜f = (3, 5)T .
Step 2. Using Algorithm 2.1, solve equations: A1 y = f,
A1 z = u1 ,
A1 w = u2 ,
then we have y = (2.2990, 1.9901, −0.8040, −0.3921, −1.2376, −0.0059, 5.0178, 6.7188)T , z = (−0.5594, −0.5941, 0.7624, 0.9752, 1.2426, 0.6436, −1.9307, −3.3713)T , w = (0.1802, −0.1980, −0.2792, 0.5584, 0.2475, 0.2812, 0.1564, −1.0238)T .
−2 1
,
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1243
Table 2 Absolute errors for Example 3.2. n
60
100
300
500
1000
∥ x − x∗ ∥
9.2222e−016
3.7304e−015
6.3273e−015
8.7052e−015
1.2851e−014
Step 3–4. Compute m1 , m2 , m3 and m4 by (2.10), then we obtain d = m1 m4 − m2 m3 = 9.0109,
d1 = −2,
d2 = 1.0000.
Step 5. Finally we obtain the solution: x = y − d1 z − d2 w = (1.0000, 1.0000, . . . , 1.0000)T . Example 3.2. Next, we consider the n × n nearly penta-diagonal Toeplitz linear systems Ax = f of the form
−1 1 −1
−1 −1 1
..
.
2 2 −1
2
..
..
.
. −1
..
−1 1
1
−1 −1 .. .
−1
..
.
..
. .
1 −1
.. ..
.
. −1 1
−1 x1 0 1 x2 2 . 0 .. .. . . = .. . 0 2 xn−2 −3 xn−1 −1 −1 xn −1
It can be verified that the solution is x∗ = (1, 1, . . . , 1)T . We used Algorithm 2.2 to compute x. Each error ∥x − x∗ ∥ is provided in Table 2. Here, ∥ · ∥ is the Euclidean vector norm. 4. Concluding remarks In this paper, we derived a numerical algorithm for solving the nearly penta-diagonal Toeplitz linear systems and showed that the computational cost is less than those of the two algorithms in [13,14] when the matrix size n ≥ 24. Numerical examples are given to illustrate the efficiency of the algorithm. Since our algorithms are very easy to implement and they have a strong potential for parallel processing, we believe that they are useful tools for solving nearly penta-diagonal Toeplitz linear systems. References [1] A.S. Householder, The Theory of Matrices in Numerical Analysis, Blaisdel, New York, 1964. [2] W.M. Yan, K.L. Chung, A fast algorithm for solving special tridiagonal systems, Computing 52 (1994) 203–211. [3] S.S. Nemani, L.E. Garey, An efficient method for solving second order boundary value problems with two point boundary conditions, Int. J. Comput. Math. 79 (9) (2002) 1001–1008. [4] E.W. Sachs, A.K. Strauss, Efficient solution of a partial integro-differential equation in finance, Appl. Numer. Math. 58 (2008) 1687–1703. [5] P.G. Martinsson, V. Rokhlin, M. Tygert, A fast algorithm for the inversion of general Toeplitz matrices, Comput. Math. Appl. 50 (2005) 741–752. [6] Y.L. Jiang, R.M.M. Chen, O. Wing, Convergence analysis of waveform relaxation for nonlinear differential–algebraic equations of index one, IEEE Trans. Circuits Syst. I 47 (11) (2000) 1639–1645. [7] L.E. Garey, R.E. Shaw, A parallel method for linear equations with tridiagonal Toeplitz coefficient matrices, Comput. Math. Appl. 42 (2001) 1–11. [8] Y.L. Jiang, A general approach to waveform relaxation solutions of nonlinear differential–algebraic equations: the continuous-time and discrete-time cases, IEEE Trans. Circuits Syst. I. Regul. Pap. 51 (9) (2004) 1770–1780. [9] T. Sogabe, New algorithms for solving periodic tridiagonal and periodic pentadiagonal linear systems, Appl. Math. Comput. 202 (2008) 850–856. [10] M. Chen, On the solution of circulant linear systems, SIAM J. Numer. Anal. 24 (1987) 668–683. [11] G. Heinig, On the reconstruction of Toeplitz matrix inverses from columns, Linear Algebra Appl. 350 (2002) 199–212. [12] F. Diele, L. Lopez, The use of the factorization of five-diagonal matrices by tridiagonal Toeplitz matrices, Appl. Math. Lett. 11 (1998) 51–59. [13] X.G. Lv, J. Le, A note on solving nearly penta-diagonal linear systems, Appl. Math. Comput. 204 (2008) 707–712. [14] S.N. Neossi Nguetchue, S. Abelman, A computational algorithm for solving nearly penta-diagonal linear systems, Appl. Math. Comput. 203 (2008) 629–634. [15] X.J. Xu, A fast algorithm for solving a Toeplitz system, Appl. Math. Comput. 217 (2010) 1944–1948. [16] J.W. Demmel, Numerical Linear Algebra, SIAM, 1997. [17] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore, London, 1996.
Contents lists available at SciVerse ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A new algorithm for solving nearly penta-diagonal Toeplitz linear systems✩ Jiteng Jia a,∗ , Qiongxiang Kong b , Tomohiro Sogabe c a
Department of Mathematical Sciences, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
b
Department of Building Environment and Services Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
c
Graduate School of Information Science and Technology, Aichi Prefectural University, 1522-3 Kumabari, Nagakute, Aichi 480-1198, Japan
article
abstract
info
Article history: Received 20 September 2011 Received in revised form 13 December 2011 Accepted 16 December 2011
A new numerical algorithm for solving nearly penta-diagonal Toeplitz linear systems is presented. The algorithm is suited for implementation using Computer Algebra Systems (CASs) such as MATLAB, MACSYMA and MAPLE. Numerical examples are given in order to illustrate the efficiency of our algorithm. © 2011 Elsevier Ltd. All rights reserved.
Keywords: Nearly penta-diagonal matrices Sherman–Morrison–Woodbury formula Toeplitz matrices Linear systems
1. Introduction Linear systems of equations with nearly penta-diagonal Toeplitz coefficient matrices appear in many applications such as differential equations, parallel computing, matrix algebra, data interpolation, boundary value problems, etc [1–12]. In this paper, we consider an n × n nearly penta-diagonal Toeplitz system of linear equations given by Ax = f,
(1.1)
where
a 0 a1 a2 0 A= . .. 0
an − 2 an − 1
a− 1 a0 a1
a− 2 a− 1 a0
0 a− 2 a− 1
···
0
0 a− 2
···
.
.
.
.
.
..
..
. ···
..
..
.
0
..
..
.
0
···
a2 0
an − 2
0
···
a2 − n 0
.
..
.
··· .. . .. .
a1 a2 0
a0 a1 a2
a− 1 a0 a1
..
..
0
..
a1−n a2−n 0
.. . , 0 a− 2
(1.2)
a− 1 a0
✩ This work was supported by the Natural Science Foundation of China (NSFC) under grant 11071192 and the International Science and Technology Cooperation Program of China under grant 2010DFA14700. ∗ Corresponding author. E-mail addresses: [email protected] (J. Jia), [email protected] (T. Sogabe).
0898-1221/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.camwa.2011.12.044
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1239
x = (x1 , x2 , . . . , xn )T , f = ( f1 , f2 , . . . , fn )T . Xiao-Guang Lv and Jiang Le [13] presented a fast computational algorithm for solving nearly penta-diagonal linear systems based on the use of any penta-diagonal linear solver. By using LU decomposition, Neossi Nguetchue and Abelman [14] gave a numerical algorithm to obtain the solution of nearly pentadiagonal linear systems. Recently, Xiangjian Xu [15] has proposed another efficient algorithm to solve symmetric Toeplitz penta-diagonal linear systems based on the Sherman–Morrison–Woodbury formula. In this study, we derive a numerical algorithm for solving the nearly penta-diagonal Toeplitz linear systems (1.1) and show that the total number of operations is less than those of the two algorithms in [13,14]. The rest of this paper is organized as follows: In Section 2, we describe a more efficient algorithm for solving nearly pentadiagonal Toeplitz linear systems. In Section 3, numerical examples are provided to show the performance and efficiency of our algorithm. Finally, we make some concluding remarks in Section 4. 2. Main results We present two algorithms for solving penta-diagonal linear systems and nearly penta-diagonal Toeplitz linear systems respectively. 2.1. An algorithm for solving penta-diagonal linear systems Consider the penta-diagonal system of linear equations A′ x = f given by
b c a2
a−1 b a1
a− 2 a−1 a0
a− 2 a− 1
a− 2
.
.
.
.
..
..
..
a2
a1 a2
..
..
a0 a1 a2
a−2 e
.
a− 1 d a1
x1 x2 x3
.. . .. .
x n −1 xn
d
f1 f2 f3 .. = . .. .
fn−1 fn
.
(2.1)
Without loss of generality, suppose that the penta-diagonal matrix A′ is nonsingular and that the element a−2 ̸= 0. We now introduce some notations for later use.
V =
b a− 1
U =
c b
··· ···
0 0
a2 a1 0 0
0 a2 a2 0
ˆf = ( f1 , f2 , . . . , fn−2 )T , xˆ = (x1 , x2 )T ,
T
···
0
0
···
a1 a2
d a1
e d
0 0
∈ R(n−2)×2 ,
∈ R2×(n−2) ,
˜f = ( fn−1 , fn )T ,
x˜ = (x3 , x4 , . . . , xn )T ,
ˆf f= ˜f , xˆ x= ˜ . x
According to the above notation, linear system (2.1) can be written in the form
V O
L U
ˆf xˆ = ˜ , x˜ f
(2.2)
where L is a lower triangular Toeplitz matrix of size (n − 2) × (n − 2). Thus (2.2) is equivalent to
V xˆ + Lx˜ = ˆf, U x˜ = ˜f.
(2.3)
Since det L = (a−2 )n−2 ̸= 0, L is invertible. Assume that UL−1 V is nonsingular, it is easy to deduce that xˆ = (UL−1 V )−1 (UL−1 ˆf − ˜f), x˜ = L−1 (ˆf − V xˆ ).
Proposition 2.1. The inverse matrix of a lower triangular Toeplitz matrix is a lower triangular Toeplitz matrix. Proof. See [16].
(2.4)
1240
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
Proposition 2.2. Let (t1 , t2 , . . . , tn )T be the first column of an n × n lower triangular Toeplitz matrix T , (ξ1 , ξ2 , . . . , ξn )T be the first column of T −1 , where T −1 is the inverse matrix of T . Then
ξ1 ξ2 ti , . . . , t2 , t1 .. = 0, . ξi
t1 ξ1 = 1,
i = 2, . . . , n.
Proof. Since (ξ1 , ξ2 , . . . , ξn )T is the first column of T −1 , we have t1 t2
t1
. . .
..
..
. ···
tn
. ···
1 ξ1 ξ2 0 . = . . . . . . ξn 0 t1
By simple algebraic manipulations, the conclusion is true. Obviously, L
−1
is a lower triangular Toeplitz matrix. Let the first column of L−1 be (η1 , η2 , . . . , ηn−2 )T . Then we have
xˆ =
br1 + cr2 + a2 r3 bs1 + cs2 + a2 s3
n −2 ri fi − fn−1 −1 a−1 r1 + br2 + a1 r3 + a2 r4 i =1 , a−1 s1 + bs2 + a1 s3 + a2 s4 n −2 si fi − fn
(2.5)
i =1
where ri = a2 ηn−i−4 + a1 ηn−i−3 + dηn−i−2 + eηn−i−1 for i = 0, 1, . . . , n − 2, sj = rj−1 − eηn−j for j = 1, 2, . . . , n − 2, and ηk = 0 for k ≤ 0 or k ≥ n − 2. After determining xˆ , set f1′ = f1 − bx1 − a−1 x2 , f2′ = f2 − cx1 − bx2 , f3′ = f3 − a2 x1 − a1 x2 , f4′ = f4 − a2 x2 and fi′ = fi for i = 5, 6, . . . , n − 2, then by using backward substitution, we can obtain x˜ as follows: x˜ = (x3 , . . . , xn−1 , xn )T , where x3 =
f1′ a−2
, xi =
1 a−2
′ ( fi− 2 −
(2.6)
i−3
j =1
ai−j−4 xj+2 ) for i = 4, 5, . . . , n, and ak = 0 for k > 2.
Remark 2.1. If ˆf is a standard unit vector, we can obtain x˜ from the second equation of (2.4) directly, and the computational costs are much less than those of (2.6). Thus, it is important to compute the first column elements ηi (i = 1, 2, . . . , n − 2) of L−1 . From Proposition 2.2, we can obtain the general recurrence formula
ηi = − where η1 =
i−1 1
a− 2 j = 1
1 a−2
ai−j−2 ηj ,
i = 2, 3, . . . , n − 2,
(2.7)
and ak = 0 for 2 < k < n − 2.
In the following, we state the algorithm for solving the linear system (2.1): Algorithm 2.1 Step 1. Input a−2 , a−1 , a0 , a1 , a2 , b, c , d, e, n, f. Step 2. Set η1 = a 1 , ak = 0 for 2 < k < n − 2, then compute η2 , η3 , . . . , ηn−2 by using (2.7). −2
Step 3. Step 4.
Step 5.
Step 6. Step 7. Step 8.
Set ηk = 0 for k ≤ 0 or k > n − 2. For i = 0, 1, . . . , n − 2, compute ri = a2 ηn−i−4 + a1 ηn−i−3 + dηn−i−2 + eηn−i−1 . End. For j = 1, 2, . . . , n − 2, compute sj = rj−1 − eηn−j . End. Compute xˆ = (x1 , x2 )T by using (2.5). Set f1′ = f1 − bx1 − a−1 x2 , f2′ = f2 − cx1 − bx2 , f3′ = f3 − a2 x1 − a1 x2 , f4′ = f4 − a2 x2 and fi′ = fi for i = 5, 6, . . . , n − 2, then compute x˜ = (x3 , x4 , . . . , xn−2 )T from the second equation of (2.4) or by using (2.6). xˆ Output the solution of system: x = . x˜
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1241
2.2. Solving nearly penta-diagonal Toeplitz linear systems Now, we consider the nearly pentadiagonal Toeplitz linear system (1.1). The matrix A of the form (1.2) can be decomposed as follows: A = A1 + U ′ V ′ ,
(2.8)
where
a − a 0 n −2 a1 − an − 1 a2 A1 =
U = (u1 , u2 ) = ′
V′ =
an − 2 an − 1
a− 1 a0 − an − 2 a1
..
a− 2 a− 1 a0
a− 2 a− 1
a− 2
.
.
.
..
.
..
a2
1 0
0 an − 2
0 1
···
0
0
···
···
0
0
···
..
a1 a2 1 0
a2−n 0
a0 a1 a2
T
0 1
a1−n a2−n
..
.
a− 1 a0 − a2 − n a1
a− 2 a− 1 − a1 − n a0 − a2 − n
,
, .
Proposition 2.3 (Sherman–Morrison–Woodbury Formula). Let A be an n × n nonsingular matrix, X and Y be n × m (n ≥ m) matrices. If Im + Y T A−1 X is nonsingular, then A + XY T is nonsingular. In addition, the inverse of the matrix A + XY T can be explicitly given by
(A + XY T )−1 = A−1 − A−1 X (Im + Y T A−1 X )−1 Y T A−1 . Proof. See [17].
By the Sherman–Morrison–Woodbury formula, the matrix A of (1.2) is invertible if and only if 1 ′ M = I2 + V ′ A − 1 U ,
is invertible, and 1 −1 ′ ′ −1 ′ −1 ′ −1 A−1 = A− V A1 . 1 − A1 U (I2 + V A1 U )
(2.9)
Let y = (y1 , y2 , . . . , yn )T , z = (z1 , z2 , . . . , zn )T , and w = (w1 , w2 , . . . , wn )T be solutions of the following equations A1 y = f,
A1 z = u1 ,
A1 w = u2 ,
respectively. Then, we have M = I2 + V (z, w) = ′
m1 m3
m2 m4
,
(2.10)
where m1 = an−2 z1 + a2−n zn−1 + a1−n zn + 1, m2 = an−2 w1 + a2−n wn−1 + a1−n wn , m3 = an−1 z1 + an−2 z2 + a2−n zn , m4 = an−1 w1 + an−2 w2 + a2−n wn + 1. On both sides of Eq. (2.9), multiply by f on the right to obtain x = y − (z, w)
1 d
m4 −m 3
−m2 m1
an−2 y1 + a2−n yn−1 + a1−n yn an−1 y1 + an−2 y2 + a2−n yn
= y − d1 z − d2 w, where d1 = d2 =
1 d 1 d
[(m4 an−2 − m2 an−1 )y1 + (m4 a1−n − m2 a2−n )yn + m4 a2−n yn−1 − m2 an−2 y2 ], [(m1 an−1 − m3 an−2 )y1 + (m1 a2−n − m3 a1−n )yn − m3 a2−n yn−1 + m1 an−2 y2 ].
1242
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243 Table 1 Total operations for solving the nearly penta-diagonal Toeplitz linear systems.
Operations
Algorithm 1 [13]
NPENTA algorithm [14]
Our algorithm
58n − 151
53n − 142
48n − 24
The resulting algorithm is summarized as follows: Algorithm 2.2 Step 1. Input A1 , U ′ , V ′ , n, f. Step 2. Using Algorithm 2.1, solve equations: A1 y = f, A1 z = u1 , A1 w = u2 , (These systems can be solved in parallel). Step 3. Compute m1 , m2 , m3 and m4 by using (2.10). Step 4. Compute d = m1 m4 − m2 m3 , d1 = 1d [(m4 an−2 − m2 an−1 )y1 + (m4 a1−n − m2 a2−n )yn + m4 a2−n yn−1 − m2 an−2 y2 ], Step 5.
d2 = 1d [(m1 an−1 − m3 an−2 )y1 + (m1 a2−n − m3 a1−n )yn − m3 a2−n yn−1 + m1 an−2 y2 ]. Output the solution of system: x = y − d1 z − d2 w.
Let us look at the Algorithm 2.2 and consequently determine its computational cost. Step1 costs four operations. In step2, obtaining y, z and w, takes 30n − 55, 7n − 15 and 7n − 15 operations respectively. Step3 and Step4 amount to 57 operations. Step5 costs 4n operations. On the whole, we need 48n − 24 operations to compute the solution of the system (1.1). The cost of our algorithm is less than those of the two algorithms, when n ≥ 24. Now, we compare computational costs for Lv and Le’s algorithm [13], NPENTA algorithm [14], and our algorithm in Table 1. 3. Examples We provide the results of numerical examples to illustrate the effectiveness of our algorithm. All tests are performed in MATLAB 7.0 using double precision arithmetic. Example 3.1. First, we consider the following 8 × 8 nearly penta-diagonal Toeplitz linear system
−1
1 1 3 0 0 0 −1 1
2 −1 1 1 3 0 0 0
1 1 3 0 0 0 −1
0 2 −1 1 1 3 0 0
0 0 2 −1 1 1 3 0
0 0 0 2 −1 1 1 3
1 0 0 0 2 −1 1 1
1 x1 −2 1 x2 4 0 x3 6 0 x4 = 6 . 0 x5 6 6 2 x6 3 x7 −1 5
x8
1
Step 1. 2 0 3 0 A1 = 0 0 0 0
U′ =
1 0
−1
2 −1 1 1 3 0 0 0
2 1 3 0 0 0 0 0 1
0 0
0 0
0 2 −1 1 1 3 0 0 0 0
ˆf = (1, 4, 6, 6, 6, 6)T ,
0 0 2 −1 1 1 3 0 0 0
1 0
0 0 0 2 −1 1 1 3
T
0 1
0 0 0 0 2 −1 0 1
,
0 0 0 0 , 0 2 1 0
V′ =
−1 1
0 −1
0 0
0 0
0 0
0 0
1 0
˜f = (3, 5)T .
Step 2. Using Algorithm 2.1, solve equations: A1 y = f,
A1 z = u1 ,
A1 w = u2 ,
then we have y = (2.2990, 1.9901, −0.8040, −0.3921, −1.2376, −0.0059, 5.0178, 6.7188)T , z = (−0.5594, −0.5941, 0.7624, 0.9752, 1.2426, 0.6436, −1.9307, −3.3713)T , w = (0.1802, −0.1980, −0.2792, 0.5584, 0.2475, 0.2812, 0.1564, −1.0238)T .
−2 1
,
J. Jia et al. / Computers and Mathematics with Applications 63 (2012) 1238–1243
1243
Table 2 Absolute errors for Example 3.2. n
60
100
300
500
1000
∥ x − x∗ ∥
9.2222e−016
3.7304e−015
6.3273e−015
8.7052e−015
1.2851e−014
Step 3–4. Compute m1 , m2 , m3 and m4 by (2.10), then we obtain d = m1 m4 − m2 m3 = 9.0109,
d1 = −2,
d2 = 1.0000.
Step 5. Finally we obtain the solution: x = y − d1 z − d2 w = (1.0000, 1.0000, . . . , 1.0000)T . Example 3.2. Next, we consider the n × n nearly penta-diagonal Toeplitz linear systems Ax = f of the form
−1 1 −1
−1 −1 1
..
.
2 2 −1
2
..
..
.
. −1
..
−1 1
1
−1 −1 .. .
−1
..
.
..
. .
1 −1
.. ..
.
. −1 1
−1 x1 0 1 x2 2 . 0 .. .. . . = .. . 0 2 xn−2 −3 xn−1 −1 −1 xn −1
It can be verified that the solution is x∗ = (1, 1, . . . , 1)T . We used Algorithm 2.2 to compute x. Each error ∥x − x∗ ∥ is provided in Table 2. Here, ∥ · ∥ is the Euclidean vector norm. 4. Concluding remarks In this paper, we derived a numerical algorithm for solving the nearly penta-diagonal Toeplitz linear systems and showed that the computational cost is less than those of the two algorithms in [13,14] when the matrix size n ≥ 24. Numerical examples are given to illustrate the efficiency of the algorithm. Since our algorithms are very easy to implement and they have a strong potential for parallel processing, we believe that they are useful tools for solving nearly penta-diagonal Toeplitz linear systems. References [1] A.S. Householder, The Theory of Matrices in Numerical Analysis, Blaisdel, New York, 1964. [2] W.M. Yan, K.L. Chung, A fast algorithm for solving special tridiagonal systems, Computing 52 (1994) 203–211. [3] S.S. Nemani, L.E. Garey, An efficient method for solving second order boundary value problems with two point boundary conditions, Int. J. Comput. Math. 79 (9) (2002) 1001–1008. [4] E.W. Sachs, A.K. Strauss, Efficient solution of a partial integro-differential equation in finance, Appl. Numer. Math. 58 (2008) 1687–1703. [5] P.G. Martinsson, V. Rokhlin, M. Tygert, A fast algorithm for the inversion of general Toeplitz matrices, Comput. Math. Appl. 50 (2005) 741–752. [6] Y.L. Jiang, R.M.M. Chen, O. Wing, Convergence analysis of waveform relaxation for nonlinear differential–algebraic equations of index one, IEEE Trans. Circuits Syst. I 47 (11) (2000) 1639–1645. [7] L.E. Garey, R.E. Shaw, A parallel method for linear equations with tridiagonal Toeplitz coefficient matrices, Comput. Math. Appl. 42 (2001) 1–11. [8] Y.L. Jiang, A general approach to waveform relaxation solutions of nonlinear differential–algebraic equations: the continuous-time and discrete-time cases, IEEE Trans. Circuits Syst. I. Regul. Pap. 51 (9) (2004) 1770–1780. [9] T. Sogabe, New algorithms for solving periodic tridiagonal and periodic pentadiagonal linear systems, Appl. Math. Comput. 202 (2008) 850–856. [10] M. Chen, On the solution of circulant linear systems, SIAM J. Numer. Anal. 24 (1987) 668–683. [11] G. Heinig, On the reconstruction of Toeplitz matrix inverses from columns, Linear Algebra Appl. 350 (2002) 199–212. [12] F. Diele, L. Lopez, The use of the factorization of five-diagonal matrices by tridiagonal Toeplitz matrices, Appl. Math. Lett. 11 (1998) 51–59. [13] X.G. Lv, J. Le, A note on solving nearly penta-diagonal linear systems, Appl. Math. Comput. 204 (2008) 707–712. [14] S.N. Neossi Nguetchue, S. Abelman, A computational algorithm for solving nearly penta-diagonal linear systems, Appl. Math. Comput. 203 (2008) 629–634. [15] X.J. Xu, A fast algorithm for solving a Toeplitz system, Appl. Math. Comput. 217 (2010) 1944–1948. [16] J.W. Demmel, Numerical Linear Algebra, SIAM, 1997. [17] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore, London, 1996.