Multidimensional Numerical Integration Using the Cuba Library
|
Calling Sequence
|
|
evalf(Int(..., method=meth))
evalf(Int(..., method=meth, methodoptions=molist))
int(..., numeric, method=meth)
int(..., numeric, method=meth, methodoptions=molist)
|
|
Parameters
|
|
meth
|
-
|
one of the Cuba methods, _CubaVegas, _CubaSuave, _CubaDivonne, or _CubaCuhre
|
molist
|
-
|
a list of equations of the form key = value
|
|
|
|
|
Description
|
|
•
|
The Cuba library provides four routines for multidimensional numerical integration. It is distributed via the webpage http://www.feynarts.de/cuba. This help page describes the usage of this library in Maple in some detail.
|
•
|
For general information about numerical integration with Maple, see the help page evalf/Int.
|
•
|
Because of the nature of these methods, you should expect relatively slow convergence; they are most useful for difficult cases, when other methods converge slowly too, or do not converge at all. In order to get sufficient convergence, you should expect to increase the error tolerance, as specified by the epsilon option to int or Int or the absepsilon method option. Examples are below.
|
•
|
The source of "random" numbers is fixed in this library: it is a Sobol quasi-random sequence.
|
|
Description of the four integration methods
|
|
•
|
The Vegas method, selected with , uses importance sampling as a variance reduction technique. It is particularly suitable if the integrand is separable; that is, if it can be written as a product of univariate functions in the integration variables. This can be viewed as the most basic of the four Cuba methods.
|
•
|
The Suave method, selected with , uses importance sampling combined with a globally adaptive subdivision strategy. This method is unique to Cuba. It is quite effective, but uses more memory than most other methods.
|
•
|
The Divonne method, selected with , uses stratified sampling for variance reduction. In particular, it attempts to find the minimum and maximum in each region using optimization methods, and attempts to partition the integration region so that all subregions, , have similar values of a quantity called the spread, , defined by:
|
|
where is the volume of . The algorithm proceeds in (up to) three phases. In phase 1, the domain of integration is partitioned and estimates for the required number of samples are made. In phase 2, these samples are taken and a test assesses whether the estimated accuracy conforms with the results of phase 1. Subregions where it does not undergo a further phase 3 for refinement.
|
•
|
The Cuhre method, selected with , is a deterministic integration method. It is most suitable for moderate dimensions. A different implementation of the same algorithm is available with .
|
|
|
Method Options
|
|
•
|
The Cuba library allows the user to specify a large variety of options. These can be passed to the routines using the argument. Some options are common to all methods, others only occur for some.
|
|
Common Options
|
|
•
|
absepsilon = positive number (default: )
|
|
This specifies an absolute tolerance. Integration stops when either the absolute tolerance, or the relative tolerance (specified by the epsilon argument in the regular evalf/int calling sequence) has been achieved.
|
•
|
minimalpoints = nonnegative integer (default: 0)
|
|
This specifies the (approximate) minimum number of integrand evaluations that should take place.
|
•
|
maximalpoints = nonnegative integer
|
|
This specifies the (approximate) maximal number of integrand evaluations that should take place. If not specified, then this setting is taken from the environment variable _EnvIntMaxPoints, if that has a positive integer value. Otherwise, it is taken as the maximum of and , where is the dimension of the region of integration.
|
•
|
allsamples = true or false (default: true)
|
|
This option determines whether all the samples taken during the integration process are used in the final evaluation, or only those in the last run.
|
|
|
Vegas-Specific Options
|
|
•
|
nstart = positive integer (default: 1000)
|
|
Number of samples to start the first round of integration with.
|
•
|
nincrease = positive integer (default: 500)
|
|
Increase in the number of samples for every following round.
|
•
|
nbatch = positive integer (default: 1000)
|
|
The Vegas algorithm samples points in batches of this size.
|
•
|
smooth = true or false (default: true)
|
|
Whether to apply smoothing to the importance function. This moderately improves convergence for relatively smooth integrands, but should be disabled if the integrand has sharp edges.
|
|
|
Suave-Specific Options
|
|
•
|
nnew = positive integer (default: 1000)
|
|
The number of new integrand evaluations for every subdivision.
|
•
|
flatness = positive number (default: 50)
|
|
This parameter determines how prominently "outliers" figure in the total fluctuation. It should be chosen large for "flat" integrands and small for "volatile" integrands. Supplying values greater than a few hundred can cause truncation of numbers in the internal code, and it is therefore not recommended.
|
•
|
smooth = true or false (default: true)
|
|
Same as for the Vegas algorithm.
|
|
|
Divonne-Specific Options
|
|
•
|
key1 = integer (default: 47)
|
|
This determines sampling in the first (partitioning) phase. Having selects the deterministic cubature rule of that degree; the rule of degree 11 is available only in dimension 3, and the rule of degree 13 in dimension 2. Otherwise, a quasi-random sample of points is used. If , then it is a Korobov quasi-random sample, and if , it is a Sobol quasi-random sample.
|
•
|
key2 = integer (default: 1)
|
|
This determines sampling in the second (regular integration) phase. Like for , having selects the corresponding deterministic cubature rule. Otherwise, the sign of determines the type of quasi-random sample. The number of samples is equal to , if that is greater than or equal to 40, or otherwise, where is the number of points needed for the prescribed accuracy as estimated by the algorithm from the results of the first phase.
|
•
|
key3 = integer (default: 1)
|
|
This determines sampling in the third (refinement) phase. If , then subregions are not treated any further in this phase. If , then suitable subregions are split up once more. Otherwise, determines the sampling parameters in the same way as .
|
•
|
safetyiterations = nonnegative integer (default: 5)
|
|
This parameter determines the thoroughness of the partitioning phase. The phase terminates when the estimated total number of integrand evaluations (partitioning plus final integration) does not decrease during consecutive iterations. Such a decrease indicates that Divonne's partitioning now accommodates structure in the integrand that it did not accommodate before, and that the partitioning is therefore more effective.
|
•
|
border = nonnegative number (default: 0)
|
|
The Divonne algorithm will search for extrema near the border of the region of integration. In some cases, the integrand is not well-defined. For such cases, the Divonne algorithm can extrapolate such values from two points further from the border. This parameter gives the relative width of the border region where this feature should be used.
|
•
|
maxchisquare = positive number (default: 10)
|
|
The maximum value for a subregion in order for it to not processed in the third (refinement) phase. If a region has a value greater than this parameter and fails the test associated to the mindeviation parameter, below, then it is processed in the third phase.
|
•
|
mindeviation = nonnegative number (default: 0.25)
|
|
The minimal difference between the two sampling averages (obtained in phases 1 and 2), expressed as a fraction of the requested error of the entire integral, that a subregion must exhibit for it to proceed to the third (refinement) phase.
|
•
|
peaks = list of list of numbers or list of numbers or Vector of numbers or Matrix of numbers (default: empty list)
|
|
If there are known peaks present in the integrand, then the Divonne algorithm can be sped up greatly by specifying these peaks. They are specified by coordinates in the domain of integration. If it is a single peak, the coordinates can be supplied as a list or a Vector. For any number of peaks, the coordinates can be supplied as a list of lists, where each inner list describes the coordinates of one peak. Alternatively, the coordinates can be supplied as a Matrix with each column containing the coordinates of a peak.
|
|
|
Cuhre-Specific Options
|
|
•
|
key = 7, 9, 11, 13, or default
|
|
Supplying a number selects the cubature rule of that degree. The rule of degree 11 is available only in dimension 3, and the rule of degree 13 in dimension 2. If default is supplied (the default), then the highest-degree rule appropriate for the dimension (9, 11, or 13) is taken.
|
|
|
|
Userinfo Messages
|
|
•
|
If you set or , you can obtain some diagnostic information.
|
|
|
|
Examples
|
|
We turn on extra information from the numeric integration code.
>
|
|
|
Simple Examples of Each Method
|
|
>
|
|
| (1) |
>
|
|
| (2) |
•
|
The default method is the NAG Cuhre algorithm.
|
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
trying DCUHRE (TOMS algorithm 698)
| |
cuhre: epsrel=.1e-4; minpts=0; maxpts=90000000
| |
cuhre: result=3.05062989428920384
| |
cuhre: abserr=.305046087432878584e-4; usedpts=7806563
| |
result=3.05062989428920384
| |
•
|
We apply the four Cuba algorithms, with standard options.
|
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 16065000
| |
cuba: estimated (absolute) error: 3.04443e-05
| |
cuba: chi-square probability that the error is not reliable: 0
| |
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 2050000
| |
cuba: estimated (absolute) error: 3.03368e-05
| |
cuba: chi-square probability that the error is not reliable: 1
| |
cuba: number of regions that the domain was divided into: 41
| |
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 103867
| |
cuba: estimated (absolute) error: 3.04855e-05
| |
cuba: chi-square probability that the error is not reliable: 0
| |
cuba: number of regions that the domain was divided into: 30
| |
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 200787
| |
cuba: estimated (absolute) error: 3.04852e-05
| |
cuba: chi-square probability that the error is not reliable: 1.11022e-16
| |
cuba: number of regions that the domain was divided into: 791
| |
|
|
The Vegas Method
|
|
•
|
We use a separable integrand.
|
>
|
|
| (8) |
>
|
|
| (9) |
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 3751000
| |
cuba: estimated (absolute) error: 5.86145e-07
| |
cuba: chi-square probability that the error is not reliable: 2.22045e-16
| |
•
|
We see that about 4000000 function evaluations were needed. If we restrict the number of function evaluations to substantially less, we will not get a floating point result.
|
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: could not attain requested accuracy
| |
cuba: # of integrand evaluations: 1007500
| |
cuba: estimated (absolute) error: 1.31319e-06
| |
cuba: chi-square probability that the error is not reliable: 0
| |
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"
| |
| (11) |
•
|
The integrand is not smooth, so we could try to turn off smoothing. However, in this example it does not make a difference.
|
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 3751000
| |
cuba: estimated (absolute) error: 5.86145e-07
| |
cuba: chi-square probability that the error is not reliable: 2.22045e-16
| |
|
|
The Suave Method
|
|
•
|
The parameter one will most typically need to adapt for the Suave method is . Lower values mean the process of integration is stopped relatively frequently to assess if further subdivisions are necessary. Also, there will be more subdivisions for a given number of integrand evaluations if is low. Hence, for low requested accuracy and highly variable integrands, lower values of are better. For higher requested accuracy and for a relatively uniform integrand, low values induce more subdivisions than necessary and higher values are more suitable.
|
•
|
In the first example above, using the default value for would be quite slow with the default number of . However, if we lower the requested accuracy, it does work.
|
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 371000
| |
cuba: estimated (absolute) error: 0.000304326
| |
cuba: chi-square probability that the error is not reliable: 1
| |
cuba: number of regions that the domain was divided into: 371
| |
•
|
Still, with set to 3000 (rather than its default of 1000), we get the result with substantially fewer integrand evaluations.
|
>
|
|
Control_multi: integrating on [-1, -1, -1] .. [1, 1, 1] the integrand
| |
cuba: transformed original integrand
| |
cuba: with lower bounds [-1., -1., -1.] and upper bounds [1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 48000
| |
cuba: estimated (absolute) error: 0.000297572
| |
cuba: chi-square probability that the error is not reliable: 1
| |
cuba: number of regions that the domain was divided into: 16
| |
•
|
The following integrand has substantial spikes. Therefore integration works a little better if we set to a low number. If we would need higher accuracy, we could determine the most effective value for at low accuracy, then do the final computation with that value of . (The calling sequence int(..., numeric, ...), used here, is the same as the one with evalf(Int(...)).)
|
>
|
|
>
|
|
| (16) |
>
|
|
| (17) |
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 4440000
| |
cuba: estimated (absolute) error: 0.00152791
| |
cuba: chi-square probability that the error is not reliable: 1
| |
cuba: number of regions that the domain was divided into: 444
| |
|
|
The Divonne Method
|
|
•
|
This method is highly customizable. Consider, for example, the previous integrand.
|
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 2066376
| |
cuba: estimated (absolute) error: 0.00134965
| |
cuba: chi-square probability that the error is not reliable: 1.11022e-16
| |
cuba: number of regions that the domain was divided into: 446
| |
•
|
For somewhat better performance, one can select Sobol, rather than Korobov, quasi-random numbers for the first phase.
|
>
|
|
Control_multi: integrating on [0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0., 0., 0., 0., 0.] and upper bounds [1., 1., 1., 1., 1., 1.], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 1662309
| |
cuba: estimated (absolute) error: 0.00137549
| |
cuba: chi-square probability that the error is not reliable: 0
| |
cuba: number of regions that the domain was divided into: 331
| |
•
|
In low dimensions, it is often advantageous to select a Cuhre integration rule for use in phases 1 and 2.
|
>
|
|
| (21) |
>
|
|
| (22) |
>
|
|
Control_multi: integrating on [0, 0] .. [10*Pi, 10*Pi] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0.] and upper bounds [31.41592653590, 31.41592653590], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 99003
| |
cuba: estimated (absolute) error: 0.120726
| |
cuba: chi-square probability that the error is not reliable: 2.22045e-16
| |
cuba: number of regions that the domain was divided into: 526
| |
|
|
The Cuhre Method
|
|
•
|
The Cuhre method is deterministic, and is therefore most suitable for dimensions lower than 6 to 10. The previous integral is computed with similar efficiency when using only the Cuhre method.
|
>
|
|
Control_multi: integrating on [0, 0] .. [10*Pi, 10*Pi] the integrand
| |
| |
cuba: transformed original integrand
| |
| |
cuba: with lower bounds [0., 0.] and upper bounds [31.41592653590, 31.41592653590], to the following integrand to be integrated over the unit n-cube:
| |
| |
cuba: integration completed successfully
| |
cuba: # of integrand evaluations: 90415
| |
cuba: estimated (absolute) error: 0.0198997
| |
cuba: chi-square probability that the error is not reliable: 0
| |
cuba: number of regions that the domain was divided into: 696
| |
•
|
The Cuhre implementation outside the Cuba library is similar.
|
>
|
|
Control_multi: integrating on [0, 0] .. [10*Pi, 10*Pi] the integrand
| |
| |
trying DCUHRE (TOMS algorithm 698)
| |
cuhre: epsrel=.1e-3; minpts=0; maxpts=195000
| |
cuhre: result=199.188283610236368
| |
cuhre: abserr=.197922435357250452e-1; usedpts=97695
| |
result=199.188283610236368
| |
|
|
|
Compatibility
|
|
•
|
The Multidimensional Numerical Integration Using the Cuba Library command was introduced in Maple 18.
|
|
|