Mercurial > hg > octave-lyh
diff scripts/general/quadl.m @ 5838:376e02b2ce70
[project @ 2006-06-01 20:23:53 by jwe]
author | jwe |
---|---|
date | Thu, 01 Jun 2006 20:23:54 +0000 |
parents | 55404f3b0da1 |
children | 93c65f2a5668 |
line wrap: on
line diff
--- a/scripts/general/quadl.m +++ b/scripts/general/quadl.m @@ -56,19 +56,19 @@ ## * 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); +function Q = quadl (f, a, b, tol, trace, varargin) + need_warning (1); if (nargin < 4) - tol=[]; + tol = []; endif if (nargin < 5) trace = []; endif - if (isempty(tol)) + if (isempty (tol)) tol = eps; endif - if (isempty(trace)) - trace=0; + if (isempty (trace)) + trace = 0; endif if (tol < eps) tol = eps; @@ -78,25 +78,38 @@ h = (b-a)/2; alpha = sqrt(2/3); beta = 1/sqrt(5); + x1 = .942882415695480; x2 = .641853342345781; x3 = .236383199662150; - x = [a,m-x1*h,m-alpha*h,m-x2*h,m-beta*h,m-x3*h,m,m+x3*h,... - m+beta*h,m+x2*h,m+alpha*h,m+x1*h,b]; - y = feval(f,x,varargin{:}); + + x = [a, m-x1*h, m-alpha*h, m-x2*h, m-beta*h, m-x3*h, m, m+x3*h, ... + m+beta*h, m+x2*h, m+alpha*h, m+x1*h, b]; + + y = feval (f, x, varargin{:}); + fa = y(1); fb = y(13); - i2 = (h/6)*(y(1)+y(13)+5*(y(5)+y(9))); - 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))+.0942738402188500 ... - *(y(2)+y(12))+.155071987336585*(y(3)+y(11))+ ... - .188821573960182*(y(4)+y(10))+.199773405226859 ... - *(y(5)+y(9))+.224926465333340*(y(6)+y(8))... - +.242611071901408*y(7)); + + i2 = (h/6)*(y(1) + y(13) + 5*(y(5)+y(9))); + + 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)) + +.0942738402188500*(y(2)+y(12)) + + .155071987336585*(y(3)+y(11)) + + .188821573960182*(y(4)+y(10)) + + .199773405226859*(y(5)+y(9)) + + .224926465333340*(y(6)+y(8)) + + .242611071901408*y(7)); + s = sign(is); - if (s == 0), - s=1; + + if (s == 0) + s = 1; endif erri1 = abs(i1-is); erri2 = abs(i2-is); @@ -105,18 +118,18 @@ R = erri1/erri2; endif if (R > 0 && R < 1) - tol=tol/R; + tol = tol/R; endif is = s*abs(is)*tol/eps; if (is == 0) is = b-a; endif - Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin{:}); + 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 +## 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 ## a string containing the name of f. The remaining @@ -124,7 +137,7 @@ ## ## Walter Gautschi, 08/03/98 -function Q = adaptlobstp(f,a,b,fa,fb,is,trace,varargin) +function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin) h = (b-a)/2; m = (a+b)/2; alpha = sqrt(2/3); @@ -133,37 +146,36 @@ ml = m-beta*h; mr = m+beta*h; mrr = m+alpha*h; - x = [mll,ml,m,mr,mrr]; - y = feval(f,x,varargin{:}); + x = [mll, ml, m, mr, mrr]; + y = feval(f, x, varargin{:}); fmll = y(1); 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()) - warning(['Interval contains no more machine number. ',... - 'Required tolerance may not be met.']); - need_warning(0); + 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 ()) + warning ("quadl: interval contains no more machine number"); + warning ("quadl: required tolerance may not be met"); + need_warning (0); endif Q = i1; if (trace) - disp([a b-a Q]); + 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{:}); + 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) +function r = need_warning (v) persistent w = []; if (nargin == 0) r = w;