Mercurial > hg > octave-nkf
annotate scripts/polynomial/pchip.m @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | e07e93c04080 |
| children | 1bf0ce0930be |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2001, 2002, 2006, 2007, 2008, 2009 Kai Habel |
| 5837 | 2 ## |
| 6440 | 3 ## This file is part of Octave. |
| 5837 | 4 ## |
| 6440 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 6440 | 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. | |
| 5837 | 14 ## |
| 15 ## You should have received a copy of the GNU General Public License | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5837 | 18 |
| 19 ## -*- texinfo -*- | |
| 7650 | 20 ## @deftypefn {Function File} {@var{pp} =} pchip (@var{x}, @var{y}) |
| 21 ## @deftypefnx {Function File} {@var{yi} =} pchip (@var{x}, @var{y}, @var{xi}) | |
| 5837 | 22 ## |
| 23 ## Piecewise Cubic Hermite interpolating polynomial. Called with two | |
| 24 ## arguments, the piece-wise polynomial @var{pp} is returned, that may | |
| 25 ## later be used with @code{ppval} to evaluate the polynomial at | |
| 26 ## specific points. | |
| 27 ## | |
| 28 ## The variable @var{x} must be a strictly monotonic vector (either | |
| 29 ## increasing or decreasing). While @var{y} can be either a vector or | |
| 30 ## array. In the case where @var{y} is a vector, it must have a length | |
| 31 ## of @var{n}. If @var{y} is an array, then the size of @var{y} must | |
| 32 ## have the form | |
| 33 ## @iftex | |
| 34 ## @tex | |
| 35 ## $$[s_1, s_2, \cdots, s_k, n]$$ | |
| 36 ## @end tex | |
| 37 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
38 ## @ifnottex |
| 5837 | 39 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
40 ## @end ifnottex |
| 7001 | 41 ## The array is then reshaped internally to a matrix where the leading |
| 5837 | 42 ## dimension is given by |
| 43 ## @iftex | |
| 44 ## @tex | |
| 45 ## $$s_1 s_2 \cdots s_k$$ | |
| 46 ## @end tex | |
| 47 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
48 ## @ifnottex |
| 5837 | 49 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
50 ## @end ifnottex |
| 7001 | 51 ## and each row in this matrix is then treated separately. Note that this |
| 5837 | 52 ## is exactly the opposite treatment than @code{interp1} and is done |
| 7001 | 53 ## for compatibility. |
| 5837 | 54 ## |
| 55 ## Called with a third input argument, @code{pchip} evaluates the | |
| 56 ## piece-wise polynomial at the points @var{xi}. There is an equivalence | |
| 57 ## between @code{ppval (pchip (@var{x}, @var{y}), @var{xi})} and | |
| 58 ## @code{pchip (@var{x}, @var{y}, @var{xi})}. | |
| 59 ## | |
| 60 ## @seealso{spline, ppval, mkpp, unmkpp} | |
| 61 ## @end deftypefn | |
| 62 | |
| 63 ## Author: Kai Habel <kai.habel@gmx.de> | |
| 64 ## Date: 9. mar 2001 | |
| 65 ## | |
| 66 ## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom) | |
| 67 ## | |
| 68 ## 4 conditions: | |
| 69 ## S_k(x_k) = y_k; | |
| 70 ## S_k(x_k+1) = y_k+1; | |
| 71 ## S_k'(x_k) = y_k'; | |
| 72 ## S_k'(x_k+1) = y_k+1'; | |
| 73 | |
| 74 function ret = pchip (x, y, xi) | |
| 75 | |
| 76 if (nargin < 2 || nargin > 3) | |
| 77 print_usage (); | |
| 78 endif | |
| 79 | |
| 80 x = x(:); | |
| 81 n = length (x); | |
| 82 | |
| 83 ## Check the size and shape of y | |
| 84 ndy = ndims (y); | |
| 85 szy = size (y); | |
| 86 if (ndy == 2 && (szy(1) == 1 || szy(2) == 1)) | |
| 87 if (szy(1) == 1) | |
| 88 y = y'; | |
| 89 else | |
| 90 szy = fliplr (szy); | |
| 91 endif | |
| 92 else | |
| 93 y = reshape (y, [prod(szy(1:end-1)), szy(end)])'; | |
| 94 endif | |
| 95 | |
| 5838 | 96 h = diff (x); |
| 97 if (all (h < 0)) | |
| 98 x = flipud (x); | |
| 99 h = diff (x); | |
| 100 y = flipud (y); | |
| 101 elseif (any (h <= 0)) | |
| 8664 | 102 error("pchip: x must be strictly monotonic"); |
| 5837 | 103 endif |
| 104 | |
| 5838 | 105 if (rows (y) != n) |
| 5837 | 106 error("pchip: size of x and y must match"); |
| 107 endif | |
| 108 | |
| 109 [ry, cy] = size (y); | |
| 110 if (cy > 1) | |
| 111 h = kron (diff (x), ones (1, cy)); | |
| 112 endif | |
| 113 | |
| 114 dy = diff (y) ./ h; | |
| 115 | |
| 116 a = y; | |
| 5838 | 117 b = __pchip_deriv__ (x, y); |
| 5837 | 118 c = - (b(2:n, :) + 2 * b(1:n - 1, :)) ./ h + 3 * diff (a) ./ h .^ 2; |
| 119 d = (b(1:n - 1, :) + b(2:n, :)) ./ h.^2 - 2 * diff (a) ./ h.^3; | |
| 120 | |
| 121 d = d(1:n - 1, :); c = c(1:n - 1, :); | |
| 122 b = b(1:n - 1, :); a = a(1:n - 1, :); | |
| 123 coeffs = [d(:), c(:), b(:), a(:)]; | |
| 124 pp = mkpp (x, coeffs, szy(1:end-1)); | |
| 125 | |
| 126 if (nargin == 2) | |
| 127 ret = pp; | |
| 128 else | |
| 129 ret = ppval (pp, xi); | |
| 130 endif | |
| 131 | |
| 132 endfunction | |
| 133 | |
| 134 %!demo | |
| 135 %! x = 0:8; | |
| 136 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; | |
| 137 %! xi = 0:0.01:8; | |
| 138 %! yspline = spline(x,y,xi); | |
| 139 %! ypchip = pchip(x,y,xi); | |
| 140 %! title("pchip and spline fit to discontinuous function"); | |
| 6746 | 141 %! plot(xi,yspline,xi,ypchip,"-",x,y,"+"); |
| 142 %! legend ("spline","pchip","data"); | |
| 5837 | 143 %! %------------------------------------------------------------------- |
| 144 %! % confirm that pchip agreed better to discontinuous data than spline | |
| 145 | |
| 146 %!shared x,y | |
| 147 %! x = 0:8; | |
| 148 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; | |
| 149 %!assert (pchip(x,y,x), y); | |
| 150 %!assert (pchip(x,y,x'), y'); | |
| 151 %!assert (pchip(x',y',x'), y'); | |
| 152 %!assert (pchip(x',y',x), y); | |
| 153 %!assert (isempty(pchip(x',y',[]))); | |
| 154 %!assert (isempty(pchip(x,y,[]))); | |
| 155 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)]) |