|
5837
|
1 ## Copyright (C) 1998 Walter Gautschi |
|
|
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{q} =} quadl (@var{f}, @var{a}, @var{b}) |
|
|
22 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}) |
|
|
23 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}) |
|
|
24 ## @deftypefnx {Function File} {@var{q} =} quadl (@var{f}, @var{a}, @var{b}, @var{tol}, @var{trace}, @var{p1}, @var{p2}, @dots{}) |
|
|
25 ## |
|
|
26 ## Numerically evaluate integral using adaptive Lobatto rule. |
|
|
27 ## @code{quadl (@var{f}, @var{a}, @var{b})} approximates the integral of |
|
|
28 ## @code{@var{f}(@var{x})} to machine precision. @var{f} is either a |
|
|
29 ## function handle, inline function or string containing the name of |
|
|
30 ## the function to evaluate. The function @var{f} must return a vector |
|
|
31 ## of output values if given a vector of input values. |
|
|
32 ## |
|
|
33 ## If defined, @var{tol} defines the relative tolerance to which to |
|
|
34 ## which to integrate @code{@var{f}(@var{x})}. While if @var{trace} is |
|
|
35 ## defined, displays the left end point of the current interval, the |
|
|
36 ## interval length, and the partial integral. |
|
|
37 ## |
|
|
38 ## Additional arguments @var{p1}, etc, are passed directly to @var{f}. |
|
|
39 ## To use default values for @var{tol} and @var{trace}, one may pass |
|
|
40 ## empty matrices. |
|
|
41 ## |
|
|
42 ## Reference: W. Gander and W. Gautschi, 'Adaptive Quadrature - |
|
|
43 ## Revisited', BIT Vol. 40, No. 1, March 2000, pp. 84--101. |
|
|
44 ## @url{http://www.inf.ethz.ch/personal/gander/} |
|
|
45 ## |
|
|
46 ## @end deftypefn |
|
|
47 |
|
|
48 ## Author: Walter Gautschi |
|
|
49 ## Date: 08/03/98 |
|
|
50 ## Reference: Gander, Computermathematik, Birkhaeuser, 1992. |
|
|
51 |
|
|
52 ## 2003-08-05 Shai Ayal |
|
|
53 ## * permission from author to release as GPL |
|
|
54 ## 2004-02-10 Paul Kienzle |
|
|
55 ## * renamed to quadl for compatibility |
|
|
56 ## * replace global variable terminate2 with local function need_warning |
|
|
57 ## * add paper ref to docs |
|
|
58 |
|
5838
|
59 function Q = quadl (f, a, b, tol, trace, varargin) |
|
|
60 need_warning (1); |
|
5837
|
61 if (nargin < 4) |
|
5838
|
62 tol = []; |
|
5837
|
63 endif |
|
|
64 if (nargin < 5) |
|
|
65 trace = []; |
|
|
66 endif |
|
5838
|
67 if (isempty (tol)) |
|
5837
|
68 tol = eps; |
|
|
69 endif |
|
5838
|
70 if (isempty (trace)) |
|
|
71 trace = 0; |
|
5837
|
72 endif |
|
|
73 if (tol < eps) |
|
|
74 tol = eps; |
|
|
75 endif |
|
|
76 |
|
|
77 m = (a+b)/2; |
|
|
78 h = (b-a)/2; |
|
|
79 alpha = sqrt(2/3); |
|
|
80 beta = 1/sqrt(5); |
|
5838
|
81 |
|
5837
|
82 x1 = .942882415695480; |
|
|
83 x2 = .641853342345781; |
|
|
84 x3 = .236383199662150; |
|
5838
|
85 |
|
|
86 x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ... |
|
|
87 m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b]; |
|
|
88 |
|
|
89 y = feval (f, x, varargin{:}); |
|
|
90 |
|
5837
|
91 fa = y(1); |
|
|
92 fb = y(13); |
|
5838
|
93 |
|
|
94 i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9))); |
|
|
95 |
|
|
96 i1 = (h/1470)*(77*(y(1)+y(13)) |
|
|
97 + 432*(y(3)+y(11)) |
|
|
98 + 625*(y(5)+y(9)) |
|
|
99 + 672*y(7)); |
|
|
100 |
|
|
101 is = h*(.0158271919734802*(y(1)+y(13)) |
|
|
102 +.0942738402188500*(y(2)+y(12)) |
|
|
103 + .155071987336585*(y(3)+y(11)) |
|
|
104 + .188821573960182*(y(4)+y(10)) |
|
|
105 + .199773405226859*(y(5)+y(9)) |
|
|
106 + .224926465333340*(y(6)+y(8)) |
|
|
107 + .242611071901408*y(7)); |
|
|
108 |
|
5837
|
109 s = sign(is); |
|
5838
|
110 |
|
|
111 if (s == 0) |
|
|
112 s = 1; |
|
5837
|
113 endif |
|
|
114 erri1 = abs(i1-is); |
|
|
115 erri2 = abs(i2-is); |
|
|
116 R = 1; |
|
|
117 if (erri2 != 0) |
|
|
118 R = erri1/erri2; |
|
|
119 endif |
|
|
120 if (R > 0 && R < 1) |
|
5838
|
121 tol = tol/R; |
|
5837
|
122 endif |
|
|
123 is = s*abs(is)*tol/eps; |
|
|
124 if (is == 0) |
|
|
125 is = b-a; |
|
|
126 endif |
|
5838
|
127 Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); |
|
5837
|
128 endfunction |
|
|
129 |
|
|
130 ## ADAPTLOBSTP Recursive function used by QUADL. |
|
|
131 ## |
|
5838
|
132 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to |
|
5837
|
133 ## approximate the integral of F(X) from A to B to |
|
|
134 ## an appropriate relative error. The argument 'F' is |
|
|
135 ## a string containing the name of f. The remaining |
|
|
136 ## arguments are generated by ADAPTLOB or by recursion. |
|
|
137 ## |
|
|
138 ## Walter Gautschi, 08/03/98 |
|
|
139 |
|
5838
|
140 function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) |
|
5837
|
141 h = (b-a)/2; |
|
|
142 m = (a+b)/2; |
|
|
143 alpha = sqrt(2/3); |
|
|
144 beta = 1/sqrt(5); |
|
|
145 mll = m-alpha*h; |
|
|
146 ml = m-beta*h; |
|
|
147 mr = m+beta*h; |
|
|
148 mrr = m+alpha*h; |
|
5838
|
149 x = [mll, ml, m, mr, mrr]; |
|
|
150 y = feval(f, x, varargin{:}); |
|
5837
|
151 fmll = y(1); |
|
|
152 fml = y(2); |
|
|
153 fm = y(3); |
|
|
154 fmr = y(4); |
|
|
155 fmrr = y(5); |
|
5838
|
156 i2 = (h/6)*(fa + fb + 5*(fml+fmr)); |
|
|
157 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm); |
|
|
158 if (is+(i1-i2) == is || mll <= a || b <= mrr) |
|
|
159 if ((m <= a || b <= m) && need_warning ()) |
|
|
160 warning ("quadl: interval contains no more machine number"); |
|
|
161 warning ("quadl: required tolerance may not be met"); |
|
|
162 need_warning (0); |
|
5837
|
163 endif |
|
|
164 Q = i1; |
|
|
165 if (trace) |
|
5838
|
166 disp ([a, b-a, Q]); |
|
5837
|
167 endif |
|
|
168 else |
|
5838
|
169 Q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) |
|
|
170 + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:}) |
|
|
171 + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:}) |
|
|
172 + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:}) |
|
|
173 + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:}) |
|
|
174 + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:})); |
|
5837
|
175 endif |
|
|
176 endfunction |
|
|
177 |
|
5838
|
178 function r = need_warning (v) |
|
5837
|
179 persistent w = []; |
|
|
180 if (nargin == 0) |
|
|
181 r = w; |
|
|
182 else |
|
|
183 w = v; |
|
|
184 endif |
|
|
185 endfunction |
