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