Mercurial > hg > octave-nkf
changeset 5824:448f9982e7fb
[project @ 2006-05-19 06:53:31 by jwe]
author | jwe |
---|---|
date | Fri, 19 May 2006 06:53:31 +0000 |
parents | 080c08b192d8 |
children | f3e37beb03aa |
files | scripts/ChangeLog scripts/polynomial/mkpp.m scripts/polynomial/ppval.m scripts/polynomial/spline.m scripts/polynomial/unmkpp.m |
diffstat | 5 files changed, 415 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2006-05-19 David Bateman <dbateman@free.fr> + + * polynomial/unmkpp.m, polynomial/mkpp.m, polynomial/spline.m, + polynomial/ppval.m: New files from Octave Forge. + 2006-05-17 John W. Eaton <jwe@octave.org> * set/intersection.m: Delete
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/mkpp.m @@ -0,0 +1,67 @@ +## Copyright (C) 2000 Paul Kienzle +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pp} = } mkpp (@var{x}, @var{p}) +## @deftypefnx {Function File} {@var{pp} = } mkpp (@var{x}, @var{p}, @var{d}) +## +## Construct a piece-wise polynomial structure from sample points +## @var{x} and coefficients @var{p}. The ith row of @var{p}, +## @code{@var{p} (@var{i},:)}, contains the coefficients for the polynomial +## over the @var{i}-th interval, ordered from highest to +## lowest. There must be one row for each interval in @var{x}, so +## @code{rows (@var{p}) == length (@var{x}) - 1}. +## +## You can concatenate multiple polynomials of the same order over the +## same set of intervals using @code{@var{p} = [ @var{p1}; @var{p2}; +## @dots{}; @var{pd} ]}. In this case, @code{rows (@var{p}) == @var{d} +## * (length (@var{x}) - 1)}. +## +## @var{d} specifies the shape of the matrix @var{p} for all except the +## last dimension. If @var{d} is not specified it will be computed as +## @code{round (rows (@var{p}) / (length (@var{x}) - 1)) instead. +## +## @seealso{unmkpp, ppval, spline} +## @end deftypefn + +function pp = mkpp (x, P, d) + if (nargin < 2 || nargin > 3) + usage ("pp = mkpp(x,P,d)"); + endif + pp.x = x(:); + pp.P = P; + pp.n = length (x) - 1; + pp.k = columns (P); + if (nargin < 3) + d = round (rows (P) / pp.n); + endif + pp.d = d; + if (pp.n*d != rows (P)) + error ("mkpp: num intervals in x doesn't match num polynomials in P"); + endif +endfunction + +%!demo # linear interpolation +%! x=linspace(0,pi,5)'; +%! t=[sin(x),cos(x)]; +%! m=diff(t)./(x(2)-x(1)); +%! b=t(1:4,:); +%! pp = mkpp(x, [m(:),b(:)]); +%! xi=linspace(0,pi,50); +%! plot(x,t,"x;control;",xi,ppval(pp,xi),";interp;");
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/ppval.m @@ -0,0 +1,58 @@ +## Copyright (C) 2000 Paul Kienzle +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp, @var{xi}) +## +## Evaluate piece-wise polynomial @var{pp} at the points @var{xi}. +## If @code{@var{pp}.d} is a scalar greater than 1, or an array, +## then the returned value @var{yi} will be an array that is +## @code{d1, d1, @dots{}, dk, length (@var{xi})]}. +## +## @seealso{mkpp, unmkpp, spline} +## @end deftypefn + +function yi = ppval (pp, xi) + + if (nargin != 2) + usage ("yi = ppval(pp, xi)") + endif + if (! isstruct (pp)) + error ("ppval: expects a pp structure"); + endif + if (isempty (xi)) + yi = []; + else + transposed = (columns (xi) == 1); + xi = xi(:); + xn = length (xi); + idx = lookup (pp.x(2:pp.n), xi) + 1; + dx = (xi - pp.x(idx))'; + dx = reshape (dx(ones(1,prod(pp.d)),:),[pp.d,xn]); + c = reshape (pp.P(:,1), pp.n, prod (pp.d)); + yi = reshape (c(idx,:)', [pp.d, xn]); + for i = 2 : pp.k; + c = reshape (pp.P(:,i), pp.n, prod (pp.d)); + yi = yi .* dx + reshape (c(idx,:)', [pp.d, xn]); + endfor + if (transposed && isscalar (pp.d) && pp.d == 1) + yi = yi.'; + endif + endif +endfunction
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/spline.m @@ -0,0 +1,224 @@ +## Copyright (C) 2000,2001 Kai Habel +## Copyright (C) 2006 David Bateman +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{pp} = } spline (@var{x}, @var{y}) +## @deftypefnx {Function File} {@var{yi} = } spline (@var{x}, @var{y}, @var{xi}) +## +## Returns the cubic spline interpolation of @var{y} at the point +## @var{x}. Called with two arguments the piece-wse polynomial @var{pp} +## that may later be used with @code{ppval} to evaluate the polynomial +## at specific points. +## +## The variable @var{x} must be a vector of length @var{n}, and @var{y} +## can be either a vector of array. In the case where @var{y} is a +## vector, it can have a length of either @var{n} or @code{@var{n} + 2}. +## If the length of @var{y} is @var{n}, then the 'not-a-knot' end +## condition is used. If the length of @var{y} is @code{@var{n} + 2}, +## then the first and last values of the vector @var{y} are the first +## derivative of the cubic spline at the end-points. +## +## If @var{y} is an array, then the size of @var{y} must have the form +## @iftex +## @tex +## $$[s_1, s_2, \cdots, s_k, n]$$ +## @end tex +## @end iftex +## @ifinfo +## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]} +## @end ifinfo +## or +## @iftex +## @tex +## $$[s_1, s_2, \cdots, s_k, n + 2]$$. +## @end tex +## @end iftex +## @ifinfo +## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n} + 2]}. +## @end ifinfo +## The array is then reshaped internally to a matrix where to leading +## dimension is given by +## @iftex +## @tex +## $$s_1 s_2 \cdots s_k$$ +## @end tex +## @end iftex +## @ifinfo +## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}} +## @end ifinfo +## and each row this matrix is then treated seperately. Note that this +## is exactly the opposite treatment than @code{interp1} and is done +## for compatiability. +## +## Called with a third input argument, @code{spline} evaluates the +## piece-wise spline at the points @var{xi}. There is an equivalence +## between @code{ppval (spline (@var{x}, @var{y}), @var{xi})} and +## @code{spline (@var{x}, @var{y}, @var{xi})}. +## +## @seealso{ppval, mkpp, unmkpp} +## @end deftypefn + +## This code is based on csape.m from octave-forge, but has been +## modified to use the sparse solver code in octave that itself allows +## special casing of tri-diagonal matrices, modified for NDArrays and +## for the treatment of vectors y 2 elements longer than x as complete +## splines. + +function ret = spline (x, y, xi) + + x = x(:); + n = length (x); + if (n < 3) + error ("spline: requires at least 3 points"); + endif + + ## Check the size and shape of y + ndy = ndims (y); + szy = size (y); + if (ndy == 2 && (szy(1) == 1 || szy(2) == 1)) + if (szy(1) == 1) + a = y'; + else + a = y; + szy = fliplr (szy); + endif + else + a = reshape (y, [prod(szy(1:end-1)), szy(end)])'; + endif + complete = false; + if (size (a, 1) == n + 2) + complete = true; + dfs = a(1,:); + dfe = a(end,:); + a = a(2:end-1,:); + endif + + b = c = zeros (size (a)); + h = diff (x); + idx = ones (columns (a),1); + + if (complete) + + if (n == 3) + dg = 1.5 * h(1) - 0.5 * h(2); + c(2:n - 1,:) = 1/dg(1); + else + dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); + dg(1) = dg(1) - 0.5 * h(1); + dg(n - 2) = dg(n-2) - 0.5 * h(n - 1); + + e = h(2:n - 2); + + size(a) + size(h) + n + + g = 3 * diff (a(2:n,:)) ./ h(2:n - 1,idx)\ + - 3 * diff (a(1:n - 1,:)) ./ h(1:n - 2,idx); + g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) \ + - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs); + g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - dfe)\ + - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2); + + c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g; + endif + + c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs + - c(2,:) * h(1)) / (2 * h(1)); + c(n,:) = - (3 / h(n - 1) * (a(n,:) - a(n - 1,:)) - 3 * dfe + + c(n - 1,:) * h(n - 1)) / (2 * h(n - 1)); + b(1:n - 1,:) = diff (a) ./ h(1:n - 1, idx)\ + - h(1:n - 1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); + d = diff (c) ./ (3 * h(1:n - 1, idx)); + + else + + g = zeros(n - 2,columns(a)); + g(1,:) = 3 / (h(1) + h(2)) * (a(3,:) - a(2,:)\ + - h(2) / h(1) * (a(2,:) - a(1,:))); + g(n - 2,:) = 3 / (h(n - 1) + h(n - 2)) *\ + (h(n - 2) / h(n - 1) * (a(n,:) - a(n - 1,:)) -\ + (a(n - 1,:) - a(n - 2,:))); + + if (n > 4) + + g(2:n - 3,:) = 3 * diff (a(3:n - 1,:)) ./ h(3:n - 2,idx)\ + - 3 * diff (a(2:n - 2,:)) ./ h(2:n - 3,idx); + + dg = 2 * (h(1:n - 2) .+ h(2:n - 1)); + dg(1) = dg(1) - h(1); + dg(n - 2) = dg(n-2) - h(n - 1); + + ldg = udg = h(2:n - 2); + udg(1) = udg(1) - h(1); + ldg(n - 3) = ldg(n-3) - h(n - 1); + c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g; + + elseif (n == 4) + + dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)]; + ldg = h(2) - h(3); + udg = h(2) - h(1); + c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g; + + else # n == 3 + + dg= [h(1) + 2 * h(2)]; + c(2:n - 1,:) = g/dg(1); + + endif + + c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:)); + c(n,:) = c(n - 1,:) + h(n - 1) / h(n - 2) * (c(n - 1,:) - c(n - 2,:)); + b = diff (a) ./ h(1:n - 1, idx)\ + - h(1:n - 1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n - 1,:)); + d = diff (c) ./ (3 * h(1:n - 1, idx)); + + endif + + 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(:)]; + ret = mkpp (x, coeffs, szy(1:end-1)); + + if (nargin == 3) + ret = ppval (ret, xi); + endif + +endfunction + +%!demo +%! x = 0:10; y = sin(x); +%! xspline = 0:0.1:10; yspline = spline(x,y,xspline); +%! title("spline fit to points from sin(x)"); +%! plot(xspline,sin(xspline),";original;",... +%! xspline,yspline,"-;interpolation;",... +%! x,y,"+;interpolation points;"); +%! %-------------------------------------------------------- +%! % confirm that interpolated function matches the original + +%!shared x,y +%! x = [0:10]; y = sin(x); +%!assert (spline(x,y,x), y); +%!assert (spline(x,y,x'), y'); +%!assert (spline(x',y',x'), y'); +%!assert (spline(x',y',x), y); +%!assert (isempty(spline(x',y',[]))); +%!assert (isempty(spline(x,y,[]))); +%!assert (spline(x,[y;y],x), [spline(x,y,x);spline(x,y,x)])
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/unmkpp.m @@ -0,0 +1,61 @@ +## Copyright (C) 2000 Paul Kienzle +## +## 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, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{x}, @var{p}, @var{n}, @var{k}, @var{d}] =} unmkpp (@var{pp}) +## +## Extract the components of a piece-wise polynomial structure @var{pp}. +## These are as follows: +## +## @table @asis +## @item @var{x} +## Samples points. +## @item @var{p} +## Polynomial coefficients for points in sample interval. @code{@var{p} +## (@var{i}, :)} contains the coefficients for the polynomial over +## interval @var{i} ordered from highest to lowest. If @code{@var{d} > +## 1}, @code{@var{p} (@var{r}, @var{i}, :)} contains the coeffients for +## the r-th polynomial defined on interval @var{i}. However, this is +## stored as a 2-D array such that @code{@var{c} = reshape (@var{p} (:, +## @var{j}), @var{n}, @var{d})} gives @code{@var{c} (@var{i}, @var{r})} +## is the j-th coefficient of the r-th polynomial over the i-th interval. +## @item @var{n} +## Number of polynomial pieces. +## @item @var{k} +## Order of the polynomial plus 1. +## @item @var{d} +## Number of polynomials defined for each interval. +## @end table +## +## @seealso{mkpp, ppval, spline} +## @end deftypefn + +function [x, P, n, k, d] = unmkpp (pp) + if (nargin == 0) + usage ("[x, P, n, k, d] = unmkpp(pp)") + endif + if (! isstruct (pp)) + error ("unmkpp: expecting piecewise polynomial structure"); + endif + x = pp.x; + P = pp.P; + n = pp.n; + k = pp.k; + d = pp.d; +endfunction