Modifying a tridiagonal matrix

# Coefficient matrix of second-order centered-difference operator (δ²u)ₙ
M           = Tridiagonal(fill(1,N₁-1), fill(-2,N₁), fill(1,N₁-1)); 
M[N₁,1]     = 1; # Periodic boundary conditions
M[1,N₁]     = 1;

results in the following error

ArgumentError: cannot set entry (5, 1) off the tridiagonal band to a nonzero value (1)

Why is such an error thrown? What can I do to perform the operator I want? That is, modifying the entries of the matrix.

Per Linear Algebra · The Julia Language it seems that Tridiagonal is a Type itself!

That error occurs because the matrix would no longer be tri-diagonal. Tridiagonal is a special matrix type that only stores tridiagonal matrices and does not handle any other sparsity pattern.

You need to use some other matrix type. e.g. you could use a dense matrix (just convert to Matrix), but of course that throws away the sparsity. Or you could use a generic sparse matrix type from SparseArrays. (I don’t think anyone has implemented a specialized library for periodic tridiagonal matrices in Julia?)

1 Like

Perhaps a Circulant matrix, or something similar? Using ToeplitzMatrices.jl:

julia> Circulant([-2,1,0,0,1])
5×5 Circulant{Int64, Vector{Int64}}:
 -2   1   0   0   1
  1  -2   1   0   0
  0   1  -2   1   0
  0   0   1  -2   1
  1   0   0   1  -2

We get condensed syntax which is nice, but what is the benefit of using the Circulant matrix versus the following?

M           = spdiagm(-1 => fill(1,N₁-1),0 => fill(-2,N₁)  ,1 => fill(1,N₁-1))
M[N₁,1]     = 1; # Periodic boundary conditions
M[1,N₁]     = 1;

A Circulant-typed matrix is likely to take advantage of the structural properties it implies. For example, many operations (e.g. multiplication, inversion, and solves) with a Circulant matrix are likely to exploit the FFT for better efficiency. A generic sparse matrix will not fully exploit this structure.

Multiplication will likely be faster with the sparse matrix (although would be better-yet via direct circular convolution). Most more-complicated operations (inverses, solves, determinant, etc.) will likely be faster with the circulant matrix. With a lot of work you could custom-make a SparseCirculant type that exploits both structures, but for now I would settle for whichever single type is best-suited to what you’re trying to do.