% Experiment 1, pg 64
clc;
clear;
%close all;

subplot(2,2,1);
a = 128; % number of subintervals = 2*a
x = (-a:a)'/a; % initialize x discretization of [1,1]
A = [x.^0 x.^1 x.^2 x.^3]; % Construct Vandermonde Matrix
[Q,R] = qr(A,0); % Find reduced QR
scale = Q(end,:); % Select last row of Q
Q = Q*diag(1./scale); % Rescale column by these numbers
plot(x,Q(:,1),'b-',x,Q(:,2),'r-',x,Q(:,3),'g-',x,Q(:,4),'k-');
axis([-1.2 1.2 -1.2 1.2])
xlabel('x axis')
ylabel('P(x)')
title('Plot of Legendre Polynomials')
legend('P_0','P_1','P_2','P_3')
grid on

%figure;
subplot(2,2,2);
plot(x,abs(Q(:,1)-1),'b-',...
     x,abs(Q(:,2)-x),'r-',...
     x,abs(Q(:,3)-3/2*x.^2+1/2),'g-',...
     x,abs(Q(:,4)-5/2*x.^3+3/2*x),'k-');
xlabel('x axis')
ylabel('Error')
title('Error of discrete Legendre Polynomials')
legend('P_0','P_1','P_2','P_3')
grid on

v3 = []; % storage variable for rates of convergence
v2 = [];

for n = 4:10
  a = 2^n;
  x = (-a:a)'/a; % initialize x discretization of [1,1]
  A = [x.^0 x.^1 x.^2 x.^3]; % Construct Vandermonde Matrix
  [Q,R] = qr(A,0); % Find reduced QR
  scale = Q(end,:); % Select last row of Q
  Q = Q*diag(1./scale); % Rescale column by these numbers

  p2_calc = Q(:,3);
  p2_exact = 3/2*x.^2 - 1/2;
  p2_error = abs(p2_calc - p2_exact);
  
  p3_calc = Q(:,4);
  p3_exact = 5/2*x.^3 - 3/2*x;
  p3_error = abs(p3_calc - p3_exact);

  v2 = [v2 p2_error(1:2^(n-2):1+2^n)];		    
  v3 = [v3 p3_error(1:2^(n-2):1+2^n)];		    
end  

n = [4:10]; 
dx = (1./2.^n);

%figure
subplot(2,2,3);
plot(log(dx),log(v2(2,:)),'ro-',...
     log(dx),log(v2(3,:)),'kx-',...
     log(dx),log(v2(4,:)),'g+-',...
     log(dx),log(v2(5,:)),'b^-');
title('rate of convergence for P_2(x) : linear convergence')
xlabel('log(grid spacing)')
ylabel('log(Error)')
legend('x = -0.75','x = -0.50','x = -0.25','x = 0')
grid on
polyfit(log(dx),log(v2(3,:)),1)

%figure
subplot(2,2,4);
plot(log(dx),log(v3(2,:)),'ro-',...
     log(dx),log(v3(3,:)),'kx-',...
     log(dx),log(v3(4,:)),'g+-');
title('rate of convergence for P_3(x) : linear convergence')
xlabel('log(grid spacing)')
ylabel('log(Error)')
legend('x = -0.75','x = -0.50','x = -0.25')
grid on
polyfit(log(dx),log(v3(3,:)),1)