Mercurial > hg > octave-nkf

comparison scripts/polynomial/pchip.m @ 5837:55404f3b0da1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-06-01 19:05:31 by jwe]
author jwe
date Thu, 01 Jun 2006 19:05:32 +0000
parents
children 376e02b2ce70
comparison
equal deleted inserted replaced
5836:ed69a3b5b3d0 5837:55404f3b0da1
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
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))
100 error("pchip: x must be strictly monotonic")
101 endif
102
103 if (rows(y) != n)
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;
115 b = __pchip_deriv__(x,y);
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)])