A fast algorithm for solving diagonally dominant symmetric quasi-pentadiagonal Toeplitz linear systems

A fast algorithm for solving diagonally dominant symmetric quasi-pentadiagonal Toeplitz linear systems Skander Belhaja,1,∗, Fahd Hcinia, Maher Moakhera, Yulin Zhangb aUniversity of Tunis El Manar, ENIT-LAMSIN, BP 37, 1002, Tunis, Tunisia bCentro de Matem ˜ A¡tica, Universidade do Minho, 4710-057 Braga, Portugal Abstract In this paper, we develop a new algorithm for solving diagonally dominant symmetric quasi-pentadiagonal Toeplitz linear systems. Numerical experiments are given in order to illustrate the validity and efficiency of our algorithm. Keywords: Quasi-pentadiagonal Toeplitz matrix, Diagonally dominant, LU decomposition. 2019 MSC: 15A23, 35L30 1. Introduction In this paper, we will focus on the problem of solving Tx =f(1) where Tis a quasi-pentadiagonal Toeplitz matrix. An n×nmatrix T= (tij) is said to be Toeplitz if ti,j =ti−j.Tis said to be banded Toeplitz if there are positive integers pand qsuch that p+q=k < n and tv= 0 if v > q or v < −p. A banded quasi-Toeplitz matrix is defined to be a banded Toeplitz matrix where there are at most paltered rows among the first5 p rows and at most qaltered rows among the last q rows. For example, when p=q= 1, and only the first row and the last row of Tare perturbed, then Tis said to be quasi-tridiagonal Toeplitz matrix, the numerical solution of Tx =ffor this kind of linear equations was studied by [1], and more general, the numerical solution of block quasi-tridiagonal Toeplitz matrix was studied by [2]. Here we will study the case when p=q= 2, and only the first two rows and the last two rows of Tare perturbed i.e., when Tis a10 quasi-pentadiagonal Toeplitz matrix. Pentadiagonal matrices and quasi-pentadiagonal matrices frequently arise in many application areas, such as computational physics, scientific and engineering computings [3,4,5,6], as well as in the wavefunction formalism [7] and density functional theory [8] in quantum chemistry. The importance of these applications motivated an extensive theoretical study of these kinds of matrices, such as determinant evalu-15 ation, eigenvalues computing and pentadiagonal linear systems solving in the last decades, see for example [9,10] and a large literature therein. In this work, we will present a fast algorithm for the numerical solution of an n×n, nonsingular, diagonally dominant, symmetric quasi-pentadiagonal Toeplitz linear system. In other words, the cofficient matrix of (1) is20 ∗Corresponding author Email address: [email protected] (Skander Belhaj) Preprint submitted to Elsevier March 25, 2021 T=               x y z p q r s c b a b c ............... ............... c b a b c t w k e g d h               . and |a|>2(|b|+|c|), c 6= 0.(2) When x=q=k=h=a,s=z=t=g=c, and d=w=e=p=r=y=b, the matrix Tbecomes a symmetric pentadiagonal Toeplitz matrix. This case was studied in [11,12,13]. For the general case of nonsymmetric pentadiagonal linear systems, algorithms have been introduced in [14,15]. In the following sections, we will introduce an algorithm for solving the diagonally dominant symmetric25 quasi-pentadiagonal Toeplitz linear systems (1). Then present the numerical results. 2. An algorithm for solving quasi-pentadiagonal Toeplitz linear systems In general, we can deal with the quasi-pentadiagonal Toeplitz linear systems (1) as an usual linear systems, and solve it by the LU decomposition without pivoting. Here we give an alterative choice, we factor the quasi-pentadiagonal Toeplitz into the following form30 T=LU +SV +PQ (3) where L=            1 l21 l1l2 ... ......... ......... l1l21            , U =           u1u2c u1u2c ......... ......c u1u2 u1           , l1, l2, u1, u2∈R, S=e1e2, P =en−1enand eiis the ith column of the identity matrix In. V=x−u1y−u2z−c0 0 . . . 0 p−l2u1q−(l2u2+u1)r−(cl2+u2)s−c0. . . 0,and Q=0··· 0t−u1l1w−(l1u2+l2u1)k−(cl1+l2u2+u1)e−(u2+cl2) 0··· 0 0 g−l1u1d−(l1u2+l2u1)h−(cl1+l2u2+u1). By (3), the system (1) becomes (LU +SV +PQ)x=f. (4) Multiplying (4) by (LU)−1, where it is assumed to be non-singular, we obtain the following system (I+ZV +W Q)x=x0(5) where Z= (LU)−1S,W= (LU)−1Pand x0= (LU)−1f. The matricees Z, W and x0can be calculated by Algorithm 1. 2 From (5), the final solution of (1) is given by x= (I+ZV +W Q)−1x0.(6) Next step, we will use Sherman-Morrison-Woodbury inversion formula to give the inverse of (I+ZV + WQ). Let G=Z Wand H=V Q, then ZV +WQ =GH, now apply Sherman-Morrison-Woodbury inversion formula directly to (I+GH)−1, we have that (I+GH)−1=I−G(I+HG)−1H=I−G(I+V Q[Z, W ])−1H=I−[Z, W]N−1V Q where N=I+V Z V W QZ I +QW  is a matrix of the order 4 ×4, which is assumed to be non-singular and its inverse is very easy to get.35 Finally we can obtain the solution xof (1) as x=x0−[Z, W ]N−1V Qx0(7) Now all we need is to determine u1,u2,l1and l2, once these values are determined, we may go to Algorithm 2 to solve our equation. 2.1. Determination of parameters u1,u2,l1and l2 In this section we discuss how to determine the parameters u1,u2,l2and l2. By (3) we have the four equations l1u1=c(8) cl2+u2=b(9) l1u2+l2u1=b(10) cl1+l2u2+u1=a. (11) By (9), u1=c l1, and by (10,) u2=b−cl2. Replacing u1and u2into equations (11) and (12), respectively, we obtain two quadratic equations (b−l2c)l2 1−bl1+cl2= 0 (12) cl1l2 2−bl1l2−(c−al1+cl2 1) = 0.(13) By solving equation (12) we have l1= 1 or l1=cl2 b−cl2 . Case 1: When l1= 1.40 In this case we have that u1=c, l2=b±pb2−4c(a−2c) 2c, and u2=b±pb2−4c(a−2c) 2. 3 These solutions are not stable in digital tests. Case 2: When l1=cl2 b−cl2. (Here we assume that b6= 0, if not l1= 1. ) First we assume that b−cl26= 0. Replacing l1=cl2 b−cl2in (13), we obtain the following quadratic equation l4 2+m1l3 2+m2l2 2+m3l2+m4= 0,(14) where m1=−2b c,m2=b2+ac+2c2 c2,m3=−2bc+ab c2,m4= (b c)2. Let l2=γ−m1 4, then (14) becomes γ4+ξγ2+η= 0,(15) with45 ξ=m2−3m2 1 8=1 2c2(4c2+ 2ac −b2), η=m4−3m4 1 256 +m2 1m2 16 −m1m3 4=1 16c2(b4+ 8b2c2−4acb2). After a simple calculation, we get the solutions of (15). Furthermore, we get the four roots of (14), they are50 l(1) 2=1 2cb−q−2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2 l(2) 2=1 2cb−q2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2 l(3) 2=1 2cb+q−2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2 l(4) 2=1 2cb+q2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2. Knowing l2, we may calculate l1,u1and u2easily. When b−cl2= 0, that is u2= 0, we may get u1=c,l1= 1, and l2=b c. 2.2. Selection of l2 In this section, we will discuss the choice of l2. There are four l2s, the selected l2must gaurantee the55 inverse of Land Uexist. Let’s look at the structure of L−1and U−1. (The inverse can be calculated by A−1=adj(A)/detA). Let L=            1 l21 l1l2 ... ......... ......... l1l21            ,then L−1=            1 π11 π2π1 ... . . .......... . . .......... πn. . . . . . π2π11            where π1=−l2, π2= l21 l1l2 ,···,πn= (−1)n−1  l21. . . 0 l1l2 ... ......1 0. . . l1l2  . 4 From here, we can see that, if |l2|>1, with the increase of n, the down-left corner of L−1will become60 larger and larger, and at the end, tends infinity. So |l2|<1 is a sufficient condition to guarantee our process going on. On the orther hand, U=             u1u2c u1u2c ......... ......... u1u2c u1u2 u1             , then U−1=               φ1φ2. . . . . . . . . . . . φn φ2φ2 .... . . .......... . . .......... . . φ1φ2 . . . φ1φ2 φ1               . where φ1=1 u1, φ2=−u2 u2 1,···,φn= (−1)n−1             u2c . . . 0 u1u2c ......... 0. . . u1u2             un 1. φn=(−1)n−1∗+···+ (n−1)cu1un−3 2+un−1 2 un 1 =(−1)n−1(∗+···+ (n−1)c(u2 u1 )n−31 u2 1 + (u2 u1 )n−11 u1 ), so if |u2 u1|<1, φnis convergent. From u2=b−cl2and l1=cl2 b−cl2=cl2 u2, then we have that u1=c l1 =u2 l2 , which gives l2=u2 u1 , so |u2 u1|<1 is equivalent to |l2|<1. Therefore |l2|<1 is sufficiently to guarantee that L−1and U−1converge.65 In the next, we take a, b, c all positive as an example to show that there exists an l2which satisfies the required conditions. For the other cases, the arguments are similar. We take l(2) 2as an example to prove that l2∈R. Theorem 1. Suppose that a,band care all positive, and let l2=l(2) 2, i.e., l2=1 2c(b−q2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2). Then l2is real. Proof. We first show 4ac +a2−4b2+ 4c2>0. By our hypothsis, the matrix Tis diagonally dominant i.e., |a|>2(|b|+|c|), or |a|−2|c|>|b|. So 4ac +a2−4b2+ 4c2= (a+ 2c)2−(2b)2>(2b)2−(2b)2= 0. 5 Next, we will show 2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2≥0. When −2ac +b2−4c2≥0, the inequality holds true. We consider only 2ac −b2+ 4c2>0. In fact,70 2c√4ac +a2−4b2+ 4c2−2ac +b2−4c2≥0 ⇔4c2(4ac +a2−4b2+ 4c2)≥(2ac −b2+ 4c2)2 ⇔4ac −8c2−b2≥0 By a > 2(b+c), we have 4ac −8c2−b2>4(2(b+c))c−8c2−b2= 8bc −b2. When c > b, then 8bc −b2>8b2−b2>0, the inequality holds true. When c < b, we consider in two cases. (i) c≤b≤2c. In this case, a > 2(b+c)≥4c.75 Then we have b2≤4c2,b2+ 8c2≤4c2+ 8c2= 12c2and 4ac > 4(4c)c= 16c2. So that 4ac −8c2−b2>16c2−(8c2+b2)>16c2−12c2>0 (ii) b > 2c. In this case, a > 2(b+c)≥6cand 2ac > 2(6c)c= 12c2. Since we consider only 2ac > b2−4c2, so we have b2+ 8c2=b2−4c2+ 12c2<2ac + 12c2 Then 4ac −8c2−b2>4ac −(2ac + 12c2)=2ac −12c2>0. So l2is real, and we end the proof.80 Theorem 2. Under the assumption of Theorem 1, we have that |l2|<1. Proof. We first prove that |l2|<1, i.e., −1< l2<1. We begin by proving the left side, that is −1< l2. −2c < b −p2c√4ac +a2−4b2+ 4c2−2ac +b2−4c2 =⇒(2c+b)2>p2c√4ac +a2−4b2+ 4c2−2ac +b2−4c22 85 =⇒2ac + 4bc + 8c2>2c√4ac +a2−4b2+ 4c2 =⇒2a+ 4b+ 8c > 2√4ac +a2−4b2+ 4c2 90 =⇒(2a+ 4b+ 8c)2>(2√4ac +a2−4b2+ 4c2)2 =⇒4a2+ 16ab + 32ac + 16b2+ 64bc + 64c2>4a2+ 16ac −16b2+ 16c2 6 =⇒4a2+ 16ab + 32ac + 16b2+ 64bc + 64c2−(4a2+ 16ac −16b2+ 16c2)>095 =⇒32b2+ 64bc + 16ab + 48c2+ 16ac > 0 Since a, b, c are positive, so the left side holds true. Now we prove the right side, that is l2<1.100 b−q2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2<2c gives b−2c < q2cp4ac +a2−4b2+ 4c2−2ac +b2−4c2. If b−2c < 0, then the inequality holds true. In the following, we suppose that b > 2c. (b−2c)2<(p2c√4ac +a2−4b2+ 4c2−2ac +b2−4c2)2 105 =⇒ −4bc + 4c2<2c√a2+ 4ac −4b2+ 4c2−2ac −4c2 =⇒(2√a2+ 4ac −4b2+ 4c2)2>(−4b+ 8c+ 2a)2 =⇒4a2+ 16ac −16b2+ 16c2>4a2−16ab + 32ac + 16b2−64bc + 64c2 110 =⇒ −2b2+ 4bc +ab −3c2−ac > 0 Since b > 2cand a > 2(b+c), so −2b2+ 4bc +ab −3c2−ac =−2b2+ 4bc −3c2+ab −ac =−2b2+ 4bc −3c2+a(b−c)>115 −2b2+ 4bc −3c2+ 2(b+c)(b−c)=4bc −5c2>4(2c)c−5c2>0. Therefore |l2|<1 is true. So the proof of this theorem is concluded. According to the signs of a, b, c, we give the following table for the selection of l2. a b c l2 +++l(2) 2 −+ + l(1) 2 +−+l(2) 2 + + −l(2) 2 − − +l(3) 2 +− − l(4) 2 −+−l(1) 2 −−−l(3) 2 2.3. Case b= 0120 In the previous section, we assume that b6= 0. Here we study what happens when b= 0. By solving equations (8)-(11), we get 6 solutions of the system. 7 (1) l1=1 c1 2a+1 2√a2−4c2, l2= 0, u1=1 2a−1 2√a2−4c2, u2= 0. (2) l1=−1 c−1 2a+1 2√a2−4c2, l2= 0, u1=1 2a+1 2√a2−4c2, u2= 0. (3) l1= 1, l2=−1 cp−c(a−2c), u1=c, u2=p−c(a−2c). (4) l1= 1, l2=1 cp−c(a−2c), u1=c, u2=−p−c(a−2c). (5) l1=−1, l2=−1 cpc(a+ 2c), u1=−c, u2=p−c(a+ 2c). (6) l1=−1, l2=1 cpc(a+ 2c), u1=−c, u2=−p−c(a+ 2c). Let’s look at solution (1). l1=1 c1 2a+1 2pa2−4c2, l2= 0, u1=1 2a−1 2pa2−4c2, u2= 0. For l1,u1to be real, we need a2−4c2≥0, since our cofficient matrix is diagonally dominant, so this125 condition is guaranteed. When l2= 0 and u2= 0, L−1=                   1 0 1 −l10... 0−l1 ...... l2 10−l1 ...... 0l2 10−l1 ...... . . ................... (−1)nl(n−1 2) 1. . . 0l2 10−l10 1                   and U−1=                  1 u10−c u2 10. . . . . . (−1)nc(n−1 2) un−1 1 1 u10−c u2 10. . . ............. . . .........0 ......−c u2 1 ...0 1 u1                  , we need |l1| ≤ 1 and |u1| ≥ 1 to guarantee the convergence of L−1and U−1. We consider first |l1| ≤ 1, that is equivalent to |a+pa2−4c2| ≤ 2|c|. When c > 0, |a+pa2−4c2| ≤ 2|c| ⇐⇒ −2c≤a+pa2−4c2≤2c. 8 When c < 0, |a+pa2−4c2| ≤ 2|c| ⇐⇒ 2c≤a+pa2−4c2≤ −2c. By a straightforward calculation, we get the solution for |l1| ≤ 1, which is a < 0. Now we consider |u1| ≥ 1. Again, by a straightforward calculation, we get the solution for |u1| ≥ 1,130 which is a≤ −2. So when a≤ −2, we have that |l1| ≤ 1 and |u1| ≥ 1. By analogous arguments, we get that when a≥2, |l1| ≤ 1 and |u1| ≥ 1 are guaranteed by solution (2), i.e., . (2) l1=−1 c−1 2a+1 2pa2−4c2, l2= 0, u1=1 2a+1 2pa2−4c2, u2= 0. So when a≤ −2, we choose solution (1) and when a≥2 we choose solution (2). And when a∈(−2,2), we may simply solve the system 2 aTx =2 af. Remark 1. The other four solutions are not suitable for the case b= 0. 2.4. The algorithm135 In this subsection we give an algorithm for solving (1). We first give the Algorithm 1 to solve LUy =f, then Algorithm 2 to solve the equation (1), i.e., Tx =f. Algorithm 1 An algorithm for solving LUy =f Input: l1,l2,u1,u2,cand f 1. (solving Lz =f)z1=f1,z2=f2−l1z1,zi=fi−l1zi−1−l2zi−2,i= 3 to n 2. (solving Uy =z)yn=zn u1,yn−1=zn−1−u2yn u1,yi=zi−u2yi+1−cyi+2 u1,i=n−2 to 1 Output: y= [y1, y2...,yn]T. Algorithm 2 : An algorithm for solving Tx =f Input: a,b,c,d,e,x,y,z,p,q,r,s,h,w,k,g,tand f; 1. Find the parameters liand ui(i= 1,2); 2. Solve linear systems LUZ =S,LUW =P, and LUx0=fby using Algorithm 1; Output: Compute xby (7). For the computational cost, when nis large this algorithm takes about ≈8n+O(1) flops. An advantage of our algorithm is that it needs less data transmission since both the subdiagonal and superdiagonal of Land Uhave constant values, respectively. It only reads one vector (the right-hand side140 vector) and writes one vector (the solution). The stability of Algorithm 2 depends on the step that solves the upper and lower pentadiagonal linear systems LDU[Z, W, f] = [S, P, x0]. More precisely, two recursive iteration steps such as zi=fi−1−l1zi−1−l2zi−2 and xn+1−i=zn+1−1−u2yn+2−i−cyn+3−ifor i= 2, . . . , n corresponding to the forward and backward substitutions as in Algorithm 1 are essential for Algorithm 2. When using finite precision arithmetic, we145 should avoid roundoff error propagation. If all the roots of their characteristic equations λ2+l1λ+l2= 0 and λ2+u2λ+c= 0 are all less than unity in magnitude, errors of ziand xn+1−iwill smaller than errors of previous values zi−1and xn+2−i, respectively, in which case the Algorithm 2 is stable. 9