home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / qtrap.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  827 b   |  39 lines
  1. #
  2. # qtrap
  3. #
  4. # Simple interface for trapzd. Returns estimate of integral of f from a to b,
  5. # using the trapezoidal rule.
  6. #
  7. # The (global) variable QTRAP_EPS can be set to desired fractional accuracy
  8. # and QTRAP_JMAX so that 2^(QTRAP_JMAX-1) is maximum allowed number of steps.
  9. #
  10. # Error 100 indicates that QTRAP_JMAX is too low to get desired accuracy.
  11. #
  12. # This routine is from "Numerical Recipes in C" by Press et. al.
  13. # It uses trapzd(), defined in trapzd.ic.
  14. #
  15. # mws, 7/92. 
  16. #
  17. silent
  18. echo "Reading trapzd..."
  19. read "trapzd.ic"
  20. echo "Defining qtrap..."
  21. #
  22. QTRAP_EPS = 1.0e-5
  23. QTRAP_JMAX = 20
  24. #
  25. func qtrap(~f,a,b) = {
  26.     local j,s,olds
  27.  
  28.     olds = -1.0e30;
  29.     for (j = 1; j <= QTRAP_JMAX; j += 1) {
  30.         s = trapzd(f,a,b,j)
  31.         if (abs(s-olds) < QTRAP_EPS*abs(olds)) return s
  32.         olds = s
  33.     }
  34.     error(99)    # too many steps
  35. }
  36.  
  37. echo "Done."
  38. verbose
  39.