%***************************************************************
% Polynomial interpolation to illustrate a "bad" numerical method
% for the following interpolation problem:
%
% Given (x_i,y_i), i=0,1,...,n, find a function f(x) such that
%     f(x_i) = y_i, i=0,1,...,n
%
% If we choose to use polynomials for f, we need to determine
% p(x) = a_n x^n + a_{n-1} x^{n-1} + ... + a_1 x + a_0
% such that p(x_i) = y_i, i=0,1,...,n
% These conditions translate to a system of (n+1) linear equations
% for the unknowns a_0, a_1, ..., a_n, for which we list as a column
% vector "asolu" in the code, where
%
%    asolu(k)=a_{k-1}, k=1,2,...,n+1
%
% Then a perturbed data set is introduced to generate a new coefficient set
% to demonstrate the unstable nature of the polynomial interpolation.
%
% For polynomial evaluations, we need to choose one from either direct evaluation
% or nested evaluation. The latter is found to be more inexpensive, more
% accurate, and easier to program.
%
clear all
n = 10; m = 50;          % degree of the polynomial, number of evaluation points
x = linspace(0,1,n+1);   % equally spaced x points

y = randn(n+1,1);        % randomly chosen y-data points
A = zeros(n+1,n+1);

% construct the matrix associated with x:
A(:,1) = ones(n+1,1);
for k=1:n
    A(:,k+1) = x'.^k;
end

% solve the matrix problem to obtain the polynomial coefficients:
asolu = A\y;

% perturb the y-values to obtain another polynomial coefficient set:
y1 = y;    % copy y to y1 first
y1(n/2) = y1(n/2) + 0.1;  % perturb only the middle component n/2

% solve the same matrix problem with a slightly different right-hand-side
asolu1 = A\y1;

% evaluate the polynomials over the interval [0,1]:

xx = linspace(0,1,m);

switch input('choose the method: (1) for direct, and (2) for nested:')

    case 1
        
    % direct evaluation:
    
    yy = zeros(1,m); yy1 = zeros(1,m);
    xxk = zeros(n+1,1);

    for kk=1:m
        for j=1:n+1
            xxk(j) = xx(kk)^(j-1);
        end
        yy(kk) = sum(asolu.*xxk);
        yy1(kk) = sum(asolu1.*xxk);
    end

    case 2
    % nested evaluation:

    yy = asolu(n+1)*ones(1,m); yy1 = asolu1(n+1)*ones(1,m);

    for k=n:-1:1
        yy = yy.*xx + asolu(k);
        yy1 = yy1.*xx + asolu1(k);
    end
  
    otherwise
        disp('method not identified (must be either 1 or 2)!');
        return;
end

ymin = min(min(yy),min(yy1)) - 1;
ymax = max(max(yy),max(yy1)) + 1;

% plot the original data points and the interpolating polynomials:

figure(1)

hold on;
plot(x,y,'gx',xx,yy,'g-');
xlabel('x'); ylabel('y'); title('Polynomial Interpolation');
axis([0 1 ymin ymax]);
pause;
plot(x,y1,'rx',xx,yy1,'r-');
hold off;

pause

% Use cubic spline to do the same interpolation job:
pp = spline(x,y); pp1 = spline(x,y1);
v = ppval(pp,xx); v1 = ppval(pp1,xx);

figure(2)
hold on;
plot(x,y,'gx',xx,v,'g-');
xlabel('x'); ylabel('y'); title('Cubic Spline Interpolation');
axis([0 1 ymin ymax]);
pause;
plot(x,y1,'rx',xx,v1,'r-');
hold off;


