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