|
5820
|
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} =} cumtrapz (@var{y}) |
|
|
22 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@var{x}, @var{y}) |
|
|
23 ## @deftypefnx {Function File} {@var{z} =} cumtrapz (@dots{}, @var{dim}) |
|
|
24 ## |
|
|
25 ## Cumulative numerical intergration using trapezodial method. |
|
|
26 ## @code{trapz (@var{y})} computes the cummulative integral of the |
|
|
27 ## @var{y} along the first non singleton dimension. If the argument |
|
|
28 ## @var{x} is omitted a equally spaced vector is assumed. @code{trapz |
|
|
29 ## (@var{x}, @var{y})} evaluates the cummulative integral with respect |
|
|
30 ## to @var{x}. |
|
|
31 ## |
|
|
32 ## @seealso{trapz,cumsum} |
|
|
33 ## @end deftypefn |
|
|
34 |
|
|
35 ## Author: Kai Habel <kai.habel@gmx.de> |
|
|
36 ## |
|
|
37 ## also: June 2000 Paul Kienzle (fixes,suggestions) |
|
|
38 ## 2006-05-12 David Bateman - Modified for NDArrays |
|
|
39 |
|
|
40 function z = cumtrapz (x, y, dim) |
|
|
41 |
|
|
42 if (nargin < 1) || (nargin > 3) |
|
6046
|
43 print_usage (); |
|
5820
|
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 * cumsum (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 * cumsum ((x(idx1{:}) - x(idx2{:})) .* |
|
|
97 (y(idx1{:}) + y(idx2{:})), dim); |
|
|
98 endif |
|
|
99 |
|
|
100 sz(dim) = 1; |
|
|
101 z = cat (dim, zeros (sz), z); |
|
|
102 |
|
|
103 endfunction |
|
|
104 |
|
|
105 %!shared x1,x2,y |
|
|
106 %! x1 = [0,0,0;2,2,2]; |
|
|
107 %! x2 = [0,2,4;0,2,4]; |
|
|
108 %! y = [1,2,3;4,5,6]; |
|
|
109 %!assert (cumtrapz(y),[0,0,0;2.5,3.5,4.5]) |
|
|
110 %!assert (cumtrapz(x1,y),[0,0,0;5,7,9]) |
|
|
111 %!assert (cumtrapz(y,1),[0,0,0;2.5,3.5,4.5]) |
|
|
112 %!assert (cumtrapz(x1,y,1),[0,0,0;5,7,9]) |
|
|
113 %!assert (cumtrapz(y,2),[0,1.5,4;0,4.5,10]) |
|
|
114 %!assert (cumtrapz(x2,y,2),[0,3,8;0,9,20]) |
