Maxima Manual
22.2 DCADRE
The following is obsolete. To make an interface to fortran
libraries in the current MAXIMA look at the examples in
"maxima/src/fortdef.lsp"
- The IMSL version of Romberg integration is now available in
Macsyma. For documentation, Do PRINTFILE(DCADRE,USAGE,IMSL1); . For
a demo, do batch("dcadre.mc");
This is a numerical integration package using cautious, adaptive
Romberg extrapolation.
The DCADRE package is written to call the IMSL fortran library routine
DCADRE. This is documentation for that program. Send bugs/comments to
KMP
To load this package, do
For a demo of this package, do
The worker function takes the following syntax:
IMSL_ROMBERG(fn,low,hi)
where fn is a function of 1 argument; low and hi should be the lower and
upper bounds of integration. fn must return floating point values.
IMSL_ROMBERG(exp,var,low,hi)
where exp should be integrated over the range var=low to hi. The result
of evaluating exp must always be a floating point number.
FAST_IMSL_ROMBERG(fn,low,hi)
This function does no error checking but may achieve a speed gain over
the IMSL_ROMBERG function. It expects that fn is a Lisp function (or
translated Macsyma function) which accepts a floating point argument
and that it always returns a floating point value.
Returns either
[SUCCESS, answer, error] where answer is the result of the integration and
error is the estimated bound on the absolute error of the output, DCADRE,
as described in PURPOSE below.
or
[WARNING, n, answer, error] where n is a warning code, answer is the answer,
and error is the estimated bound on the absolute error of the output, DCADRE,
as described in PURPOSE below. The following warnings may occur:
65 = One or more singularities were successfully handled.
66 = In some subinterval(s), the estimate of the integral was accepted
merely because the estimated error was small, even though no regular
behavior was recognized.
or
[ERROR, errorcode] where error code is the IMSL-generated
error code. The following error codes may occur:
131 = Failure due to insufficient internal working storage.
132 = Failure. This may be due to too much noise in function
(relative to the given error requirements) or due to an
ill-behaved integrand.
133 = RERR is greater than 0.1 or less than 0.0 or is too small
for the precision of the machine.
The following flags have an influence upon the operation of IMSL_ROMBERG --
ROMBERG_AERR [Default 1.0E-5] -- Desired absolute error in answer.
ROMBERG_RERR [Default 0.0] -- Desired relative error in the answer.
Note: If IMSL signals an error, a message will be printed on the user's
console stating the nature of the error. (This error message
may be supressed by setting IMSLVERBOSE to FALSE.)
Note: Because this uses a translated Fortran routine, it may not be
recursively invoked. It does not call itself, but the user should
be aware that he may not type ^A in the middle of an IMSL_ROMBERG
computation, begin another calculation using the same package,
and expect to win -- IMSL_ROMBERG will complain if it was already
doing one project when you invoke it. This should cause minimal
problems.
Purpose (modified version of the IMSL documentation)
----------------------------------------------------
DCADRE attempts to solve the following problem: Given a real-valued
function F of one argument, two real numbers A and B, find a number
DCADRE such that:
| | | / B | [ | / B | ]
| [ | [ | [ | ]
| I F(x)dx - DCADRE | <= max [ ROMBERG_AERR, ROMBERG_RERR * | I F(x)dx | ]
| ] | [ | ] | ]
| / A | [ | / A | ]
|
Algorithm (modified version of the IMSL documentation)
This routine uses a scheme whereby DCADRE is computed as the sum of
estimates for the integral of F(x) over suitably chosen subintervals of
the given interval of integration. Starting with the interval of
integration itself as the first such subinterval, cautious Romberg
extrapolation is used to find an acceptable estimate on a given
subinterval. If this attempt fails, the subinterval is divided into two
subintervals of equal length, each of which is considered separately.
Programming Notes (modified version of the IMSL documentation)
-
1. DCADRE (the translated-Fortran base for IMSL_ROMBERG) can, in many cases,
handle jump discontinuities and certain algebraic discontinuities. See
reference for full details.
-
2. The relative error parameter ROMBERG_RERR must be in the interval [0.0,0.1].
For example, ROMBERG_RERR=0.1 indicates that the estimate of the intergral
is to be correct to one digit, where as ROMBERG_RERR=1.0E-4 calls for four
digits of accuracy. If DCADRE determines that the relative accuracy
requirement cannot be satisfied, IER is set to 133 (ROMBERG_RERR should be
large enough that, when added to 100.0, the result is a number greater than
100.0 (this will not be true of very tiny floating point numbers due to
the nature of machine arithmetic)).
-
3. The absolute error parameter, ROMBERG_AERR, should be nonnegative. In
order to give a reasonable value for ROMBERG_AERR, the user must know
the approximate magnitude of the integral being computed. In many cases,
it is satisfactory to use AERR=0.0. In this case, only the relative error
requirement is satisfied in the compuatation.
-
4. We quote from the reference, "A very cautious man would accept DCADRE
only if IER [the warning or error code] is 0 or 65. The merely reasonable
man would keep the faith even if IER is 66. The adventurous man is quite
often right in accepting DCADRE even if the IER is 131 or 132." Even when
IER is not 0, DCADRE returns the best estimate that has been computed.
For references on this technique, see
de Boor, Calr, "CADRE: An Algorithm for Numerical Quadrature,"
Mathematical Software (John R. Rice, Ed.), New York, Academic Press,
1971, Chapter 7.
