%
%  Math 551 - 3/12/13 demos for interpolation, 
%  ill-conditioning of Vandermonde and Hilbert matrices, etc.
%


help vander
 VANDER Vandermonde matrix.
    A = VANDER(V) returns the Vandermonde matrix whose columns
    are powers of the vector V, that is A(i,j) = v(i)^(n-j).
 
    Class support for input V:
       float: double, single

    Reference page in Help browser
       <a href="matlab:doc vander">doc vander</a>

n=3;
v=1:n

v =

     1     2     3

A=vander(v)

A =

     1     1     1
     4     2     1
     9     3     1

cond(A)

ans =

  70.923133806678663

n=5;
v=1:n

v =

     1     2     3     4     5

A=vander(v)

A =

     1     1     1     1     1
    16     8     4     2     1
    81    27     9     3     1
   256    64    16     4     1
   625   125    25     5     1

cond(A)

ans =

     2.616968797063434e+04

n=10;
v=1:n

v =

     1     2     3     4     5     6     7     8     9    10

A=vander(v)

A =

   1.0e+09 *

  Columns 1 through 3

   0.000000001000000   0.000000001000000   0.000000001000000
   0.000000512000000   0.000000256000000   0.000000128000000
   0.000019683000000   0.000006561000000   0.000002187000000
   0.000262144000000   0.000065536000000   0.000016384000000
   0.001953125000000   0.000390625000000   0.000078125000000
   0.010077696000000   0.001679616000000   0.000279936000000
   0.040353607000000   0.005764801000000   0.000823543000000
   0.134217728000000   0.016777216000000   0.002097152000000
   0.387420489000000   0.043046721000000   0.004782969000000
   1.000000000000000   0.100000000000000   0.010000000000000

  Columns 4 through 6

   0.000000001000000   0.000000001000000   0.000000001000000
   0.000000064000000   0.000000032000000   0.000000016000000
   0.000000729000000   0.000000243000000   0.000000081000000
   0.000004096000000   0.000001024000000   0.000000256000000
   0.000015625000000   0.000003125000000   0.000000625000000
   0.000046656000000   0.000007776000000   0.000001296000000
   0.000117649000000   0.000016807000000   0.000002401000000
   0.000262144000000   0.000032768000000   0.000004096000000
   0.000531441000000   0.000059049000000   0.000006561000000
   0.001000000000000   0.000100000000000   0.000010000000000

  Columns 7 through 9

   0.000000001000000   0.000000001000000   0.000000001000000
   0.000000008000000   0.000000004000000   0.000000002000000
   0.000000027000000   0.000000009000000   0.000000003000000
   0.000000064000000   0.000000016000000   0.000000004000000
   0.000000125000000   0.000000025000000   0.000000005000000
   0.000000216000000   0.000000036000000   0.000000006000000
   0.000000343000000   0.000000049000000   0.000000007000000
   0.000000512000000   0.000000064000000   0.000000008000000
   0.000000729000000   0.000000081000000   0.000000009000000
   0.000001000000000   0.000000100000000   0.000000010000000

  Column 10

   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000
   0.000000001000000

cond(A)

ans =

     2.106257537043968e+12

help hilb
 HILB   Hilbert matrix.
    HILB(N) is the N by N matrix with elements 1/(i+j-1),
    which is a famous example of a badly conditioned matrix.
    See INVHILB for the exact inverse.
 
    HILB(N,CLASSNAME) produces a matrix of class CLASSNAME.
    CLASSNAME must be either 'single' or 'double' (the default).
 
    This is also a good example of efficient MATLAB programming
    style where conventional FOR or DO loops are replaced by
    vectorized statements.  This approach is faster, but uses
    more storage.
 
    See also <a href="matlab:help invhilb">invhilb</a>.

    Reference page in Help browser
       <a href="matlab:doc hilb">doc hilb</a>

hilb(2)

ans =

   1.000000000000000   0.500000000000000
   0.500000000000000   0.333333333333333

H=hilb(2)

H =

   1.000000000000000   0.500000000000000
   0.500000000000000   0.333333333333333

cond(H)

ans =

  19.281470067903971

H=hilb(3)

H =

   1.000000000000000   0.500000000000000   0.333333333333333
   0.500000000000000   0.333333333333333   0.250000000000000
   0.333333333333333   0.250000000000000   0.200000000000000

cond(H)

ans =

     5.240567775860644e+02

H=hilb(4);
cond(H)

ans =

     1.551373873892814e+04

H=hilb(5);
cond(H)

ans =

     4.766072502417230e+05

H=hilb(10);
cond(H)

ans =

     1.602466539329798e+13

% look at ncm code polyinterp.m

help polyval
 POLYVAL Evaluate polynomial.
    Y = POLYVAL(P,X) returns the value of a polynomial P evaluated at X. P
    is a vector of length N+1 whose elements are the coefficients of the
    polynomial in descending powers.
 
        Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)
 
    If X is a matrix or vector, the polynomial is evaluated at all
    points in X.  See POLYVALM for evaluation in a matrix sense.
 
    [Y,DELTA] = POLYVAL(P,X,S) uses the optional output structure S created
    by POLYFIT to generate prediction error estimates DELTA.  DELTA is an
    estimate of the standard deviation of the error in predicting a future
    observation at X by P(X).
 
    If the coefficients in P are least squares estimates computed by
    POLYFIT, and the errors in the data input to POLYFIT are independent,
    normal, with constant variance, then Y +/- DELTA will contain at least
    50% of future observations at X.
 
    Y = POLYVAL(P,X,[],MU) or [Y,DELTA] = POLYVAL(P,X,S,MU) uses XHAT =
    (X-MU(1))/MU(2) in place of X. The centering and scaling parameters MU
    are optional output computed by POLYFIT.
 
    Class support for inputs P,X,S,MU:
       float: double, single
 
    See also <a href="matlab:help polyfit">polyfit</a>, <a href="matlab:help polyvalm">polyvalm</a>.

    Reference page in Help browser
       <a href="matlab:doc polyval">doc polyval</a>

