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