Mercurial > hg > octave-lyh
changeset 13008:85dac13a911b
quadl.m: Fix integration with large error tolerances (Bug #33792)
quadl.m: Force recursion to occur at least once through global state variable.
Miscellaneous spacing changes to improve code appearance.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Thu, 25 Aug 2011 09:51:27 -0700 |
parents | 61be447052c3 |
children | 38b52a073cfa |
files | scripts/general/quadl.m |
diffstat | 1 files changed, 53 insertions(+), 54 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/general/quadl.m +++ b/scripts/general/quadl.m @@ -62,14 +62,12 @@ ## * replace global variable terminate2 with local function need_warning ## * add paper ref to docs -function q = quadl (f, a, b, tol, trace, varargin) - need_warning (1); - if (nargin < 4) - tol = []; +function q = quadl (f, a, b, tol = [], trace = false, varargin) + + if (nargin < 3) + print_usage (); endif - if (nargin < 5) - trace = []; - endif + if (isa (a, "single") || isa (b, "single")) myeps = eps ("single"); else @@ -79,16 +77,23 @@ tol = myeps; endif if (isempty (trace)) - trace = 0; + trace = false; endif if (tol < myeps) tol = myeps; endif + ## Track whether recursion has occurred + global __quadl_recurse_done__; + __quadl_recurse_done__ = false; + ## Track whether warning about machine precision has been issued + global __quadl_need_warning__; + __quadl_need_warning__ = true; + m = (a+b)/2; h = (b-a)/2; - alpha = sqrt(2/3); - beta = 1/sqrt(5); + alpha = sqrt (2/3); + beta = 1/sqrt (5); x1 = .942882415695480; x2 = .641853342345781; @@ -104,12 +109,12 @@ i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9))); - i1 = (h/1470)*(77*(y(1)+y(13)) + i1 = (h/1470)*( 77*(y(1)+y(13)) + 432*(y(3)+y(11)) + 625*(y(5)+y(9)) + 672*y(7)); - is = h*(.0158271919734802*(y(1)+y(13)) + is = h*( .0158271919734802*(y(1)+y(13)) +.0942738402188500*(y(2)+y(12)) + .155071987336585*(y(3)+y(11)) + .188821573960182*(y(4)+y(10)) @@ -117,102 +122,96 @@ + .224926465333340*(y(6)+y(8)) + .242611071901408*y(7)); - s = sign(is); - + s = sign (is); if (s == 0) s = 1; endif - erri1 = abs(i1-is); - erri2 = abs(i2-is); - R = 1; + erri1 = abs (i1-is); + erri2 = abs (i2-is); if (erri2 != 0) R = erri1/erri2; + else + R = 1; endif if (R > 0 && R < 1) tol = tol/R; endif - is = s*abs(is)*tol/myeps; + is = s * abs(is) * tol/myeps; if (is == 0) is = b-a; endif + q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin{:}); + endfunction ## ADAPTLOBSTP Recursive function used by QUADL. ## ## Q = ADAPTLOBSTP('F', A, B, FA, FB, IS, TRACE) tries to ## approximate the integral of F(X) from A to B to -## an appropriate relative error. The argument 'F' is +## an appropriate relative error. The argument 'F' is ## a string containing the name of f. The remaining ## arguments are generated by ADAPTLOB or by recursion. ## ## Walter Gautschi, 08/03/98 function q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) + global __quadl_recurse_done__; + global __quadl_need_warning__; + h = (b-a)/2; m = (a+b)/2; - alpha = sqrt(2/3); - beta = 1/sqrt(5); + alpha = sqrt (2/3); + beta = 1 / sqrt(5); mll = m-alpha*h; - ml = m-beta*h; - mr = m+beta*h; + ml = m-beta*h; + mr = m+beta*h; mrr = m+alpha*h; x = [mll, ml, m, mr, mrr]; - y = feval(f, x, varargin{:}); + y = feval (f, x, varargin{:}); fmll = y(1); - fml = y(2); - fm = y(3); - fmr = y(4); + fml = y(2); + fm = y(3); + fmr = y(4); fmrr = y(5); i2 = (h/6)*(fa + fb + 5*(fml+fmr)); i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm); - if (is+(i1-i2) == is || mll <= a || b <= mrr) - if ((m <= a || b <= m) && need_warning ()) + if ((is+(i1-i2) == is || mll <= a || b <= mrr) && __quadl_recurse_done__) + if ((m <= a || b <= m) && __quadl_need_warning__) warning ("quadl: interval contains no more machine number"); warning ("quadl: required tolerance may not be met"); - need_warning (0); + __quadl_need_warning__ = false; endif q = i1; if (trace) disp ([a, b-a, q]); endif else - q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:}) - + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:}) - + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:}) - + adaptlobstp (f, m, mr, fm, fmr, is, trace, varargin{:}) - + adaptlobstp (f, mr, mrr, fmr, fmrr, is, trace, varargin{:}) - + adaptlobstp (f, mrr, b, fmrr, fb, is, trace, varargin{:})); - endif -endfunction - -function r = need_warning (v) - persistent w = []; - if (nargin == 0) - r = w; - else - w = v; + __quadl_recurse_done__ = true; + q = ( adaptlobstp (f, a , mll, fa , fmll, is, trace, varargin{:}) + + adaptlobstp (f, mll, ml , fmll, fml , is, trace, varargin{:}) + + adaptlobstp (f, ml , m , fml , fm , is, trace, varargin{:}) + + adaptlobstp (f, m , mr , fm , fmr , is, trace, varargin{:}) + + adaptlobstp (f, mr , mrr, fmr , fmrr, is, trace, varargin{:}) + + adaptlobstp (f, mrr, b , fmrr, fb , is, trace, varargin{:})); endif endfunction ## basic functionality -%!assert( quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16) +%!assert (quadl (@(x) sin (x), 0, pi, [], []), 2, -3e-16) ## the values here are very high so it may be unavoidable that this fails -%!assert ( quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15), +%!assert (quadl (@(x) sin (3*x).*cosh (x).*sinh (x),10,15), %! 2.588424538641647e+10, -9e-15) ## extra parameters %!assert (quadl (@(x,a,b) sin (a + b*x), 0, 1, [], [], 2, 3), %! cos(2)/3 - cos(5)/3, - 3e-16) -## test different tolerances. This test currently fails for a very high -## tolerances. -%!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []), +## test different tolerances. +%!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, []), %! (60 + sin(4) - sin(64))/12, -0.3) +%!assert (quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []), +%! (60 + sin(4) - sin(64))/12, -0.1) - -## for lower tolerances the test passes. -%!assert ( quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.1, []), -%! (60 + sin(4) - sin(64))/12, -0.1) \ No newline at end of file