Mercurial > hg > octave-lyh

comparison scripts/general/quadv.m @ 11587:c792872f8942

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
all script files: untabify and strip trailing whitespace
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:35:29 -0500
parents 3c6e8aaa9555
children 16cca721117b
comparison
equal deleted inserted replaced
11586:12df7854fa7c 11587:c792872f8942
21 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}) 21 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol})
22 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) 22 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
23 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) 23 ## @deftypefnx {Function File} {@var{q} =} quadv (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
24 ## @deftypefnx {Function File} {[@var{q}, @var{fcnt}] =} quadv (@dots{}) 24 ## @deftypefnx {Function File} {[@var{q}, @var{fcnt}] =} quadv (@dots{})
25 ## 25 ##
26 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} 26 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b}
27 ## using adaptive Simpson's rule. 27 ## using adaptive Simpson's rule.
28 ## @var{f} is either a function handle, inline function or string 28 ## @var{f} is either a function handle, inline function or string
29 ## containing the name of the function to evaluate. 29 ## containing the name of the function to evaluate.
30 ## The function defined by @var{f} may be a scalar, vector or array-valued. 30 ## The function defined by @var{f} may be a scalar, vector or array-valued.
31 ## 31 ##
32 ## If a value for @var{tol} is given, it defines the tolerance used to stop 32 ## If a value for @var{tol} is given, it defines the tolerance used to stop
33 ## the adaptation procedure, otherwise the default value of 1e-6 is used. 33 ## the adaptation procedure, otherwise the default value of 1e-6 is used.
34 ## 34 ##
35 ## The algorithm used by @code{quadv}, involves recursively subdividing the 35 ## The algorithm used by @code{quadv}, involves recursively subdividing the
36 ## integration interval and applying Simpson's rule on each sub-interval. 36 ## integration interval and applying Simpson's rule on each sub-interval.
37 ## If @var{trace} is @var{true}, after computing each of these partial 37 ## If @var{trace} is @var{true}, after computing each of these partial
49 function [q, fcnt] = quadv (f, a, b, tol, trace, varargin) 49 function [q, fcnt] = quadv (f, a, b, tol, trace, varargin)
50 if (nargin < 3) 50 if (nargin < 3)
51 print_usage (); 51 print_usage ();
52 endif 52 endif
53 if (nargin < 4) 53 if (nargin < 4)
54 tol = []; 54 tol = [];
55 endif 55 endif
56 if (nargin < 5) 56 if (nargin < 5)
57 trace = []; 57 trace = [];
58 endif 58 endif
59 if (isa (a, "single") || isa (b, "single")) 59 if (isa (a, "single") || isa (b, "single"))
60 myeps = eps ("single"); 60 myeps = eps ("single");
61 else 61 else
62 myeps = eps; 62 myeps = eps;
63 endif 63 endif
64 if (isempty (tol)) 64 if (isempty (tol))
65 tol = 1e-6; 65 tol = 1e-6;
66 endif 66 endif
67 if (isempty (trace)) 67 if (isempty (trace))
68 trace = 0; 68 trace = 0;
69 endif 69 endif
70 70
71 ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule 71 ## Split the interval into 3 abscissa, and apply a 3 point Simpson's rule
72 c = (a + b) / 2; 72 c = (a + b) / 2;
73 fa = feval (f, a, varargin{:}); 73 fa = feval (f, a, varargin{:});
84 fb = feval (f, b - myeps * (b-a), varargin{:}); 84 fb = feval (f, b - myeps * (b-a), varargin{:});
85 endif 85 endif
86 86
87 h = (b - a); 87 h = (b - a);
88 q = (b - a) / 6 * (fa + 4 * fc + fb); 88 q = (b - a) / 6 * (fa + 4 * fc + fb);
89 89
90 [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, fcnt, abs (h), 90 [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q, fcnt, abs (h),
91 tol, trace, varargin{:}); 91 tol, trace, varargin{:});
92 92
93 if (fcnt > 10000) 93 if (fcnt > 10000)
94 warning ("maximum iteration count reached"); 94 warning ("maximum iteration count reached");
95 elseif (isnan (q) || isinf (q)) 95 elseif (isnan (q) || isinf (q))
97 elseif (hmin < (b - a) * myeps) 97 elseif (hmin < (b - a) * myeps)
98 warning ("minimum step size reached -- possibly singular integral"); 98 warning ("minimum step size reached -- possibly singular integral");
99 endif 99 endif
100 endfunction 100 endfunction
101 101
102 function [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0, 102 function [q, fcnt, hmin] = simpsonstp (f, a, b, c, fa, fb, fc, q0,
103 fcnt, hmin, tol, trace, varargin) 103 fcnt, hmin, tol, trace, varargin)
104 if (fcnt > 10000) 104 if (fcnt > 10000)
105 q = q0; 105 q = q0;
106 else 106 else
107 d = (a + c) / 2; 107 d = (a + c) / 2;