Symbolic Algorithm for Solving SLAEs with Heptadiagonal Coefficient Matrices

TL;DR

This paper introduces a symbolic algorithm specifically designed for solving systems of linear equations with heptadiagonal matrices, including a stability condition theorem, advancing computational methods for such structured matrices.

Contribution

The paper presents a novel symbolic algorithm for heptadiagonal SLAEs and provides a proven stability condition, enhancing solution techniques for structured band matrices.

Findings

Algorithm successfully solves heptadiagonal SLAEs

Stability condition theorem established and proven

Potential for improved computational efficiency

Abstract

This paper presents a symbolic algorithm for solving band matrix systems of linear algebraic equations with heptadiagonal coefficient matrices. The algorithm is given in pseudocode. A theorem which gives the condition for the algorithm to be stable is formulated and proven.

Equations

$$\left(\begin{array}{cccccccccc}{d}_0& a_0& b_0& c_0& 0& \dots & \dots & \dots & \dots & 0\\ a_1^{\ast }& d_1& a_1& b_1& c_1& 0& \dots & \dots & \dots & 0\\ b_2^{\star }& a_2^{\star }& d_2& a_2& b_2& c_2& 0& \dots & \dots & 0\\ c_3^{\star }& b_3^{\star }& a_3^{\star }& d_3& a_3& b_3& c_3& 0& \dots & 0\\ 0& c_4^{\star }& b_4^{\star }& a_4^{\star }& d_4& a_4& b_4& c_4& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& c_{N-4}^{\star }& b_{N-4}^{\star }& a_{N-4}^{\star }& d_{N-4}& a_{N-4}& b_{N-4}& c_{N-4}\\ 0& \dots & \dots & 0& c_{N-3}^{\star }& b_{N-3}^{\star }& a_{N-3}^{\star }& d_{N-3}& a_{N-3}& b_{N-3}\\ 0& \dots & \dots & \dots & 0& c_{N-2}^{\star }& b_{N-2}^{\star }& a_{N-2}^{\star }& d_{N-2}& a_{N-2}\\ 0& \dots & \dots & \dots & \dots & 0& c_{N-1}^{\star }& b_{N-1}^{\star }& a_{N-1}^{\star }& d_{N-1}\end{array}\right)\left(\begin{array}{c}{x}_0\\ x_1\\ x_2\\ x_3\\ x_4\\ ⋮\\ x_{N-4}\\ x_{N-3}\\ x_{N-2}\\ x_{N-1}\end{array}\right)=\left(\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ y_4\\ ⋮\\ y_{N-4}\\ y_{N-3}\\ y_{N-2}\\ y_{N-1}\end{array}\right).$$

$$\text{det}(A)=\prod _{i=0}^{N-1}{\mu }_i{|}_{\text{symb}=0}.$$

$$\text{det}(A)$$

$$=M_0\frac{M_1}{{\mu }_0}\frac{M_2}{{\mu }_0{\mu }_1}\dots \frac{M_{N-2}}{{\mu }_0{\mu }_1\dots {\mu }_{N-3}}\frac{M_{N-1}}{{\mu }_0{\mu }_1\dots {\mu }_{N-2}}=$$

$$=\frac{{\prod }_{i=0}^{N-1}{M}_i}{{\mu }_0^{N-1}{\mu }_1^{N-2}{\mu }_2^{N-3}\dots {\mu }_{N-3}^2{\mu }_{N-2}^1}=$$

$$=\frac{{\prod }_{i=0}^{N-1}{M}_i}{M_0^{N-1}\frac{M_1^{N-2}}{{\mu }_0^{N-2}}\frac{M_2^{N-3}}{{\mu }_0^{N-3}{\mu }_1^{N-3}}\dots \frac{M_{N-3}^2}{{\mu }_0^2{\mu }_1^2\dots {\mu }_{N-4}^2}\frac{M_{N-2}^1}{{\mu }_0^1{\mu }_1^1\dots {\mu }_{N-3}^1}}=$$

$$=\frac{{\prod }_{i=0}^{N-1}{M}_i}{M_0^{N-1}\frac{M_1^{N-2}}{{\mu }_0^{N-2}}\frac{M_2^{N-3}}{{\mu }_0^{N-3}{\left(\frac{M_1}{{\mu }_0}\right)}^{N-3}}\dots \frac{M_{N-3}^2}{{\mu }_0^2{\left(\frac{M_1}{{\mu }_0}\right)}^2\dots {\left(\frac{M_{N-4}}{{\mu }_{N-3}}\right)}^2}\frac{M_{N-2}^1}{{\mu }_0^1{\left(\frac{M_1}{{\mu }_0}\right)}^1\dots {\left(\frac{M_{N-3}}{{\mu }_{N-4}}\right)}^1}}=$$

