Mercurial > hg > octave-nkf
annotate scripts/polynomial/pchip.m @ 9978:13a85d3e13bb
Fix typo and distribute new ppder.m script
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sun, 13 Dec 2009 19:35:24 -0800 |
| parents | 31900e17b5f5 |
| children | be55736a0783 |
| 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 ## |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
23 ## Piecewise Cubic Hermite interpolating polynomial. Called with two |
| 5837 | 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 | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
29 ## increasing or decreasing). While @var{y} can be either a vector or |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
30 ## array. In the case where @var{y} is a vector, it must have a length |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
31 ## of @var{n}. If @var{y} is an array, then the size of @var{y} must |
| 5837 | 32 ## have the form |
| 33 ## @tex | |
| 34 ## $$[s_1, s_2, \cdots, s_k, n]$$ | |
| 35 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
36 ## @ifnottex |
| 5837 | 37 ## @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
|
38 ## @end ifnottex |
| 7001 | 39 ## The array is then reshaped internally to a matrix where the leading |
| 5837 | 40 ## dimension is given by |
| 41 ## @tex | |
| 42 ## $$s_1 s_2 \cdots s_k$$ | |
| 43 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7650
diff
changeset
|
44 ## @ifnottex |
| 5837 | 45 ## @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
|
46 ## @end ifnottex |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
47 ## and each row in this matrix is then treated separately. Note that this |
| 5837 | 48 ## is exactly the opposite treatment than @code{interp1} and is done |
| 7001 | 49 ## for compatibility. |
| 5837 | 50 ## |
| 51 ## Called with a third input argument, @code{pchip} evaluates the | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
52 ## piece-wise polynomial at the points @var{xi}. There is an equivalence |
| 5837 | 53 ## between @code{ppval (pchip (@var{x}, @var{y}), @var{xi})} and |
| 54 ## @code{pchip (@var{x}, @var{y}, @var{xi})}. | |
| 55 ## | |
| 56 ## @seealso{spline, ppval, mkpp, unmkpp} | |
| 57 ## @end deftypefn | |
| 58 | |
| 59 ## Author: Kai Habel <kai.habel@gmx.de> | |
| 60 ## Date: 9. mar 2001 | |
| 61 ## | |
| 62 ## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom) | |
| 63 ## | |
| 64 ## 4 conditions: | |
| 65 ## S_k(x_k) = y_k; | |
| 66 ## S_k(x_k+1) = y_k+1; | |
| 67 ## S_k'(x_k) = y_k'; | |
| 68 ## S_k'(x_k+1) = y_k+1'; | |
| 69 | |
| 70 function ret = pchip (x, y, xi) | |
| 71 | |
| 72 if (nargin < 2 || nargin > 3) | |
| 73 print_usage (); | |
| 74 endif | |
| 75 | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
76 x = x(:).'; |
| 5837 | 77 n = length (x); |
| 78 | |
| 79 ## Check the size and shape of y | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
80 if (isvector (y)) |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
81 y = y(:).'; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
82 szy = size (y); |
| 5837 | 83 else |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
84 szy = size (y); |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
85 y = reshape (y, [prod(szy(1:end-1)), szy(end)]); |
| 5837 | 86 endif |
| 87 | |
| 5838 | 88 h = diff (x); |
| 89 if (all (h < 0)) | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
90 x = fliplr (x); |
| 5838 | 91 h = diff (x); |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
92 y = fliplr (y); |
| 5838 | 93 elseif (any (h <= 0)) |
| 8664 | 94 error("pchip: x must be strictly monotonic"); |
| 5837 | 95 endif |
| 96 | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
97 if (columns (y) != n) |
| 5837 | 98 error("pchip: size of x and y must match"); |
| 99 endif | |
| 100 | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
101 f1 = y(:,1:n-1); |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
102 |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
103 ## Compute derivatives. |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
104 d = __pchip_deriv__ (x, y, 2); |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
105 d1 = d(:,1:n-1); |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
106 d2 = d(:,2:n); |
| 5837 | 107 |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
108 ## This is taken from SLATEC. |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
109 h = diag (h); |
| 5837 | 110 |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
111 delta = diff (y, 1, 2) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
112 del1 = (d1 - delta) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
113 del2 = (d2 - delta) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
114 c3 = del1 + del2; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
115 c2 = -c3 - del1; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
116 c3 = c3 / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
117 |
|
9768
31900e17b5f5
improve Matlab compatibility & performance of ppval/mkpp and some associated funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
9754
diff
changeset
|
118 coeffs = cat (3, c3, c2, d1, f1); |
| 5837 | 119 pp = mkpp (x, coeffs, szy(1:end-1)); |
| 120 | |
| 121 if (nargin == 2) | |
| 122 ret = pp; | |
| 123 else | |
| 124 ret = ppval (pp, xi); | |
| 125 endif | |
| 126 | |
| 127 endfunction | |
| 128 | |
| 129 %!demo | |
| 130 %! x = 0:8; | |
| 131 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; | |
| 132 %! xi = 0:0.01:8; | |
| 133 %! yspline = spline(x,y,xi); | |
| 134 %! ypchip = pchip(x,y,xi); | |
| 135 %! title("pchip and spline fit to discontinuous function"); | |
| 6746 | 136 %! plot(xi,yspline,xi,ypchip,"-",x,y,"+"); |
| 137 %! legend ("spline","pchip","data"); | |
| 5837 | 138 %! %------------------------------------------------------------------- |
| 139 %! % confirm that pchip agreed better to discontinuous data than spline | |
| 140 | |
| 141 %!shared x,y | |
| 142 %! x = 0:8; | |
| 143 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; | |
| 144 %!assert (pchip(x,y,x), y); | |
| 145 %!assert (pchip(x,y,x'), y'); | |
| 146 %!assert (pchip(x',y',x'), y'); | |
| 147 %!assert (pchip(x',y',x), y); | |
| 148 %!assert (isempty(pchip(x',y',[]))); | |
| 149 %!assert (isempty(pchip(x,y,[]))); | |
| 150 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)]) |