octave_packages/m/polynomial/ppjumps.m

   1 ## Copyright (C) 2008-2012 VZLU Prague, a.s., Czech Republic
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this software; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {@var{jumps} =} ppjumps (@var{pp})
  21 ## Evaluate the boundary jumps of a piecewise polynomial.
  22 ## If there are @math{n} intervals, and the dimensionality of @var{pp} is
  23 ## @math{d}, the resulting array has dimensions @code{[d, n-1]}.
  24 ## @seealso{mkpp}
  25 ## @end deftypefn
  26
  27 function jumps = ppjumps (pp)
  28   if (nargin != 1)
  29     print_usage ();
  30   endif
  31
  32   if (! (isstruct (pp) && strcmp (pp.form, "pp")))
  33     error ("ppjumps: PP must be a structure");
  34   endif
  35
  36   ## Extract info.
  37   [x, P, n, k, d] = unmkpp(pp);
  38   nd = length (d) + 1;
  39
  40   ## Offsets.
  41   dx = diff(x(1:n));
  42   dx = repmat (dx, [prod(d), 1]);
  43   dx = reshape (dx, [d, n-1]);
  44   dx = shiftdim (dx, nd - 1);
  45
  46   ## Use Horner scheme.
  47   if (k>1)
  48     llim = shiftdim (reshape (P(1:(n-1) * prod(d), 1), [d, n-1]), nd - 1);
  49   endif
  50
  51   for i = 2 : k;
  52     llim .*= dx;
  53     llim += shiftdim (reshape (P(1:(n-1) * prod (d), i), [d, n-1]), nd - 1);
  54   endfor
  55
  56   rlim = shiftdim (ppval (pp, x(2:end-1)), nd - 1);
  57   jumps = shiftdim (rlim - llim, 1);
  58 endfunction
  59
  60
  61 %!test
  62 %! p = [1 6 11 6];
  63 %! x = linspace (5, 6, 4);
  64 %! y = polyval (p, x);
  65 %! pp = spline (x, y);
  66 %! jj = ppjumps (pp);
  67 %! assert (jj, [0 0], eps)
  68
  69 %!test
  70 %!
  71 %! breaks = [0 1 2];
  72 %! pp1 = poly (-[1 2 3]);
  73 %! pp2 = poly (-([1 2 3]+1));
  74 %! pp = mkpp (breaks, [pp1;pp2]);
  75 %! assert (ppjumps (pp), 0, eps)
  76
  77 %!test
  78 %!
  79 %! breaks = [0 1 2];
  80 %! pp1 = poly (-[1 2 3]);
  81 %! pp2 = poly (([1 2 3]+1));
  82 %! pp = mkpp (breaks, [pp1;pp2]);
  83 %! j  = - 2 * polyval (pp1, 1);
  84 %! assert (ppjumps (pp), j, eps)