home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / qsimp.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  818 b   |  41 lines
  1. #
  2. # qsimp
  3. #
  4. # Simple interface for trapzd. Returns estimate of integral of f from a to b,
  5. # using Simpson's rule.
  6. #
  7. # The (global) variable QSIMP_EPS can be set to desired fractional accuracy
  8. # and QSIMP_JMAX so that 2^(QSIMP_JMAX-1) is maximum allowed number of steps.
  9. #
  10. # Error 101 indicates that QSIMP_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 qsimp..."
  21. #
  22. QSIMP_EPS = 1e-6
  23. QSIMP_JMAX = 20
  24. #
  25. func qsimp(~f,a,b) = {
  26.     local j,s,st,ost,os
  27.  
  28.     ost = os = -1e30
  29.     for (j=1;j<=QSIMP_JMAX;j+=1) {
  30.         st=trapzd(f,a,b,j)
  31.         s=(4*st-ost)/3
  32.         if (abs(s-os) < QSIMP_EPS*abs(os)) return s
  33.         os=s
  34.         ost=st
  35.     }
  36.     error(101)
  37. }
  38.  
  39. echo "Done."
  40. verbose
  41.