How do I extract the weights used by NIntegrate?

I'd like to use NIntegrate to adaptively integrate an expression, then use the same points and weights to integrate that expression with some parameters changed. I know I can extract the points used by NIntegrate using EvaluationMonitor. Is there an analogous way to get the weights?

• Have you seen this? – J. M.'s technical difficulties Mar 3 '16 at 15:22
• I hadn't seen that, no. IntegrationMonitor doesn't seem to contain the relevant information, but it might be possible to get it out of the tracing code provided in that answer. Thanks! – Matt Mar 3 '16 at 19:00
• @Matt I thought that the application of IntegrationMonitor to this question is not that straightforward so that is why I provided an answer. – Anton Antonov Mar 3 '16 at 19:08
• – george2079 Mar 3 '16 at 19:46

Procedure outline

The concrete steps depend on which strategy and integration rules are used. Let us assume that:

1. "GlobalAdaptive" is used with one of the deterministic rules (e.g. "GaussKronrodRule" or "ClenshawCurtisRule"),
2. no singularity handling is applied (using the option "SingularityHandler"->None).

With these assumptions we can do the re-integration in the following way.

1. Find the last collection of regions used by "GlobalAdaptive". IntegrationMonitor can be used to get the iteration step and the integration regions at that step. See the explantions in the second update of my answer to "Determining which rule NIntegrate selects automatically".

2. For each region object robj change the integrand with the new integrand and apply the integration rule of robj.

3. Sum the results of the previous step.

There are some problems applying step 2 directly, so instead we can use the function IRuleEstimate defined in NIntegrate[]'s advanced documentation (see NIntegrate Integration Strategies ).

Code

Preliminary definitions:

Clear[IRuleEstimate]
IRuleEstimate[f_, {a_, b_}, absc_, weights_, errweights_] :=
Module[{integral,
error}, {integral,
error} = (b - a) Total@
MapThread[{f[#1] #2, f[#1] #3} &, {Rescale[
absc, {0, 1}, {a, b}], weights, errweights}];
{integral, Abs[error]}];

nPoints = 5;
{absc, weights, errweights} =
NIntegrateGaussKronrodRuleData[nPoints, MachinePrecision];


First integration:

F[x_] := 1/Sqrt[Abs[x]]

Block[{k = 0}, {val, {iregions}} =
Reap@NIntegrate[F[x], {x, -1, 10}, PrecisionGoal -> 6,
Method -> {"GaussKronrodRule", "Points" -> nPoints},
"SingularityHandler" -> None},
IntegrationMonitor :> (Sow[{k++, #}] &)]
];


Re-integration:

newF[x_] := 1/Sqrt[1 + x + x^2]

Total@
Map[IRuleEstimate[newF, #["Boundaries"[]][[1]], absc, weights,
errweights] &, iregions[[-1, 2]]]

(* {3.73937, 7.95498*10^-9} *)


Check:

Integrate[newF[x], {x, -1, 10}]
% // N

(* ArcSinh[1/Sqrt[3]] + ArcSinh[7 Sqrt[3]] *)

(* 3.73937 *)

• Thanks! I'm assuming the generalization to multiple dimensions is fairly straightforward, depending on how IntegrationMonitor records boundaries? Also, if I wanted to use Singularity Handling (say if I knew that the singularities would be similar in the functions I'm looking at), how would the answer change? – Matt Mar 3 '16 at 19:38
• @Matt For multi-dimensional integrals the function IRuleEstimate` (i) has to be given abscissas and weights corresponding to the dimension, and (ii) the rescaling has to be done for all axes. As for the singularity handling, I would say it is not that straightforward because of the related variable transformations. (Too lengthy to explain further here.) I would consider posting an update for multi-dimensional integrals... – Anton Antonov Mar 3 '16 at 20:15
• @Matt, actually, I would suggest asking the multidimensional case as a separate question instead of as an update to this one; it's sufficiently different and complicated to warrant a separate thread, I reckon. – J. M.'s technical difficulties Mar 3 '16 at 22:21
• Anton: if memory serves, the Iri-Moriguti-Takasawa handler for singularities is already sufficiently complicated... :) – J. M.'s technical difficulties Mar 3 '16 at 22:23
• @J.M. Agreed -- "IMT" is complicated and "DuffyCoordinates" is not that simple either. – Anton Antonov Mar 3 '16 at 22:31