Matlab example 1: Legendre polynomials

I recommend you try to write a script that, given a weight function $w(x)$ and an interval $I=[a,b]$ generates the sequence of orthogonal monic polynomials using the three term recurrence we discussed. Then, you can verify that for $w(x) =1$, $I=[-1,1]$ and requiring $p_k(1) = 1$, we obtain exactly the Lagendre polynomials.

close all
clear all


%% Legendre Polynomials
a = -1;
b = 1;
w = @(x) x./x;

k = 10;
l = legendre_poly(k);
n = 100;
x = linspace(a,b,n);

subplot(131)
hold on
yline(0,'--')

for i = 1 : length(l)
    L = legendre_poly(i-1);
    plot(x,L(x))
    ylim([-1 1])
    xlim([-1.1 1.1])
    axis square
    pause(.01)
end
hold off

%%% Compute the scalar products between each polynomial and put it into
%%% the matrix 'orth'
for i = 1 : k
    for j = 1 : k
        orth(i,j) = scalarp(legendre_poly(i-1),legendre_poly(j-1),w,a,b);
    end
end

%%% Color plots the values of the matrix orth
subplot(132)
imagesc(orth)
axis square

%%% Replace all small-enough entries in orth with zeros
orth(orth<1e-16)=0;

%%% Sparsity plot of the new matrix orth
subplot(133)
spy(orth)




function sp = scalarp(f,g,w,a,b)
    fgw = @(x) f(x).*g(x).*w(x);
    sp = integral(fgw,a,b);
end

function p = legendre_poly(k)
    if k == 0
        p = @(x) x./x;
    elseif k == 1
        p = @(x) x;
    else
        L1 = legendre_poly(k-1);
        L2 = legendre_poly(k-2);
        p = @(x) ( (2*(k-1)+1).*x.*L1(x)-(k-1).*L2(x) ) ./ k;
    end
end