export/document 3-term recurrence coefficients

[stevengj]

This package includes functions to compute 3-term recurrence coefficients and related integrals (e.g. jacobi_jacobimatrix and jacobimoment), for many families of orthogonal polynomials. It would be nice if this functionality could be exported and documented.

(For example, it could be used in conjunction with QuadGK.jl for #120.)

(Or should this be in https://jverzani.github.io/SpecialPolynomials.jl? Looks like some of these are currently in both places.)

[dlfivefifty]

There is also functionality in ClassicalOrthogonalPolynomials.jl: there the recurrences are lazy so can be used without allocation (though that makes everything a bit "heavy" in terms of compiler load etc.). But the syntax is pretty straightforward:

julia> using ClassicalOrthogonalPolynomials

julia> jacobimatrix(Legendre())[1:10,1:10]
10×10 BandedMatrix{Float64} with bandwidths (1, 1):
 0.0  0.333333   ⋅    ⋅         ⋅        …   ⋅         ⋅         ⋅ 
 1.0  0.0       0.4   ⋅         ⋅            ⋅         ⋅         ⋅ 
  ⋅   0.666667  0.0  0.428571   ⋅            ⋅         ⋅         ⋅ 
  ⋅    ⋅        0.6  0.0       0.444444      ⋅         ⋅         ⋅ 
  ⋅    ⋅         ⋅   0.571429  0.0           ⋅         ⋅         ⋅ 
  ⋅    ⋅         ⋅    ⋅        0.555556  …   ⋅         ⋅         ⋅ 
  ⋅    ⋅         ⋅    ⋅         ⋅           0.466667   ⋅         ⋅ 
  ⋅    ⋅         ⋅    ⋅         ⋅           0.0       0.470588   ⋅ 
  ⋅    ⋅         ⋅    ⋅         ⋅           0.533333  0.0       0.473684
  ⋅    ⋅         ⋅    ⋅         ⋅            ⋅        0.529412  0.0

julia> A,B,C = ClassicalOrthogonalPolynomials.recurrencecoefficients(Legendre()); A[1:10]
10-element Vector{Float64}:
 1.0
 1.5
 1.6666666666666667
 1.75
 1.8
 1.8333333333333333
 1.8571428571428572
 1.875
 1.8888888888888888
 1.9

FastTransforms.jl also has support for Clenshaw: https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/src/clenshaw.jl

These could be moved into a more basic package for orthogonal polynomials.

[stevengj]

How come it uses a BandedMatrix instead of Tridiagonal? (I thought that the Jacobi matrix is usually defined as the symmetric tridiagonal version? That is just a rescaling, I guess, but would be nice to have a function that does this for you.)

[dlfivefifty]

jacobimatrix(T) is an infinite Tridiagonal:

julia> jacobimatrix(Legendre())
ℵ₀×ℵ₀ LazyBandedMatrices.Tridiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}} with indices OneToInf()×OneToInf():
 0.0  0.333333   ⋅    ⋅         ⋅        …  
 1.0  0.0       0.4   ⋅         ⋅           
  ⋅   0.666667  0.0  0.428571   ⋅           
  ⋅    ⋅        0.6  0.0       0.444444     
  ⋅    ⋅         ⋅   0.571429  0.0          
  ⋅    ⋅         ⋅    ⋅        0.555556  …  
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅        …  
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅        …  
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅        …  
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅        …  
  ⋅    ⋅         ⋅    ⋅         ⋅           
  ⋅    ⋅         ⋅    ⋅         ⋅           
 ⋮                                       ⋱  

But when you take a slice it can't tell its square so just returns a banded matrix.

For the symmetric version one can use:

julia> jacobimatrix(Normalized(Legendre()))
ℵ₀×ℵ₀ LazyBandedMatrices.SymTridiagonal{Float64, FillArrays.Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, LazyArrays.BroadcastVector{Float64, typeof(sqrt), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}, LazyArrays.BroadcastVector{Float64, typeof(/), Tuple{InfiniteArrays.InfStepRange{Float64, Float64}, InfiniteArrays.InfStepRange{Int64, Int64}}}}}}}} with indices OneToInf()×OneToInf():
 0.0      0.57735    ⋅         ⋅        …  
 0.57735  0.0       0.516398   ⋅           
  ⋅       0.516398  0.0       0.507093     
  ⋅        ⋅        0.507093  0.0          
  ⋅        ⋅         ⋅        0.503953     
  ⋅        ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅        …  
  ⋅        ⋅         ⋅         ⋅           
  ⋅        ⋅         ⋅         ⋅           
 ⋮                                      ⋱  

I agree its a confusing abuse of terminology.