archive/share/trash/qq.usg

   1
   2 The file SHARE1;QQ FASL contains a function QUANC8 which can take either
   3 3 or 4 arguments. The 3 arg version computes the integral of the function
   4 specified as the first argument over the interval from lo to hi as in
   5 QUANC8('function name,lo,hi).  The function name should be quoted.
   6 The 4 arg version will compute the integral of the function or expression
   7 (first arg) with respect to the variable (second arg) over the interval
   8 from lo to hi as in QUANC8(<f(x) or expression in x>,x,lo,hi).
   9
  10 The method used is the Newton-Cotes 8th order polynomial quadrature, and
  11 the routine is adaptive. It will thus spend time dividing the interval only
  12 when necessary to achieve the error conditions specified by the global
  13 variables QUANC8_RELERR (default value=1.0e-4) and QUANC8_ABSERR
  14 (default value=1.0e-8) which give
  15 the relative error test
  16 |integral(function)-computed value|< quanc8_relerr*|integral(function)|
  17 and the absolute error test
  18 |integral(function)-computed value|<quanc8_abserr.
  19
  20 The error from each subinterval is estimated and the contribution from a
  21 subinterval is accepted only when the integral over the subinterval satisfies
  22 the error test over the subinterval. The total estimated error of the integral
  23 is contained in the global variable QUANC8_ERREST (default value=0.0).
  24
  25 The global variable QUANC8_FLAG (default value=0.0) will contain valuable
  26 information if the computation fails to satisfy the error conditions. The
  27 integer part will tell you how many subintervals failed to converge and the
  28 fractional part will tell you where the singular behavior is, as follows:
  29 singular point=lo+(1.-frac part)*(hi-lo).
  30 Thus quanc8(tan(x),x,1.57,1.6); gives frac=.97 so trouble is at 1.57+.03*.03=
  31 1.5709 (=half pi). If QUANC8_FLAG is not 0.0, you should be cautious in using
  32 the return value, and should try ROMBERG or a Simpson method and see if the
  33 result checks. Analysis of possible singular behavior might be advisable. You
  34 may get QUANC8_FLAG=<integer>.0 and an error message (such as division by 0)
  35 when a singular point is hit in the interval. You will have to find the
  36 singularity and eliminate it before QUANC8 will get an answer. Functions
  37 which have very large derivatives may throw the error estimate way off and
  38 cause the wrong points to be used, and a wrong answer returned. Try
  39 romberg(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with the default tolerance, and
  40 quanc8(exp(-.002*x^2)*cos(x)^2,x,0.,100.); with quanc8_relerr=1.e-7 and 1.e-8.
  41 The last result is consistent with romberg while the previous one is off by
  42 a factor of 2! This is due to the bad behavior of the derivatives near
  43 x=10.0 which cause the adaptive routine to have trouble. If you use
  44 quanc8('f,a,c)+quanc8('f,c,b) where a<c<b, you will do better in such cases.
  45
  46 Typing a control-right-bracket (^]) will give a printout of where the
  47 computation is at the moment, and how much time has been used.
  48
  49 You can do batch("qq demo"); for some comparisons with the ROMBERG numerical
  50 integrator (which is not adaptive). Note that ROMBERG usually gives more
  51 accurate answers for comparable tolerances, while QUANC8 will get the same
  52 answer faster even with a smaller tolerance, because ROMBERG subdivides the
  53 whole interval if the total result is not within error tolerance, while QUANC8
  54 improves only where needed, thus saving many function calls. ROMBERG will also
  55 fail to converge when oscillatory behavior is overwhelming, while QUANC8 will
  56 adapt in the regions as it sees fit. (The global variable ROMBERGMIN is
  57 designed to allow you a minimum number of function calls in such cases, so
  58 that exp(-x)*sin(12*x) can be integrated from 0 to 4*%pi without erroneously
  59 giving 0. from the first few function calls.)
  60
  61 To make your macsyma user function callable in the fastest way, you must use
  62 MODE_DECLARE and then translate and compile the function. Read TRANSL;TRANSL
  63 NEWS for info on this, or ask me for more info. The speed of the computation
  64 may be increased by well over an order of magnitude when compilation is used.
  65 If you do multiple integrals, it is really necessary to compile the function
  66 in order to avoid the time spent on function calls. A sample use of QUANC8
  67 for a double integral is in the batch("qq demo"); and compilation is nearly
  68 a hundred times faster in doing the work!
  69
  70 Comments, bugs, etc. should be sent to LPH@MIT-MC.