Mercurial > hg > octave-lyh

comparison scripts/general/quadv.m @ 7771:680631e787aa

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add quadv, quadgk, dblquad and triplequad functions
author David Bateman <dbateman@free.fr>
date Sat, 10 May 2008 21:55:40 +0200
parents
children df9519e9990c
comparison
equal deleted inserted replaced
7770:c6a1a217ac3c 7771:680631e787aa
1 ## Copyright (C) 2008 David Bateman
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} =} quadv (@var{f}, @var{a}, @var{b})
21 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
22 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
23 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
24 ## @deftypefnx {Function File} {[@var{q}, @var{fcnt}] =} quadl (@dots{})
25 ##
26 ## Numerically evaluate integral using adaptive Simpson's rule.
27 ## @code{quadv (@var{f}, @var{a}, @var{b})} approximates the integral of
28 ## @code{@var{f}(@var{x})} to the default absolute tolerance of @code{1e-6}.
29 ## @var{f} is either a function handle, inline function or string
30 ## containing the name of the function to evaluate. The function @var{f}
31 ## must accept a string, and can return a vector representing the
32 ## approximation to @var{n} different sub-functions.
33 ##
34 ## If defined, @var{tol} defines the absolute tolerance to which to
35 ## which to integrate each sub-interval of @code{@var{f}(@var{x})}.
36 ## While if @var{trace} is defined, displays the left end point of the
37 ## current interval, the interval length, and the partial integral.
38 ##
39 ## Additional arguments @var{p1}, etc, are passed directly to @var{f}.
40 ## To use default values for @var{tol} and @var{trace}, one may pass
41 ## empty matrices.
42 ## @seealso{triplequad, dblquad, quad, quadl, quadgk, trapz}
43 ## @end deftypefn
44
45 function [Q, fcnt] = quadv (f, a, b, tol, trace, varargin)
46 if (nargin < 3)
47 print_usage ();
48 endif
49 if (nargin < 4)
50 tol = [];
51 endif
52 if (nargin < 5)
53 trace = [];
54 endif
55 if (isempty (tol))
56 tol = 1e-6;
57 endif
58 if (isempty (trace))
59 trace = 0;
60 endif
61
62 ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule
63 c = (a + b) / 2;
64 fa = feval (f, a, varargin{:});
65 fc = feval (f, c, varargin{:});
66 fb = feval (f, b, varargin{:});
67 fcnt = 3;
68
69 ## If have edge singularities, move edge point by eps*(b-a) as
70 ## discussed in Shampine paper used to implement quadgk
71 if (isinf (fa))
72 fa = feval (f, a + eps * (b-a), varargin{:});
73 endif
74 if (isinf (fb))
75 fb = feval (f, b - eps * (b-a), varargin{:});
76 endif
77
78 h = (b - a) / 2;
79 Q = (b - a) / 6 * (fa + 4 * fc + fb);
80
81 [Q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, Q, fcnt, abs (b - a),
82 tol, trace, varargin{:});
83
84 if (fcnt > 10000)
85 warning("Maximum iteration count reached");
86 elseif (isnan(Q) || isinf (Q))
87 warning ("Infinite or NaN function evaluations were returned");
88 elseif (hmin < (b - a) * eps)
89 warning ("Minimum step size reached. Possibly singular integral");
90 endif
91 endfunction
92
93 function [Q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, Q0,
94 fcnt, hmin, tol, trace, varargin)
95 if (fcnt > 10000)
96 Q = Q0;
97 else
98 d = (a + c) / 2;
99 e = (c + b) / 2;
100 fd = feval (f, d, varargin{:});
101 fe = feval (f, e, varargin{:});
102 fcnt += 2;
103 Q1 = (c - a) / 6 * (fa + 4 * fd + fc);
104 Q2 = (b - c) / 6 * (fc + 4 * fe + fb);
105 Q = Q1 + Q2;
106
107 if (abs(a - c) < hmin)
108 hmin = abs (a - c);
109 endif
110
111 if (trace)
112 disp ([fcnt, a, b-a, Q]);
113 endif
114
115 ## Force at least one adpative step
116 if (fcnt == 5 || abs (Q - Q0) > tol)
117 [Q1, fcnt, hmin] = simpsonstp (f, a, c, d, fa, fc, fd, Q1, fcnt, hmin,
118 tol, trace, varargin{:});
119 [Q2, fcnt, hmin] = simpsonstp (f, c, b, e, fc, fb, fe, Q2, fcnt, hmin,
120 tol, trace, varargin{:});
121 Q = Q1 + Q2;
122 endif
123 endif
124 endfunction
125
126 %!assert (quadv (@sin, 0, 2 * pi), 0, 1e-5)
127 %!assert (quadv (@sin, 0, pi), 2, 1e-5)
128
129 %% Handles weak singularities at the edge
130 %!assert (quadv (@(x) 1 ./ sqrt(x), 0, 1), 2, 1e-5)
131