In [38]:
using Plots
using Interact
gadfly()
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]:
t -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
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]:
t -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 y1 y2 y3 y4 y5 y6 y7 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
In [42]:
""" Generate the n-point Gauss-Legendre quadrature rule. """
function gauss_legendre_quadrature(n)
    a,b = legendre_coeffs(n)
    Jsym = jacobi_matrix_sym(a,b)
    x,V = eig(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(3)
Out[42]:
([-0.7745966692414822,7.771561172376096e-16,0.7745966692414833],[0.555555555555557,0.8888888888888867,0.5555555555555558])