% script: PPfcns.m
% last rev: 11/05/09
% description: this package contains a discussion of structures and
% routines for implementing a piecewise polynomial package suitable
% for Octave and Matlab. Because of incompatibilities, the data
% structures are kept simple and do not use all the features of
% either language.  One could add some of these features to the package.
% Matlab users will have to unpackage functions of this file into 
% separate function files for each of the routines below, while
% Octave users can simply read in this single file. 
% The basic structure, whether formally declared or not, is a piecewise
% polynomial (p.p.), represented by a PP structure. An instance of this
% structure named pp has the following properties (components): 
%	pp.order: the order of the p.p. approximation (order = degree+1)
%	pp.knum: the number of knots for this p.p.
%	pp.knots: a row vector of size pp.order, consisting of knots in
%                 strictly ascending order.
%	pp.coefs: a matrix of pp.knum - 1 rows, pp.order columns, each row
%                 of which is array of coefficients of a degree m-1
%                 polynomial in the form [a_1, a_2, ... , a_m], where
%                 si(x) = a_1(x-xi)^(m-1) + ... + a_{m-1}(x-xi) + a_m.
%                 (This is the standard form for a "polynomial" vector in
%                 Matlab and Octave.  Thus, their builtin routines like
%                 polyeval are available for use.  However, the polynomial
%                 variable for these routines is y = x-xi, not x.)
% 
% Matlab and Octave handle objects rather differently.  Matlab has a full
% featured OOP environment with classes, objects with properties,methods,
% inheritance, operator overloading, etc. However, Octave has only 
% C-style structures (last time I checked -- this may have changed).
% Since our packages are on the simple side, we won't worry about these 
% details and simply use structures with parts, with our self-enforced 
% convention that PP structures have the properties given above. 
% Furthermore, we only do a little bit of bulletproofing.  It is advisable
% to use the PP constructor to ensure a correct PP structure.  Outside the 
% constructor one could, for example, define a PP structure with
% non-strictly increasing knots vector, but results would be unpredictable.

%
% not a function file
1;
%

function pp = PP(order,knum,knots,coefs)
% usage: pp = PP(order,knum,knots,coefs)
% description: given input of at least positive integers order and
% knum, this function constructs a PP object with this order and
% knum knots. Last two arguments are optional, but if entered,
% user must define knots before coefs. If these are not supplied,
% the  the coef matrix is initialized to zeros and knots to 0:knum-1.

% local variables:
% j: loop variable

if (nargin < 2 | nargin >4)
  error('Incorrect paramater list');
end;
pp.order = order;
pp.knum = knum;
pp.knots = 0:(knum-1);
pp.coefs = zeros(knum-1,order);
if (nargin >= 3)
  if (sum([1,knum] ~=size(knots)))
    clear pp;
    error('Incorrect knots vector')
  end;
  for j = 2:knum
    if (knots(j) <= knots(j-1))
      clear pp;
      error('Knots vector is not strictly increasing')
    end;
  end;
  pp.knots = knots;
end;
if (nargin == 4)
  if (sum([knum-1,order] ~=size(coefs)))
    clear pp;
    error('Incorrect coefs matrix')
    return;
  end;
  pp.coefs = coefs;
end;

endfunction % Matlab users must cut above this line for a function file

%

function y = PPeval(pp,x)
% usage: y = PPeval(pp,x)
% description:  given input pp a PP structure and x a vector of
% absicssas, this routine returns the vector of values of the PP
% structure pp at the coordinates of x.
 
% local variables: 
% j,k,m:  index variables
% n: size of node array
% klim: number of knots
% ord: order of pp
% xtmp: temporary variable

n = length(x);

y = zeros(size(x));
klim = pp.knum;
ord = pp.order;
for j = 1:n % loop over entries of x 
  xtmp = x(j); 
  % find index of interval containing entry with binary search
  low = 1;
  high = klim;
  while (low + 1 < high)
    mid = fix((low+high)/2);
    if (xtmp >= pp.knots(mid))
      low = mid;
    else
      high = mid;
    end;
  end;
  dx = xtmp - pp.knots(low);
  % use Horner's method to evaluate the polynomial
  y(j) = pp.coefs(low,1);
  for m = 2:ord
    y(j) = dx*y(j)+pp.coefs(low,m);
  end;
