Mercurial > octave

comparison scripts/general/quadl.m @ 30893:e1788b1a315f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use "fcn" as preferred abbreviation for "function" in m-files. * accumarray.m, accumdim.m, quadl.m, quadv.m, randi.m, structfun.m, __is_function__.m, uigetfile.m, uimenu.m, uiputfile.m, doc_cache_create.m, colorspace_conversion_input_check.m, imageIO.m, argnames.m, vectorize.m, vectorize.m, normest1.m, inputname.m, nthargout.m, display_info_file.m, decic.m, ode15i.m, ode15s.m, ode23.m, ode23s.m, ode45.m, odeset.m, check_default_input.m, integrate_adaptive.m, ode_event_handler.m, runge_kutta_23.m, runge_kutta_23s.m, runge_kutta_45_dorpri.m, runge_kutta_interpolate.m, starting_stepsize.m, __all_opts__.m, fminbnd.m, fminsearch.m, fminunc.m, fsolve.m, fzero.m, sqp.m, fplot.m, plotyy.m, __bar__.m, __ezplot__.m, flat_entry.html, profexport.m, movfun.m, bicg.m, bicgstab.m, cgs.m, eigs.m, gmres.m, pcg.m, __alltohandles__.m, __sprand__.m, qmr.m, tfqmr.m, dump_demos.m: Replace "func", "fun", "fn" in documentation and variable names with "fcn".
author Rik <rik@octave.org>
date Mon, 04 Apr 2022 18:14:56 -0700
parents 796f54d4ddbf
children 597f3ee61a48
comparison
equal deleted inserted replaced
30892:1a3cc2811090 30893:e1788b1a315f
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}) 27 ## @deftypefn {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b})
28 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}) 28 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol})
29 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) 29 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace})
30 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) 30 ## @deftypefnx {} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{})
31 ## @deftypefnx {} {[@var{q}, @var{nfun}] =} quadl (@dots{}) 31 ## @deftypefnx {} {[@var{q}, @var{nfev}] =} quadl (@dots{})
32 ## 32 ##
33 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using 33 ## Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
34 ## an adaptive @nospell{Lobatto} rule. 34 ## an adaptive @nospell{Lobatto} rule.
35 ## 35 ##
36 ## @var{f} is a function handle, inline function, or string containing the name 36 ## @var{f} is a function handle, inline function, or string containing the name
53 ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass 53 ## @var{f}. To use default values for @var{tol} and @var{trace}, one may pass
54 ## empty matrices ([]). 54 ## empty matrices ([]).
55 ## 55 ##
56 ## The result of the integration is returned in @var{q}. 56 ## The result of the integration is returned in @var{q}.
57 ## 57 ##
58 ## The optional output @var{nfun} indicates the total number of function 58 ## The optional output @var{nfev} indicates the total number of function
59 ## evaluations performed. 59 ## evaluations performed.
60 ## 60 ##
61 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature - 61 ## Reference: @nospell{W. Gander and W. Gautschi}, @cite{Adaptive Quadrature -
62 ## Revisited}, BIT Vol.@: 40, No.@: 1, March 2000, pp.@: 84--101. 62 ## Revisited}, BIT Vol.@: 40, No.@: 1, March 2000, pp.@: 84--101.
63 ## @url{https://www.inf.ethz.ch/personal/gander/} 63 ## @url{https://www.inf.ethz.ch/personal/gander/}
70 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992. 70 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992.
71 71
72 ## 2003-08-05 Shai Ayal 72 ## 2003-08-05 Shai Ayal
73 ## * permission from author to release as GPL 73 ## * permission from author to release as GPL
74 74
75 function [q, nfun] = quadl (f, a, b, tol = [], trace = false, varargin) 75 function [q, nfev] = quadl (f, a, b, tol = [], trace = false, varargin)
76 76
77 if (nargin < 3) 77 if (nargin < 3)
78 print_usage (); 78 print_usage ();
79 endif 79 endif
80 80
95 if (isempty (trace)) 95 if (isempty (trace))
96 trace = false; 96 trace = false;
97 endif 97 endif
98 98
99 y = feval (f, [a, b], varargin{:}); 99 y = feval (f, [a, b], varargin{:});
100 nfun = 1; 100 nfev = 1;
101 101
102 fa = y(1); 102 fa = y(1);
103 fb = y(2); 103 fb = y(2);
104 104
105 h = b - a; 105 h = b - a;
106 106
107 [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, Inf, nfun, abs (h), 107 [q, nfev, hmin] = adaptlobstp (f, a, b, fa, fb, Inf, nfev, abs (h),
108 tol, trace, varargin{:}); 108 tol, trace, varargin{:});
109 109
110 if (nfun > 10_000) 110 if (nfev > 10_000)
111 warning ("quadl: maximum iteration count reached -- possible singular integral"); 111 warning ("quadl: maximum iteration count reached -- possible singular integral");
112 elseif (any (! isfinite (q(:)))) 112 elseif (any (! isfinite (q(:))))
113 warning ("quadl: infinite or NaN function evaluations were returned"); 113 warning ("quadl: infinite or NaN function evaluations were returned");
114 elseif (hmin < (b - a) * eps) 114 elseif (hmin < (b - a) * eps)
115 warning ("quadl: minimum step size reached -- possible singular integral"); 115 warning ("quadl: minimum step size reached -- possible singular integral");
116 endif 116 endif
117 117
118 endfunction 118 endfunction
119 119
120 function [q, nfun, hmin] = adaptlobstp (f, a, b, fa, fb, q0, nfun, hmin, 120 function [q, nfev, hmin] = adaptlobstp (f, a, b, fa, fb, q0, nfev, hmin,
121 tol, trace, varargin) 121 tol, trace, varargin)
122 122
123 persistent alpha = sqrt (2/3); 123 persistent alpha = sqrt (2/3);
124 persistent beta = 1 / sqrt (5); 124 persistent beta = 1 / sqrt (5);
125 125
126 if (nfun > 10_000) 126 if (nfev > 10_000)
127 q = q0; 127 q = q0;
128 return; 128 return;
129 endif 129 endif
130 130
131 h = (b - a) / 2; 131 h = (b - a) / 2;
134 ml = m - beta*h; 134 ml = m - beta*h;
135 mr = m + beta*h; 135 mr = m + beta*h;
136 mrr = m + alpha*h; 136 mrr = m + alpha*h;
137 x = [mll, ml, m, mr, mrr]; 137 x = [mll, ml, m, mr, mrr];
138 y = feval (f, x, varargin{:}); 138 y = feval (f, x, varargin{:});
139 nfun += 1; 139 nfev += 1;
140 fmll = y(1); 140 fmll = y(1);
141 fml = y(2); 141 fml = y(2);
142 fm = y(3); 142 fm = y(3);
143 fmr = y(4); 143 fmr = y(4);
144 fmrr = y(5); 144 fmrr = y(5);
148 if (abs (b - a) < hmin) 148 if (abs (b - a) < hmin)
149 hmin = abs (b - a); 149 hmin = abs (b - a);
150 endif 150 endif
151 151
152 if (trace) 152 if (trace)
153 disp ([nfun, a, b-a, i1]); 153 disp ([nfev, a, b-a, i1]);
154 endif 154 endif
155 155
156 ## Force at least one adaptive step (nfun > 2 test). 156 ## Force at least one adaptive step (nfev > 2 test).
157 if ((abs (i1-i2) < tol || mll <= a || b <= mrr) && nfun > 2) 157 if ((abs (i1-i2) < tol || mll <= a || b <= mrr) && nfev > 2)
158 q = i1; 158 q = i1;
159 else 159 else
160 q = zeros (6, 1, class (x)); 160 q = zeros (6, 1, class (x));
161 [q(1), nfun, hmin] = adaptlobstp (f, a , mll, fa , fmll, q0/6, nfun, hmin, 161 [q(1), nfev, hmin] = adaptlobstp (f, a , mll, fa , fmll, q0/6, nfev, hmin,
162 tol, trace, varargin{:}); 162 tol, trace, varargin{:});
163 [q(2), nfun, hmin] = adaptlobstp (f, mll, ml , fmll, fml , q0/6, nfun, hmin, 163 [q(2), nfev, hmin] = adaptlobstp (f, mll, ml , fmll, fml , q0/6, nfev, hmin,
164 tol, trace, varargin{:}); 164 tol, trace, varargin{:});
165 [q(3), nfun, hmin] = adaptlobstp (f, ml , m , fml , fm , q0/6, nfun, hmin, 165 [q(3), nfev, hmin] = adaptlobstp (f, ml , m , fml , fm , q0/6, nfev, hmin,
166 tol, trace, varargin{:}); 166 tol, trace, varargin{:});
167 [q(4), nfun, hmin] = adaptlobstp (f, m , mr , fm , fmr , q0/6, nfun, hmin, 167 [q(4), nfev, hmin] = adaptlobstp (f, m , mr , fm , fmr , q0/6, nfev, hmin,
168 tol, trace, varargin{:}); 168 tol, trace, varargin{:});
169 [q(5), nfun, hmin] = adaptlobstp (f, mr , mrr, fmr , fmrr, q0/6, nfun, hmin, 169 [q(5), nfev, hmin] = adaptlobstp (f, mr , mrr, fmr , fmrr, q0/6, nfev, hmin,
170 tol, trace, varargin{:}); 170 tol, trace, varargin{:});
171 [q(6), nfun, hmin] = adaptlobstp (f, mrr, b , fmrr, fb , q0/6, nfun, hmin, 171 [q(6), nfev, hmin] = adaptlobstp (f, mrr, b , fmrr, fb , q0/6, nfev, hmin,
172 tol, trace, varargin{:}); 172 tol, trace, varargin{:});
173 q = sum (q); 173 q = sum (q);
174 endif 174 endif
175 175
176 endfunction 176 endfunction
187 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3), 187 %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3),
188 %! cos (2)/3 - cos (5)/3, 1e-6) 188 %! cos (2)/3 - cos (5)/3, 1e-6)
189 189
190 ## test different tolerances. 190 ## test different tolerances.
191 %!test 191 %!test
192 %! [q, nfun1] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.5, []); 192 %! [q, nfev1] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.5, []);
193 %! assert (q, (60 + sin (4) - sin (64))/12, 0.5); 193 %! assert (q, (60 + sin (4) - sin (64))/12, 0.5);
194 %! [q, nfun2] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []); 194 %! [q, nfev2] = quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []);
195 %! assert (q, (60 + sin (4) - sin (64))/12, 0.1); 195 %! assert (q, (60 + sin (4) - sin (64))/12, 0.1);
196 %! assert (nfun2 > nfun1); 196 %! assert (nfev2 > nfev1);
197 197
198 %!test # test single input/output 198 %!test # test single input/output
199 %! assert (class (quadl (@sin, 0, 1)), "double"); 199 %! assert (class (quadl (@sin, 0, 1)), "double");
200 %! assert (class (quadl (@sin, single (0), 1)), "single"); 200 %! assert (class (quadl (@sin, single (0), 1)), "single");
201 %! assert (class (quadl (@sin, 0, single (1))), "single"); 201 %! assert (class (quadl (@sin, 0, single (1))), "single");