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