type polyval

function [y, delta] = polyval(p,x,S,mu)
%POLYVAL Evaluate polynomial.
%   Y = POLYVAL(P,X) returns the value of a polynomial P evaluated at X. P
%   is a vector of length N+1 whose elements are the coefficients of the
%   polynomial in descending powers.
%
%       Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)
%
%   If X is a matrix or vector, the polynomial is evaluated at all
%   points in X.  See POLYVALM for evaluation in a matrix sense.
%
%   [Y,DELTA] = POLYVAL(P,X,S) uses the optional output structure S created
%   by POLYFIT to generate prediction error estimates DELTA.  DELTA is an
%   estimate of the standard deviation of the error in predicting a future
%   observation at X by P(X).
%
%   If the coefficients in P are least squares estimates computed by
%   POLYFIT, and the errors in the data input to POLYFIT are independent,
%   normal, with constant variance, then Y +/- DELTA will contain at least
%   50% of future observations at X.
%
%   Y = POLYVAL(P,X,[],MU) or [Y,DELTA] = POLYVAL(P,X,S,MU) uses XHAT =
%   (X-MU(1))/MU(2) in place of X. The centering and scaling parameters MU
%   are optional output computed by POLYFIT.
%
%   Class support for inputs P,X,S,MU:
%      float: double, single
%
%   See also POLYFIT, POLYVALM.

%   Copyright 1984-2007 The MathWorks, Inc.
%   $Revision: 5.16.4.7 $  $Date: 2007/10/15 22:54:57 $

%   DELTA can be used to compute a 100(1-alpha)% prediction interval
%   for future observations at X, as Y +/- DELTA*t(alpha/2,df), where
%   t(alpha/2,df) is the upper (alpha/2) quantile of the Student's t
%   distribution with df degrees of freedom.  Since t(.25,df) < 1 for any
%   degrees of freedom, Y +/- DELTA is at least a 50% prediction interval
%   in all cases.  For large degrees of freedom, the confidence level
%   approaches approximately 68%.

% Check input is a vector
if ~(isvector(p) || isempty(p))
    error('MATLAB:polyval:InvalidP',...
            'P must be a vector.');
end

nc = length(p);
if isscalar(x) && (nargin < 3) && nc>0 && isfinite(x) && all(isfinite(p(:)))
    % Make it scream for scalar x.  Polynomial evaluation can be
    % implemented as a recursive digital filter.
    y = filter(1,[1 -x],p);
    y = y(nc);
    return
end

siz_x = size(x);
if nargin == 4
   x = (x - mu(1))/mu(2);
end

% Use Horner's method for general case where X is an array.
y = zeros(siz_x, superiorfloat(x,p));
if nc>0, y(:) = p(1); end
for i=2:nc
    y = x .* y + p(i);
end

if nargout > 1
    if nargin < 3 || isempty(S)
        error('MATLAB:polyval:RequiresS',...
                'S is required to compute error estimates.');
    end
    
    % Extract parameters from S
    if isstruct(S),  % Use output structure from polyfit.
      R = S.R;
      df = S.df;
      normr = S.normr;
    else             % Use output matrix from previous versions of polyfit.
      [ms,ns] = size(S);
      if (ms ~= ns+2) || (nc ~= ns)
          error('MATLAB:polyval:SizeS',...
                'S matrix must be n+2-by-n where n = length(p).');
      end
      R = S(1:nc,1:nc);
      df = S(nc+1,1);
      normr = S(nc+2,1);
    end

    % Construct Vandermonde matrix for the new X.
    x = x(:);
    V(:,nc) = ones(length(x),1,class(x));
    for j = nc-1:-1:1
        V(:,j) = x.*V(:,j+1);
    end

    % S is a structure containing three elements: the triangular factor of
    % the Vandermonde matrix for the original X, the degrees of freedom,
    % and the norm of the residuals.
    E = V/R;
    e = sqrt(1+sum(E.*E,2));
    if df == 0
        warning('MATLAB:polyval:ZeroDOF',['Zero degrees of freedom implies ' ...
                'infinite error bounds.']);
        delta = repmat(Inf,size(e));
    else
        delta = normr/sqrt(df)*e;
    end
    delta = reshape(delta,siz_x);
end


% rungeinterp
% Modify expression to add input arguments.
% Example:
%   a = [1 2 3; 4 5 6]; 
%   foo(a);

rungeinterp
% see m551_2_28_12_Weier

uiopen('/home/user33/Documents/DeLillo/Math551_S2010/ncm/piecelin.m', true);
x=[1 2 3 4]

x =

     1     2     3     4

y=[-1 3 2 5]

y =

    -1     3     2     5

u=0:.1:5;
v=piecelin(x,y,u);
plot(x,y,'o',u,v)
plot(x,y,'o',u,v)
plot(x,y)

% we looked at ncm code pchiptx - more next time

diary off