In [1]:
## Orthogonal polynomials and tridiagonal matrices

In [2]:
## The Legendre polynomials are orthogonal polynomials over [-1,1]
# From Wikipedia, http://en.wikipedia.org/wiki/Legendre_polynomials
ps = Vector{Function}()
push!(ps, t -> 1.0)
push!(ps, t -> t)
push!(ps, t -> 1/2*(3t^2 - 1.0))
push!(ps, t -> 1/2*(5t^3 - 3t))
push!(ps, t -> 1/8*(35t^4 - 30t^2 + 3))

Out[2]:
5-element Array{Function,1}:
#1
#3
#5
#7
#9
In [3]:
##
using Plots
gr(dpi = 200)
plot() # create a new plot
for i=1:length(ps)
plot!(-1:0.01:1, ps[i], label="\$(i-1)", lw=3.0)
end
plot!()

WARNING: Couldn't initialize gr.  (might need to install it?)
INFO: 
In [4]:
## Now let's look at the three term recurrences
## Three-term recurrence for the
# (n+1) Pn+1(t) = (2n+1)*t*Pn(t) - n*Pn-1(t)
n = 5
B = zeros(Int64,n,n)
b = zeros(Int64,n)
t = 0.0
B[1,1] = 1
b[1] = 1
for i = 1:n-1
B[i+1,i] = -t*(2*(i-1) + 1)
B[i+1,i+1] = i
if i>1
B[i+1,i-1] = i-1
end
end

In [5]:
##
B\b

Out[5]:
5-element Array{Float64,1}:
1.0
0.0
-0.5
0.0
0.375
In [6]:
## We can code this into a function!
function legendre_eval(k,t)
n = k
B = zeros(n,n)
b = zeros(n)
B[1,1] = 1
b[1] = 1
for i = 1:n-1
B[i+1,i] = -t*(2*(i-1) + 1)
B[i+1,i+1] = i
if i>1
B[i+1,i-1] = i-1
end
end
p = B\b
return p
end

Out[6]:
legendre_eval (generic function with 1 method)
In [7]:
##
legendre_eval(10,1/sqrt(2))

Out[7]:
10-element Array{Float64,1}:
1.0
0.707107
0.25
-0.176777
-0.40625
-0.37565
-0.148437
0.127058
0.29834
0.285536
In [8]:
##
k = 55
pts = -1:0.01:1
P = zeros(k,length(pts))
for i = 1:length(pts)
P[:,i] = legendre_eval(k, pts[i])
end
f = plot(pts,P',legend=false, lw=0.5, alpha=0.5)

WARNING: Couldn't initialize gr.  (might need to install it?)
INFO: