From octave-sources-request at bevo dot che dot wisc dot edu Fri Apr 6 14:52:17 2001 Subject: Re: spline, mkpp and ppval From: Paul Kienzle To: Mats Jansson Cc: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 6 Apr 2001 20:51:51 +0100 You may want to compare your implementation with Kai Habel's from http://www.sourceforge.net/projects/octave also included in http://users.powernet.co.uk/kienzle/octave/matcompat His is modelled on the spline toolbox, and so it allows more options. Particularly, you can select the end conditions. I wrote an Octave interface to LAPACK's tridiagonal solver (with an extension to handle the cyclic tridiagonal system required for the periodic spline), so sparse matrices are not necessary. Paul Kienzle pkienzle at kienzle dot powernet dot co dot uk On Fri, Apr 06, 2001 at 04:05:07PM +0200, Mats Jansson wrote: > 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 > ------------------------------------------------------------- > > ------------------------------------------------------------- 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 -------------------------------------------------------------