$$=\frac{{\prod }_{i=0}^{N-1}{M}_i}{{\prod }_{i=0}^{N-2}{M}_i}=M_{N-1}$$

Full text

1 Introduction

Systems of linear algebraic equations (SLAEs) with heptadiagonal coefficient matrices may arise after many different scientific and engineering problems, as well as problems of the computational linear algebra where finding the solution of a SLAE is considered to be one of the most important problems. For instance, a semi-implicit formulation for the discretization of the transient terms of the system of partial differential equations (PDEs) which models a multiphase fluid flow in porous media yields to a heptadiagonal system of pressure equations for each time step (see [1]). On the other hand, in [2] the 3D problem, simulating the incompressible blood flow in arteries with a structured mesh domain leads to a heptadiagonal SLAE.

A whole branch of symbolic algorithms for solving systems of linear algebraic equations with different coefficient matrices exists in the literature. [3] considers a tridiagonal matrix and a symbolic version of the Thomas method [4] is formulated. The authors of [5] build an algorithm in the case of a general bordered tridiagonal SLAE, while in [6] the coefficient matrix taken into consideration is a general opposite-bordered tridiagonal one. A pentadiagonal coefficient matrix is of interest in [7], while a cyclic pentadiagonal coefficient matrix is considered in [8]. The latter algorithm can be applied to periodic tridiagonal and periodic pentadiagonal SLAE either by setting the corresponding matrix terms to be zero.

A performance analysis of effective methods (both numerical and symbolic) for solving band matrix SLAEs (with three and five diagonals) being implemented in C++ and run on modern (as of 2018) computer systems is made by us in [9]. Different strategies (symbolic included) for solving band matrix SLAEs (with three and five diagonals) are explored by us in [11]. A performance analysis of effective symbolic algorithms for solving band matrix SLAEs with coefficient matrices with three, five and seven diagonals being implemented in both C++ and Python and run on modern (as of 2018) computer systems is made by us in [10]. Note that the algorithm for solving a SLAE with a heptadiagonal coefficient matrix considered in [10] is the one that is going to be introduced in the next Section.

After obtaining the algorithm independently, it has been found in the article [12] where it is applied for cyclic heptadiagonal SLAEs. Thus, we do not claim out priority to this algorithm. However, the novelties of this work are as follows: pure heptadiagonality, proved necessary requirements, classical Thomas expressions i. e. the algorithm’s formalism follows the form of expressions that are usually used in the Thomas algorithm for tridiagonal SLAE [4], that is, the solution is searched in the form: $y_{i}=\alpha_{i+1}y_{i+1}+\beta_{i+1}$.

2 The Algorithm

Let us consider a SLAE $Ax=y$, where $A$ is a $N\times N$ heptadiagonal matrix,

$A=\textrm{heptadiag}(\mathbf{c^{*}},\mathbf{b^{*}},\mathbf{a^{*}},\mathbf{d},\mathbf{a},\mathbf{b},\mathbf{c})$, $x$ and $y$ are vectors of length $N$:

[TABLE]

A symbolic algorithm for solving SLAEs with a heptadiagonal coefficient matrix is considered. It is based on LU decomposition in which the system $Ax=y$ is rewritten as $LUx=y$, where $L$ and $U$ are a lower triangular and an upper triangular matrices, respectively. The algorithm consists of two steps – the $LU$ decomposition together with the downwards sweep $Lz=y$ happen during the first step, leading us from $Ax=y$ to $Ux=z$, while the upwards sweep (solving $Ux=z$ for $x$) is done during the second step.

Remark: seems that the expression for $g_{i-2}$ is missing in the for loop on p. 435 of [12]. Also, the expression for $k_{1}$ on p. 435 of the same paper is not used anywhere, so it is probably a leftover from a previous algorithm.

Now we shall formulate a symbolic algorithm for solving such a SLAE.

 

Algorithm 1. Symbolic algorithm for solving a SLAE $Ax=y$.

1:$\mathbf{c^{*}},\mathbf{b^{*}},\mathbf{a^{*}},\mathbf{d},\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{y},\varepsilon$

2:$\mathbf{x}$

3:if $det(A)==0$ then

4:     Exit.

5:end if

6:bool flag = False

7:$\mu_{0}:=d_{0}$ $\triangleright$ Step 1.(0)

8:if $|\mu_{0}|<\varepsilon$ then

9:     $\mu_{0}:=$ symb; flag = True

10:end if

11:$\alpha_{0}:=\dfrac{a_{0}}{\mu_{0}};\quad\beta_{0}:=\dfrac{b_{0}}{\mu_{0}};\quad\gamma_{0}:=\dfrac{c_{0}}{\mu_{0}};\quad\delta_{0}:=0.0$

12:$\delta_{1}:=a_{1}^{*};\quad\mu_{1}:=d_{1}-\alpha_{0}\delta_{1}$ $\triangleright$ (1)

