spline, mkpp and ppval

[Mats Jansson]

I wrote the functions spline.m, mkpp.m and ppval.m. I don't know exactly
how Matlab's functions with the same names behave, please tell me if you
find any differences.

/Mats


spline.m
========

## Copyright (C) 1996, 1997 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden =} spline (@var{x}, @var{y})
## @deftypefnx {Function File} address@hidden =} spline (@var{x}, @var{y}, 
@var{xi})
## Piece-wise interpolation with cubic splines.
##
## The argument @var{x} should be a vector and @var{y} should either be
## a vector of the same length as @var{x} or a matrix with length
## (@var{x}) number of columns. The spline-function uses cubic spline
## interpolation to find @var{yi} corresponding to @var{xi}. If the
## third argument is omitted, spline returns a piece-wise polynom on
## pp-form.
##
## Example:
## 
## @example
## x  = -1:0.25:1;
## xi = -1:0.01:1;
## y  = [exp(x)./(1 + 9*x.^2); sin(4*x)];
## y1 = [exp(xi)./(1 + 9*xi.^2); sin(4*xi)];
## yi = spline (x, y, xi);
## plot (xi, y1, ["g-;calculated;";"g-;;"], \
##       xi, yi, ["b-;interpolated;";"b-;;"], \
##       x, y, ["ro;interpolation-points;";"ro;;"])
## @end example
## @end deftypefn
## @seealso{ppval}

## Author: Mats Jansson <address@hidden>
## Created: March 5, 2000 - April 6, 2001
## Keywords: piecewise interpolation cubic spline

function yi = spline (x, y, xi)

  if (nargin < 2),
    usage ("spline: expecting 2 or 3 arguments.");
  endif
  
  x = x(:); ## Make sure x is a column vector
  n = size (x, 1);
  if (n < 2),
    usage ("yi = spline (x, y, xi): x should have at least 2 elements.");
  endif
  
  if (is_vector (y)),
    y = y(:);    ## Make sure y is a columnvector
    [nry, ncy] = size (y);
    if (nry != n),
      error ("yi = spline (x, y, xi): x and y should have the same number of 
elements.");
    endif
  elseif (is_matrix (y)),
    [nry, ncy] = size (y);
    if (ncy != n),
      ## The number of columns of y should match the length of x (I
      ## think this is the way Matlabs' spline-function behaves).
      error ("yi = spline (x, y, xi): The number of columns of y should match 
the length of x.");
    else
      y = y.';
      [nry, ncy] = size (y);
    endif
  endif
  
  
  [x, ind] = sort (x);
  if (! all (diff (x))),
    usage ("yi = spline (x, y, xi): x should be distinct.");
  endif
  
  y = y(ind,:);
  
  h = diff (x);
  d = diff (y);
  
  if (n == 2),
    ## Linear interpolation
    coef = [d./h(:,ones (1, ncy)); y(1,:)].';
  else
    ## We want the differentiate and the second differentiate to be
    ## continuous in each point. To find the differentiate, k, in each
    ## point, we have to solve the equation system:
    ##
    ## H * k = 3b, where H is the n x n tri-diagonal matrix
    ##   _                                                           _
    ##  |  2h(1)       h(1)          0         0      ...        0    |
    ##  |   h(2)   2h((2)+h(1))     h(1)       0      ...        0    |
    ##  |    0         h(3)     2(h(3)+h(2))  h(2)    ...        0    |
    ##  |    .          .               .                        .    |
    ##  |    .          .                  .                     .    |
    ##  |    .          .         h(n-1)  2(h(n-1)+h(n-2)    h(n-2)   |
    ##  |_   0          0           0          h(n-1)       2h(n-1)  _|
    ##
    ## and b is the n-by-ncy matrix
    ##   _                                              _
    ##  |                     d(1,:)                     |
    ##  |        h(1)/h(2)*d(2,:)+h(2)/h(1)*d(1,:)       |
    ##  |        h(2)/h(3)*d(3,:)+h(3)/h(2)*d(2,:)       |
    ##  |                      ...                       |
    ##  |  h(n-2)/h(n-1)*d(n-1,:)+h(n-1)/h(n-2)*d(n-2,:) |
    ##  |_                   d(n-1,:)                   _|
    ##
    
    H = diag (2 * [h(1); h(2:n-1)+h(1:n-2); h(n-1)]) + \
        diag ([h(1); h(1:n-2)], 1) + \
        diag ([h(2:n-1); h(n-1)], -1);

    b = zeros (n, ncy);
    
    ## Repeat the column-vector h ncy times to form a matrix of the
    ## same size as d.
    h = h(:,ones (1, ncy));

    b(1,:) = d(1,:);
    b(n,:) = d(n-1,:);
    b(2:n-1,:) = h(1:n-2,:)./h(2:n-1,:).*d(2:n-1,:) + \
        h(2:n-1,:)./h(1:n-2,:).*d(1:n-2,:);
    
    ## H is a sparse matrix, this system could be solved more efficient.
    k = H \ (3*b);
    
    ## We should have n-1 3:rd degree interpolation polynoms for each
    ## column of y.
    coef = \
        [((h(1:n-1,:).*(k(1:n-1,:)+k(2:n,:))-2*d(1:n-1,:)) ./ \
          h(1:n-1,:).^3).'(:); \
         ((3*d(1:n-1,:)-2*h(1:n-1,:).*k(1:n-1,:)-h(1:n-1,:).*k(2:n,:)) ./ \
          h(1:n-1,:).^2).'(:); \
         k(1:n-1,:).'(:); \
         y(1:n-1,:).'(:)];
  endif
  
  yi = mkpp (x, coef, ncy);

  if (nargin >= 3),
    yi = ppval (yi, xi);
  endif
endfunction


mkpp.m
======

## Copyright (C) 1996, 1997 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden =} mkpp (@var{breaks}, 
@var{coefficients})
## @deftypefnx {Function File} address@hidden =} mkpp (@var{breaks}, 
@var{coefficients}, @var{number_of_sets})
## Define piecewise polynomials.
## @end deftypefn
## @seealso{ppval and spline}

