home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / Integration.notes < prev    next >
Encoding:
Text File  |  1992-08-03  |  4.8 KB  |  153 lines
  1. =============================================================================
  2.  
  3.     icalc - a complex-number expression language
  4.  
  5.     (C) 1991, 1992 Martin W. Scott. All Rights Reserved.
  6.  
  7. =============================================================================
  8.  
  9.     Integration Notes
  10.  
  11. =============================================================================
  12.  
  13. There are a few icalc scripts that deal with the problem of real
  14. integration. These scripts provide 4 basic routines which vary in their
  15. sophistication, and also in their speed of calculation. They are:
  16.  
  17.     gl.ic        - contains definition of function gl, which
  18.               calculates integrals using a 10-point
  19.               Gauss-Legendre quadrature.
  20.  
  21.     trapzd.ic    - contains definition of function trapzd, used
  22.               by qtrap, qsimp and qromb routines.
  23.  
  24.     polint.ic    - contains definition of function polint, used
  25.               by qromb, but useful in its own right.
  26.  
  27.     qtrap.ic    - contains definition of function qtrap, which
  28.               performs numerical integration using the
  29.               Trapezium Rule.
  30.  
  31.     qsimp.ic    - contains definition of function qtrap, which
  32.               performs numerical integration using the
  33.               Simpson's Rule.
  34.  
  35.     qromb.ic    - contains definition of function qtrap, which
  36.               performs Romberg numerical integration 
  37.  
  38. All the routines assume that the array base is 1, i.e. given an array a,
  39. its first element is a[1]. If you have changed the array base (with a
  40. call to arraybase()) you should restore its value to 1 before using
  41. these routines.
  42.  
  43.  
  44.  
  45. How to use the integration routines
  46. -----------------------------------
  47. Start icalc with the desired routine, from the Workbench or a Shell.
  48. Here, xxx denotes gl, qromb, qtrap or qsimp.
  49.  
  50.     From Shell:    type "icalc xxx.ic -" at shell prompt
  51.  
  52.     From WB:    click once on icalc icon, hold down shift
  53.             and double-click on xxx.ic icon.
  54.  
  55. All the routines gl, qtrap, qsimp and qromb take 3 arguments - an
  56. expression in x, and the interval endpoints a and b.
  57.  
  58. The expression in x can be any expression at all; here are some
  59. examples:
  60.  
  61.     > gl(x*x,0,1)
  62.         0.33333333328
  63.  
  64.     > gl(sqr(x),0,1)
  65.         0.33333333328
  66.  
  67.     > gl(sin(ln(x)),1,2)
  68.         0.36972237492
  69.  
  70.     > qtrap(sin(ln(x)),1,2)
  71.         0.36972159245
  72.  
  73.     > qsimp(log(cos(x)),0,PI/4)
  74.         -0.03752900436
  75.  
  76.     > 2*sqr(qromb(exp(-sqr(x)/2), 0, 100))
  77.         3.14159314113
  78.  
  79. and so on. The question remains, when to use which routine?
  80.  
  81. The qtrap routine may be used for integrating functions which are not
  82. very smooth; since most expression-based functions are, you should
  83. hardly ever need to use it. It is the slowest of the methods, and can
  84. take many seconds to produce a result.
  85.  
  86. The qsimp routine is usually more efficient than qtrap when the
  87. integrand has a continuous 3rd derivative.
  88.  
  89. The qromb routine is best for sufficiently smooth functions (e.g.
  90. analytic functions) and is usually fastest at producing a result.
  91.  
  92. The gl routine is useful for quick results, although usually qromb is
  93. quick enough (about a second or two) for most integrations. gl is useful
  94. when used in a routine for multi-dimensional integration.
  95.  
  96. For more information on all routines, consult "Numerical Recipes in C",
  97. by Press et al. That is the source of these algorithms.
  98.  
  99.  
  100.  
  101. Improving Accuracy
  102. ------------------
  103. You can improve the accuracy of results by setting the global variables
  104. controlling whichever routine you are using (this does not apply to gl).
  105. The variables, with their default values, are as follows:
  106.  
  107.     Routine        Variable    Value
  108.     -------        --------    -----
  109.     qtrap        QTRAP_EPS    1e-5
  110.     qsimp        QSIMP_EPS    1e-6
  111.     qromb        QROMB_EPS    1e-6
  112.  
  113. The xxx_EPS variables denote the fractional accuracy desired, i.e. if
  114. the (error estimate ds) < xxx_EPS*(integral estimate s) then the routine
  115. returns the integral estimate, s.
  116.  
  117. In some cases, desired accuracy cannot be achieved within the maximum
  118. allowed iterations, causing a user error. In such cases, try another
  119. routine.
  120.  
  121. Finally, a warning: because of the type of termination condition
  122. (fractional accuracy) integrations with a result of zero cause problems.
  123. If your routine is taking an excessively long time (more than 10 seconds
  124. for qromb, say) it may be because the result is zero. Try splitting the
  125. integration limits up, and doing two integrations. If this fails, graph
  126. the function to see what is going on. The routines may also have
  127. problems with functions containing singularities.
  128.  
  129. Examples (all using best routine, qromb):
  130.  
  131.     > qromb(sin(x),0,PI)
  132.         1.99999999459
  133.     > QROMB_EPS=1e-8
  134.         0.00000001
  135.     > qromb(sin(x),0,PI)
  136.         2
  137.  
  138.     > qromb(sin(x),0,2*PI)
  139.  
  140. No result appears for some time. Press control-c to abort.
  141.  
  142.     <CONTROL>-C
  143.     *** break
  144.     > qromb(sin(x),0,PI) + qromb(sin(x),PI,2*PI)
  145.         -4.21884749358e-15
  146.  
  147. The last integration was very close to zero, which is why the
  148. qromb(sin(x),0,2*PI) call was taking so long.
  149.  
  150. That's about all I want to say about these routines. The scripts
  151. containing them are partially-commented, but you'd best get a hold of
  152. Numerical Recipes if you need more information.
  153.