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.
##Introduction
A rule of the type ``NIintegrate`GeneralRule[{abscissas, weights, errweights}]``
computes the integral and error estimates via dot products in more or less the
following way:
fvals = f /@ abscissas;
integral = fvals . weights;
error = fvals . errweights;
Many 1D integration rules are of this form (hence the name `GeneralRule`).
There is a level above this. Recall that in adaptive integration `NIntegrate`
subdivides regions (intervals in 1D), and applies the integration rule to each
subregion. If the error is too great in a subregion, `NIntegrate` divides it again. When
`NIntegrate` applies a `GeneralRule` to a subregion of integration, the result
returned to `NIntegrate` has the form
{integral value, error estimate, axis to be subdivided}
A `GeneralRule` is always 1D, so the axis is always `1`. In multidimensional
rules the number will indicate along which axis the region is to be split in
half, if the error is too great. Ultimately, the triple above is what an
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 at least two ways, one of which I have figured out. The call has the form
chebRule[data___]["ApproximateIntegral"[region]]
where `chebRule[data___]` is the structure returned by
``NIntegrate`InitializeIntegrationRule `` and `region` is a
``NIntegrate`IntegrationRegion[]`` that contains data (including the integrand)
for computing the integral over the region. (An ``NIntegrate`IntegrationRegion[]``
is also the argument `#1` in the option `IntegrationMonitor`, an undocumented option which has been discussed
[in these Q&A](https://mathematica.stackexchange.com/search?q=IntegrationMonitor).)
An `IntegrationRegion` has [many
methods](https://mathematica.stackexchange.com/questions/26401/determining-which-rule-nintegrate-selects-automatically/96663#96663).
We will need two:
region@"Boundaries"
region["EvaluateTransformedIntegrand"[{##}]] &
which give the boundaries of the region and a function for evaluating the
integrand, respectively. You evaluate the function on the rule's sampling
points taken from the interval `{0, 1}`. The method takes care of the rescaling
for you. All you need to do inside
`chebRule[data___]["ApproximateIntegral"[region]]` is calculate the three things
above, the integral, error, and, if implementing a multidimensional rule, which
axis to subdivide next.
##Example: Computing iterated integrals via Chebyshev series
The initialization code processes the options, of which currently there is only
one, and gets the Clenshaw-Curtis sampling points for use in the rule. The only
option is `"Points"`, which works like the
[`"ClenshawCurtisRule"`](http://reference.wolfram.com/language/tutorial/NIntegrateIntegrationRules.html#486402291)
option and generates `2 n - 1` points for a setting of `n`. Unlike
`"ClenshawCurtisRule"`, it will take a list of point settings for
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 the convergence theory for Chebyshev series (see Boyd, reference below).
What really matters is that it computes a number that estimates the error of the
quadrature to pass back to `NIntegrate`. The code below is used to come up with
an error for each axis (from the iterated integral), which is then used in the
main integrator below to determine how the subregion should be divided next.
(* self["localError"[coeff]] --
* estimates error of Chebyshev series `coeff`,
* provided truncated series is long enough *)
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*r2/(1 - Sqrt[r2]), (* extrapolate error *)
e2] (* use last two terms as upper bound *)
]]
];
(* called 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`) *)
chebGE = Compile[{{errors, _Real, 1}, {prevabsc, _Real, 1}, {volume, _Real}},
Module[{error, ints},
(* integrate errors along next axis up (trapezoidal rule) *)
error = Partition[Flatten[errors], Length@prevabsc];
ints =
error.(Join[{First@#}, Rest@# + Most@#, {Last@#}] &@
Differences@prevabsc);
0.5 volume*Max@ints (*
max integral will be the error bound *)
]
];
chebRule[absc_]["globalError"[errors_, volume_, 1]] := volume*Max[errors];
(* compiling won't help much, except maybe in high dimensional integrals *)
chebRule[absc_]["globalError"[errors_, volume_, axis_]] :=
chebGE[errors, absc[[axis - 1]], volume];
The computation of the Chebyshev series follows Boyd, which I first encountered
[here](https://mathematica.stackexchange.com/questions/5663/about-multi-root-search-in-mathematica-for-transcendental-equations/89692#89692).
Finding the integral from a Chebyshev series is
[well-known](https://en.wikipedia.org/wiki/Chebyshev_polynomials#Differentiation_and_integration)
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, integral, widths, error, axis},
absc = self@"Abscissas"[];
ranges = region@"Boundaries";
dim = Length@ranges;
widths = Flatten[Differences /@ ranges];
error = Table[{}, {dim}];
fvals = (* evaluate integrand *)
Outer[region["EvaluateTransformedIntegrand"[{##}]] &,
Sequence @@ absc];
axis = dim; (* convenient iterator with Nest[] *)
integral = Nest[
(* replaces each bottom level list of function/integral values by its integral *)
Function[{vals},
(*first compute Chebyshev series *)
coeff = Map[self["series"[#, widths[[axis]]]] &, vals, {-2}];
(* estimate the errors *)
error[[axis]] = (* first of the Chebyshev approximation *)
Map[self["localError"[#]] &, coeff, {-2}];
error[[axis]] = (* next for the whole level (axis n) *)
self["globalError"[error[[axis]], Times @@ widths[[;; axis]],
axis]];
(* caveat: decrement axis one step too early *)
axis--;
Map[self["integrate"[#, widths[[axis + 1]]]] &, coeff, {-2}]
],
fvals,
dim];
axis = First@Ordering[error, -1];
dbPrint[ranges -> {integral, error, axis}];
error = (* convert from MachinePrecision to working precision *)
SetPrecision[Max@error, region@"WorkingPrecision"];
{integral, error, axis}
];
##Examples
**1.**
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. (There are four iterations; they can be grouped as the first
printed rule and then groups of two; the last two pairs show the subdivision of the first two pairs, in reversed order according to their error estimates, `1.93 x 10^-6` vs. `3.27 * 10^-6`. The errors that `IntegrationMonitor` sows are the sums at each step of the errors over the current subintervals.)
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}]
**2.**
A 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 the rules `"GaussKronrod"`, `"ClenshawCurtisRule"`, except for the `{9, 17}` setting, which they do not support. For comparison, the exact value to 16 digits is
exact = 0.1998439750348762
and the default method `"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 for `chebRule` is closer to the goal, which should 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 over `"ClenshawCurtisRule"`.
**3.**
An example where `chebRule` can beat `"ClenshawCurtisRule"` is shown below. One thing to realize is that one does not take advantage of the convergence rate of the Chebyshev series if a small value for `"Points"` is used (this is true for the Clenshaw-Curtis rule, too). Another is that if we do not use `MachinePrecision`, then we cannot take advantage of current hardware that optimizes matrix multiplication, and `"ClenshawCurtisRule"` loses some of its advantage.
exact = Integrate[Cos[x]^2 Cos[10 x]^2, {x, 0, 50}];
NIntegrate[Cos[x]^2 Cos[10 x]^2, {x, 0, 50},
Method -> {chebRule, "Points" -> 1 + 2^6},
WorkingPrecision -> 32] - exact // AbsoluteTiming
NIntegrate[Cos[x]^2 Cos[10 x]^2, {x, 0, 50},
Method -> {"ClenshawCurtisRule", "Points" -> 1 + 2^6},
WorkingPrecision -> 32] - exact // AbsoluteTiming
(*
{0.137261, -2.*10^-30}
{0.317895, 0.*10^-31}
*)
One can improve `chebRule` quite a bit by writing one's own "`chebPoints"` routine that calculates just the abscissas (and not the weights and error weights).
----
##Background (TBD)
I decided to put the back story in the back (at the end of the post).
The error weights in a `GeneralRule` are typically the difference of the weights
of two related integration rules which is known to approximate (or at least
bound) the error. In my case, I wanted to compute the error in a different way,
basically for exploring my own interests and not to create a better integrator.
From Anton's examples and a little spelunking, I seem to have figured how to do
this. Since it is an answer to this question and it extends the (currently
only) other answer, I thought I should share it. I doubt my code is
bullet-proof. It *can* perform quite well; it can also do just okay; and I
haven't tested enough pathological cases to find one where it breaks, but that is not my main interest. Hopefully, it will be accepted as an example of how such a integration rule can be incorporated in `NIntegrate[]`.
The rule `chebRule` is closely related to the Clenshaw-Curtis rule, which itself
follows the form of a `GeneralRule`. Instead of the usual error estimate, I
wanted to estimate the error from the Chebyshev series of integrand over the
integration region. If the function is analytic in a (complex) neighborhood of
the interval of integration, then the coefficients converge to zero
geometrically and one can estimate the error of approximation from the tail of
such a truncated series. (See, for instance, [Boyd, *Chebyshev and Fourier
Spectral Methods*, Dover, New York,
2001](https://books.google.com/books?id=b4TCAgAAQBAJ), ch. 2). Since one can
construct a series in which any number of coefficients are zero, this is not a
fool-proof way to estimate the error. Except for (nearly) even or odd
functions, which have series in which every other coefficient is (nearly) zero,
functions in practice tend to have only a few (finite number of) terms that do
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] &