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

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.


Introduction
In this paper, we will focus on the problem of solving where T is a quasi-pentadiagonal Toeplitz matrix.
An n × n matrix T = (t ij ) is said to be Toeplitz if t i,j = t i−j . T is said to be banded Toeplitz if there are positive integers p and q such that p + q = k < n and t v = 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 p altered rows among the first p rows and at most q altered rows among the last q rows. For example, when p = q = 1 , and only the first row and the last row of T are perturbed, then T is said to be quasitridiagonal Toeplitz matrix, the numerical solution of Tx = f for 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 (1) Tx = f 1 3 case when p = q = 2 , and only the first two rows and the last two rows of T are perturbed i.e., when T is a 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 evaluation, 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) is and When x = q = k = h = a , s = z = t = g = c , and d = w = e = p = r = y = b , the matrix T becomes 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 symmetric quasi-pentadiagonal Toeplitz linear systems (1). Then present the numerical results.

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 form (2) |a| > 2(|b| + |c|), c ≠ 0.
(3) T = LU + SV + PQ Journal of Mathematical Chemistry (2021) 59:757-774 S = e 1 e 2 , P = e n−1 e n and e i is the ith column of the identity matrix I n .
From (5), the final solution of (1) is given by 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 where is a matrix of the order 4 × 4 , which is assumed to be non-singular and its inverse is very easy to get. Finally we can obtain the solution x of (1) as Now all we need is to determine u 1 , u 2 , l 1 and l 2 , once these values are determined, we may go to Algorithm 2 to solve our equation.
2.1 Determination of parameters u 1 , u 2 , l 1 and l 2 In this section we discuss how to determine the parameters u 1 , u 2 , l 2 and l 2 . By (3) we have the four equations , and by (10), u 2 = b − cl 2 . Replacing u 1 and u 2 into equations (11) and (12), respectively, we obtain two quadratic equations By solving equation (12) we have in (13), we obtain the following quadratic equation where After a simple calculation, we get the solutions of (15). Furthermore, we get the four roots of (14), they are Knowing l 2 , we may calculate l 1 , u 1 and u 2 easily.

Selection of l 2
In this section, we will discuss the choice of l 2 . There are four l 2 s, the selected l 2 must gaurantee the inverse of L and U exist. Let's look at the structure of L −1 and From here, we can see that, if |l 2 | > 1 , with the increase of n, the down-left corner of L −1 will become larger and larger, and at the end, tends infinity. So |l 2 | < 1 is a sufficient condition to guarantee our process going on. On the orther hand, u2 , then we have that which gives so Therefore |l 2 | < 1 is sufficiently to guarantee that L −1 and U −1 converge.
In the next, we take a, b, c all positive as an example to show that there exists an l 2 which satisfies the required conditions. For the other cases, the arguments are similar.
We take l (2) 2 as an example to prove that l 2 ∈ ℝ.
So l 2 is real, and we end the proof. ◻ Theorem 2 Under the assumption of Theorem 1, we have that |l 2 | < 1.
Since a, b, c are positive, so the left side holds true. Now we prove the right side, that is l 2 < 1.
gives If b − 2c < 0 , then the inequality holds true. In the following, we suppose that b > 2c.
Therefore |l 2 | < 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

Case b = 0
In the previous section, we assume that b ≠ 0 . Here we study what happens when b = 0 . By solving equations (8)-By analogous arguments,(11), we get 6 solutions of the system.

Let's look at solution (1).
For l 1 , u 1 to be real, we need a 2 − 4c 2 ≥ 0 , since our cofficient matrix is diagonally dominant, so this condition is guaranteed. When l 2 = 0 and u 2 = 0, and we need |l 1 | ≤ 1 and |u 1 | ≥ 1 to guarantee the convergence of L −1 and U −1 . We consider first |l 1 | ≤ 1 , that is equivalent to When c > 0, By a straightforward calculation, we get the solution for |l 1 | ≤ 1 , which is a < 0. Now we consider |u 1 | ≥ 1 . Again, by a straightforward calculation, we get the solution for |u 1 | ≥ 1 , which is a ≤ −2 . So when a ≤ −2 , we have that |l 1 | ≤ 1 and |u 1 | ≥ 1.

The algorithm
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 .
For the computational cost, when n is 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 L and U have constant values, respectively. It only reads one vector (the right-hand side 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, x � ] . More precisely, two recursive iteration steps such as , n corresponding to the forward and backward substitutions as in Algorithm 1 are essential for Algorithm 2. When using finite precision arithmetic, we should avoid roundoff error propagation. If all the roots of their characteristic equations 2 + l 1 + l 2 = 0 and 2 + u 2 + c = 0 are all less than unity in magnitude, errors of z i and x n+1−i will smaller than errors of previous values z i−1 and x n+2−i , respectively, in which case the Algorithm 2 is stable.

3 3 Numerical examples
In this section, numerical results are presented to confirm the effectiveness of our algorithm. All algorithms are implemented in MATLAB R2018a and the computations are done on an Intel PENTIUM computer, (2.2 GHz), 6 GB memory. We fix the exact solution to be x * = [1, 1 … , 1] T and the right-hand side vector is set to be f = Ax * .

Experiment 1: Kuramoto Sivashinsky equation
In this experiment, we will take the pentadiagonal matrices which appear in the numerical solution of Kuramoto Sivashinsky (KS) equation as an example. KS equation is a nonlinear partial differential equation first derived for the study of chemical reaction system, see [16,17].
In paper [17], the initial vector is calculated by solving the following linear equations By appling von-Neumann boundary conditions [17], the coefficient matrix becomes By appling second order mixed boundary conditions [16], the coefficient matrix becomes

3
Journal of Mathematical Chemistry (2021) 59:757-774 It can be seen that our proposed algorithm takes less CPU time than the LU method.

Experiment 2
Some artificial examples were used in this experiment. The values of a, b, c, d, e, x, y, z, p, q, r, s, t, w, k, h, and g corresponding to each examples are presented in Table 3. From Tables 4, 5 Journal of Mathematical Chemistry (2021) 59:757-774    Thus, comparing with the initial LU method for a sparse matrix , our algorithm improves the computational cost of the numerical solution remarkably. For the accuracy, our algorithm is similar to the other well-known existing methods.

Conclusion
In this paper, we have proposed a new algorithm for solving diagonally dominant symmetric quasi-pentadiagonal Toeplitz linear systems. We discussed possible choices for the parametres for each situation. We implemented our method in Matlab with respect to computational costs. The numerical results show the robustness of our method. The required memory and the computational time of our algorithm are lower than those of other well-known existing methods. The effectiveness of our algorithm is confirmed by numerical experiments.