In [38]:
using Plots
using Interact

Out[38]:
Plots.GadflyPackage()
In [39]:
"""
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(linspace(-1,1,1000))
@manipulate for n = 2:25
a,b = legendre_coeffs(n)
plot(ts,ortho_eval(a,b,ts))
xlabel!("t").o
end

Out[39]:
In [16]:
sqrt(1/3)

Out[16]:
0.5773502691896257
In [40]:
# Use the eigenvalues to get the zeros!

"""
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)
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)
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)))

5x5 Array{Float64,2}:
0.0       1.0       0.0       0.0       0.0
0.333333  0.0       1.0       0.0       0.0
0.0       0.266667  0.0       1.0       0.0
0.0       0.0       0.257143  0.0       1.0
0.0       0.0       0.0       0.253968  0.0
5x5 Array{Float64,2}:
0.0      0.57735   0.0       0.0       0.0
0.57735  0.0       0.516398  0.0       0.0
0.0      0.516398  0.0       0.507093  0.0
0.0      0.0       0.507093  0.0       0.503953
0.0      0.0       0.0       0.503953  0.0     
In [41]:
ts = collect(linspace(-1,1,1000))
@manipulate for n = 2:10
a,b = legendre_coeffs(n)
plot(ts,ortho_eval(a,b,ts))
lams = eigvals(jacobi_matrix_sym(a[1:n-1],b[1:n-1]))
scatter!(lams, zeros(n-1))
xlabel!("t").o
end

Out[41]:
In [42]:
""" Generate the n-point Gauss-Legendre quadrature rule. """

([-0.7745966692414822,7.771561172376096e-16,0.7745966692414833],[0.555555555555557,0.8888888888888867,0.5555555555555558])