home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / QUAD.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.8 KB  |  64 lines
  1. function [Q,cnt] = quad(funfcn,a,b,tol,trace)
  2. %QUAD    Numerical evaluation of an integral, low order method.
  3. %    Q = QUAD('f',a,b) approximates the integral of f(x) from a to b
  4. %    to within a relative error of 1e-3.  'f' is a string
  5. %    containing the name of the function.  Function f must return a
  6. %    vector of output values if given a vector of input values.
  7. %    Q = QUAD(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 = QUAD(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. %    QUAD uses an adaptive recursive Simpson's rule.
  14. %
  15. %    See also QUAD8.
  16.  
  17. %    C.B. Moler, 3-22-87.
  18. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  19.  
  20. % [Q,cnt] = quad(F,a,b,tol) also returns a function evaluation count.
  21.  
  22. if nargin < 4, tol = 1.e-3; trace = 0; end
  23. if nargin < 5, trace = 0; end
  24. c = (a + b)/2;
  25.  
  26. % Top level initialization
  27. x = [a b c a:(b-a)/10:b];
  28. y = feval(funfcn,x);
  29. fa = y(1);
  30. fb = y(2);
  31. fc = y(3);
  32. if trace
  33.     lims = [min(x) max(x) min(y) max(y)];
  34.     if any(imag(y))
  35.         lims(3) = min(min(real(y)),min(imag(y)));
  36.         lims(4) = max(max(real(y)),max(imag(y)));
  37.     end
  38.     ind = find(~finite(lims));
  39.     if ~isempty(ind)
  40.         [mind,nind] = size(ind);
  41.         lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2));
  42.     end
  43.     axis(lims);
  44.     plot([c b],[fc fb],'.')
  45.     hold on
  46.     if any(imag(y))
  47.         plot([c b],imag([fc fb]),'+')
  48.     end
  49. end
  50. lev = 1;
  51.  
  52. % Adaptive, recursive Simpson's quadrature
  53. if any(imag(y))
  54.     Q0 = 1e30;
  55. else
  56.     Q0 = inf;
  57. end
  58. [Q,cnt] = quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace);
  59. cnt = cnt + 3;
  60. if trace
  61.     hold off
  62.     axis('auto');
  63. end
  64.