home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / QUAD8.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.9 KB  |  60 lines
  1. function [Q,cnt] = quad8(funfcn,a,b,tol,trace)
  2. %QUAD8    Numerical evaluation of an integral, higher order method.
  3. %    Q = QUAD8('f',a,b) approximates the integral of f(x) from a
  4. %    to b to within a relative error of 1e-3.  'f' is a string
  5. %    containing the name of the function.  The function must return
  6. %    a vector of output values if given a vector of input values.
  7. %    Q = QUAD8(F,a,b,tol) integrates to a relative error of tol.
  8. %    Q = Inf is returned if an excessive recursion level is reached,
  9. %    indicating a possibly singular integral.
  10. %    Q = QUAD8(F,a,b,tol,trace) integrates to a relative error of tol and
  11. %    traces the function evaluations with a point plot of the integrand.
  12. %
  13. %    QUAD8 uses an adaptive recursive Newton Cotes 8 panel rule.  See
  14. %    also QUAD.
  15.  
  16. %    Cleve Moler, 5-08-88.
  17. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  18.  
  19. % [Q,cnt] = quad8(F,a,b,tol) also returns a function evaluation count.
  20.  
  21. if nargin < 4, tol = 1.e-3; trace = 0; end
  22. if nargin < 5, trace = 0; end
  23. % QUAD8 usually does better than the default 1e-3.
  24. h = b - a;
  25.  
  26. % Top level initialization, Newton-Cotes weights
  27. w = [3956 23552 -3712 41984 -18160 41984 -3712 23552 3956]/14175;
  28. x = a + (0:8)*(b-a)/8;
  29. y = feval(funfcn,x);
  30. yllow = min([min(real(y)) min(imag(y))]);
  31. ylhi = max([max(real(y)) max(imag(y))]);
  32. lims = [min(x) max(x) yllow ylhi];
  33. ind = find(~finite(lims));
  34. if ~isempty(ind)
  35.     [mind,nind] = size(ind);
  36.     lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
  37. end
  38. if trace
  39.     axis(lims);
  40. % doesn't take care of complex case
  41.     plot([a b],[real(y(1)) real(y(9))],'.'), hold on
  42.     if any(imag(y))
  43.         plot([a b],[imag(y(1)) imag(y(9))],'+')
  44.     end
  45. end
  46. lev = 1;
  47.  
  48. % Adaptive, recursive Newton-Cotes 8 panel quadrature
  49. if any(imag(y))
  50.     Q0 = 1e30;
  51. else
  52.     Q0 = inf;
  53. end
  54. [Q,cnt] = quad8stp(funfcn,a,b,tol,lev,w,x,y,Q0,trace);
  55. cnt = cnt + 9;
  56. if trace
  57.     hold off
  58.     axis('auto');
  59. end
  60.