# Modifying a tridiagonal matrix

``````N₁=5
# 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
``````
2 Likes

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.

2 Likes