Mercurial > hg > octave-lyh

changeset 9976:702b998698ea

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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