13:if !flag then

14:     if $|\mu_{1}|<\varepsilon$ then

15:         $\mu_{1}:=$ symb; flag = True

16:     end if

17:end if

18:$\alpha_{1}:=\dfrac{a_{1}-\beta_{0}\delta_{1}}{\mu_{1}};\quad\beta_{1}:=\dfrac{b_{1}-\gamma_{0}\delta_{1}}{\mu_{1}};\quad\gamma_{1}:=\dfrac{c_{1}}{\mu_{1}}$

19:$\delta_{2}:=a_{2}^{*}-\alpha_{0}b_{2}^{*};\quad\mu_{2}:=d_{2}-\alpha_{1}\delta_{2}-\beta_{0}b_{2}^{*}$ $\triangleright$ (2)

20:if !flag then

21:     if $|\mu_{2}|<\varepsilon$ then

22:         $\mu_{2}:=$ symb; flag = True

23:     end if

24:end if

25:$\alpha_{2}:=\dfrac{a_{2}-\beta_{1}\delta_{2}-\gamma_{0}b_{2}^{*}}{\mu_{2}};\quad\beta_{2}:=\dfrac{b_{2}-\gamma_{1}\delta_{2}}{\mu_{2}};\quad\gamma_{2}:=\dfrac{c_{2}}{\mu_{2}}$

26:$z_{0}:=\dfrac{y_{0}}{\mu_{0}};\quad z_{1}:=\dfrac{y_{1}-\delta_{1}z_{0}}{\mu_{1}};\quad z_{2}:=\dfrac{y_{2}-\delta_{2}z_{1}-b_{2}^{*}z_{0}}{\mu_{2}}$

27:$\xi_{0}:=0.0;\quad\xi_{1}:=0.0;\quad\xi_{2}:=0.0$

28:for $\overline{i=3,\ldots N-4}$ do $\triangleright$ ($i$)

29:     $\delta_{i}:=a_{i}^{*}-\alpha_{i-2}b_{i}^{*}-c_{i}^{*}(\beta_{i-3}-\alpha_{i-3}\alpha_{i-2})$

30:     $\xi_{i}:=b_{i}^{*}-\alpha_{i-3}c_{i}^{*}$

31:     $\mu_{i}:=d_{i}-\alpha_{i-1}\delta_{i}-\beta_{i-2}\xi_{i}-\gamma_{i-3}c_{i}^{*}$

32:     if !flag then

33:         if $|\mu_{i}|<\varepsilon$ then

34:              $\mu_{i}:=$ symb; flag = True

35:         end if

36:     end if

37:     $\alpha_{i}:=\dfrac{a_{i}-\beta_{i-1}\delta_{i}-\gamma_{i-2}\xi_{i}}{\mu_{i}};\quad\beta_{i}:=\dfrac{b_{i}-\gamma_{i-1}\delta_{i}}{\mu_{i}};\quad\gamma_{i}:=\dfrac{c_{i}}{\mu_{i}}$

38:     $z_{i}:=\dfrac{y_{i}-\delta_{i}z_{i-1}-\xi_{i}z_{i-2}-c_{i}^{*}z_{i-3}}{\mu_{i}}$

39:end for

40:$\delta_{N-3}:=a_{N-3}^{*}-\alpha_{N-5}b_{N-3}^{*}-c_{N-3}^{*}(\beta_{N-6}-\alpha_{N-6}\alpha_{N-5})$ $\triangleright$ ($N-3$)

41:$\xi_{N-3}:=b_{N-3}^{*}-\alpha_{N-6}c_{N-3}^{*}$

42:$\mu_{N-3}:=d_{N-3}-\alpha_{N-4}\delta_{N-3}-\beta_{N-5}\xi_{N-3}-\gamma_{N-6}c_{N-3}^{*}$

43:if !flag then

44:     if $|\mu_{N-3}|<\varepsilon$ then

45:         $\mu_{N-3}:=$ symb; flag = True

46:     end if

47:end if

48:$\alpha_{N-3}:=\dfrac{a_{N-3}-\beta_{N-4}\delta_{N-3}-\gamma_{N-5}\xi_{N-3}}{\mu_{N-3}}$

49:$\beta_{N-3}:=\dfrac{b_{N-3}-\gamma_{N-4}\delta_{N-3}}{\mu_{N-3}}$

50:$z_{N-3}:=\dfrac{y_{N-3}-\delta_{N-3}z_{N-4}-\xi_{N-3}z_{N-5}-c_{N-3}^{*}z_{N-6}}{\mu_{N-3}}$

51:$\delta_{N-2}:=a_{N-2}^{*}-\alpha_{N-4}b_{N-2}^{*}-c_{N-2}^{*}(\beta_{N-5}-\alpha_{N-5}\alpha_{N-4})$ $\triangleright$ ($N-2$)

