% Piecewise polynomial package
% programer: t. shores		date: 10/7/00
% 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
%	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 an 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.)
% 
% Matlab and Octave handle objects rather differently.  Matlab has a full
% featured OOP environment with classes, objects with properties and methods,
% inheritance, operator overloading, etc. Octave has only C-style structures.
% 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. Results would then 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)
  disp('Incorrect paramater list')
  return;
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)))
    disp('Incorrect knots vector')
    clear pp;
    return;
  end;
  for j = 2:knum
    if (knots(j) <= knots(j-1))
      disp('Knots vector is not strictly increasing')
      clear pp;
      return;
    end;
  end;
  pp.knots = knots;
end;
if (nargin == 4)
  if (sum([knum-1,order] ~=size(coefs)))
    disp('Incorrect coefs matrix')
    clear pp;
    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
  for k = 1:klim-1
    if (xtmp <  pp.knots(k+1))
      break;
    end;
  end;
  dx = xtmp - pp.knots(k);
  % use Horner's method to evaluate the polynomial
  y(j) = pp.coefs(k,1);
  for m = 2:ord
    y(j) = dx*y(j)+pp.coefs(k,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+1 where  
%     k = number of knots -1.

% 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
  disp('Size mismatch of abscissa and function ordinate');
  return;
end;
m = length(fp);
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and derivative ordinate');
  return;
end;
pp = PP(4, n, x);
% determine the poly coefs in each knot subinterval
for j = 1:(n-1)
  % set up and solve the system on this subinterval
  h = x(j+1) - x(j);
  pp.coefs(j,:) = ([0, 0, 0, 1; h^3, h^2, h, 1; 0 0 1 0;
3*h^2, 2*h, 1, 0] \ [f(j); f(j+1); fp(j); fp(j+1)] )';  
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
  disp('Size mismatch of abscissa and function ordinate');
  return;
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 Atkinson, p. 169, 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');
% evaulation gives di = yi and ci = (y(i+1) - ai*hi^3-bi*hi^2-di)/hi
coefs(:,4) = f(1:n-1)';
coefs(:,3) = (f(2:n)' - coefs(:,1).*h'.^3 - coefs(:,2).*h'.^2 - coefs(:,4)) ./ 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
  disp('Size mismatch of abscissa and function ordinate or size problem');
  return;
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 Atkinson, p. 169, row by row
Coef = zeros(n,n); 
Rhs = zeros(n,1);
% first row (use equation of exercise 3.34)
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 (use equation of exercise 3.34)
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');
% evaulation gives di = yi and ci = (y(i+1) - ai*hi^3-bi*hi^2-di)/hi
coefs(:,4) = f(1:n-1)';
coefs(:,3) = (f(2:n)' - coefs(:,1).*h'.^3 - coefs(:,2).*h'.^2 - coefs(:,4)) ./ 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
  disp('Size mismatch of abscissa and function ordinate');
  return;
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 Atkinson, p. 169, 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');
% evaulation gives di = yi and ci = (y(i+1) - ai*hi^3-bi*hi^2-di)/hi
coefs(:,4) = f(1:n-1)';
coefs(:,3) = (f(2:n)' - coefs(:,1).*h'.^3 - coefs(:,2).*h'.^2 - coefs(:,4)) ./ h';
pp = PP(4,n,x,coefs);

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

function pp = RNCSpp(x,a,alp)
% usage: pp = RNCSpp(x,a,alp)
% description: this function returns the regularized natural 
% cubic spline with nodes at the entries of the vector
% of abscissas x and corresponding ordinates in the 
% vector a. Here alp is  the regularization parameter.
% Given that the unregularized system is R*c=Q*a, where c is
% the column of moments (2nd deriv values), with 
% both R and Q tridiagonal and diagonally dominant, 
% the regularized system is 
%             (R + alp*Q*Q^{T}*R^{-1})*c = Q*a
% which is equivalent to the interpolating natural cubic 
% spline when alp=0. 

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

% columnize the input data
[m,n]=size(x);
if (n > 1)
    x = x';
end;
[m,n] = size(a);
if (n > 1)
    a = a';
end;
n = length(x);
m= length(a);
if (m ~= n) % error check for size and order
  disp('Size mismatch of abscissa and function ordinate');
  return;
end;
% construct the stepsize and divided difference vectors
h = x(2:n) - x(1:n-1);
% now set up the system as in my notes along lines of de Boor
R = zeros(n,n);
Q = R;
% first row
R(1,1) = 1;
Q(1,1) = 1;
for j = 2:(n-1) % interior rows
  R(j,(j-1):(j+1)) = [h(j-1), (h(j-1)+h(j))*2, h(j)];
  Q(j,(j-1):(j+1)) = [3/h(j-1), -3/h(j-1)-3/h(j), 3/h(j)];
end;
% last row
R(n,n) = 1;
Q(n,n) = 1;
% fix R to be symmetric
R(2,1) = 0;
R(n-1,n) = 0;
% solve for c
c = (R+alp*Q*Q'*inv(R))\(Q*a);
% find new nodes for spline
a = Q\(R*c);
% fix first and last entries, currently spline values at end nodes
c(1) = 0;
c(n) = 0;
% finally construct the PP structure to be returned
coefs = zeros(n-1,4);
% given si(x)=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 
% we deduce from si''(x)=6*di*(x-xi)+2*ci that
% ai = (ci(+1)-ci)/(3*hi)
coefs(:,2) = c(1:n-1);
coefs(:,1) = (c(2:n) - c(1:n-1))./(3*h);
% evaulation gives ai = yi and bi = (y(i+1) - di*hi^3-ci*hi^2-ai)/hi
coefs(:,4) = a(1:n-1);
coefs(:,3) = (a(2:n) - coefs(:,1).*h .^3 - coefs(:,2).*h .^2 - ...
    coefs(:,4)) ./ h;
pp = PP(4,n,x',coefs);

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

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

function [y,x] = MaxNorm(f_f,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.

% 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(f_f,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)
% description:  inputs a function of one variable fcn and ordered
% vector nodes 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=3.

% 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

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
    disp('This node number is not supported');
    return;
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.

%