Effective Methods for Solving Band SLEs after Parabolic Nonlinear PDEs

A class of models of heat transfer processes in a multilayer domain is considered. The governing equation is a nonlinear heat-transfer equation with different temperature-dependent densities and thermal coefficients in each layer. Homogeneous Neumann boundary conditions and ideal contact ones are applied. A finite difference scheme on a special uneven mesh with a second-order approximation in the case of a piecewise constant spatial step is built. This discretization leads to a pentadiagonal system of linear equations (SLEs) with a matrix which is neither diagonally dominant, nor positive definite. Two different methods for solving such a SLE are developed -- diagonal dominantization and symbolic algorithms.


Introduction and Mathematical Model
In this note, we focus on solving pentadiagonal (PD) and tridiagonal (TD) systems of linear algebraic equations (SLEs) which are obtained after the discretization of parabolic nonlinear partial differential equations (PDEs), using finite difference methods (FDM) of second-order approximation. Such a problem was solved in [1]. There, a finite difference scheme of first-order approximation was built that leads to a TD SLE with a diagonally dominant coefficient matrix. The system was solved using the Thomas method ( [2]). The following nonlinear model of a cylindrical multilayer structure is considered: ∂u ∂r = 0 ∀ r ∈ {r min , r max }; where (r, z) ∈ Ω ∪ ∂Ω, t ≥ 0; m -index of the subdomain. Equation (1)   Focusing on the radial term of the GE with the assumption that the other terms will be moved to the right-hand side (RHS), we consider the following special mesh with grid points on the inner boundaries: A finite difference scheme with a second-order approximation has the following form (three-point stencils are taken for the GE and the outer boundary conditions (BC), and five-point stencil for the inner The matrix form of the considered system is: A û = ϕ(û), where A is a PD sparse matrix which does not have any special properties, e.g. diagonally dominance or positive definiteness. In order to preserve the band structure of the matrix, we cannot use the Gaussian elimination with pivoting. Within this, we could obtain division by zero at some point of the procedure, which is going to make our algorithm unstable. For that reason, we alter the initial PD matrix by adding the minimum values to the diagonal elements so as to transform the matrix into a weakly diagonally dominant one: .
The Gaussian elimination with pivoting (the procedure we use to transform the PD matrix into a TD one) does not preserve the diagonal dominance of the matrix. The use of the Gaussian elimination to the initial matrix A (not A DD ) yields a transformed matrixÃ. Then, the diagonal dominantization method is used in order to transform the obtained TD matrixÃ into a diagonally dominant one. To that purpose, the nondiagonal elements are added to the diagonal ones:

Numerical and Symbolic Algorithms
Two different approaches for solving the SLE are considered -numerical and symbolic. The complexity of all the suggested numerical algorithms is O(N ). Since it is unknown what stands behind the symbolic library, evaluating the complexity of the symbolic algorithms is a very complicated task. Numerical Algorithms. Two different numerical algorithms are applied to the system with a PD matrix. Both of them are based on LU decomposition. The first one ( [4]) is intended for a dense PD matrix. In the case of the considered problem, the matrix is sparse. For that reason, a modified algorithm is built. The main idea is that after the mesh was defined, the indexes of the discontinuity points are known. Since these indexes coincide with those of the matrix rows which are not sparse, we can reference them to the algorithm and conduct the full calculation only for them. For the rest of the rows, the algorithm is reduced to a problem similar to solving a system with TD matrix. This way, the complexity of the algorithm is decreased but at the cost of additional N + 2 check-ups for the non-sparse rows. In the case of the TD matrix, the system is solved using the Thomas method.
Symbolic Algorithms. The symbolic algorithm in the case of a PD matrix is also based on LU decomposition ( [4]). For the TD matrix, a symbolic version of Thomas method is created. As it is known, Thomas method is not suitable for nondiagonally dominant matrices ( [2]). In order to cope with this problem, in the case of a zero quotient of two subsequent leading principal minors, a symbolic zero is assigned instead and the calculations are continued. At the end of the algorithm, this symbolic zero is substituted with zero. The same approach is suggested in [5].

Implementation and Results
All the numerical algorithms are implemented using C++. The matrix needs to be nonsingular and diagonally dominant so as the methods to be stable. Two symbolic algorithms are implemented, using the GiNaC library (version 1.7.2) ( [6]) of C++. The symbolic algorithms require the matrix to be nonsingular only. In Table 1 one can find the wall-clock time results from the conducted experiments. Since the largest supported precision in the GiNaC library is double, during all the experiments double data type is used. The notation is as follows: NPDM stands for numerical PD method, MNPDM -modified numerical PD method, SPDM -symbolic PD method, NTDM -numerical TD method, STDM -symbolic TD method. The achieved accuracy is summarized, using infinity norm. On the penultimate row of the table, one can find the complexity of all the considered methods. On the last row, the characteristics of the computer which is used are described. The number of considered discontinuity points is K = 11.

Discussion and Conclusions
A nonlinear heat transfer equation in a multilayer domain was considered. The suggested discretization scheme always has a first-order approximation. In the case of piecewise constant thermal conductivities or when λ i+1/2 −λ i−1/2 ∞ = O( h i ∞ ), the order of approximation is going to be second. Focusing on the radial term, a SLE with a PD matrix was obtained. Then, applying Gaussian elimination, a TD matrix was derived. For both these matrices, a diagonal dominantization procedure was suggested. This approach ensures the stability of the suggested methods. A modified version of the numerical method for solving a SLE with a PD matrix was built. Since the complexity of this method is lower than the complexity of the general algorithm (usually K N ), better computational time was achieved. The fastest numerical algorithm was found to come from the Thomas method. All the experiments gave an accuracy of an order of magnitude of 10 −16 . As a next step symbolic algorithms were used. They do not require the matrices to be of a special form and are exact. However, they are not comparable with the numerical algorithms with respect to the required time in the case of a numerical solving of the heat equation when one needs to solve the SLE many times. On the other hand, these symbolic methods are not as restrictive as the numerical ones when it comes to the matrix properties. Another upside of the symbolic algorithms is that in the case of a piecewise linear equation, they do not add nonlinearity to the RHS of the system and hence, there is no need of iterations for the time step to be executed. In future, the approach suggested in this note will be investigated in detail.