end;

endfunction % Matlab users must cut above this line for a function file

%

function dd = PPder(pp)
% usage: dd = PPder(pp)
% description: given input pp a PP structure, this function
% returns the PP structure representing the derivative of the p.p.
% function represented by pp.

% local variables: 
% j,k:  index variables
% klim: number of knots
% ord: order of pp

dd.knum = pp.knum;
dd.knots = pp.knots;
dd.order =pp.order -1;
dd.coefs = zeros(dd.knum -1, dd.order);
deg = pp.order - 1;
for j = 1:dd.knum-1
  for k = 1:dd.order
    dd.coefs(j,k) = (deg-k+1)*pp.coefs(j,k);
  end;
end;

endfunction % Matlab users must cut above this line for a function file

%

function pp = LCpp(x,f) 
% usage: pp = LCpp(x,f) 
% description:  this function returns a PP structure representing
% the Lagrange cubic p.p. with knot row vector consisting of 
% the first entry and every subsequent third entry of x, counting
% from the second entry of x.  It is understood that the ordered
% node row vector x includes the intermediate nodes needed for
% Lagrange interpolation and that these are not counted as knots.
% Therefore, x and f must have length 3*k-2 where
%     k = number of knots.

% local variables: 
% n,m: size of node and value arrays
% j,knotindex:  index variables

n = length(x);
m= size(f);
m = m(1)+m(2)-1;
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and ordinate');
  return;
end;
if (rem(n,3)~=1)
  disp('Abscissa size is not 3*k+1');
  return;
end;
for j = 2:n
  if (x(j) <= x(j-1))
    disp('Nodes vector is not strictly increasing')
    return;
  end;
