Computing exponential of tridiagonal Toeplitz matrices with applications to numerical solution of the heat equation

TL;DR

This paper develops a stable, efficient method for computing the exponential of tridiagonal Toeplitz matrices, with applications to solving the heat equation, including error analysis and extensions to 2D problems.

Contribution

It extends existing methods to general tridiagonal Toeplitz matrices, stabilizes the computation via banded approximations, and applies the approach to heat equation solutions with proven stability.

Findings

The method achieves size-independent complexity.

Error bounds for the exponential approximation are established.

Numerical results confirm efficiency and stability.

Abstract

The computation of the exponential of a tridiagonal matrix and its applications have always been of interest. One application considered here is when the method of lines is used to solve the heat equation, where the equation is transformed into a system of ordinary differential equations (ODEs), and this system has a solution that depends on the exponential of a tridiagonal Toeplitz matrix. Strang and MacNamara presented an approximate method for computing the exponential of a symmetric tridiagonal Toeplitz matrix that appears in the solution of ODEs. Their method is based on approximating the entries of the exponential matrix with modified Bessel functions of the first kind at certain values, and accordingly, the exponential matrix is decomposed as the difference of a Toeplitz matrix and a Hankel matrix. Here, we aim to extend this idea to the general case of tridiagonal Toeplitz matrices and stabilize the method by approximating the matrix exponential with a banded matrix, which makes the complexity of computing the exponential matrix independent of the matrix size. Additionally, we provide an error analysis for these methods and a bound for the entries of the exponential of the tridiagonal Toeplitz matrices. As a main contribution of this work, the idea is implemented to solve the heat equation, and the uniform stability of the method is proved. By using a splitting approach, the method is generalized for two-dimensional problems. Numerical illustrations demonstrate the efficiency of the new methods and bounds.