|
5837
|
1 ## Copyright (C) 2001,2002 Kai Habel |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## -*- texinfo -*- |
|
|
18 ## @deftypefn {Function File} {@var{pp} = } pchip (@var{x}, @var{y}) |
|
|
19 ## @deftypefnx {Function File} {@var{yi} = } pchip (@var{x}, @var{y}, @var{xi}) |
|
|
20 ## |
|
|
21 ## Piecewise Cubic Hermite interpolating polynomial. Called with two |
|
|
22 ## arguments, the piece-wise polynomial @var{pp} is returned, that may |
|
|
23 ## later be used with @code{ppval} to evaluate the polynomial at |
|
|
24 ## specific points. |
|
|
25 ## |
|
|
26 ## The variable @var{x} must be a strictly monotonic vector (either |
|
|
27 ## increasing or decreasing). While @var{y} can be either a vector or |
|
|
28 ## array. In the case where @var{y} is a vector, it must have a length |
|
|
29 ## of @var{n}. If @var{y} is an array, then the size of @var{y} must |
|
|
30 ## have the form |
|
|
31 ## @iftex |
|
|
32 ## @tex |
|
|
33 ## $$[s_1, s_2, \cdots, s_k, n]$$ |
|
|
34 ## @end tex |
|
|
35 ## @end iftex |
|
|
36 ## @ifinfo |
|
|
37 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} |
|
|
38 ## @end ifinfo |
|
|
39 ## The array is then reshaped internally to a matrix where to leading |
|
|
40 ## dimension is given by |
|
|
41 ## @iftex |
|
|
42 ## @tex |
|
|
43 ## $$s_1 s_2 \cdots s_k$$ |
|
|
44 ## @end tex |
|
|
45 ## @end iftex |
|
|
46 ## @ifinfo |
|
|
47 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} |
|
|
48 ## @end ifinfo |
|
|
49 ## and each row this matrix is then treated seperately. Note that this |
|
|
50 ## is exactly the opposite treatment than @code{interp1} and is done |
|
|
51 ## for compatiability. |
|
|
52 ## |
|
|
53 ## Called with a third input argument, @code{pchip} evaluates the |
|
|
54 ## piece-wise polynomial at the points @var{xi}. There is an equivalence |
|
|
55 ## between @code{ppval (pchip (@var{x}, @var{y}), @var{xi})} and |
|
|
56 ## @code{pchip (@var{x}, @var{y}, @var{xi})}. |
|
|
57 ## |
|
|
58 ## @seealso{spline, ppval, mkpp, unmkpp} |
|
|
59 ## @end deftypefn |
|
|
60 |
|
|
61 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
62 ## Date: 9. mar 2001 |
|
|
63 ## |
|
|
64 ## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom) |
|
|
65 ## |
|
|
66 ## 4 conditions: |
|
|
67 ## S_k(x_k) = y_k; |
|
|
68 ## S_k(x_k+1) = y_k+1; |
|
|
69 ## S_k'(x_k) = y_k'; |
|
|
70 ## S_k'(x_k+1) = y_k+1'; |
|
|
71 |
|
|
72 function ret = pchip (x, y, xi) |
|
|
73 |
|
|
74 if (nargin < 2 || nargin > 3) |
|
|
75 print_usage (); |
|
|
76 endif |
|
|
77 |
|
|
78 x = x(:); |
|
|
79 n = length (x); |
|
|
80 |
|
|
81 ## Check the size and shape of y |
|
|
82 ndy = ndims (y); |
|
|
83 szy = size (y); |
|
|
84 if (ndy == 2 && (szy(1) == 1 || szy(2) == 1)) |
|
|
85 if (szy(1) == 1) |
|
|
86 y = y'; |
|
|
87 else |
|
|
88 szy = fliplr (szy); |
|
|
89 endif |
|
|
90 else |
|
|
91 y = reshape (y, [prod(szy(1:end-1)), szy(end)])'; |
|
|
92 endif |
|
|
93 |
|
5838
|
94 h = diff (x); |
|
|
95 if (all (h < 0)) |
|
|
96 x = flipud (x); |
|
|
97 h = diff (x); |
|
|
98 y = flipud (y); |
|
|
99 elseif (any (h <= 0)) |
|
5837
|
100 error("pchip: x must be strictly monotonic") |
|
|
101 endif |
|
|
102 |
|
5838
|
103 if (rows (y) != n) |
|
5837
|
104 error("pchip: size of x and y must match"); |
|
|
105 endif |
|
|
106 |
|
|
107 [ry, cy] = size (y); |
|
|
108 if (cy > 1) |
|
|
109 h = kron (diff (x), ones (1, cy)); |
|
|
110 endif |
|
|
111 |
|
|
112 dy = diff (y) ./ h; |
|
|
113 |
|
|
114 a = y; |
|
5838
|
115 b = __pchip_deriv__ (x, y); |
|
5837
|
116 c = - (b(2:n, :) + 2 * b(1:n - 1, :)) ./ h + 3 * diff (a) ./ h .^ 2; |
|
|
117 d = (b(1:n - 1, :) + b(2:n, :)) ./ h.^2 - 2 * diff (a) ./ h.^3; |
|
|
118 |
|
|
119 d = d(1:n - 1, :); c = c(1:n - 1, :); |
|
|
120 b = b(1:n - 1, :); a = a(1:n - 1, :); |
|
|
121 coeffs = [d(:), c(:), b(:), a(:)]; |
|
|
122 pp = mkpp (x, coeffs, szy(1:end-1)); |
|
|
123 |
|
|
124 if (nargin == 2) |
|
|
125 ret = pp; |
|
|
126 else |
|
|
127 ret = ppval (pp, xi); |
|
|
128 endif |
|
|
129 |
|
|
130 endfunction |
|
|
131 |
|
|
132 %!demo |
|
|
133 %! x = 0:8; |
|
|
134 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; |
|
|
135 %! xi = 0:0.01:8; |
|
|
136 %! yspline = spline(x,y,xi); |
|
|
137 %! ypchip = pchip(x,y,xi); |
|
|
138 %! title("pchip and spline fit to discontinuous function"); |
|
|
139 %! plot(xi,yspline,";spline;",... |
|
|
140 %! xi,ypchip,"-;pchip;",... |
|
|
141 %! x,y,"+;data;"); |
|
|
142 %! %------------------------------------------------------------------- |
|
|
143 %! % confirm that pchip agreed better to discontinuous data than spline |
|
|
144 |
|
|
145 %!shared x,y |
|
|
146 %! x = 0:8; |
|
|
147 %! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0]; |
|
|
148 %!assert (pchip(x,y,x), y); |
|
|
149 %!assert (pchip(x,y,x'), y'); |
|
|
150 %!assert (pchip(x',y',x'), y'); |
|
|
151 %!assert (pchip(x',y',x), y); |
|
|
152 %!assert (isempty(pchip(x',y',[]))); |
|
|
153 %!assert (isempty(pchip(x,y,[]))); |
|
|
154 %!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)]) |
