|
5837
|
1 ## Copyright (C) 2001,2002 Kai Habel |
|
|
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 -*- |
|
|
20 ## @deftypefn {Function File} {@var{pp} = } pchip (@var{x}, @var{y}) |
|
|
21 ## @deftypefnx {Function File} {@var{yi} = } pchip (@var{x}, @var{y}, @var{xi}) |
|
|
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 |
|
|
38 ## @ifinfo |
|
|
39 ## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} |
|
|
40 ## @end ifinfo |
|
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 |
|
|
48 ## @ifinfo |
|
|
49 ## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} |
|
|
50 ## @end ifinfo |
|
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)) |
|
5837
|
102 error("pchip: x must be strictly monotonic") |
|
|
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)]) |
