Mercurial > hg > octave-nkf
annotate scripts/general/quadl.m @ 9245:16f53d29049f
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 22 May 2009 10:46:00 -0400 |
| parents | 923c7cb7f13f |
| children | 95c3e38098bf |
| rev | line source |
|---|---|
| 9245 | 1 ## Copyright (C) 1998, 2006, 2007, 2008, 2009 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 | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
27 ## @code{@var{f}(@var{x})} to machine precision. @var{f} is either a |
| 5837 | 28 ## function handle, inline function or string containing the name of |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
29 ## the function to evaluate. The function @var{f} must return a vector |
| 5837 | 30 ## of output values if given a vector of input values. |
| 31 ## | |
| 32 ## If defined, @var{tol} defines the relative tolerance to which to | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
33 ## which to integrate @code{@var{f}(@var{x})}. While if @var{trace} is |
| 5837 | 34 ## defined, displays the left end point of the current interval, the |
| 35 ## interval length, and the partial integral. | |
| 36 ## | |
|
9209
923c7cb7f13f
Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction.
Rik <rdrider0-list@yahoo.com>
parents:
9051
diff
changeset
|
37 ## Additional arguments @var{p1}, etc., are passed directly to @var{f}. |
| 5837 | 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 | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
66 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
67 myeps = eps ("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
68 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
69 myeps = eps; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
70 endif |
| 5838 | 71 if (isempty (tol)) |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
72 tol = myeps; |
| 5837 | 73 endif |
| 5838 | 74 if (isempty (trace)) |
| 75 trace = 0; | |
| 5837 | 76 endif |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
77 if (tol < myeps) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
78 tol = myeps; |
| 5837 | 79 endif |
| 80 | |
| 81 m = (a+b)/2; | |
| 82 h = (b-a)/2; | |
| 83 alpha = sqrt(2/3); | |
| 84 beta = 1/sqrt(5); | |
| 5838 | 85 |
| 5837 | 86 x1 = .942882415695480; |
| 87 x2 = .641853342345781; | |
| 88 x3 = .236383199662150; | |
| 5838 | 89 |
| 90 x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ... | |
| 91 m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b]; | |
| 92 | |
| 93 y = feval (f, x, varargin{:}); | |
| 94 | |
| 5837 | 95 fa = y(1); |
| 96 fb = y(13); | |
| 5838 | 97 |
| 98 i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9))); | |
| 99 | |
| 100 i1 = (h/1470)*(77*(y(1)+y(13)) | |
| 101 + 432*(y(3)+y(11)) | |
| 102 + 625*(y(5)+y(9)) | |
| 103 + 672*y(7)); | |
| 104 | |
| 105 is = h*(.0158271919734802*(y(1)+y(13)) | |
| 106 +.0942738402188500*(y(2)+y(12)) | |
| 107 + .155071987336585*(y(3)+y(11)) | |
| 108 + .188821573960182*(y(4)+y(10)) | |
| 109 + .199773405226859*(y(5)+y(9)) | |
| 110 + .224926465333340*(y(6)+y(8)) | |
| 111 + .242611071901408*y(7)); | |
| 112 | |
| 5837 | 113 s = sign(is); |
| 5838 | 114 |
| 115 if (s == 0) | |
| 116 s = 1; | |
| 5837 | 117 endif |
| 118 erri1 = abs(i1-is); | |
| 119 erri2 = abs(i2-is); | |
| 120 R = 1; | |
| 121 if (erri2 != 0) | |
| 122 R = erri1/erri2; | |
| 123 endif | |
| 124 if (R > 0 && R < 1) | |
| 5838 | 125 tol = tol/R; |
| 5837 | 126 endif |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
127 is = s*abs(is)*tol/myeps; |
| 5837 | 128 if (is == 0) |
| 129 is = b-a; | |
| 130 endif | |
| 5838 | 131 Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); |
| 5837 | 132 endfunction |
| 133 | |
| 134 ## ADAPTLOBSTP Recursive function used by QUADL. | |
| 135 ## | |
| 5838 | 136 ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to |
| 5837 | 137 ## approximate the integral of F(X) from A to B to |
| 138 ## an appropriate relative error. The argument 'F' is | |
| 139 ## a string containing the name of f. The remaining | |
| 140 ## arguments are generated by ADAPTLOB or by recursion. | |
| 141 ## | |
| 142 ## Walter Gautschi, 08/03/98 | |
| 143 | |
| 5838 | 144 function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) |
| 5837 | 145 h = (b-a)/2; |
| 146 m = (a+b)/2; | |
| 147 alpha = sqrt(2/3); | |
| 148 beta = 1/sqrt(5); | |
| 149 mll = m-alpha*h; | |
| 150 ml = m-beta*h; | |
| 151 mr = m+beta*h; | |
| 152 mrr = m+alpha*h; | |
| 5838 | 153 x = [mll, ml, m, mr, mrr]; |
| 154 y = feval(f, x, varargin{:}); | |
| 5837 | 155 fmll = y(1); |
| 156 fml = y(2); | |
| 157 fm = y(3); | |
| 158 fmr = y(4); | |
| 159 fmrr = y(5); | |
| 5838 | 160 i2 = (h/6)*(fa + fb + 5*(fml+fmr)); |
| 161 i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm); | |
| 162 if (is+(i1-i2) == is || mll <= a || b <= mrr) | |
| 163 if ((m <= a || b <= m) && need_warning ()) | |
| 164 warning ("quadl: interval contains no more machine number"); | |
| 165 warning ("quadl: required tolerance may not be met"); | |
| 166 need_warning (0); | |
| 5837 | 167 endif |
| 168 Q = i1; | |
| 169 if (trace) | |
| 5838 | 170 disp ([a, b-a, Q]); |
| 5837 | 171 endif |
| 172 else | |
| 5838 | 173 Q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) |
| 174 + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:}) | |
| 175 + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:}) | |
| 176 + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:}) | |
| 177 + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:}) | |
| 178 + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:})); | |
| 5837 | 179 endif |
| 180 endfunction | |
| 181 | |
| 5838 | 182 function r = need_warning (v) |
| 5837 | 183 persistent w = []; |
| 184 if (nargin == 0) | |
| 185 r = w; | |
| 186 else | |
| 187 w = v; | |
| 188 endif | |
| 189 endfunction |