octave_packages/m/general/cumtrapz.m

   1 ## Copyright (C) 2000-2012 Kai Habel
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {@var{q} =} cumtrapz (@var{y})
  21 ## @deftypefnx {Function File} {@var{q} =} cumtrapz (@var{x}, @var{y})
  22 ## @deftypefnx {Function File} {@var{q} =} cumtrapz (@dots{}, @var{dim})
  23 ##
  24 ## Cumulative numerical integration of points @var{y} using the trapezoidal
  25 ## method.
  26 ## @w{@code{cumtrapz (@var{y})}} computes the cumulative integral of @var{y}
  27 ## along the first non-singleton dimension.  Where @code{trapz} reports
  28 ## only the overall integral sum, @code{cumtrapz} reports the current partial
  29 ## sum value at each point of @var{y}.  When the argument @var{x} is omitted
  30 ## an equally spaced @var{x} vector with unit spacing (1) is assumed.
  31 ## @code{cumtrapz (@var{x}, @var{y})} evaluates the integral with respect to
  32 ## the spacing in @var{x} and the values in @var{y}.  This is useful if the
  33 ## points in @var{y} have been sampled unevenly.  If the optional @var{dim}
  34 ## argument is given, operate along this dimension.
  35 ##
  36 ## If @var{x} is not specified then unit spacing will be used.  To scale
  37 ## the integral to the correct value you must multiply by the actual spacing
  38 ## value (deltaX).
  39 ## @seealso{trapz, cumsum}
  40 ## @end deftypefn
  41
  42 ## Author:      Kai Habel <kai.habel@gmx.de>
  43 ##
  44 ## also: June 2000 Paul Kienzle (fixes,suggestions)
  45 ## 2006-05-12 David Bateman - Modified for NDArrays
  46
  47 function z = cumtrapz (x, y, dim)
  48
  49   if (nargin < 1) || (nargin > 3)
  50     print_usage ();
  51   endif
  52
  53   have_xy = have_dim = false;
  54
  55   if (nargin == 3)
  56     have_xy = true;
  57     have_dim = true;
  58   elseif (nargin == 2)
  59     if (! size_equal (x, y) && isscalar (y))
  60       dim = y;
  61       have_dim = true;
  62     else
  63       have_xy = true;
  64     endif
  65   endif
  66
  67   if (have_xy)
  68     nd = ndims (y);
  69     sz = size (y);
  70   else
  71     nd = ndims (x);
  72     sz = size (x);
  73   endif
  74
  75   if (! have_dim)
  76     ## Find the first non-singleton dimension.
  77     (dim = find (sz > 1, 1)) || (dim = 1);
  78   else
  79     if (!(isscalar (dim) && dim == fix (dim))
  80         || !(1 <= dim && dim <= nd))
  81       error ("trapz: DIM must be an integer and a valid dimension");
  82     endif
  83   endif
  84
  85   n = sz(dim);
  86   idx1 = idx2 = repmat ({':'}, [nd, 1]);
  87   idx1{dim} = 2 : n;
  88   idx2{dim} = 1 : (n - 1);
  89
  90   if (! have_xy)
  91     z = 0.5 * cumsum (x(idx1{:}) + x(idx2{:}), dim);
  92   else
  93     if (isvector (x) && !isvector (y))
  94       if (length (x) != sz(dim))
  95         error ("cumtrapz: length of X and length of Y along DIM must match");
  96       endif
  97       ## Reshape vector to point along dimension DIM
  98       shape = ones (nd, 1);
  99       shape(dim) = sz(dim);
 100       x = reshape (x, shape);
 101       z = 0.5 * cumsum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim);
 102     else
 103       if (! size_equal (x, y))
 104         error ("cumtrapz: X and Y must have same shape");
 105       endif
 106       z = 0.5 * cumsum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim);
 107     endif
 108   endif
 109
 110   sz(dim) = 1;
 111   z = cat (dim, zeros (sz), z);
 112
 113 endfunction
 114
 115
 116 %!shared x1,x2,y
 117 %! x1 = [0,0,0;2,2,2];
 118 %! x2 = [0,2,4;0,2,4];
 119 %! y = [1,2,3;4,5,6];
 120 %!assert (cumtrapz(y), [0,0,0;2.5,3.5,4.5])
 121 %!assert (cumtrapz(x1,y), [0,0,0;5,7,9])
 122 %!assert (cumtrapz(y,1), [0,0,0;2.5,3.5,4.5])
 123 %!assert (cumtrapz(x1,y,1), [0,0,0;5,7,9])
 124 %!assert (cumtrapz(y,2), [0,1.5,4;0,4.5,10])
 125 %!assert (cumtrapz(x2,y,2), [0,3,8;0,9,20])
 126 %% Test ND-array implementation
 127 %!shared x1,x2,y
 128 %! x1 = 1:3;
 129 %! x2 = reshape ([0,2,4;0,2,4], [1 2 3]);
 130 %! y = reshape ([1,2,3;4,5,6], [1 2 3]);
 131 %!assert (cumtrapz(y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3]))
 132 %!assert (cumtrapz(x1,y,3), reshape([0,1.5,4;0,4.5,10],[1 2 3]))
 133 %!assert (cumtrapz(x2,y,3), reshape([0,3,8;0,9,20],[1 2 3]))
 134