easy way to solve problems with tridiagonal matrix

General Usage
1

I’m learning JuliaLang as an alternative to python for its speed. I would like to know if there exist functions in julia for tridiagonal matrix problem typically encountered in quantum kinetic energy operator in coordinate grid, like eigen value and vectors of the tridiagonal matrix and solution: Tridg*x = b where b’s are known without full form of Tridg matrix.

2

Julia has Tridiagonal and SymTriDiagonal types in the standard library package LinearAlgebra.

For example (Julia 1.0):

julia> using LinearAlgebra

julia> n = 10000

julia> S = SymTridiagonal(randn(n), randn(n-1));

julia> @time evals = eigvals(S);
  2.216645 seconds (13 allocations: 313.047 KiB)
4 Likes
3

Dear @dpsanders, I would like to ask you something. We can write a symmetric tridiagonal matrix by using the Tridiagonal function:

S = Tridiagonal( randn(n-1),randn(n), randn(n-1));

However when I pass this S matrix to eigvals function I obtain:

MethodError: no method matching eigvals!(::Tridiagonal{Float64,Array{Float64,1}})
Closest candidates are: eigvals!(!Matched::SymTridiagonal{#s662,V}

So, what is going on here?
Thanks in advance,
Crist.

4

@ccmejia That is because your S is in general not symmetric, hence its eigenvalues are in general complex. I do not know exactly why, but Julia forbids using eigvals with non-symmetric real matrices, but I guess it uses specialized fast algorithms, which make use of the real entries and the symmetry of the matrix. You have at least 2 options:

  1. Make the matrix symmetric
  2. Do not specialize the tridimensional structure and cast it to an array
S = Array( Tridiagonal( randn(n-1),randn(n), randn(n-1)) )

if it comes in that type or avoid it at all.

1 Like
5

Not every specialized matrix type has solvers for every factorization (eigensolvers, SVD, inv, etcetera), so you might have to convert it to Matrix first.

The Tridiagonal type was initially designed mainly for solving tridiagonal systems (e.g. with \), so it doesn’t have a specialized eigensolver. It would be reasonable to implement one, however: No method `eigen(::Tridiagonal)` · Issue #29396 · JuliaLang/julia · GitHub

3 Likes