Here is another type of integration rule. It has a customized error estimator,
incorporating which was my main interest. In fact, the rule below computes an integral of arbitrary dimension, using its own internal formula for the integral as well.
A rule of the type NIintegrate`GeneralRule[{abscissas, weights, errweights}] computes
computes the integral and error estimates via dot products in more or less the following
following way:
There is a level above this. Recall that in adaptive integration NIntegrate subdivides
subdivides regions (intervals in 1D), and applies the integration rule to each subregion
subregion. If the error is too great in a subregion, it divides again. When
NIntegrate applies a GeneralRule to a subregion of integration, the result returned
returned to NIntegrate has the form
A GeneralRule is always 1D, so the axis is always 1. In multidimensional rules
rules the number will indicate along which axis the region is to be split in half
half, if the error is too great. Ultimately, the triple above is what an integration
integration rule needs to return when it is called by NIntegrate.
A rule, let's call it chebRule as we will below, can be called by NIntegrate in
in at least two ways, one of which I have figured out. The call has the form
where chebRule[data___] is the structure returned by
NIntegrate`InitializeIntegrationRule and region is a
NIntegrate`IntegrationRegion[] that contains data (including the integrand) for
for computing the integral over the region. An ``NIntegrateIntegrationRegion[] is
is the argument #1 in the option
IntegrationMonitor has
has many methods
methods. We
We will need two:
which give the boundaries of the region and a function for evaluating the integrand
integrand, respectively. You evaluate the function on the rule's sampling points
points taken from the interval {0, 1}. The method takes care of the rescaling for
for you.
All All you need to do inside
chebRule[data___]["ApproximateIntegral"[region]] is calculate the three things above
above, the integral, error, and, if implementing a multidimensional rule, which axis
axis to subdivide next.
The initialization code processes the options, of which currently there is only one
one, and gets the Clenshaw-Curtis sampling points for use in the rule.
The The only option
option is "Points", which works like the
"ClenshawCurtisRule" option
option and generates 2 n +- 1 points for a setting of n. Unlike
"ClenshawCurtisRule", it will take a list of point settings for multidimensional
multidimensional integrals.
ClearAll[chebRule];
Options[chebRule] = {"Points" -> Automatic};
chebRuleProperties = Part[Options[chebRule], All, 1];
chebRule::dim = "The value of option `` should be a list of length ``";
(* initialization of rule *)
chebRule /: NIntegrate`InitializeIntegrationRule[
chebRule, nfs_, ranges_, ruleOpts_, allOpts_] :=
Module[{absc, points, minpts = 3, dim, wprec, t},
wprec = WorkingPrecision /. allOpts;
t = NIntegrate`GetMethodOptionValues[chebRule, chebRuleProperties, ruleOpts];
If[t === $Failed, Return[$Failed]];
{points} = t;
points = points /. Automatic -> 5;
dim = Length@ranges;
If[! AllTrue[Flatten@{points}, # >= minpts &],
Message[NIntegrate::mintmin, First@FilterRules[ruleOpts, "Points"], minpts];
Return[$Failed]
];
(* "Points" -> list of length dim *)
If[ListQ@points && Length@points != dim,
Message[chebRule::dim, First@FilterRules[ruleOpts, "Points"], dim];
Return[$Failed]
];
(* absc = list of ascissas for each dimension *)
If[ListQ@points,
absc = First@NIntegrate`ClenshawCurtisRuleData[#, wprec] & /@ points,
absc = First@NIntegrate`ClenshawCurtisRuleData[points, wprec];
absc = Table[absc, {dim}]
];
chebRule[absc]
]
Here is the integration code. First there is the error estimator. It comes from
from the convergence theory for Chebyshev series (see Boyd, reference below). What
What really matters is that it computes a number that estimates the error of the quadrature
quadrature to pass back to NIntegrate. The code below is used to come up with an
an error for each axis (from the iterated integral), which is then used in the main
main integrator below to determine how the subregion should be divided next.
(* call: self["localError"[coeff]] --
* estimates error of Chebyshev series `coeff`,
* provided truncated series is long enough *)
chebRule[_]["localError"[coeff_]]chebLE := Compile[{{coeff, _Real, 1}},
With[{e2 = Total[Abs@coeff[[{-2, ;;]]]-1}]]], e4 = Total[Abs@coeff[[{-4, ;;]]]-3, -2, -1}]]]},
With[{r2 = e2/e4},
If[Length@coeff > 8 && e4 != 0 && r2 < 1,
1.21 e2^2e2*r2/e4(1 - Sqrt[r2]), (* extrapolate error *)
e2] (* use last two terms as upper bound *)
]; ]]
];
(* callcalled many times, MachinePrecision OK(?), hence compiled *)
chebRule[_]["localError"[coeff_]] := chebLE[coeff];
(* self["globalError"[error[[n]],Times@@widths[[;;n]],n]] --
* estimates "global" error of integrating along axis `n`
* by integrating local error bound along axis `(n-1)`
* and multiplying by the other dimensions (`volume`) *)
chebRule[absc_]["globalError"[errors_chebGE = Compile[{{errors, volume__Real, dim_]]1}, :={prevabsc,
_Real, 1}, Module[{error =volume, errors_Real}},
If[dim >Module[{error, 1ints},
(* integrate errors along next axis up (trapezoidal rule) *)
error = Partition[Flatten[error]Partition[Flatten[errors], Length@absc[[dim - 1]]];Length@prevabsc];
errorints = volume *
error.(PadRight[#Join[{First@#}, Length@absc[[dimRest@# -+ 1]]Most@#, 0.{Last@#}] + &@
Differences@prevabsc);
PadLeft[#, Length@absc[[dim0.5 -volume*Max@ints 1]],(*
0.] &@ Differences@absc[[dim -max 1]]integral will be the error bound *);
error]
= Max@error/2, ];
chebRule[absc_]["globalError"[errors_, volume_, 1]] := volume*Max[errors];
(* topcompiling axis:won't volumehelp ==much, lengthexcept ofmaybe intervalin high dimensional integrals *)
chebRule[absc_]["globalError"[errors_, volume_, erroraxis_]] := volume*Max[error]
chebGE[errors, ]
absc[[axis - 1]], ];volume];
The computation of the Chebyshev series follows Boyd, which I first encountered
here. Finding
Finding the integral from a Chebyshev series is
well-known and
and easy to derive.
chebRule[absc_]["Abscissas"[]] := absc;
(* Subroutines for computing and integrating Chebyshev series *)
chebRule[_]["series"[vals_, width_]] := Module[{coeff},
coeff =
Sqrt[2/(Length@vals - 1)]*FourierDCT[Reverse@vals, 1]/width;
coeff[[{1, -1}]] /= 2;
coeff
];
chebRule[_]["integrate"[coeff_, width_]] :=
Map[ (* integration rule for a Chebyshev series *)
Total[width/(1 - Range[0, Length@# - 1, 2]^2)*#[[1 ;; ;; 2]]] &,
coeff,
{-2}];
(*
* Main integrator
*)
(self : chebRule[___])["ApproximateIntegral"[region_]] :=
Module[{absc, ranges, dim, fvals, coeff, intintegral, widths, error, naxis},
absc = self@ "Abscissas"[];self@"Abscissas"[];
ranges = region@ "Boundaries";region@"Boundaries";
dim = Length@ ranges;Length@ranges;
widths = Flatten[Differences /@ ranges];
error = Table[{}, {dim}];
fvals = (* evaluate integrand *)
Outer[region["EvaluateTransformedIntegrand"[{##}]] &,
Sequence @@ absc];
naxis = dim; (* convenient iterator with Nest[] *)
intintegral = Nest[
(* replacereplaces each {-2}bottom level list of function/integral values by its integral *)
Function[{vals},
(* first*first compute Chebyshev series and approximation error *)
coeff = Map[
(* compute cosine transformsMap[self["series"[#, onwidths[[axis]]]] bottom&, levelvals, *){-2}];
(coeff =* Sqrt[2/(Length@#estimate -the 1)]errors * FourierDCT[Reverse@#, 1]/widths[[n]]; )
coeff[[{1, -1}]]error[[axis]] /= 2;
(* first of error[[n]]the =Chebyshev {error[[n]],approximation self["localError"[coeff]]};*)
coeff)Map[self["localError"[#]] &,
valscoeff, {-2}];
(*error[[axis]] estimate= an error(* next for the whole level (axis n) *)
error[[n]] = self["globalError"[error[[n]]self["globalError"[error[[axis]], Times @@ widths[[;; n]]axis]], n]];
axis]];
(* caveat: decrement n axis one step too early *)
naxis--;
(* integration rule for a Chebyshev series *)
Map[
Total[widths[[n + 1]]/(1 - Range[0, Length@# - 1Map[self["integrate"[#, 2]^2) #[[1 ;;widths[[axis ;;+ 2]]]1]]]] &,
coeff,
{-2}]
],
fvals,
dim];
(* define dbPrintaxis = Print for verbose outputFirst@Ordering[error, *)-1];
dbPrint[ranges -> {intintegral, error, First@Ordering[error, -1]axis}];
error = (* convert from MachinePrecision to working precision *)
{int, Max@error, First@Ordering[errorSetPrecision[Max@error, -1]}region@"WorkingPrecision"];
(* {integral, error, axis *)}
];
##Examples (TBD)
Here you can see the recursive subdivision and the convergence error estimates.
If you follow carefully, you can see how the error and the axis determine the subdivisions
subdivisions. (There are three iterations; they can be grouped as the first printed
printed rule, the next two, and the last four.)
Block[{dbPrint = Print},
Reap@NIntegrate[Cos[7 x] Sin[5 y], {x, 0, 2}, {y, 0, 3},
Method -> {chebRule, "Points" -> 9},
IntegrationMonitor :> ((Sow@Total@Through[#1@"Error"]) &)]
]
First@% - Integrate[Cos[7 x] Sin[5 y], {x, 0, 2}, {y, 0, 3}]
MoreA relatively innocuous integrand is
f1 = Sin[2 x] Sin[y^3] Sin[x y];
We will integrate it with calls of the form
NIntegrate[f1, {x, 0, 2}, {y, 0, 3}, Method -> ...]
with Method settings like
Method -> chebRule
Method -> {chebRule, "Points" -> 9}
Method -> {chebRule, "Points" -> 17}
Method -> {chebRule, "Points" -> {9, 17}}
We compare it with similar calls to comethe rules "GaussKronrod"``, ``"ClenshawCurtisRule", except for the last. By comparison the exact value to 16 digits is
exact = 0.1998439750348762
and the default "MultidimensionalRule" takes 2.38428 seconds, 8325 function evaluations and has an error of 3.75823*10^-10.
Comparison of methods (code at end)
First of all, the error is closer to the goal, which be around 10^-8 for this integral. The function f1 evaluates quickly so one can see that chebRule appears to be slower per function evaluation than "ClenshawCurtisRule", but it seems also to use fewer function calls. The larger error and fewer calls are related to the difference in the error estimators of the two methods. If the integrand is expensive to calculate, then chebRule may have an advantage.
The error weights in a GeneralRule are typically the difference of the weights of
of two related integration rules which is known to approximate (or at least bound
bound) the error.
In In my case, I wanted to compute the error in a different way, basically
basically for exploring my own interests and not to create a better integrator. From
From Anton's examples and a little spelunking, I seem to have figured how to do this
this. Since it is an answer to this question and it extends the (currently only
only) other answer, I thought I should share it. I doubt my code is bullet
bullet-proof. It can perform quite well; it can also do just okay; and I haven't
haven't tested enough pathological cases to find one where it breaks, but I haven't
haven't been too cruel to it, either.
The rule below is closely related to the Clenshaw-Curtis rule, which itself follows
follows the form of a GeneralRule. Instead of the usual error estimate, I wanted
wanted to estimate the error from the Chebyshev series of integrand over the integration
integration region. If the function is analytic in a (complex) neighborhood of the
the interval of integration, then the coefficients converge to zero geometrically
geometrically and one can estimate the error of approximation from the tail of such
such a truncated series.
(See, for instance, Boyd, Chebyshev and Fourier Spectral
Spectral Methods, Dover, New York, 2001
2001, ch. 2). Since one can construct
construct a series in which any number of coefficients are zero, this is not a fool
fool-proof way to estimate the error. Except for (nearly) even or odd functions
functions, which have series in which every other coefficient is (nearly) zero, functions
functions in practice tend to have only a few (finite number of) terms that do not
not behave.
##Extra code
Code for the comparison of methods:
f1 = Sin[2 x] Sin[y^3] Sin[x y];
exact = 0.19984397503487617461507951200812381171`16.;
Needs["Integration`NIntegrateUtilities`"];
ClearAll[ni];
ni[{"ClenshawCurtisRule" | "GaussKronrod", "Points" -> {__}}] := "";
ni[meth_] := {"IntegralEstimate", "IntegralEstimate" - exact,
"Evaluations", "Timing"} /.
NIntegrateProfile@NIntegrate[f1, {x, 0, 2}, {y, 0, 3},
Method -> meth] /.
Plus[InputForm[a_?NumericQ], b_?NumericQ] :> a + b // Column;
$rules = {Automatic; "GaussKronrod", "ClenshawCurtisRule", chebRule};
tbl = Table[ni[{rule, pts}],
{pts, {"Points" -> Automatic, "Points" -> 9, "Points" -> 17,
"Points" -> {9, 17}}},
{rule, $rules}];
Join[{{"Points\\Method", SpanFromLeft, Sequence @@ $rules}},
Transpose@
Join[{{"Auto.", 9, 17, {9, 17}},
Table[Column[{"Est", "Error", "Evals", "Time"}], {4}]},
Transpose[tbl]]
] // Grid[#, Dividers -> All, BaseStyle -> Smaller] &
