solveTridiag {sdetorus}

Thomas algorithm for solving tridiagonal matrix systems, with optional presaving of LU decomposition

Description

Implementation of the Thomas algorithm to solve efficiently the tridiagonal matrix system

b_1 x_1 + c_1 x_2 + a_1x_n = d_1

a_2 x_1 + b_2 x_2 + c_2x_3 = d_2

\vdots \vdots \vdots

a_{n-1} x_{n-2} + b_{n-1} x_{n-1} + c_{n-1}x_{n} = d_{n-1}

c_n x_1 + a_{n} x_{n-1} + b_nx_n = d_n

with a_1=c_n=0 (usual tridiagonal matrix). If a_1\neq0 or c_n\neq0 (circulant tridiagonal matrix), then the Sherman–Morrison formula is employed.

Usage

solveTridiag(a, b, c, d, LU = 0L)

solveTridiagMatConsts(a, b, c, d, LU = 0L)

solvePeriodicTridiag(a, b, c, d, LU = 0L)

forwardSweepTridiag(a, b, c)

forwardSweepPeriodicTridiag(a, b, c)

Arguments

a, b, c

subdiagonal (below main diagonal), diagonal and superdiagonal (above main diagonal), respectively. They all are vectors of length n.

d

vector of constant terms, of length n. For solveTridiagMatConsts, it can be a matrix with n rows.

LU

flag denoting if the forward sweep encoding the LU decomposition is supplied in vectors b and c. See details and examples.

Details

The Thomas algorithm is stable if the matrix is diagonally dominant.

For the periodic case, two non-periodic tridiagonal systems with different constant terms (but same coefficients) are solved using solveTridiagMatConsts. These two solutions are combined by the Sherman–Morrison formula to obtain the solution to the periodic system.

Note that the output of solveTridiag and solveTridiagMatConsts are independent from the values of a[1] and c[n], but solvePeriodicTridiag is not.

If LU is TRUE, then b and c must be precomputed with forwardSweepTridiag or
forwardSweepPeriodicTridiag for its use in the call of the appropriate solver, which will be slightly faster.

Value

• solve* functions: the solution, a vector of length n and a matrix with n rows for
  solveTridiagMatConsts.

• forward* functions: the matrix cbind(b, c) creating the suitable b and c arguments for calling solve* when LU is TRUE.

References

Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Springer, New York. doi:10.1007/978-1-4899-7278-1

Conte, S. D. and de Boor, C. (1980). Elementary Numerical Analysis: An Algorithmic Approach. Third edition. McGraw-Hill, New York. doi:10.1137/1.9781611975208

Examples

# Tridiagonal matrix
n <- 10
a <- rnorm(n, 3, 1)
b <- rnorm(n, 10, 1)
c <- rnorm(n, 0, 1)
d <- rnorm(n, 0, 1)
A <- matrix(0, nrow = n, ncol = n)
diag(A) <- b
for (i in 1:(n - 1)) {
  A[i + 1, i] <- a[i + 1]
  A[i, i + 1] <- c[i]
}
A

# Same solutions
drop(solveTridiag(a = a, b = b, c = c, d = d))
solve(a = A, b = d)

# Presaving the forward sweep (encodes the LU factorization)
LU <- forwardSweepTridiag(a = a, b = b, c = c)
drop(solveTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1))

# With equal coefficient matrix
solveTridiagMatConsts(a = a, b = b, c = c, d = cbind(d, d + 1))
cbind(solve(a = A, b = d), solve(a = A, b = d + 1))
LU <- forwardSweepTridiag(a = a, b = b, c = c)
solveTridiagMatConsts(a = a, b = LU[, 1], c = LU[, 2], d = cbind(d, d + 1), LU = 1)

# Periodic matrix
A[1, n] <- a[1]
A[n, 1] <- c[n]
A

# Same solutions
drop(solvePeriodicTridiag(a = a, b = b, c = c, d = d))
solve(a = A, b = d)

# Presaving the forward sweep (encodes the LU factorization)
LU <- forwardSweepPeriodicTridiag(a = a, b = b, c = c)
drop(solvePeriodicTridiag(a = a, b = LU[, 1], c = LU[, 2], d = d, LU = 1))