[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
spline, mkpp and ppval
From: |
Mats Jansson |
Subject: |
spline, mkpp and ppval |
Date: |
Fri, 06 Apr 2001 13:54:51 +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 -*-
## @deftypefn {Function File} address@hidden =} spline (@var{x}, @var{y})
## @deftypefnx {Function File} address@hidden =} spline (@var{x}, @var{y},
@var{xi})
## Piece-wise interpolation with cubic splines.
##
## The argument @var{x} should be a vector and @var{y} should either be
## a vector of the same length as @var{x} or a matrix with length
## (@var{x}) number of columns. The spline-function uses cubic spline
## interpolation to find @var{yi} corresponding to @var{xi}. If the
## third argument is omitted, spline returns a piece-wise polynom on
## pp-form.
##
## Example:
##
## @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;;"])
## @end example
## @end deftypefn
## @seealso{ppval}
## Author: Mats Jansson <address@hidden>
## 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 -*-
## @deftypefn {Function File} address@hidden =} mkpp (@var{breaks},
@var{coefficients})
## @deftypefnx {Function File} address@hidden =} mkpp (@var{breaks},
@var{coefficients}, @var{number_of_sets})
## Define piecewise polynomials.
## @end deftypefn
## @seealso{ppval and spline}
## Author: Mats Jansson <address@hidden>
## 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 -*-
## @deftypefn {Function File} address@hidden =} ppval (@var{pp}, @var{x})
## Return the values of the piecewise polynomials @var{pp} at the points
## @var{x}.
##
## If @var{pp} contains one set of piecewise polynomials, the result has
## the same size as @var{x}. 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 @var{pp} in the points
## @var{x}, were @var{x} is taken in column-major order.
## @end deftypefn
## @seealso{mkpp and spline}
## Author: Mats Jansson <address@hidden>
## 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
-------------------------------------------------------------
- spline, mkpp and ppval,
Mats Jansson <=