%  quick discussion of polyomial interpolation, showing slightly edited
%  versions of the book's m-files for it.

>> type interpn2

  function p = InterpN2(x,y)
% p = InterpN2(x,y)
% The Newton polynomial interpolant.
% x is a column n-vector with distinct components and y is
% a column n-vector. c = p.c is  a column n-vector with the property that if
%
%      p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1))
% then
%      p(x(i)) = y(i), i=1:n.

n = length(x);

for k = 1:n-1
   y(k+1:n) = (y(k+1:n)-y(k:n-1)) ./ (x(k+1:n) - x(1:n-k));  
end

p.c = y; p.x = x;


>> type hornern

function pVal = HornerN(p,z)

% pVal = HornerN(p,z)
% Evaluates the Newton interpolant on z where
% c = p.c and x = p.x are n-vectors and z is an m-vector.
%
% pVal is a vector the same size as z with the property that if
%
%         p(x) = c(1) +  c(2)(x-x(1))+ ... + c(n)(x-x(1))...(x-x(n-1))
% then 
%         pval(i) = p(z(i)) , i=1:m.

pVal = repmat(p.c(end),size(z));

for k=length(p.c)-1:-1:1
   pVal = (z-p.x(k)).*pVal + p.c(k);
end


% Note that I made InterpN2 return a structure, with the two fields c and x,
% and that I made, correspondingly, HornerN, to have just two input arguments
% rather than 3, with the first one the structure p, expected to contain both
% the coefficients and the corresponding points x .

% (To be honest, in class I only talked about editing these two m-files in this
% fashion. Here, I actually made those changes, for the record.)

% Here is an example

>> x = rand(1,10)*(2*pi); y = sin(x);

>> p = interpn2(x,y)
p = 
    c: [1x10 double]
    x: [1x10 double]
>> xx = linspace(0,2*pi);
>> plot(x,y,'o',xx,hornern(p,xx))

% note that this looks very good, even though it uses a polynomial of degree
% <= 9. Compare it to the function being interpolated:

>> hold on, plot(xx,sin(xx),'k'), hold off

%  .. it's the same graphically!

% On to piecewise polynomials, i.e., chapter 3. Will do it slightly
% differently; the notes for this are on the web already.

%  matlab command  diff
%  What is diff(x) in case x has only one entry?
>> diff(1)

ans = Empty matrix: 1-by-0


%  Now a version of the book's piecewise linear interpolation example, but
%  using matlab's commands  mkpp  and  ppval , instead of the book's homemade
%  ones.
>> type myshowpl

% script: myshowpl
n = 10;
x=linspace(0,2*pi,n);
y = sin(x);
pl = mkpp(x,[diff(y)./diff(x) y(1:n-1)]);
plot(x,y,'o')
xx = linspace(-2,2*pi+2,100);
hold on, plot(xx,ppval(pl,xx)), hold off
title('broken line interpolant to sin at linspace(0,2*pi,10)')


%% Now run it:

>> myshowpl

%% discussed at length the use of matlab's ppval ...
%% specifically how to vectorize the locating of a given z within a given
%% increasing sequence x
% stressed the matlab command  [yo,I] = sort(y) 

>> [yo,I] = sort([3 1 2])
yo = 1     2     3
I =  2     3     1

%% ... i.e., the second output tells you how to reorder the input vector in
%% order to achieve the ordered output vector

%  Looked at I  in case  [ignored, I] = sort([x z])  in case  x  is increasing

%% ran the book's script for comparing static with adaptive piecewise linear
%  interpolation

>> showpwl2
diary off