using Plots
using Blink,Interact
"""
Generate the coefficients of the Legendre polynomials in
three-term recurrence form.
From Gautschi (Numerical Analysis) Table 3.1 (p.178)
"""
function legendre_coeffs(n::Int)
a = zeros(n)
b = zeros(n)
b[1] = 2.
b[2:end] = 1 ./ (4 .- (1.:(n-1)).^(-2))
return a, b
end
"""
Evaluate a set of orthogonal polynomials given
the coefficients of the 3-term recurrence
and a set of points. The result
is a matrix where column j is the
polynomial pi_{j-1} evaluated at points ts.
"""
function ortho_eval(a, b, ts)
n = length(a)
N = length(ts)
P = zeros(N,n)
P[:,1] .= 1.
for j=2:n
if j==2 # annoying start condition
P[:,j] = (ts .- a[j-1]).*P[:,j-1]
else
P[:,j] = (ts .- a[j-1]).*P[:,j-1] - b[j-1].*P[:,j-2]
end
end
return P
end
ts = collect(range(-1,stop=1,length=1000))
ui = @manipulate for n = 2:25
a,b = legendre_coeffs(n)
plot(ts,ortho_eval(a,b,ts))
xlabel!("t")
end
w = Window()
body!(w,ui)
sqrt(1/3)
# Use the eigenvalues to get the zeros!
using LinearAlgebra
"""
Build the Jacobi matrix for a given set of orthogonal polynomials
"""
function jacobi_matrix(a,b)
n = length(a)
# J = diagm(a) +diagm(b[2:end],-1) + diagm(ones(n-1),1)
J = diagm(0=>a) +diagm(-1=>b[2:end]) + diagm(1=>ones(n-1))
end
jacobi_matrix(ab) = jacobi_matrix(ab[1],ab[2]) # little Julia helper
function jacobi_matrix_sym(a,b)
n = length(a)
#J = diagm(a) +diagm(sqrt(b[2:end]),-1) + diagm(sqrt(b[2:end]),1)
J = diagm(0=>a) +diagm(-1=>sqrt.(b[2:end])) + diagm(1=>sqrt.(b[2:end]))
end
jacobi_matrix_sym(ab) = jacobi_matrix_sym(ab[1],ab[2]) # little Julia helper
display(jacobi_matrix(legendre_coeffs(5)))
display(jacobi_matrix_sym(legendre_coeffs(5)))
using Blink,Interact
ts = collect(range(-1,stop=1,length=1000))
ui = @manipulate for n = 2:10
a,b = legendre_coeffs(n)
plot(ts,ortho_eval(a,b,ts),framestyle=:origin)
lams = eigvals(jacobi_matrix_sym(a[1:n-1],b[1:n-1]))
scatter!(lams, zeros(n-1))
xlabel!("t")
end
w = Window()
body!(w,ui)
""" Generate the n-point Gauss-Legendre quadrature rule. """
using LinearAlgebra
function gauss_legendre_quadrature(n)
a,b = legendre_coeffs(n)
Jsym = jacobi_matrix_sym(a,b)
x,V = eigen(Symmetric(Jsym)) # force the symmetric solver
w = V[1,:].^2*b[1]
p = sortperm(x) # get them in sorted order
return x[p],w[p]
end
gauss_legendre_quadrature(100)