From octave-sources-request at bevo dot che dot wisc dot edu Fri Apr 6 09:05:39 2001 Subject: spline, mkpp and ppval From: Mats Jansson To: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 06 Apr 2001 16:05:07 +0200 I wrote the functions spline.m, mkpp.m and ppval.m. I don't know exactly how Matlab's functions with the same names behave, please tell me if you find any differences. /Mats spline.m ======== ## Copyright (C) 1996, 1997 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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, or (at your option) ## any later version. ## ## Octave 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 Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## -*- texinfo -*- ## at deftypefn {Function File} {@var{pp} =} spline (@var{x}, @var{y}) ## at deftypefnx {Function File} {@var{yi} =} spline (@var{x}, @var{y}, @var{xi}) ## Piece-wise interpolation with cubic splines. ## ## The argument at var{x} should be a vector and @var{y} should either be ## a vector of the same length as at var{x} or a matrix with length ## ( at var{x}) number of columns. The spline-function uses cubic spline ## interpolation to find at var{yi} corresponding to @var{xi}. If the ## third argument is omitted, spline returns a piece-wise polynom on ## pp-form. ## ## Example: ## ## at example ## x = -1:0.25:1; ## xi = -1:0.01:1; ## y = [exp(x)./(1 + 9*x.^2); sin(4*x)]; ## y1 = [exp(xi)./(1 + 9*xi.^2); sin(4*xi)]; ## yi = spline (x, y, xi); ## plot (xi, y1, ["g-;calculated;";"g-;;"], \ ## xi, yi, ["b-;interpolated;";"b-;;"], \ ## x, y, ["ro;interpolation-points;";"ro;;"]) ## at end example ## at end deftypefn ## at seealso{ppval} ## Author: Mats Jansson ## Created: March 5, 2000 - April 6, 2001 ## Keywords: piecewise interpolation cubic spline function yi = spline (x, y, xi) if (nargin < 2), usage ("spline: expecting 2 or 3 arguments."); endif x = x(:); ## Make sure x is a column vector n = size (x, 1); if (n < 2), usage ("yi = spline (x, y, xi): x should have at least 2 elements."); endif if (is_vector (y)), y = y(:); ## Make sure y is a columnvector [nry, ncy] = size (y); if (nry != n), error ("yi = spline (x, y, xi): x and y should have the same number of elements."); endif elseif (is_matrix (y)), [nry, ncy] = size (y); if (ncy != n), ## The number of columns of y should match the length of x (I ## think this is the way Matlabs' spline-function behaves). error ("yi = spline (x, y, xi): The number of columns of y should match the length of x."); else y = y.'; [nry, ncy] = size (y); endif endif [x, ind] = sort (x); if (! all (diff (x))), usage ("yi = spline (x, y, xi): x should be distinct."); endif y = y(ind,:); h = diff (x); d = diff (y); if (n == 2), ## Linear interpolation coef = [d./h(:,ones (1, ncy)); y(1,:)].'; else ## We want the differentiate and the second differentiate to be ## continuous in each point. To find the differentiate, k, in each ## point, we have to solve the equation system: ## ## H * k = 3b, where H is the n x n tri-diagonal matrix ## _ _ ## | 2h(1) h(1) 0 0 ... 0 | ## | h(2) 2h((2)+h(1)) h(1) 0 ... 0 | ## | 0 h(3) 2(h(3)+h(2)) h(2) ... 0 | ## | . . . . | ## | . . . . | ## | . . h(n-1) 2(h(n-1)+h(n-2) h(n-2) | ## |_ 0 0 0 h(n-1) 2h(n-1) _| ## ## and b is the n-by-ncy matrix ## _ _ ## | d(1,:) | ## | h(1)/h(2)*d(2,:)+h(2)/h(1)*d(1,:) | ## | h(2)/h(3)*d(3,:)+h(3)/h(2)*d(2,:) | ## | ... | ## | h(n-2)/h(n-1)*d(n-1,:)+h(n-1)/h(n-2)*d(n-2,:) | ## |_ d(n-1,:) _| ## H = diag (2 * [h(1); h(2:n-1)+h(1:n-2); h(n-1)]) + \ diag ([h(1); h(1:n-2)], 1) + \ diag ([h(2:n-1); h(n-1)], -1); b = zeros (n, ncy); ## Repeat the column-vector h ncy times to form a matrix of the ## same size as d. h = h(:,ones (1, ncy)); b(1,:) = d(1,:); b(n,:) = d(n-1,:); b(2:n-1,:) = h(1:n-2,:)./h(2:n-1,:).*d(2:n-1,:) + \ h(2:n-1,:)./h(1:n-2,:).*d(1:n-2,:); ## H is a sparse matrix, this system could be solved more efficient. k = H \ (3*b); ## We should have n-1 3:rd degree interpolation polynoms for each ## column of y. coef = \ [((h(1:n-1,:).*(k(1:n-1,:)+k(2:n,:))-2*d(1:n-1,:)) ./ \ h(1:n-1,:).^3).'(:); \ ((3*d(1:n-1,:)-2*h(1:n-1,:).*k(1:n-1,:)-h(1:n-1,:).*k(2:n,:)) ./ \ h(1:n-1,:).^2).'(:); \ k(1:n-1,:).'(:); \ y(1:n-1,:).'(:)]; endif yi = mkpp (x, coef, ncy); if (nargin >= 3), yi = ppval (yi, xi); endif endfunction mkpp.m ====== ## Copyright (C) 1996, 1997 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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, or (at your option) ## any later version. ## ## Octave 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 Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## -*- texinfo -*- ## at deftypefn {Function File} {@var{pp} =} mkpp (@var{breaks}, @var{coefficients}) ## at deftypefnx {Function File} {@var{pp} =} mkpp (@var{breaks}, @var{coefficients}, @var{number_of_sets}) ## Define piecewise polynomials. ## at end deftypefn ## at seealso{ppval and spline} ## Author: Mats Jansson ## Created: April 4, 2001 - April 6, 2001 ## Keywords: piecewise polynomial function pp = mkpp (b, c, s) ## b: breakpoints ## c: polynomial coefficients ## s: number of sets if (nargin == 2), s = 1; elseif (nargin < 2 || nargin > 3), usage ("mkpp (breaks, coefs) or mkpp (breaks, coefs, number_of_sets)"); endif s = fix (s); b = b(:); ## Assume n breakpoints (b-values), n-1 polynomials of degree p for ## each set. ## ## c(1): The highest coefficient of the polynomial valid in ## b(1) < x <= b(2) for the first set. ## c(2): The highest coefficient of the polynomial valid in ## b(1) < x <= b(2) for the second set. ## ... ## ## c(s): The highest coefficient of the polynomial valid in ## b(1) < x <= b(2) for the set number s. ## c(s+1): ## The highest coefficient of the polynomial valid in ## b(2) < x <= b(3) for the first set. ## c(s+2): ## The highest coefficient of the polynomial valid in ## b(2) < x <= b(3) for the second set. ## ... ## c(2*s): ## The highest coefficient of the polynomial valid in ## b(2) < x <= b(3) for set number s. ## ... ## c((n-1)*s): ## The highest coefficient of the polynomial valid in ## b(n-1) < x <= b(n) for set number s. ## c((n-1)*s+1): ## The second highest coefficient of the polynomial valid in ## b(1) < x <= b(2) ## ... ## c(2*(n-1)*s): ## The second highest coefficient of the polynomial valid in ## b(n-1) < x <= b(n) ## ... ## c((p+1)*(n-1)*s): ## The last (lowest) coefficient of the polynomial valid in ## b(n-1) < x <= b(n) c = c(:); n = length (b); p = fix (length (c)/(n-1)/s + 100*eps) - 1; pp.breaks = b; pp.coefs = reshape (c, (n-1)*s, p+1); pp.degree = p; pp.sets = s; endfunction ppval.m ======= ## Copyright (C) 1996, 1997 John W. Eaton ## ## This file is part of Octave. ## ## Octave 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, or (at your option) ## any later version. ## ## Octave 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 Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## -*- texinfo -*- ## at deftypefn {Function File} {@var{y} =} ppval (@var{pp}, @var{x}) ## Return the values of the piecewise polynomials at var{pp} at the points ## at var{x} dot ## ## If at var{pp} contains one set of piecewise polynomials, the result has ## the same size as at var{x} dot If @var{pp} contains more than one set of ## piecewise polynomials, each row of the result represents one set. In ## this case each column contains the values of at var{pp} in the points ## at var{x}, were @var{x} is taken in column-major order. ## at end deftypefn ## at seealso{mkpp and spline} ## Author: Mats Jansson ## Created: April 4, 2001 - April 6, 2001 ## Keywords: piecewise polynomial function y = ppval (pp, x) if (nargin != 2), usage ("Y = ppval (PP, X)"); endif if (!isnumeric (x)), error ("ppval (PP, X): Expecting X to be numeric."); endif if (!is_struct (pp)), error ("ppval (PP, X): Expecting PP to be a structure."); endif if (!struct_contains (pp, "breaks") || \ !struct_contains (pp, "coefs") || \ !struct_contains (pp, "degree") || \ !struct_contains (pp, "sets")), error ("ppval (PP, X): PP doesn't seem to be a piece-wise polynomial created with mkpp."); endif [nrb, ncb] = size (pp.breaks); if (ncb != 1 || nrb < 2), error ("ppval (PP, X): Expecting PP.breaks to be a column-vector."); endif [nrx, ncx] = size (x); if (nrx*ncx == 0), y = x; return endif x = x(:).'; ## x is now a row-vector [x, indx] = sort (x); ## How many x-values do we have in each interval? if (nrb == 2), ## Only one interval number_of_x_in_interval_n = nrx*ncx; else xx = x(ones (1, nrb-1),:); bb = pp.breaks(1:nrb-1,ones (1, nrx*ncx)); bb2 = [bb(2:nrb-1,:); Inf*ones(1, nrx*ncx)]; bb = [-Inf*ones(1, nrx*ncx); bb(2:nrb-1,:)]; if (nrx*ncx == 1), number_of_x_in_interval_n = (xx >= bb & xx < bb2)'; else number_of_x_in_interval_n = sum ((xx >= bb & xx < bb2)'); endif endif y = zeros (pp.sets, nrx*ncx); pow = (pp.degree:-1:0).'; ind = find (number_of_x_in_interval_n); x_pointer = 0; ## Step through the intervals with x-values. for nn = 1:length (ind) n = ind(nn); ## Interval number ## The coefficients valid in this interval. The first column is the ## coefficients for the first set, the second column is the ## coefficients for the second set and so on. The first row is the ## coefficients for the highest power, the second row is the ## coefficients for the second highest power... c = pp.coefs((n-1)*pp.sets+1:n*pp.sets,:).'; ## Repeat the coefficients so that we have one set of coefficients ## for each x-value in the interval. cind = (1:pp.sets)'; cind = cind(:, ones (1, number_of_x_in_interval_n(n))); c = c((1:pp.degree+1)', cind); ## A row-vector with the x-values in the interval (local coordinate ## system) xx = x(x_pointer+1:x_pointer+number_of_x_in_interval_n(n)) - \ pp.breaks(n); ## Repeat each x-value as many times as we have sets xx = xx(ones (1, pp.sets),:); xx = xx(:).'; ## Repeat the row-vector of x so that we have one row for each ## coefficient. xx = xx(ones (1, pp.degree+1),:); if (pp.degree == 0), yn = c; else yn = sum (c.*xx.^pow(:, \ ones (1, \ number_of_x_in_interval_n(n)*pp.sets))); endif yn = reshape (yn, pp.sets, number_of_x_in_interval_n(n)); y(:, x_pointer+1:x_pointer+number_of_x_in_interval_n(n)) = yn; x_pointer += number_of_x_in_interval_n(n); endfor y(:,indx) = y; if (pp.sets == 1 && nrx > 1), ## y has the same number of elements as x. Return y with the ## same number of rows and columns. y = reshape (y, nrx, ncx); endif endfunction ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------