% CP7p6.m
% Interpolating the square root function using 
% monomial basis functions and cubic spline.
% Computer problem 7.6 on p. 337 of Heath
% Ming Leung, 04/02/04

clear all; hold off;
tV = [0 1 4 9 16 25 36 49 64]';
yV = [0 1 2 3 4 5 6 7 8]';

xPlot = 0:0.02:64;  % x points for plotting

% Interpolation using monomials
n = length(tV); AS = ones(n,n);
for i=2:n, AS(:,i) = tV.*AS(:,i-1); end,
xV = AS\yV, 
% Evaluate the polynomial using Matlab's own 
% polynomial evaluation function polyval:
% Note ordering of coefficients of polynomial in Matlab:
%  Y = a(1)*X^N + a(2)*X^(N-1) + ... + a(N)*X + a(N+1)
yMonomial1 = polyval(xV(n:-1:1),xPlot);
% We can also use our own version of Horner's method:
yMonomial2 = Horner(xV,xPlot);
% The 2 result are the same as you can check:
diffMonomial = max(abs(yMonomial1-yMonomial2)),
plot(xPlot,yMonomial2,'b-','lineWidth',2);
axis([0 2 0 1.5]);
hold on;

% plot the exact square root function
plot(xPlot,sqrt(xPlot),'k-','lineWidth',2); 
xlabel('t','FontSize',18); ylabel('y','FontSize',18);
title('Interpolating the square root from 1 to 64','FontSize',18);

% plot the points to be interpolated as red open circles
plot(tV,yV,'ro');

% Interpolation using cubic spline
[xPlotCSpline, yPlotCSpline] = cubicSplineNaturalF(tV,yV,10);
plot(xPlotCSpline,yPlotCSpline,'g','lineWidth',2);

legend('monomial basis','exact','interpolation points','cubic spline',2);
[xPlotCSpline' (yPlotCSpline-sqrt(xPlotCSpline))' ]