|
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 |
|
|
59 function Q = quadl(f,a,b,tol,trace,varargin) |
|
|
60 need_warning(1); |
|
|
61 if (nargin < 4) |
|
|
62 tol=[]; |
|
|
63 endif |
|
|
64 if (nargin < 5) |
|
|
65 trace = []; |
|
|
66 endif |
|
|
67 if (isempty(tol)) |
|
|
68 tol = eps; |
|
|
69 endif |
|
|
70 if (isempty(trace)) |
|
|
71 trace=0; |
|
|
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); |
|
|
81 x1 = .942882415695480; |
|
|
82 x2 = .641853342345781; |
|
|
83 x3 = .236383199662150; |
|
|
84 x = [a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,... |
|
|
85 m+beta*h,m+x2*h,m+alpha*h,m+x1*h,b]; |
|
|
86 y = feval(f,x,varargin{:}); |
|
|
87 fa = y(1); |
|
|
88 fb = y(13); |
|
|
89 i2 = (h/6)*(y(1)+y(13)+5*(y(5)+y(9))); |
|
|
90 i1 = (h/1470)*(77*(y(1)+y(13))+432*(y(3)+y(11))+ ... |
|
|
91 625*(y(5)+y(9))+672*y(7)); |
|
|
92 is = h*(.0158271919734802*(y(1)+y(13))+.0942738402188500 ... |
|
|
93 *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ... |
|
|
94 .188821573960182*(y(4)+y(10))+.199773405226859 ... |
|
|
95 *(y(5)+y(9))+.224926465333340*(y(6)+y(8))... |
|
|
96 +.242611071901408*y(7)); |
|
|
97 s = sign(is); |
|
|
98 if (s == 0), |
|
|
99 s=1; |
|
|
100 endif |
|
|
101 erri1 = abs(i1-is); |
|
|
102 erri2 = abs(i2-is); |
|
|
103 R = 1; |
|
|
104 if (erri2 != 0) |
|
|
105 R = erri1/erri2; |
|
|
106 endif |
|
|
107 if (R > 0 && R < 1) |
|
|
108 tol=tol/R; |
|
|
109 endif |
|
|
110 is = s*abs(is)*tol/eps; |
|
|
111 if (is == 0) |
|
|
112 is = b-a; |
|
|
113 endif |
|
|
114 Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:}); |
|
|
115 endfunction |
|
|
116 |
|
|
117 ## ADAPTLOBSTP Recursive function used by QUADL. |
|
|
118 ## |
|
|
119 ## Q = ADAPTLOBSTP('F',A,B,FA,FB,IS,TRACE) tries to |
|
|
120 ## approximate the integral of F(X) from A to B to |
|
|
121 ## an appropriate relative error. The argument 'F' is |
|
|
122 ## a string containing the name of f. The remaining |
|
|
123 ## arguments are generated by ADAPTLOB or by recursion. |
|
|
124 ## |
|
|
125 ## Walter Gautschi, 08/03/98 |
|
|
126 |
|
|
127 function Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin) |
|
|
128 h = (b-a)/2; |
|
|
129 m = (a+b)/2; |
|
|
130 alpha = sqrt(2/3); |
|
|
131 beta = 1/sqrt(5); |
|
|
132 mll = m-alpha*h; |
|
|
133 ml = m-beta*h; |
|
|
134 mr = m+beta*h; |
|
|
135 mrr = m+alpha*h; |
|
|
136 x = [mll,ml,m,mr,mrr]; |
|
|
137 y = feval(f,x,varargin{:}); |
|
|
138 fmll = y(1); |
|
|
139 fml = y(2); |
|
|
140 fm = y(3); |
|
|
141 fmr = y(4); |
|
|
142 fmrr = y(5); |
|
|
143 i2 = (h/6)*(fa+fb+5*(fml+fmr)); |
|
|
144 i1 = (h/1470)*(77*(fa+fb)+432*(fmll+fmrr)+625*(fml+fmr) ... |
|
|
145 +672*fm); |
|
|
146 if ((is+(i1-i2) == is) || (mll <= a) || (b <= mrr)) |
|
|
147 if (((m <= a) || (b <= m)) && need_warning()) |
|
|
148 warning(['Interval contains no more machine number. ',... |
|
|
149 'Required tolerance may not be met.']); |
|
|
150 need_warning(0); |
|
|
151 endif |
|
|
152 Q = i1; |
|
|
153 if (trace) |
|
|
154 disp([a b-a Q]); |
|
|
155 endif |
|
|
156 else |
|
|
157 Q = adaptlobstp(f,a,mll,fa,fmll,is,trace,varargin{:})+... |
|
|
158 adaptlobstp(f,mll,ml,fmll,fml,is,trace,varargin{:})+... |
|
|
159 adaptlobstp(f,ml,m,fml,fm,is,trace,varargin{:})+... |
|
|
160 adaptlobstp(f,m,mr,fm,fmr,is,trace,varargin{:})+... |
|
|
161 adaptlobstp(f,mr,mrr,fmr,fmrr,is,trace,varargin{:})+... |
|
|
162 adaptlobstp(f,mrr,b,fmrr,fb,is,trace,varargin{:}); |
|
|
163 endif |
|
|
164 endfunction |
|
|
165 |
|
|
166 function r = need_warning(v) |
|
|
167 persistent w = []; |
|
|
168 if (nargin == 0) |
|
|
169 r = w; |
|
|
170 else |
|
|
171 w = v; |
|
|
172 endif |
|
|
173 endfunction |
