Mercurial > hg > octave-shane
changeset 9976:702b998698ea
implement ppder, ppint, ppjmups
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Sun, 13 Dec 2009 13:18:27 +0100 |
parents | 14ed68363284 |
children | 711aa22ff83d |
files | scripts/ChangeLog scripts/polynomial/module.mk scripts/polynomial/ppder.m scripts/polynomial/ppint.m scripts/polynomial/ppjumps.m |
diffstat | 5 files changed, 162 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,9 @@ +2009-12-13 Jaroslav Hajek <highegg@gmail.com> + + * ppder.m: New function. + * ppint.m: New function. + * ppjumps.m: New function. + 2009-12-09 Rik <octave@nomad.inbox5.com> * Makefile.am: remove install-images target and use automake
--- a/scripts/polynomial/module.mk +++ b/scripts/polynomial/module.mk @@ -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 \
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/ppder.m @@ -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 +## <http://www.gnu.org/licenses/>. + +## -*- 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 +
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/ppint.m @@ -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 +## <http://www.gnu.org/licenses/>. + +## -*- 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 + +
new file mode 100644 --- /dev/null +++ b/scripts/polynomial/ppjumps.m @@ -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 +## <http://www.gnu.org/licenses/>. + +## -*- 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