# How to implement custom NIntegrate integration strategies?

How can new integration strategies algorithms be used with NIntegrate?

This is a different type of extension than the extensions with new integration rules, as described in the answer for the question "How to implement custom integration rules for use by NIntegrate?". (Integration strategies perform and guide the core integration process, using or leveraging different integration rules and/or preprocessing algorithms.)

## Motivation (for a new semi-symbolic integration strategy)

Consider the following integral, which cannot be done neigther by Integrate:

Integrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1}]

(* Integrate[If[Re[y] > -(1/3), Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]),
Integrate[BesselJ[y, x^3], {x, 0, Infinity},
Assumptions -> Re[y] <= -(1/3)]], {y, 0, 1}] *)


nor NIntegrate:

NIntegrate[BesselJ[y, x^3], {x, 0, ∞}, {y, 0, 1},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 2000}]


NIntegrate::slwcon: Numerical integration converging too slowly;

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has ...

(* 0.524338 *)


Here is a plot of the integrand function over a much smaller domain:

Plot3D[BesselJ[y, x^3], {x, 0, 10}, {y, 0, 1}, PlotPoints -> {100, 10},
MaxRecursion -> 5, PlotRange -> All, BoxRatios -> {10, 3}] Because of the oscillatory nature of the integrand we can see why NIntegrate has difficulties.

On the other hand, Integrate can find the value of the integral integrating over $x$:

In:= Integrate[BesselJ[y, x^3], {x, 0, ∞}]

Out= ConditionalExpression[Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]),
Re[y] > -(1/3)]


but it has problems integrating over $y$:

In:= Integrate[BesselJ[y, x^3], {y, 0, 1}]

Out= Integrate[BesselJ[y, x^3], {y, 0, 1}]


Since Integrate can do partially the integral along one of the axes, we can just take that symbolic expression and give it to NIntegrate for integration over the other axis.

## Semi-symbolic NIntegrate implementation

Here we make an integration strategy that combines Integrate and NIntegrate -- it uses Integrate over some of the integration range(s) and then NIntegrate for the rest of the range(s) with the symbolic expressions obtained by Integrate.

The following defintion is for the initialization of the integration strategy SemiSymbolic.

Clear[SemiSymbolic];
Options[SemiSymbolic] = {"AnalyticalVariables" -> {}};
SemiSymbolicProperties = Options[SemiSymbolic][[All, 1]];
SemiSymbolic /:
NIntegrateInitializeIntegrationStrategy[SemiSymbolic, nfs_, ranges_,
strOpts_, allOpts_] :=
Module[{t, anVars},
t = NIntegrateGetMethodOptionValues[SemiSymbolic, SemiSymbolicProperties,
strOpts];
If[t === $Failed, Return[$Failed]];
{anVars} = t;
SemiSymbolic[First /@ ranges, anVars]
];


This is the implementation of the integaration strategy SemiSymbolic:

SemiSymbolic[vars_, anVars_]["Algorithm"[regions_, opts___]] :=
Module[{ranges, anRanges, funcs, t},

ranges = Map[Flatten /@ Transpose[{vars, #@"OriginalBoundaries"}] &, regions];
ranges = Map[Flatten, ranges, {-2}];
anRanges = Map[Select[#, MemberQ[anVars, #[]] &] &, ranges];
ranges = Map[Select[#, ! MemberQ[anVars, #[]] &] &, ranges];
funcs = (#@"Integrand"[])@"FunctionExpression"[] & /@ regions;

Integrate[#1, Sequence @@ #2,
Assumptions -> (#[] <= #[] <= #[] & /@ #3)] &, {funcs,
anRanges, ranges}];
Print["SemiSymbolic::Integrate's result:", t];

If[! FreeQ[t, Integrate], Return[\$Failed]];

NIntegrate[#1, Sequence @@ #2 // Evaluate,
Sequence @@ DeleteCases[opts, Method -> _] // Evaluate] &, {t, ranges}]]
];


(Note the implementation prints the intermediate result obtained by Integrate.)

## Signatures

### Initialization

We can see that the new rule SemiSymbolic is defined through TagSetDelayed for SemiSymbolic and NIntegrateInitializeIntegrationStrategy. The rest of the arguments are:

nfs -- numerical function objects; several might be given depending on the integrand and ranges;

ranges -- a list of ranges for the integration variables;

strOpts -- the options given to the strategy;

allOpts -- all options given to NIntegrate.

### Algorithm

StrategySymbol[strategyData___]["Algorithm"[regions_, opts___]] := ...


The algorithm can use regions objects as described in this answer of "Determining which rule NIntegrate selects automatically".

## Testing SemiSymbolic

The strategy works without (observable) problems for the motivational integral:

In:= NIntegrate[BesselJ[y, x^3], {x, 0, Infinity}, {y, 0, 1},
Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}]

During evaluation of In:= SemiSymbolic::Integrate's result:{(2^-y HypergeometricPFQ[{1/6+y/2},{7/6+y/2,1+y},-(1/4)])/((1+3 y) Gamma[1+y])}

Out= 0.524448


Note the printout for the intermediate result by Integrate.

Since SemiSymbolic passes inside its body the non-method NIntegrate options it was invoked with we can also see the sampling points used by SemiSymbolic through EvaluationMonitor.

res =
Reap@NIntegrate[BesselJ[y, x^3], {x, 0, Infinity}, {y, 0, 1},
Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}},
EvaluationMonitor :> Sow[{x, y}]]

During evaluation of In:= SemiSymbolic::Integrate's result:{(2^-y HypergeometricPFQ[{1/6+y/2},{7/6+y/2,1+y},-(1/4)])/((1+3 y) Gamma[1+y])}

(* {0.524448, {{{x, 0.00795732}, {x, 0.0469101}, {x, 0.122917}, {x,
0.230765}, {x, 0.360185}, {x, 0.5}, {x, 0.639815}, {x,
0.769235}, {x, 0.877083}, {x, 0.95309}, {x, 0.992043}, {x,
0.00397866}, {x, 0.023455}, {x, 0.0614583}, {x, 0.115383}, {x,
0.180092}, {x, 0.25}, {x, 0.319908}, {x, 0.384617}, {x,
0.438542}, {x, 0.476545}, {x, 0.496021}, {x, 0.503979}, {x,
0.523455}, {x, 0.561458}, {x, 0.615383}, {x, 0.680092}, {x,
0.75}, {x, 0.819908}, {x, 0.884617}, {x, 0.938542}, {x,
0.976545}, {x, 0.996021}}}} *)

ListPlot[res[[2, 1, All, 2]], Frame -> True] ## Further tests

Below are some other tests / examples.

In:= NIntegrate[x^2 + y^2 + z^2, {x, 0, 1}, {y, 0, 1}, {z, 0, 1},
Method -> {SemiSymbolic, "AnalyticalVariables" -> {x, y}}]

During evaluation of In:= SemiSymbolic::Integrate's result:{2/3+z^2}

Out= 1.


Note that the symbolic integration was done over two variables.

Let us use the same integrand but with different range boundaries for the different variables in order to evaluate better the variable correspondence in the 2D sampling points pattern.

In:= res =
Reap@NIntegrate[x^2 + y^2 + z^2, {x, 0, 1}, {y, 0, 2}, {z, 0, 10},
Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}},
EvaluationMonitor :> Sow[{y, z}]]

During evaluation of In:= SemiSymbolic::Integrate's result:{1/3+y^2+z^2}

Out= {700., {{{1., 5.}, {1.35857, 5.}, {0.641431, 5.}, {1.94868,
5.}, {0.0513167, 5.}, {1., 6.79284}, {1., 3.20716}, {1.,
9.74342}, {1., 0.256584}, {1.94868, 9.74342}, {1.94868,
0.256584}, {0.0513167, 0.256584}, {0.0513167, 9.74342}, {0.311753,
1.55876}, {0.311753, 8.44124}, {1.68825, 1.55876}, {1.68825,
8.44124}}}}

In:= ListPlot[res[[2, 1]], Frame -> True] ## Another example

This Lebesgue integration implementation, AdaptiveNumericalLebesgueIntegration.m -- discussed in detail in "Adaptive numerical Lebesgue integration by set measure estimates" -- has implementations of integration strategy (and rules) with the complete signatures for the plug-in mechanism.

• +1. However, I've found a damaging bug of Boundaries in IntegrationMonitor which gives wrong ranges {0,1} for {x,0,Infinity},(my guess is that some transformation work is make internally before returning the correct boundaries) so your SemiSymbolic actually integrate over region {x,0,1},{y,0,1}. (BTW, the correct result is 0.524448.) Do you know how to fix this? Since you've helped to designed the numerical integration, did you know when WR will add these instructions of user defined rules to the docs? – luyuwuli Nov 28 '16 at 14:40
• Not a all. I'm very grateful for your reply. I've already post a new question related to this issue here mathematica.stackexchange.com/q/132355/6468, could you explain it there with some details? Michael E2 is also interested in this issue. – luyuwuli Dec 1 '16 at 1:30
• There was a bug in SemiSymbolic` -- I fixed it. I deleted my comment-answer to @luyuwuli since I posted a longer explanation here. – Anton Antonov Dec 1 '16 at 16:51