# HG changeset patch # User Jaroslav Hajek # Date 1260706707 -3600 # Node ID 702b998698ea4f38f4dc60dec4ad47debd1f2e64 # Parent 14ed68363284206e40a4c74d4429f6a495b47959 implement ppder, ppint, ppjmups diff -r 14ed68363284 -r 702b998698ea scripts/ChangeLog --- a/scripts/ChangeLog Sat Dec 12 06:18:31 2009 +0100 +++ b/scripts/ChangeLog Sun Dec 13 13:18:27 2009 +0100 @@ -1,3 +1,9 @@ +2009-12-13 Jaroslav Hajek + + * ppder.m: New function. + * ppint.m: New function. + * ppjumps.m: New function. + 2009-12-09 Rik * Makefile.am: remove install-images target and use automake diff -r 14ed68363284 -r 702b998698ea scripts/polynomial/module.mk --- a/scripts/polynomial/module.mk Sat Dec 12 06:18:31 2009 +0100 +++ b/scripts/polynomial/module.mk Sun Dec 13 13:18:27 2009 +0100 @@ -20,6 +20,9 @@ polynomial/polyval.m \ polynomial/polyvalm.m \ polynomial/ppval.m \ + polynomial/ppval.m \ + polynomial/ppint.m \ + polynomial/ppjumps.m \ polynomial/residue.m \ polynomial/roots.m \ polynomial/spline.m \ diff -r 14ed68363284 -r 702b998698ea scripts/polynomial/ppder.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/ppder.m Sun Dec 13 13:18:27 2009 +0100 @@ -0,0 +1,44 @@ +## Copyright (C) 2008, 2009 VZLU Prague, a.s., Czech Republic +## +## 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 3 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 software; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn{Function File} {ppd =} ppder (pp) +## Computes the piecewise derivative of a piecewise polynomial struct @var{pp}. +## @seealso{mkpp,ppval} +## @end deftypefn + +function ppd = ppder (pp) + if (nargin != 1) + print_usage (); + endif + if (! isstruct (pp)) + error ("ppder: expects a pp structure"); + endif + + [x, p, n, k, d] = unmkpp (pp); + p = reshape (p, [], k); + if (k <= 1) + pd = zeros (rows (p), 1); + k = 1; + else + k -= 1; + pd = p(:,1:k) * diag (k:-1:1); + endif + ppd = mkpp (x, pd, d); +endfunction + diff -r 14ed68363284 -r 702b998698ea scripts/polynomial/ppint.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/ppint.m Sun Dec 13 13:18:27 2009 +0100 @@ -0,0 +1,55 @@ +## Copyright (C) 2008, 2009 VZLU Prague, a.s., Czech Republic +## +## 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 3 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 software; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn{Function File} {ppd =} ppint (pp, c) +## Computes the antiderivative of a piecewise polynomial struct @var{pp}. +## @var{c}, if given, is the constant of integration. +## @seealso{mkpp,ppval} +## @end deftypefn + +function ppi = ppint (pp, c) + if (nargin < 1 || nargin > 2) + print_usage (); + endif + if (! isstruct (pp)) + error ("ppder: expects a pp structure"); + endif + + [x, p, n, k, d] = unmkpp (pp); + p = reshape (p, [], k); + + ## Get piecewise antiderivatives + pi = p / diag (k:-1:1); + k += 1; + if (nargin == 1) + pi(:,k) = 0; + else + pi(:,k) = repmat (c(:), n, 1); + endif + + ppi = mkpp (x, pi, d); + + ## Adjust constants so the the result is continuous. + + jumps = reshape (ppjumps (ppi), prod (d), n-1); + ppi.P(:,2:n,k) -= cumsum (jumps, 2); + +endfunction + + diff -r 14ed68363284 -r 702b998698ea scripts/polynomial/ppjumps.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/ppjumps.m Sun Dec 13 13:18:27 2009 +0100 @@ -0,0 +1,54 @@ +## Copyright (C) 2008, 2009 VZLU Prague, a.s., Czech Republic +## +## 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 3 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 software; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn{Function File} {ppd =} ppjumps (pp) +## Evaluates the boundary jumps of a piecewise polynomial. +## If there are n intervals, and the dimensionality of pp is d, +## the resulting array has dimensions @code{[d, n-1]}. +## @end deftypefn + +function jumps = ppjumps (pp) + if (nargin != 1) + print_usage (); + endif + if (! isstruct (pp)) + error ("ppjumps: expects a pp structure"); + endif + + ## Extract info. + x = pp.x; + P = pp.P; + d = pp.d; + [nd, n, k] = size (P); + + ## Offsets. + dx = diff (x(1:n)).'; + dx = dx(ones (1, nd), :); # spread (do nothing in 1D) + + ## Use Horner scheme to get limits from the left. + llim = P(:,1:n-1,1); + for i = 2:k; + llim .*= dx; + llim += P(:,1:n-1,i); + endfor + + rlim = P(:,2:n,k); # limits from the right + jumps = reshape (rlim - llim, [d, n-1]); + +endfunction