% Problem interpolating a noisy straight line
% Comparing results of the following 3 methods:
% polynomial interpolation
% cubic splines
% linear least squares
% K. Ming Leung, 04/10/03

clear all; hold off;
ti=[0; 0.2; 0.8; 1.0; 1.2; 1.9; 2.0; 2.1; 2.95; 3.0];
yi=[0.01; 0.22; 0.76; 1.03; 1.18; 1.94; 2.01; 2.08; 2.9; 2.95];

% plot the original data
subplot(2,2,1)
plot(ti,yi,'r.','markerSize',14); % original data
axis([0 3 -1 4]);
title('Original Data: A Noisy Straight Line');

% using finite differences
subplot(2,2,2);
x = dividedDifference1(ti,yi);
t=0:0.001:3;
N=length(x); % N=n+1
p = x(N);
for k=N-1:-1:1
p = x(k)+(t-ti(k)).*p;
end
plot(t,p,'-b','lineWidth',2); hold on;
plot(ti,yi,'r.','markerSize',14); % original data
title('Polynomial Interpolation: Finite Differences');
axis([0 3 -1 4]);

% using a cubic spline
subplot(2,2,3);
[xPlot, yPlot]=cubicSplineNaturalF(ti,yi,10);
plot(xPlot,yPlot,'b','lineWidth',2); hold on;
plot(ti,yi,'r.','markerSize',14); % original data
title('Cubic Splines');
axis([0 3 -1 4]);

% using linear least squares
subplot(2,2,4);
m=length(ti);
A=[ones(m,1) ti]; % Set up matrix A
n=size(A,2);
b=yi;
AT = A.';
ATA = AT*A; %
ATb = AT*b;
% Solve the system of normal equations
L=Cholesky(ATA);
y=forwardSubstitution(L,ATb,n);
x=backSubstitution(L.',y,n);
% Plot the resulting cubic
tPlot=0:0.01:3;
yPlot=x(2).*tPlot+x(1);
plot(tPlot,yPlot,'b-','lineWidth',1.5);
hold on;
plot(ti,yi,'r.','markerSize',14); % original data
title('Linear Least Squares');
axis([0 3 -1 4]);