Mercurial > hg > octave-nkf
annotate scripts/polynomial/pchip.m @ 13767:2b98014771b4
fileread.m: Add functional test.
* fileread.m: Add functional test.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Thu, 27 Oct 2011 22:17:03 -0700 |
| parents | 59e2460acae1 |
| children | 614505385171 |
| rev | line source |
|---|---|
| 11523 | 1 ## Copyright (C) 2001-2011 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 -*- | |
|
10793
be55736a0783
Grammarcheck the documentation from m-files.
Rik <octave@nomad.inbox5.com>
parents:
9768
diff
changeset
|
20 ## @deftypefn {Function File} {@var{pp} =} pchip (@var{x}, @var{y}) |
| 7650 | 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 |
|
11536
702dbd0c53f5
Add undocumented ppder, ppint, ppjumps functions to documentation.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
24 ## arguments, the piecewise polynomial @var{pp} is returned, that may |
| 5837 | 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 |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
30 ## an array. In the case where @var{y} is a vector, it must have the |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
31 ## length @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 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
40 ## dimension is given by |
| 5837 | 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 ## |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
51 ## Called with a third input argument, @code{pchip} evaluates the |
|
11536
702dbd0c53f5
Add undocumented ppder, ppint, ppjumps functions to documentation.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
52 ## piecewise 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 | |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
76 ## make row vector |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
77 x = x(:).'; |
| 5837 | 78 n = length (x); |
| 79 | |
| 80 ## Check the size and shape of y | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
81 if (isvector (y)) |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
82 y = y(:).'; ##row vector |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
83 szy = size (y); |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
84 if !(size_equal (x, y)) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
85 error ("pchip: length of X and Y must match") |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
86 endif |
| 5837 | 87 else |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
88 szy = size (y); |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
89 if (n != szy(end)) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
90 error ("pchip: length of X and last dimension of Y must match") |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
91 endif |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
92 y = reshape (y, [prod(szy(1:end-1)), szy(end)]); |
| 5837 | 93 endif |
| 94 | |
| 5838 | 95 h = diff (x); |
| 96 if (all (h < 0)) | |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
97 x = fliplr (x); |
| 5838 | 98 h = diff (x); |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
99 y = fliplr (y); |
| 5838 | 100 elseif (any (h <= 0)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
10793
diff
changeset
|
101 error("pchip: X must be strictly monotonic"); |
| 5837 | 102 endif |
| 103 | |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
104 f1 = y(:, 1:n-1); |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
105 |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
106 ## Compute derivatives. |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
107 d = __pchip_deriv__ (x, y, 2); |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
108 d1 = d(:, 1:n-1); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
109 d2 = d(:, 2:n); |
| 5837 | 110 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
111 ## This is taken from SLATEC. |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
112 h = diag (h); |
| 5837 | 113 |
|
9754
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
114 delta = diff (y, 1, 2) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
115 del1 = (d1 - delta) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
116 del2 = (d2 - delta) / h; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
117 c3 = del1 + del2; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
118 c2 = -c3 - del1; |
|
4219e5cf773d
improve interp1 and pchip
Jaroslav Hajek <highegg@gmail.com>
parents:
9211
diff
changeset
|
119 c3 = c3 / h; |
|
9768
31900e17b5f5
improve Matlab compatibility & performance of ppval/mkpp and some associated funcs
Jaroslav Hajek <highegg@gmail.com>
parents:
9754
diff
changeset
|
120 coeffs = cat (3, c3, c2, d1, f1); |
| 5837 | 121 |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
122 ret = mkpp (x, coeffs, szy(1:end-1)); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
123 |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
124 if (nargin == 3) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
125 ret = ppval (ret, xi); |
| 5837 | 126 endif |
| 127 | |
| 128 endfunction | |
| 129 | |
| 130 %!demo | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
131 %! x = 0:8; |
| 5837 | 132 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
133 %! xi = 0:0.01:8; |
| 5837 | 134 %! yspline = spline(x,y,xi); |
| 135 %! ypchip = pchip(x,y,xi); | |
| 136 %! title("pchip and spline fit to discontinuous function"); | |
| 6746 | 137 %! plot(xi,yspline,xi,ypchip,"-",x,y,"+"); |
| 138 %! legend ("spline","pchip","data"); | |
| 5837 | 139 %! %------------------------------------------------------------------- |
| 140 %! % confirm that pchip agreed better to discontinuous data than spline | |
| 141 | |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
142 %!shared x,y,y2,pp,yi1,yi2,yi3 |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11536
diff
changeset
|
143 %! x = 0:8; |
| 5837 | 144 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; |
| 145 %!assert (pchip(x,y,x), y); | |
| 146 %!assert (pchip(x,y,x'), y'); | |
| 147 %!assert (pchip(x',y',x'), y'); | |
| 148 %!assert (pchip(x',y',x), y); | |
| 149 %!assert (isempty(pchip(x',y',[]))); | |
| 150 %!assert (isempty(pchip(x,y,[]))); | |
| 151 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)]) | |
|
12608
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
152 %!assert (pchip(x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)]) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
153 %!assert (pchip(x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)]) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
154 %!assert (pchip(x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)]) |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
155 %!test |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
156 %! x=(0:8)*pi/4;y=[sin(x);cos(x)]; |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
157 %! y2(:,:,1)=y;y2(:,:,2)=y+1;y2(:,:,3)=y-1; |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
158 %! pp=pchip(x,shiftdim(y2,2)); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
159 %! yi1=ppval(pp,(1:4)*pi/4); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
160 %! yi2=ppval(pp,repmat((1:4)*pi/4,[5,1])); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
161 %! yi3=ppval(pp,[pi/2,pi]); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
162 %!assert(size(pp.coefs),[48,4]); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
163 %!assert(pp.pieces,8); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
164 %!assert(pp.order,4); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
165 %!assert(pp.dim,[3,2]); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
166 %!assert(ppval(pp,pi),[0,-1;1,0;-1,-2],1e-14); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
167 %!assert(yi3(:,:,2),ppval(pp,pi),1e-14); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
168 %!assert(yi3(:,:,1),[1,0;2,1;0,-1],1e-14); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
169 %!assert(squeeze(yi1(1,2,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
170 %!assert(size(yi2),[3,2,5,4]); |
|
59e2460acae1
make piecewise polynomial (pp) functions more compatible
Kai Habel <kai.habel@gmx.de>
parents:
11587
diff
changeset
|
171 %!assert(squeeze(yi2(1,2,3,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14); |