## Author: Mats Jansson <address@hidden>
## Created: April 4, 2001 - April 6, 2001
## Keywords: piecewise polynomial

function pp = mkpp (b, c, s)

  ## b: breakpoints
  ## c: polynomial coefficients
  ## s: number of sets

  if (nargin == 2),
    s = 1;
  elseif (nargin < 2 || nargin > 3),
    usage ("mkpp (breaks, coefs) or mkpp (breaks, coefs, number_of_sets)");
  endif
  s = fix (s);
  b = b(:);
  ## Assume n breakpoints (b-values), n-1 polynomials of degree p for
  ## each set. 
  ##
  ## c(1): The highest coefficient of the polynomial valid in 
  ##       b(1) < x <= b(2) for the first set.
  ## c(2): The highest coefficient of the polynomial valid in
  ##       b(1) < x <= b(2) for the second set.
  ## ...
  ##
  ## c(s): The highest coefficient of the polynomial valid in
  ##       b(1) < x <= b(2) for the set number s.
  ## c(s+1): 
  ##       The highest coefficient of the polynomial valid in
  ##       b(2) < x <= b(3) for the first set.
  ## c(s+2):
  ##       The highest coefficient of the polynomial valid in
  ##       b(2) < x <= b(3) for the second set.
  ## ...
  ## c(2*s):
  ##       The highest coefficient of the polynomial valid in
  ##       b(2) < x <= b(3) for set number s.
  ## ...
  ## c((n-1)*s):
  ##       The highest coefficient of the polynomial valid in
  ##       b(n-1) < x <= b(n) for set number s.
  ## c((n-1)*s+1):
  ##       The second highest coefficient of the polynomial valid in
  ##       b(1) < x <= b(2)
  ## ...
  ## c(2*(n-1)*s):
  ##       The second highest coefficient of the polynomial valid in
  ##       b(n-1) < x <= b(n)
  ## ...
  ## c((p+1)*(n-1)*s):
  ##       The last (lowest) coefficient of the polynomial valid in
  ##       b(n-1) < x <= b(n)
  c = c(:);
  n = length (b);
  p = fix (length (c)/(n-1)/s + 100*eps) - 1;
  
  pp.breaks = b;
  pp.coefs = reshape (c, (n-1)*s, p+1);
  pp.degree = p;
  pp.sets = s;
endfunction


ppval.m
=======

