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	Thu Jun 01 19:05:32 2006 +0000
+++ b/scripts/general/quadl.m	Thu Jun 01 20:23:54 2006 +0000
@@ -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;