Mercurial > forge
diff main/splines/pchip.m @ 0:6b33357c7561 octave-forge
Initial revision
author | pkienzle |
---|---|
date | Wed, 10 Oct 2001 19:54:49 +0000 |
parents | |
children | 4890c0365bc5 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/splines/pchip.m Wed Oct 10 19:54:49 2001 +0000 @@ -0,0 +1,87 @@ +## Copyright (C) 2001 Kai Habel +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pp} = } pchip (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{yi} = } pchip (@var{x}, @var{y}, @var{xi}) +## piecewise cubic hermite interpolating polynom. +## pchip preserves the monotonicity of (x,y) +## +## @seealso{ppval, spline, csape} +## @end deftypefn + +## Author: Kai Habel <kai.habel@gmx.de> +## Date: 9. mar 2001 +## 2001-04-03 Paul Kienzle +## * move (:) from definition of l,r to use of l,r so it works with 2.0 + +## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom) +## +## 4 conditions: +## S_k(x_k) = y_k; +## S_k(x_k+1) = y_k+1; +## S_k'(x_k) = y_k'; +## S_k'(x_k+1) = y_k+1'; +## +## 22. april 2001: Hmm, something is wrong, it seems pchip doesn't +## preserve the monotonicity, as expected. So if _you_ know how to +## do it right, please contact me. Kai + +function ret = pchip (x, y, xi) + + if (nargin < 2 || nargin > 3) + usage ("pchip (x, y, [xi])"); + endif + + x = x(:); + n = length (x); + + if (columns(y) == n) y = y'; endif + + [ry,cy] = size (y); + if (cy > 1) + h = kron (diff (x), ones (1, cy)); + else + h = diff (x); + endif + + a = y; + + dy = diff (y) ./ h; + t = diff (sign (dy)) == 0; + l = dy(1:n - 2, :); + r = dy(2:n - 1, :); + b = zeros (size (y)); + s = reshape (0.5 * (l(:) + r(:)), ry - 2, cy); + b(2:ry - 1,:) = s .* t; + + c = - (b(2:n, :) + 2 * b(1:n - 1, :)) ./ h + 3 * diff (a) ./ h .^ 2; + + d = (b(1:n - 1, :) + b(2:n, :)) ./ h.^2 - 2 * diff (a) ./ h.^3; + + d = d(1:n - 1, :); c = c(1:n - 1, :); + b = b(1:n - 1, :); a = a(1:n - 1, :); + + coeffs = [d(:), c(:), b(:), a(:)]; + pp = mkpp (x, coeffs); + + if (nargin == 2) + ret = pp; + else + ret = ppval(pp,xi); + endif + +endfunction