octave_packages/m/polynomial/pchip.m

   1 ## Copyright (C) 2001-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{pp} =} pchip (@var{x}, @var{y})
  21 ## @deftypefnx {Function File} {@var{yi} =} pchip (@var{x}, @var{y}, @var{xi})
  22 ## Return the Piecewise Cubic Hermite Interpolating Polynomial (pchip) of
  23 ## points @var{x} and @var{y}.
  24 ##
  25 ## If called with two arguments, return the piecewise polynomial @var{pp}
  26 ## that may be used with @code{ppval} to evaluate the polynomial at specific
  27 ## points.  When called with a third input argument, @code{pchip} evaluates
  28 ## the pchip polynomial at the points @var{xi}.  The third calling form is
  29 ## equivalent to @code{ppval (pchip (@var{x}, @var{y}), @var{xi})}.
  30 ##
  31 ## The variable @var{x} must be a strictly monotonic vector (either
  32 ## increasing or decreasing) of length @var{n}.  @var{y} can be either a
  33 ## vector or array.  If @var{y} is a vector then it must be the same length
  34 ## @var{n} as @var{x}.  If @var{y} is an array then the size of @var{y} must
  35 ## have the form
  36 ## @tex
  37 ## $$[s_1, s_2, \cdots, s_k, n]$$
  38 ## @end tex
  39 ## @ifnottex
  40 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
  41 ## @end ifnottex
  42 ## The array is reshaped internally to a matrix where the leading
  43 ## dimension is given by
  44 ## @tex
  45 ## $$s_1 s_2 \cdots s_k$$
  46 ## @end tex
  47 ## @ifnottex
  48 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
  49 ## @end ifnottex
  50 ## and each row of this matrix is then treated separately.  Note that this
  51 ## is exactly opposite to @code{interp1} but is done for @sc{matlab}
  52 ## compatibility.
  53 ##
  54 ## @seealso{spline, ppval, mkpp, unmkpp}
  55 ## @end deftypefn
  56
  57 ## Author:  Kai Habel <kai.habel@gmx.de>
  58 ## Date: 9. mar 2001
  59 ##
  60 ## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom)
  61 ##
  62 ## 4 conditions:
  63 ## S_k(x_k) = y_k;
  64 ## S_k(x_k+1) = y_k+1;
  65 ## S_k'(x_k) = y_k';
  66 ## S_k'(x_k+1) = y_k+1';
  67
  68 function ret = pchip (x, y, xi)
  69
  70   if (nargin < 2 || nargin > 3)
  71     print_usage ();
  72   endif
  73
  74   ## make row vector
  75   x = x(:).';
  76   n = length (x);
  77
  78   ## Check the size and shape of y
  79   if (isvector (y))
  80     y = y(:).'; ##row vector
  81     szy = size (y);
  82     if !(size_equal (x, y))
  83       error ("pchip: length of X and Y must match")
  84     endif
  85   else
  86     szy = size (y);
  87     if (n != szy(end))
  88       error ("pchip: length of X and last dimension of Y must match")
  89     endif
  90     y = reshape (y, [prod(szy(1:end-1)), szy(end)]);
  91   endif
  92
  93   h = diff (x);
  94   if (all (h < 0))
  95     x = fliplr (x);
  96     h = diff (x);
  97     y = fliplr (y);
  98   elseif (any (h <= 0))
  99     error("pchip: X must be strictly monotonic");
 100   endif
 101
 102   f1 = y(:, 1:n-1);
 103
 104   ## Compute derivatives.
 105   d = __pchip_deriv__ (x, y, 2);
 106   d1 = d(:, 1:n-1);
 107   d2 = d(:, 2:n);
 108
 109   ## This is taken from SLATEC.
 110   h = diag (h);
 111
 112   delta = diff (y, 1, 2) / h;
 113   del1 = (d1 - delta) / h;
 114   del2 = (d2 - delta) / h;
 115   c3 = del1 + del2;
 116   c2 = -c3 - del1;
 117   c3 = c3 / h;
 118   coeffs = cat (3, c3, c2, d1, f1);
 119
 120   ret = mkpp (x, coeffs, szy(1:end-1));
 121
 122   if (nargin == 3)
 123     ret = ppval (ret, xi);
 124   endif
 125
 126 endfunction
 127
 128 %!demo
 129 %! x = 0:8;
 130 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
 131 %! xi = 0:0.01:8;
 132 %! yspline = spline(x,y,xi);
 133 %! ypchip = pchip(x,y,xi);
 134 %! title("pchip and spline fit to discontinuous function");
 135 %! plot(xi,yspline,xi,ypchip,"-",x,y,"+");
 136 %! legend ("spline","pchip","data");
 137 %! %-------------------------------------------------------------------
 138 %! % confirm that pchip agreed better to discontinuous data than spline
 139
 140 %!shared x,y,y2,pp,yi1,yi2,yi3
 141 %! x = 0:8;
 142 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
 143 %!assert (pchip(x,y,x), y);
 144 %!assert (pchip(x,y,x'), y');
 145 %!assert (pchip(x',y',x'), y');
 146 %!assert (pchip(x',y',x), y);
 147 %!assert (isempty(pchip(x',y',[])));
 148 %!assert (isempty(pchip(x,y,[])));
 149 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
 150 %!assert (pchip(x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
 151 %!assert (pchip(x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
 152 %!assert (pchip(x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
 153 %!test
 154 %! x=(0:8)*pi/4;y=[sin(x);cos(x)];
 155 %! y2(:,:,1)=y;y2(:,:,2)=y+1;y2(:,:,3)=y-1;
 156 %! pp=pchip(x,shiftdim(y2,2));
 157 %! yi1=ppval(pp,(1:4)*pi/4);
 158 %! yi2=ppval(pp,repmat((1:4)*pi/4,[5,1]));
 159 %! yi3=ppval(pp,[pi/2,pi]);
 160 %!assert(size(pp.coefs),[48,4]);
 161 %!assert(pp.pieces,8);
 162 %!assert(pp.order,4);
 163 %!assert(pp.dim,[3,2]);
 164 %!assert(ppval(pp,pi),[0,-1;1,0;-1,-2],1e-14);
 165 %!assert(yi3(:,:,2),ppval(pp,pi),1e-14);
 166 %!assert(yi3(:,:,1),[1,0;2,1;0,-1],1e-14);
 167 %!assert(squeeze(yi1(1,2,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
 168 %!assert(size(yi2),[3,2,5,4]);
 169 %!assert(squeeze(yi2(1,2,3,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
 170
 171 %!error (pchip (1,2));
 172 %!error (pchip (1,2,3));