## Copyright (C) 1996, 1997 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden =} ppval (@var{pp}, @var{x})
## Return the values of the piecewise polynomials @var{pp} at the points
## @var{x}.
##
## If @var{pp} contains one set of piecewise polynomials, the result has
## the same size as @var{x}. If @var{pp} contains more than one set of
## piecewise polynomials, each row of the result represents one set. In
## this case each column contains the values of @var{pp} in the points
## @var{x}, were @var{x} is taken in column-major order.
## @end deftypefn
## @seealso{mkpp and spline}

## Author: Mats Jansson <address@hidden>
## Created: April 4, 2001 - April 6, 2001
## Keywords: piecewise polynomial

function y = ppval (pp, x)
  if (nargin != 2),
    usage ("Y = ppval (PP, X)");
  endif
  if (!isnumeric (x)),
    error ("ppval (PP, X): Expecting X to be numeric.");
  endif
  if (!is_struct (pp)),
    error ("ppval (PP, X): Expecting PP to be a structure.");
  endif

  if (!struct_contains (pp, "breaks") || \
      !struct_contains (pp, "coefs") || \
      !struct_contains (pp, "degree") || \
      !struct_contains (pp, "sets")),
    error ("ppval (PP, X): PP doesn't seem to be a piece-wise polynomial 
created with mkpp.");
  endif

  [nrb, ncb] = size (pp.breaks);
  if (ncb != 1 || nrb < 2),
    error ("ppval (PP, X): Expecting PP.breaks to be a column-vector.");
  endif
  [nrx, ncx] = size (x);
  if (nrx*ncx == 0),
    y = x;
    return
  endif
  x = x(:).'; ## x is now a row-vector
  [x, indx] = sort (x);

  ## How many x-values do we have in each interval?
  if (nrb == 2),
    ## Only one interval
    number_of_x_in_interval_n = nrx*ncx;
  else
    xx = x(ones (1, nrb-1),:);
    bb = pp.breaks(1:nrb-1,ones (1, nrx*ncx));
    bb2 = [bb(2:nrb-1,:); Inf*ones(1, nrx*ncx)];
    bb  = [-Inf*ones(1, nrx*ncx); bb(2:nrb-1,:)];
    if (nrx*ncx == 1),
      number_of_x_in_interval_n = (xx >= bb & xx < bb2)';
    else
      number_of_x_in_interval_n = sum ((xx >= bb & xx < bb2)');
    endif
  endif

  y = zeros (pp.sets, nrx*ncx);
  pow = (pp.degree:-1:0).';
  ind = find (number_of_x_in_interval_n);
  x_pointer = 0;

  ## Step through the intervals with x-values.
  for nn = 1:length (ind)
    n = ind(nn); ## Interval number
    
    ## The coefficients valid in this interval. The first column is the
    ## coefficients for the first set, the second column is the
    ## coefficients for the second set and so on. The first row is the
    ## coefficients for the highest power, the second row is the
    ## coefficients for the second highest power...
    c = pp.coefs((n-1)*pp.sets+1:n*pp.sets,:).';
    
    ## Repeat the coefficients so that we have one set of coefficients
    ## for each x-value in the interval.
    cind = (1:pp.sets)';
    cind = cind(:, ones (1, number_of_x_in_interval_n(n)));
    c = c((1:pp.degree+1)', cind);

    ## A row-vector with the x-values in the interval (local coordinate
    ## system)
    xx = x(x_pointer+1:x_pointer+number_of_x_in_interval_n(n)) - \
        pp.breaks(n);
    
    ## Repeat each x-value as many times as we have sets
    xx = xx(ones (1, pp.sets),:);
    xx = xx(:).';
    
    ## Repeat the row-vector of x so that we have one row for each
    ## coefficient.
    xx = xx(ones (1, pp.degree+1),:);
    
    if (pp.degree == 0),
      yn = c;
    else
      yn = sum (c.*xx.^pow(:, \
                           ones (1, \
                                 number_of_x_in_interval_n(n)*pp.sets)));
    endif

    yn = reshape (yn, pp.sets, number_of_x_in_interval_n(n));
    y(:, x_pointer+1:x_pointer+number_of_x_in_interval_n(n)) = yn;

    x_pointer += number_of_x_in_interval_n(n);
  endfor

  y(:,indx) = y;

  if (pp.sets == 1 && nrx > 1),
    ## y has the same number of elements as x. Return y with the
    ## same number of rows and columns.
    y = reshape (y, nrx, ncx);
  endif

endfunction



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------