home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / FUNFUN.DI$ / QUAD8STP.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.5 KB  |  47 lines
  1. function [Q,cnt] = quad8stp(FunFcn,a,b,tol,lev,w,x0,f0,Q0,trace)
  2. %QUAD8STP Recursive function used by QUAD8.
  3. %    [Q,cnt] = quad8stp(F,a,b,tol,lev,w,f,Q0) tries to approximate
  4. %    the integral of f(x) from a to b to within a relative error of tol.  
  5. %    F is a string containing the name of f.  The remaining arguments are 
  6. %    generated by quad8 or by the recursion.  lev is the recursion level.
  7. %    w is the weights in the 8 panel Newton Cotes formula.
  8. %    x0 is a vector of 9 equally spaced abscissa is the interval.
  9. %    f0 is a vector of the 9 function values at x.
  10. %    Q0 is an approximate value of the integral.
  11. %    See also QUAD8 and QUAD.
  12.  
  13. %    Cleve Moler, 5-08-88.
  14. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  15.  
  16. LEVMAX = 10;
  17.  
  18. % Evaluate function at midpoints of left and right half intervals.
  19. x = zeros(1,17);
  20. f = zeros(1,17);
  21. x(1:2:17) = x0;
  22. f(1:2:17) = f0;
  23. x(2:2:16) = (x0(1:8) + x0(2:9))/2;
  24. f(2:2:16) = feval(FunFcn,x(2:2:16));
  25. if trace
  26.     plot(x(2:2:16),f(2:2:16),'.');
  27.     if any(imag(f))
  28.         plot(x(2:2:16),imag(f(2:2:16)),'+');
  29.     end
  30. end
  31. cnt = 8;
  32.  
  33. % Integrate over half intervals.
  34. h = (b-a)/16;
  35. Q1 = h*w*f(1:9).';
  36. Q2 = h*w*f(9:17).';
  37. Q = Q1 + Q2;
  38.  
  39. % Recursively refine approximations.
  40. if abs(Q - Q0) > tol*abs(Q) & lev <= LEVMAX
  41.    c = (a+b)/2;
  42.    [Q1,cnt1] = quad8stp(FunFcn,a,c,tol/2,lev+1,w,x(1:9),f(1:9),Q1,trace);
  43.    [Q2,cnt2] = quad8stp(FunFcn,c,b,tol/2,lev+1,w,x(9:17),f(9:17),Q2,trace);
  44.    Q = Q1 + Q2;
  45.    cnt = cnt + cnt1 + cnt2;
  46. end
  47.