Multidimensional Numerical Integration Using the Cuba Library
evalf(Int(..., method=meth, methodoptions=molist))
int(..., numeric, method=meth)
int(..., numeric, method=meth, methodoptions=molist)
one of the Cuba methods, _CubaVegas, _CubaSuave, _CubaDivonne, or _CubaCuhre
a list of equations of the form key = value
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 method=_CubaVegas, 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 method=_CubaSuave, 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 method=_CubaDivonne, 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, r, have similar values of a quantity called the spread, s⁡r, defined by:
where V⁡r is the volume of r. 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 χ2 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 method=_CubaCuhre, is a deterministic integration method. It is most suitable for moderate dimensions. A different implementation of the same algorithm is available with method=_cuhre.
The Cuba library allows the user to specify a large variety of options. These can be passed to the routines using the methodoptions=molist argument. Some options are common to all methods, others only occur for some.
absepsilon = positive number (default: 10−Digits)
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 4000000⁢d+4000000 and 50000000, where d 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.
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.
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.
Same as for the Vegas algorithm.
key1 = integer (default: 47)
This determines sampling in the first (partitioning) phase. Having key1∈7,9,11,13 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 key1 points is used. If 0<key1, then it is a Korobov quasi-random sample, and if key1<0, it is a Sobol quasi-random sample.
key2 = integer (default: 1)
This determines sampling in the second (regular integration) phase. Like for key1, having key2∈7,9,11,13 selects the corresponding deterministic cubature rule. Otherwise, the sign of key2 determines the type of quasi-random sample. The number of samples is equal to key2, if that is greater than or equal to 40, or key2⁢nneed otherwise, where nneed 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 key3=0, then subregions are not treated any further in this phase. If key3=1, then suitable subregions are split up once more. Otherwise, key3 determines the sampling parameters in the same way as key2.
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 safetyiterations 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 χ2 value for a subregion in order for it to not processed in the third (refinement) phase. If a region has a χ2 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.
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.
If you set infolevelevalf or infolevelevalf/int, you can obtain some diagnostic information.
We turn on extra information from the numeric integration code.
infolevel`evalf/int` ≔ 4:
Simple Examples of Each Method
integrand1 ≔ 11+x−0.2+y−0.7+z+0.3
region1 ≔ x=−1..1,y=−1..1,z=−1..1
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: abserr=.305046087432878584e-4; usedpts=7806563
We apply the four Cuba algorithms, with standard options.
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
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
cuba: # of integrand evaluations: 103867
cuba: estimated (absolute) error: 3.04855e-05
cuba: number of regions that the domain was divided into: 30
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.
integrand2 ≔ ⅇadd⁡−i⁢xi−i2010,i=5..15
region2 ≔ seq⁡xi=0..1,i=5..15
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: 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: # 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.
cuba: could not attain requested accuracy
cuba: # of integrand evaluations: 1007500
cuba: estimated (absolute) error: 1.31319e-06
evalf/int: error from Control_multi was:
"could not attain requested accuracy; try increasing epsilon or absepsilon or maximalpoints"
The integrand is not smooth, so we could try to turn off smoothing. However, in this example it does not make a difference.
The Suave Method
The parameter one will most typically need to adapt for the Suave method is nnew. 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 nnew is low. Hence, for low requested accuracy and highly variable integrands, lower values of nnew 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 nnew would be quite slow with the default number of maximalpoints. However, if we lower the requested accuracy, it does work.
cuba: # of integrand evaluations: 371000
cuba: estimated (absolute) error: 0.000304326
cuba: number of regions that the domain was divided into: 371
Still, with nnew set to 3000 (rather than its default of 1000), we get the result with substantially fewer integrand evaluations.
cuba: # of integrand evaluations: 48000
cuba: estimated (absolute) error: 0.000297572
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 flatness to a low number. If we would need higher accuracy, we could determine the most effective value for flatness at low accuracy, then do the final computation with that value of flatness. (The calling sequence int(..., numeric, ...), used here, is the same as the one with evalf(Int(...)).)
spikes ≔ Statistics:-Sample⁡Uniform⁡0,1,6,4
integrand3 ≔ add⁡mul⁡ln⁡1.7⁢xi−spikesi,j,i=1..6,j=1..4
region3 ≔ seq⁡xi=0..1,i=1..6
Control_multi: integrating on [0, 0, 0, 0, 0, 0] .. [1, 1, 1, 1, 1, 1] the 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: # of integrand evaluations: 4440000
cuba: estimated (absolute) error: 0.00152791
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.
cuba: # of integrand evaluations: 2066376
cuba: estimated (absolute) error: 0.00134965
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.
cuba: # of integrand evaluations: 1662309
cuba: estimated (absolute) error: 0.00137549
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.
integrand4 ≔ cos⁡x+sin⁡y36⁢sin⁡y−cos⁡x32
region4 ≔ x=0..10⁢Pi,y=0..10⁢Pi
Control_multi: integrating on [0, 0] .. [10*Pi, 10*Pi] the 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: # of integrand evaluations: 99003
cuba: estimated (absolute) error: 0.120726
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.
cuba: # of integrand evaluations: 90415
cuba: estimated (absolute) error: 0.0198997
cuba: number of regions that the domain was divided into: 696
The Cuhre implementation outside the Cuba library is similar.
cuhre: epsrel=.1e-3; minpts=0; maxpts=195000
cuhre: abserr=.197922435357250452e-1; usedpts=97695
The Multidimensional Numerical Integration Using the Cuba Library command was introduced in Maple 18.
For more information on Maple 18 changes, see Updates in Maple 18.
Download Help Document