share/numeric/dblint.usg

   1 The file SHARE1;DBLINT FASL contains a double-integral routine which was
   2 written in top-level macsyma and then translated and compiled to
   3 machine code. It uses the Simpson's Rule method (see note below) in
   4 both the x and y directions to calculate
   5
   6 /B /S(X)
   7 |  |
   8 |  |    F(X,Y) DY DX .
   9 |  |
  10 /A /R(X)
  11
  12 The function F(X,Y) must be a translated or compiled function of two
  13 variables, and R(X) and S(X) must each be a translated or compiled
  14 function of one variable, while A and B must be floating point
  15 numbers. The routine has two global variables which determine the
  16 number of divisions of the x and y intervals: DBLINT_X and DBLINT_Y,
  17 both of which are initially 10, and can be changed independently to
  18 other integer values [there are 2*DBLINT_X+1 points computed in the x
  19 direction, and 2*DBLINT_Y+1 in the y direction].
  20
  21 The routine subdivides the x axis and then for each value of x it
  22 first computes r(x) and s(x); then the y axis between r(x) and s(x) is
  23 subdivided and the integral along the y axis is performed using
  24 Simpson's Rule; then the integral along the x axis is done using
  25 Simpson's Rule with the function values being the y-integrals. This
  26 procedure may be numerically unstable for a great variety of reasons,
  27 but is reasonably fast: avoid using it on highly oscillatory functions
  28 and functions with singularities (poles or branch points in the
  29 region).  The y integrals depend on how far apart r(x) and s(x) are,
  30 so if the distance s(x)-r(x) varies rapidly with x, there may be
  31 substantial errors arising from truncation with different step-sizes
  32 in the various y integrals. One can increase dblint_x and dblint_y in
  33 an effort to improve the coverage of the region, at the expense of
  34 computation time. The function values are not saved, so if the
  35 function is very time-consuming, you will have to wait for
  36 re-computation if you change anything (sorry).
  37
  38 It is required that the functions F, R, and S be either translated or
  39 compiled prior to calling DBLINT. This will result in orders of
  40 magnitude speed improvement over interpreted code in many cases! Ask
  41 LPH (or GJC) about using these numerical aids.  The file SHARE1;DBLINT
  42 DEMO can be run in batch or demo mode to illustrate the usage on a
  43 sample problem.
  44 Please send all bug notes and questions to LPH@MIT-MC.
  45
  46 Note: Simpson's Rule specifies that
  47
  48 /X[2*N]
  49 |
  50 |       F(X) DX = H/3* (F(X[0]) +
  51 |
  52 /X[0]                   4*(F(X[1])+F(X[3])+...+F(X[2*N-1])) +
  53
  54                         2*(F(X[2])+F(X[4])+...+F(X[2*N-2])) +
  55
  56                         F(X[2*N]))
  57
  58 in one dimension, where H is the distance between the equally spaced
  59 X[N]'s, and DBLINT_X=N. The error in this formulation is of order
  60 H^5*N*DIFF(F(X),X,4) for some X in (X[0],X[2*N]).