Mercurial > octave-nkf

comparison scripts/general/trapz.m @ 5820:27c966e4b2dc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-05-17 21:00:54 by jwe]
author jwe
date Wed, 17 May 2006 21:00:54 +0000
parents
children 34f96dd5441b
comparison
equal deleted inserted replaced
5819:e54c11df0524 5820:27c966e4b2dc
1 ## Copyright (C) 2000 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 2, or (at your option)
8 ## 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, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {@var{z} =} trapz (@var{y})
22 ## @deftypefnx {Function File} {@var{z} =} trapz (@var{x}, @var{y})
23 ## @deftypefnx {Function File} {@var{z} =} trapz (@dots{}, @var{dim})
24 ##
25 ## Numerical intergration using trapezodial method. @code{trapz
26 ## (@var{y})} computes the integral of the @var{y} along the first
27 ## non singleton dimension. If the argument @var{x} is omitted a
28 ## equally spaced vector is assumed. @code{trapz (@var{x}, @var{y})}
29 ## evaluates the integral with respect to @var{x}.
30 ##
31 ## @seealso{cumtrapz}
32 ## @end deftypefn
33
34 ## Author: Kai Habel <kai.habel@gmx.de>
35 ##
36 ## also: June 2000 - Paul Kienzle (fixes,suggestions)
37 ## 2006-05-12 David Bateman - Modified for NDArrays
38
39 function z = trapz (x, y, dim)
40
41
42 if (nargin < 1) || (nargin > 3)
43 usage ("trapz (x, y, dim)");
44 endif
45
46 nd = ndims (x);
47 sz = size (x);
48
49 have_x = false;
50 have_dim = false;
51 if (nargin == 3)
52 have_x = true;
53 have_dim = true;
54 endif
55 if (nargin == 2)
56 if (size (x) != size (y) && isscalar (y))
57 dim = y;
58 have_dim = true;
59 else
60 have_x = true;
61 endif
62 endif
63
64 if (! have_dim)
65 ## Find the first singleton dimension.
66 dim = 0;
67 while (dim < nd && sz(dim+1) == 1)
68 dim++;
69 endwhile
70 dim++;
71 if (dim > nd)
72 dim = 1;
73 endif
74 else
75 dim = floor (dim);
76 if (dim < 1 || dim > nd)
77 error ("cumtrapz: invalid dimension along which to sort");
78 endif
79 endif
80
81 n = sz(dim);
82 idx1 = cell ();
83 for i = 1:nd
84 idx1{i} = 1:sz(i);
85 endfor
86 idx2 = idx1;
87 idx1{dim} = 2 : n;
88 idx2{dim} = 1 : (n - 1);
89
90 if (! have_x)
91 z = 0.5 * sum (x(idx1{:}) + x(idx2{:}), dim);
92 else
93 if (size (x) != size (y))
94 error ("cumtrapz: x and y must have same shape");
95 endif
96 z = 0.5 * sum ((x(idx1{:}) - x(idx2{:})) .*
97 (y(idx1{:}) + y(idx2{:})), dim);
98 endif
99 endfunction
100
101 %!assert (trapz(1:5), 12)
102 %!assert (trapz(0:0.5:2,1:5), 6)
103 %!assert (trapz([1:5;1:5],2),[12;12])
104 %!assert (trapz([1:5;1:5].',1),[12,12])
105 %!assert (trapz([0:0.5:2;0:0.5:2],[1:5;1:5],2),[6;6])
106 %!assert (trapz([0:0.5:2;0:0.5:2].',[1:5;1:5].',1),[6,6])