52:$\xi_{N-2}:=b_{N-2}^{*}-\alpha_{N-5}c_{N-2}^{*}$

53:$\mu_{N-2}:=d_{N-2}-\alpha_{N-3}\delta_{N-2}-\beta_{N-4}\xi_{N-2}-\gamma_{N-5}c_{N-2}^{*}$

54:if !flag then

55:     if $|\mu_{N-2}|<\varepsilon$ then

56:         $\mu_{N-2}:=$ symb; flag = True

57:     end if

58:end if

59:$\alpha_{N-2}:=\dfrac{a_{N-2}-\beta_{N-3}\delta_{N-2}-\gamma_{N-4}\xi_{N-2}}{\mu_{N-2}}$

60:$z_{N-2}:=\dfrac{y_{N-2}-\delta_{N-2}z_{N-3}-\xi_{N-2}z_{N-4}-c_{N-2}^{*}z_{N-5}}{\mu_{N-2}}$

61:$\delta_{N-1}:=a_{N-1}^{*}-\alpha_{N-3}b_{N-1}^{*}-c_{N-1}^{*}(\beta_{N-4}-\alpha_{N-4}\alpha_{N-3})$ $\triangleright$ ($N-1$)

62:$\xi_{N-1}:=b_{N-1}^{*}-\alpha_{N-4}c_{N-1}^{*}$

63:$\mu_{N-1}:=d_{N-1}-\alpha_{N-2}\delta_{N-1}-\beta_{N-3}\xi_{N-1}-\gamma_{N-4}c_{N-1}^{*}$

64:if !flag then

65:     if $|\mu_{N-1}|<\varepsilon$ then

66:         $\mu_{N-1}:=$ symb; flag = True

67:     end if

68:end if

69:$z_{N-1}:=\dfrac{y_{N-1}-\delta_{N-1}z_{N-2}-\xi_{N-1}z_{N-3}-c_{N-1}^{*}z_{N-4}}{\mu_{N-1}}$

70:$x_{N-1}:=z_{N-1}$ $\triangleright$ Step 2. Solution

71:$x_{N-2}:=z_{N-2}-\alpha_{N-2}x_{N-1}$

72:$x_{N-3}:=z_{N-3}-\alpha_{N-3}x_{N-2}-\beta_{N-3}x_{N-1}$

73:for $\overline{j=N-4,\ldots 0}$ do

74:     $x_{j}:=z_{j}-\alpha_{j}x_{j+1}-\beta_{j}x_{j+2}-\gamma_{j}x_{x+3}$

75:end for

76:Cancel the common factors in the numerators and denominators of $\mathbf{x}$, making them coprime. Substitute $\textrm{symb}:=0$ in $\mathbf{x}$ and simplify.

Remark: If any $\mu_{i}$ expression has been evaluated to be zero or numerically zero, then it is assigned to be a symbolic variable. We cannot compare any of the next $\mu$ expressions with $\varepsilon$, because any further $\mu$ is going to be a symbolic expression. To that reason, we use a boolean flag which tells us if any previous $\mu$ is a symbolic expression. In that case, comparison with $\varepsilon$ is not conducted as being not needed.

3 Stability of the Algorithm

Some observations on the stability of the proposed algorithm can be made. Firstly, assigning $\mu_{i},i=\overline{0,N-1}$ to be equal to a symbolic variable in case it is zero or numerically zero, ensures correctness of the formulae for computing the solution of the considered SLAE. This action does not add any additional requirements to the coefficient matrix, except:

Theorem 1**.**

The only requirement to the coefficient matrix so as the algorithm to be stable is nonsingularity.

Proof.

As a direct consequence of the transformations done so as the matrix $A$ to be factorized and then the downwards sweep to be conducted, it follows that the determinant of the matrix $A$ in the terms of the introduced notation is:

[TABLE]

(This formula could be used so as the nonsingularity of the coefficient matrix to be checked.) If $\mu_{i}$ for any $i$ is assigned to be equal to a symbolic variable, then it is going to appear in both the numerator and the denominator of the expression for the determinant and so it can be cancelled:

[TABLE]

where $M_{i}$ is the $i$-th leading principal minor, and $\mu_{0}=M_{0}$. This means that the only constraint on the coefficient matrix is $M_{N-1}\neq 0$. ∎

The requirement on the coefficient matrix to be nonsingular is not limiting at all since this is a standard requirement so as the SLAE to have one solution only.

Acknowledgements

The work is partially supported by the Russian Foundation for Basic Research under project #18-51-18005.

4 Conclusions

A symbolic algorithm for solving band matrix SLAEs with heptadiagonal coefficient matrices was presented in pseudocode. Some notes on the stability of the algorithm were made.