end;
% size ok, construct PP structure
pp = PP(4, round((n-1)/3)+1, x(1:3:n));
% determine the poly coefs in each knot subinterval
knotindex = 1;
for j = 1:3:(n-1)
  % set up and solve the system on this subinterval
  pp.coefs(knotindex,:) = (vander(x(j:j+3) - x(j))\f(j:j+3)')'; 
  knotindex = knotindex+1;
end;

endfunction % Matlab users must cut above this line for a function file.

%

function pp = HCpp(x,f,fp) 
% usage: pp = HCpp(x,f,fp) 
% description: this function returns a PP structure representing
% the Hermite cubic p.p. with node row vector x, corresponding
% function values in the row vector f and derivative values in the
% row vector fp.


% local variables: 
% n,m: size of node and value arrays
% j: index variable
% h: stepsize

n = length(x);
m= size(f);
m = m(1)+m(2)-1;
if (m ~= n) % error check for size and order
  error('Size mismatch of abscissa and function ordinate');
end;
m = length(fp);
if (m ~= n) % error check for size and order
  error('Size mismatch of abscissa and derivative ordinate');
end;
pp = PP(4, n, x);
% determine the poly coefs in each knot subinterval
for j = 1:(n-1)
  % use differences to determine the coefficients
  h = x(j+1) - x(j);
  fab = (f(j+1)-f(j))/h;
  pp.coefs(j,:) = [(fp(j)+fp(j+1)-2*fab)/(h*h),(3*fab-2*fp(j)-fp(j+1))/h,fp(j), f(j)];
end;

endfunction % Matlab users must cut above this line for a function file.

%

function pp = CCSpp(x,f,fpl,fpr)
% usage: pp = CCSpp(x,f,fpl,fpr)
% description:  this function returns the clamped cubic spline
% structure with knots at the entries of the row vector of abscissas
% x and corresponding ordinates in the row vector f. The scalars
% fpl and fpr are values of f'(x) at left and right endpoints of x.

% local variables:
% n,m: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

n = length(x);
m= length(f);
if (m ~= n) % error check for size and order
  error('Size mismatch of abscissa and function ordinate');
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
% now set up the system as in ClassroomNotes, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row
Coef(1, 1:2) = [h(1)/3, h(1)/6];
Rhs(1) = d(1) - fpl;
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row
Coef(n, (n-1):n) = [h(n-1)/6, h(n-1)/3]; 
Rhs(n) = fpr - d(n-1);
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaluation gives di = yi
coefs(:,4) = f(1:n-1)';
% evaluation of s' gives formula for ci
coefs(:,3) = d - ((M(2:n) + 2*M(1:n-1))'/6).*h;
pp = PP(4,n,x,coefs);

endfunction % Matlab users must cut above this line for a function file.

%

function pp = NKCSpp(xx,ff)
% usage: pp = NKCSpp(x,f)
% description: this function returns the not-a-knot cubic spline
% with nodes at the entries of the row vector of abscissas x and
% corresponding ordinates in the row vector f, excluding the 
% second and penultimate entries of each.

% local variables:
% x,f: vectors with not-a-knot entries removed
% nn,m,n: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% h0, h1, hnm2, hnm1: extra step sizes
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

nn = length(xx);
m= length(ff);
if (m ~= nn | nn<4)  % error check for size and order
  error('Size mismatch of abscissa and function ordinate or size problem');
end;
x = [xx(1), xx(3:nn-2), xx(nn)];
f = [ff(1), ff(3:nn-2), ff(nn)];
n = nn - 2; % adjust size to reflect size of x
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
hh = xx(2:nn) - xx(1:nn-1);
% now set up the system as in ClassNotes, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row 
Coef(1, 1:2) = [hh(1)+2*hh(2), 2*hh(1)+hh(2)];
Rhs(1) = 6*(ff(1)/hh(1) - h(1)*ff(2)/(hh(1)*hh(2)) + ff(3)/hh(2));
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row 
Coef(n, (n-1):n) = [hh(nn-2)+2*hh(nn-1), 2*hh(nn-2)+hh(nn-1)]; 
Rhs(n) = 6*(ff(nn-2)/hh(nn-2) - h(n-1)*ff(nn-1)/(hh(nn-2)*hh(nn-1)) + ff(nn)/hh(nn-1));
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaluation gives di = yi
coefs(:,4) = f(1:n-1)';
% evaluation of s' gives formula for ci
coefs(:,3) = d - ((M(2:n) + 2*M(1:n-1))'/6).*h;
pp = PP(4,n,x,coefs);


endfunction % Matlab users must cut above this line for a function file.

%

function pp = NCSpp(x,f)
% usage: pp = NCSpp(x,f)
% description: this function returns the natural cubic spline
% with nodes at the entries of the row vector of abscissas x and
% corresponding ordinates in the row vector f. 

% local variables:
% n,m: size of node and value arrays
% j: index variable
% h: stepsize vector
% d: divided difference vector
% Coef: coefficient matrix for the system
% Rhs: right hand side vector
% M: vector of second derivative values
% coefs: coefficient matrix of PP structure to be returned

n = length(x);
m= length(f);
if (m ~= n) % error check for size and order
  error('Size mismatch of abscissa and function ordinate');
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
d = (f(2:n) - f(1:n-1))./h;
% now set up the system as in ClassroomNotes, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row
Coef(1, 1) = 1;
Rhs(1) = 0;
for j = 2:(n-1) % interior rows
  Coef(j, (j-1):(j+1)) = [h(j-1)/6, (h(j-1)+h(j))/3, h(j)/6];
  Rhs(j) = d(j) - d(j-1);
end;
% last row
Coef(n, n) = 1; 
Rhs(n) = 0;
% solve for M
M = Coef\Rhs;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai*(x-xi)^3+bi*(x-xi)^2+ci*(x-xi)+di 
% we deduce from si''(x)=6*ai*(x-xi)+2*bi that
% bi = Mi/2 and ai = (M(i+1)-Mi)/(6*hi)
coefs(:,2) = M(1:n-1)/2;
coefs(:,1) = (M(2:n) - M(1:n-1))./(6*h');
% evaluation gives di = yi
coefs(:,4) = f(1:n-1)';
% evaluation of s' gives formula for ci
coefs(:,3) = d - ((M(2:n) + 2*M(1:n-1))'/6).*h;
pp = PP(4,n,x,coefs);

endfunction % Matlab users must cut above this line for a function file.

%
% the following useful functions are not really specific to p.p. fcns
% 

function [y,x] = MaxNorm(fn_,a,b,n)
% usage: [y,x] = MaxNorm(fcn,a,b,n) or y=MaxNorm(fcn,a,b,n)
% description: returns the max norm y of the function f(x) on the
% interval a<=x<=b based on a uniform sample of n points and
% returns the first point x where the max occurs. Observe that the
% function name fcn is a string, so must be quoted, or use a handle.

% local variables:
% h: stepsize for sampling
% xc: current sample point
% fxc: abs of function value at current sample point
% j: index variable

xc = a;
x = a;
y = 0;
h = (b-a)/(n-1);
for j=1:n
  fxc = abs(feval(fn_,xc));
  if (y < fxc)
    y = fxc;
    x = xc;
  end;
  xc = xc + h;
end;

endfunction % Matlab users must cut above this line for a function file.

%

function retval = GaussInt(fn_,nodes,n)
% usage retval = GaussInt(fcn,nodes,n)
% description:  inputs a function of one variable fcn, ordered
% vector nodes and degree n, and returns an approximate value for
% the definite integral of fcn(x) on the interval enclosing the
% coordinates of nodes, calculated by using a n-point Gaussian
% quadrature method between adjacent nodes. The function fcn
% must handle vector arguments correctly: evaluate each coordinate.
% Only 2<=n<=8 are implemented.  This routine is useful for integrating 
% pp functions exactly. E.g., for cubic splines use spline nodes and n=2.
% Function name fcn is a string, so must be quoted, or use a handle.
% Degree n is optional and default is n=2.

% local variables:
% xabsc: Gaussian abscissas 
% t: translates of Gaussian abscissas and vals
% wts: Gaussian weights
% h: width of subinterval of Gauss integration
% k: index variables

if (nargin ~= 3)
  n = 2;
end
switch n %decide on size of interpolation
  case 2
    xabsc = [-0.577350269189626, 0.577350269189626];
    wts = [1.0, 1.0 ];
  case 3
    xabsc = [-0.774596669241483 0 0.774596669241483];
    wts = [0.555555555555556 0.888888888888889 0.555555555555556];
  case 4
    xabsc = [-0.861136311594052 -0.339981043584856 0.339981043584856 0.861136311594052]; 
wts = [0.347854845137454 0.652145154862546 0.652145154862546 0.347854845137454];
  case 5
    xabsc = [-0.9061798459386639928, -0.5384693101056830910, 0, 0.5384693101056830910, 0.9061798459386639928];      
    wts = [0.236926885056189, 0.4786286704993667, 0.5688888888888887, 0.4786286704993664, 0.2369268850561892];
  case 6
    xabsc = [-0.932469514203151 -0.661209386466264 -0.238619186083197 0.238619186083197 0.661209386466264 0.932469514203151];
    wts = [0.17132449237917 0.360761573048138 0.467913934572691 0.467913934572691 0.360761573048138 0.17132449237917];
  case 7
    xabsc = [-0.949107912342758 -0.741531185599394 -0.405845151377397 0 0.405845151377397 0.741531185599394 0.949107912342758];
    wts = [0.129484966168879 0.279705391489274 0.381830050505119 0.417959183673469 0.381830050505119 0.279705391489274 0.129484966168879];
  case 8
    xabsc = [-0.960289856497535 -0.796666477413627 -0.52553240991633 -0.183434642495649 0.183434642495649 0.52553240991633 0.796666477413627 0.960289856497535];
    wts = [0.101228536290382 0.22238103445337 0.313706645877889 0.362683783378362 0.362683783378362 0.313706645877889 0.22238103445337 0.101228536290382];
  otherwise
    error('This node number is not supported');
end
retval = 0;
nn = length(nodes);
for k = 1:nn-1
  t = nodes(k) + (nodes(k+1)-nodes(k))/2*(xabsc+1);
  t = feval(fn_,t);
  retval = retval + (nodes(k+1)-nodes(k))/2*wts*t';
end;

endfunction % Matlab users must cut above this line for a function file.