home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / QUADSTP.M < prev   
Encoding:
Text File  |  1993-03-07  |  1.6 KB  |  55 lines
  1. function [Q,cnt] = quadstp(FunFcn,a,b,tol,lev,fa,fc,fb,Q0,trace)
  2. %QUADSTP Recursive function used by QUAD.
  3. %    [Q,cnt] = quadstp(F,a,b,tol,lev,fa,fc,fb,Q0) tries to
  4. %    approximate the integral of f(x) from a to b to within a
  5. %    relative error of tol.  F is a string containing the name
  6. %    of f.  The remaining arguments are generated by quad or by
  7. %    the recursion.  lev is the recursion level.
  8. %    fa = f(a). fc = f((a+b)/2). fb = f(b).
  9. %    Q0 is an approximate value of the integral.
  10. %    See also QUAD and QUAD8.
  11.  
  12. %    C.B. Moler, 3-22-87.
  13. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  14.  
  15. LEVMAX = 10;
  16. if lev > LEVMAX
  17.     disp('Recursion level limit reached in quad.  Singularity likely.')
  18.     Q = Q0;
  19.     cnt = 0;
  20.     c = (a + b)/2;
  21.     if trace
  22.         yc = feval(FunFcn,c);
  23.         plot(c,yc,'*');    
  24.         if any(imag(yc))
  25.             plot(c,imag(yc),'*');
  26.         end
  27.     end
  28. else
  29.     % Evaluate function at midpoints of left and right half intervals.
  30.     h = b - a;
  31.     c = (a + b)/2;
  32.     x = [a+h/4, b-h/4];
  33.     f = feval(FunFcn,x);
  34.     if trace
  35.         plot(x,f,'.');
  36.         if any(imag(f))
  37.             plot(x,imag(f),'+')
  38.         end
  39.     end
  40.     cnt = 2;
  41.     
  42.     % Simpson's rule for half intervals.
  43.     Q1 = h*(fa + 4*f(1) + fc)/12;
  44.     Q2 = h*(fc + 4*f(2) + fb)/12;
  45.     Q = Q1 + Q2;
  46.     
  47.     % Recursively refine approximations.
  48.     if abs(Q - Q0) > tol*abs(Q)
  49.         [Q1,cnt1] = quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace);
  50.         [Q2,cnt2] = quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace);
  51.         Q = Q1 + Q2;
  52.         cnt = cnt + cnt1 + cnt2;
  53.     end
  54. end
  55.