octave_packages/m/general/trapz.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} =} trapz (@var{y})
  21 ## @deftypefnx {Function File} {@var{q} =} trapz (@var{x}, @var{y})
  22 ## @deftypefnx {Function File} {@var{q} =} trapz (@dots{}, @var{dim})
  23 ##
  24 ## Numerically evaluate the integral of points @var{y} using the trapezoidal
  25 ## method.
  26 ## @w{@code{trapz (@var{y})}} computes the integral of @var{y} along the first
  27 ## non-singleton dimension.  When the argument @var{x} is omitted an
  28 ## equally spaced @var{x} vector with unit spacing (1) is assumed.
  29 ## @code{trapz (@var{x}, @var{y})} evaluates the integral with respect
  30 ## to the spacing in @var{x} and the values in @var{y}.  This is useful if
  31 ## the points in @var{y} have been sampled unevenly.
  32 ## If the optional @var{dim} argument is given, operate along this dimension.
  33 ##
  34 ## If @var{x} is not specified then unit spacing will be used.  To scale
  35 ## the integral to the correct value you must multiply by the actual spacing
  36 ## value (deltaX).  As an example, the integral of @math{x^3} over the range
  37 ## [0, 1] is @math{x^4/4} or 0.25.  The following code uses @code{trapz} to
  38 ## calculate the integral in three different ways.
  39 ##
  40 ## @example
  41 ## @group
  42 ## x = 0:0.1:1;
  43 ## y = x.^3;
  44 ## q = trapz (y)
  45 ##   @result{} q = 2.525   # No scaling
  46 ## q * 0.1
  47 ##   @result{} q = 0.2525  # Approximation to integral by scaling
  48 ## trapz (x, y)
  49 ##   @result{} q = 0.2525  # Same result by specifying @var{x}
  50 ## @end group
  51 ## @end example
  52 ##
  53 ## @seealso{cumtrapz}
  54 ## @end deftypefn
  55
  56 ## Author:      Kai Habel <kai.habel@gmx.de>
  57 ##
  58 ## also: June 2000 - Paul Kienzle (fixes,suggestions)
  59 ## 2006-05-12 David Bateman - Modified for NDArrays
  60
  61 function z = trapz (x, y, dim)
  62
  63   if (nargin < 1) || (nargin > 3)
  64     print_usage ();
  65   endif
  66
  67   have_xy = have_dim = false;
  68
  69   if (nargin == 3)
  70     have_xy = true;
  71     have_dim = true;
  72   elseif (nargin == 2)
  73     if (! size_equal (x, y) && isscalar (y))
  74       dim = y;
  75       have_dim = true;
  76     else
  77       have_xy = true;
  78     endif
  79   endif
  80
  81   if (have_xy)
  82     nd = ndims (y);
  83     sz = size (y);
  84   else
  85     nd = ndims (x);
  86     sz = size (x);
  87   endif
  88
  89   if (! have_dim)
  90     ## Find the first non-singleton dimension.
  91     (dim = find (sz > 1, 1)) || (dim = 1);
  92   else
  93     if (!(isscalar (dim) && dim == fix (dim))
  94         || !(1 <= dim && dim <= nd))
  95       error ("trapz: DIM must be an integer and a valid dimension");
  96     endif
  97   endif
  98
  99   n = sz(dim);
 100   idx1 = idx2 = repmat ({':'}, [nd, 1]);
 101   idx1{dim} = 2 : n;
 102   idx2{dim} = 1 : (n - 1);
 103
 104   if (! have_xy)
 105     z = 0.5 * sum (x(idx1{:}) + x(idx2{:}), dim);
 106   else
 107     if (isvector (x) && !isvector (y))
 108       if (length (x) != sz(dim))
 109         error ("trapz: length of X and length of Y along DIM must match");
 110       endif
 111       ## Reshape vector to point along dimension DIM
 112       shape = ones (nd, 1);
 113       shape(dim) = sz(dim);
 114       x = reshape (x, shape);
 115       z = 0.5 * sum (bsxfun (@times, diff (x), y(idx1{:}) + y(idx2{:})), dim);
 116     else
 117       if (! size_equal (x, y))
 118         error ("trapz: X and Y must have same shape");
 119       endif
 120       z = 0.5 * sum (diff (x, 1, dim) .* (y(idx1{:}) + y(idx2{:})), dim);
 121     endif
 122   endif
 123 endfunction
 124
 125
 126 %!assert (trapz(1:5), 12)
 127 %!assert (trapz(0:0.5:2,1:5), 6)
 128 %!assert (trapz([1:5;1:5].',1), [12,12])
 129 %!assert (trapz([1:5;1:5],2), [12;12])
 130 %!assert (trapz(repmat(reshape(1:5,1,1,5),2,2), 3), [12 12; 12 12])
 131 %!assert (trapz([0:0.5:2;0:0.5:2].',[1:5;1:5].',1), [6, 6])
 132 %!assert (trapz([0:0.5:2;0:0.5:2],[1:5;1:5],2), [6; 6])
 133 %!assert (trapz(repmat(reshape([0:0.5:2],1,1,5),2,2), ...
 134 %!              repmat(reshape(1:5,1,1,5),2,2), 3), [6 6; 6 6])
 135 %!assert (trapz(0:0.5:2,[(1:5)',(1:5)']), [6, 6])
 136 %!assert (trapz(0:0.5:2,[(1:5);(1:5)],2), [6; 6])
 137 %!assert (trapz(0:0.5:2,repmat(reshape(1:5,1,1,5),2,2),3), [6 6; 6 6])
 138