divonne(ndim, ncomp, integrand, ..., lower=rep(0,ndim), upper=rep(1,ndim), rel.tol= 0.001, abs.tol = 0, flags=list(verbose=1, final=1, pseudo.random=0, mersenne.seed=NULL), min.eval=0, max.eval=50000, key1=47, key2=1, key3=1, max.pass=5, border=0, max.chisq=10, min.deviation=0.25, xgiven=NULL, nextra=0, peakfinder=NULL)cuhrecuhrecuhre.
But, here, the input argument phw indicates the integration
phase:
0: sampling of the points in xgiven,
1: partitioning phase,
2: final integration phase,
3: refinement phase.
This information might be useful if the integrand takes long to compute and a sufficiently
accurate approximation of the integrand is available. The actual value of the integrand is only
of minor importance in the partitioning phase, which is instead much more dependent on
the peak structure of the integrand to find an appropriate tessellation. An approximation
which reproduces the peak structure while leaving out the fine details might hence be a
perfectly viable and much faster substitute when phw < 2.In all other instances, phw can be ignored.
cuhrecuhrecuhrecuhrecuhrecuhre
pseudo.random and mersenne.seed are only taken into
account when the argument key1 is negative.
cuhrecuhrekey1 = 7, 9, 11, 13 selects the cubature rule of degree key1. Note that the degree-11
rule is available only in 3 dimensions, the degree-13 rule only in 2 dimensions.
For other values of key1, a quasi-random sample of
$\code{n=|key1|}$ points is used, where
the sign of key1 determines the type of sample,key1 = 0, use the default rule.
key1 > 0, use a Korobov quasi-random sample,
key1 < 0, use a “standard” sample (a Mersenne Twister pseudo-random sample
if flags$pseudo.random=1, otherwise a Sobol quasi-random sample).
key1, but here
$\code{n = |key2|}$
determines the number of
points, $\code{n > 39}$, sample $n$ points,
$\code{n < 40}$, sample $\code{n}$
nneed points, where nneed is the number of points needed to
reach the prescribed accuracy, as estimated by Divonne from the results of the
partitioning phase.key3 = 0, do not treat the subregion any further.
key3 = 1, split the subregion up once more.
Otherwise, the subregion is sampled a third time with key3 specifying the sampling
parameters exactly as key2 above.max.pass successive
iterations.A decrease in points generally indicates that Divonne discovered new structures of
the integrand and was able to find a more effective partitioning.
max.pass can be
understood as the number of “safety” iterations that are performed before the partition is accepted as final and counting consequently restarts at zero whenever new
structures are found.
border if the integrand
subroutine cannot produce values directly on the integration
boundary. The relative width of the border
is identical in all the dimensions.
For example, set border=0.1 for a border of width equal
to 10% of the width of the integration region.min.deviation
move on to the refinement phase.ndim, ngiven).
A list of ngiven points where the
integrand might have peaks.Divonne will consider these points when partitioning the integration region. The idea here is to help the integrator find the extrema of the integrand in the presence of very narrow peaks. Even if only the approximate location of such peaks is known, this can considerably speed up convergence.
nextra is zero, peakfinder is
not called and an arbitrary object may be passed in its place, e.g. just
0.xgiven. It is expected to be declared aspeakfinder <- function(bounds)
where bounds is a
matrix of dimension (ndim, 2) which contains
the upper and lower bounds of the subregion. The names of the columns
are c("lower", "upper").
The returned value should be a matrix (ndim, nx)
where nx is the actual number of
points (should be less or equal to
nextra).
cuhre.
Here ifail may be >1 when
the accuracy goal was not met within the allowed maximum number of
integrand evaluations. Divonne
can estimate the number of points by which
maxeval needs to be increased to
reach the desired accuracy and returns this value.
See details in the documentation.
J. H. Friedman, M. H. Wright (1981) User's guide for DIVONNE. SLAC Report CGTM-193-REV, CGTM-193, Stanford University.
T. Hahn (2005) CUBA-a library for multidimensional numerical integration. Computer Physics Communications, 168, 78-95.
cuhre, suave, vegasNDIM <- 3
NCOMP <- 1
integrand <- function(arg, phase) {
x <- arg[1]
y <- arg[2]
z <- arg[3]
ff <- sin(x)*cos(y)*exp(z);
return(ff)
}
divonne(NDIM, NCOMP, integrand, rel.tol=1e-3, abs.tol=1e-12,
flags=list(verbose=2), key1= 47)
# Example with a peak-finder function
NMAX <- 4
peakf <- function(bounds) {
# print(bounds) # matrix (ndim,2)
x <- matrix(0, ncol=NMAX, nrow=NDIM)
pas <- 1/(NMAX-1)
# 1ier point
x[,1] <- rep(0, NDIM)
# Les autres points
for (i in 2:NMAX) {
x[,i] <- x[,(i-1)] + pas
}
return(x)
} #end peakf
divonne(NDIM, NCOMP, integrand,
flags=list(verbose=0) ,
peakfinder=peakf, nextra=NMAX)
Run the code above in your browser using DataLab