# Different Methods in NIntegrate

I am trying to numerically calculate a multidimensional integral which involves Jacobi elliptic theta functions. The integrand is the following:

integrand[d_, x_, y_, xp_, x0_, T_] :=
T^(-(d + 1)/2) (d-1 - y^2/(4T)) Exp[-y^2/(8T)] *
( EllipticTheta[3, 1/2 Pi (-x-x0), Exp[-Pi^2 T] ] + EllipticTheta[3, 1/2 Pi (-x+x0), Exp[-Pi^2 T] ] ) *
( EllipticTheta[3, 1/2 Pi (-xp-x0), Exp[-Pi^2 T] ] + EllipticTheta[3, 1/2 Pi (-xp+x0), Exp[-Pi^2 T] ])


My goal is to integrate this expression with respect to x0 and T for d=3, and get a 3d plot of the result as a function of x and xp (both variables between 0 and 1) - while manipulating y. After this, I need to take the derivative of the integrated result with respect to both x and xp.

For the integration, I have tried 3 different strategies. In the first one, I do not specify the Method:

integral[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
PrecisionGoal -> 10, MinRecursion -> 10, MaxRecursion -> 20]


I've found that increasing the MinRecursion changes the results, and 10 seems to work well (higher values do not seem to improve the results). Since generating the full 3D plot takes somewhat long, I then generated the following table:

Table[integral[3, x, 1, 0], {x, 0.05, 1, 0.05}]


with the outcome

{-43.386, -38.7746, -34.1253, -31.4359, -26.9778, -22.7969, -19.8602, -20.2972, -13.8984, -6.49645, -3.3476, -3.31147, 6.20662, 8.2472, 12.0905, 13.7228, 14.896, 15.814, 16.3162, 16.463}


In a second attempt, I tried Method->"LocalAdaptive" for the integration:

adaptintegral[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
PrecisionGoal -> 10, MinRecursion -> 10, MaxRecursion -> 20,


which produces the following numbers for the same table:

{-20.7877, -19.7131, -17.9935, -15.7272, -13.0363, -10.0544, -6.91493, -3.74124, -0.63984, 2.30356, 5.02495, 7.48073, 9.64493, 11.5056, 13.061, 14.316, 15.2788, 15.9584, 16.3626, 16.4967}


The outcome is very different compared to the first table, and since I did not get any error messages, I wonder if there is a way to tell which gives a more accurate estimate of the actual result.

I also tried the Method->"MonteCarlo":

mcintegral[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
PrecisionGoal -> 10, MinRecursion -> 10, MaxRecursion -> 20,
Method -> "MonteCarlo"]


which gives the following values for the same table

{-21.2913, -19.2249, -18.663, -16.2671, -13.3218, -9.81518, -4.44489, -3.11635, -0.264413, 2.72884, 4.44556, 8.09827, 9.49501, 11.4452, 13.0165, 14.0828, 15.279, 16.3008, 16.6255, 16.5606}


This one works much faster, but I also get a few error messages like this one

NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained -21.2913 and 1.3762138095540868 for the integral and error estimates.


## Questions

1. Is there a good way to compare these methods and make sure which results are reliable? I suspect this is due to a singularity in (part of) the integrated - since as T->0 the EllipticTheta function approaches to a sum of Dirac delta functions. Analytically, this does not seem to be a problem since the Exp[-y^2/(8T)] factor causes the integrand to become identical to zero. However, I imagine things are not as straight in numerics, but I also do not know how to overcome this hurdle.

2. What can I do to speed up these computations? Especially, for generating and Manipulating the Plot3D of integral (or different variants of it) with {x,0,1},{xp,0,1}

3. How to (numerically) take derivatives from integral w.r.t both x and xp? I both need to plot this derivative, as well as to integrate it against another kernel.

## Edit

The Table I mentioned above is just an example to show these methods give different results, two of them without generating any errors. Eventually, I need to get things like

Manipulate[ Plot3D[integral[3,x,y,xp],{x,0,1},{y,-5,5}],{xp,0,1}]


or

Manipulate[ Plot3D[integral[3,x,y,xp],{x,0,1},{xp,0,1}],{y,-5,5}]


This is also the case for taking derivatives, i.e. I need things like

Manipulate[ Plot3D[ D[integral[3,x,y,xp],x,xp], {x,0,1}, {y,-5,5} ], {xp,0,1}]


## Update

I realized that for the derivatives, I can use the built-in function EllipticThetaPrime. For the plots, I eventually had to use a ListPlot3D for which entries are calculated using the answer by @Michael E2.

• I fixed your code, which was leading to us getting different results. Jul 1, 2020 at 12:52

### New answer (b/c the old one was based on a mistake in the posted code):

First, multidimensional integrals can be hard to compute. Both easy and hard ones are common in dimension 2. The proportion of hard ones seems to increase with the dimension. Integrating over infinite domains can be hard if the integrand is oscillatory, which is not the case here. Integrands with singularities can be hard, too, which also is not the case here. Each of these problems is sufficiently common to have methods to address them.

The Monte Carlo methods are modestly useful when all else fails. They give a rough approximation somewhat quickly. They converge very slowly and using them to pursue high precision is usually futile.

This seems a moderately difficult integral. The integrand does not seem pathological, but the default rule, a medium-order "MultidimensionalRule", seems to struggle. In fact, it seems to get the wrong answer with the global-adaptive strategy. It turns out that the local-adaptive strategy in the OP is accurate. How to verify that?

Generally, a cartesian-product rule based on the Gauss-Kronrod or Clenshaw-Curtis rule will be effective on a smooth integrand. The main draw back is that they tend to be slow in high dimensional integrals because of excessive sampling. We can use them to verify the local-adaptive result.

In fact, though, my usual first step with a smooth integrand is to raise the order of the multidimensional rule with the suboption "Generators" - > 9. This turns out to be a good method here, too.

There is no need to use MinRecursion or other options. I'll use both a medium- and high-order Gauss-Kronrod rules to check consistency. (Another way to check consistency is to double the working precision to WorkingPrecision -> 32, but I'll omit that.)

(* high-order multidimensional rule *)
i1[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
Method -> {"MultidimensionalRule", "Generators" -> 9}];

(* Gauss-Kronrod cartesian product rule *)
i2[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
Method -> "GaussKronrodRule"];

(* High-order Gauss-Kronrod cartesian product rule: a double check *)
i3[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
integrand[d, x, y, xp, x0, T], {T, 0, ∞}, {x0, 0, 1},
Method -> {"GaussKronrodRule", "Points" -> 11}];


The OP's table with these methods agrees with each:

Table[i1[3, x, 1, 0], {x, 0.05, 1, 0.05}] // AbsoluteTiming
(*
{4.46711, {-20.7877, -19.7131, -17.9935, -15.7272, -13.0363,
-10.0544, -6.91493, -3.74124, -0.63984, 2.30356, 5.02495, 7.48073,
9.64493, 11.5056, 13.061, 14.316, 15.2788, 15.9584, 16.3626,
16.4967}}
*)

Table[i2[3, x, 1, 0], {x, 0.05, 1, 0.05}] // AbsoluteTiming
(*
{4.37294, {-20.7877, < same as above >, 16.4967}}
*)

Table[i3[3, x, 1, 0], {x, 0.05, 1, 0.05}] // AbsoluteTiming
(*
{7.19945, {-20.7877, < same as above>, 16.4967}}
*)


The derivative with respect to y

One way is to differentiate under the integral sign:

i2dy[d_?NumericQ, x_?NumericQ, y_?NumericQ, xp_?NumericQ] :=
NIntegrate[
D[integrand[d, x, \[FormalY], xp, x0, T], \[FormalY]] /. \[FormalY] -> y,
{T, 0, ∞}, {x0, 0, 1},
Method -> "GaussKronrodRule"];


Another is to use complex-step differentiation. A third way is to use the central difference formula. Below is an example of each:

i2dy[3, 0.1, 1, 0]
i2[3, 0.1, 1 + Sqrt@$$MachineEpsilon*I, 0]/Sqrt@$$MachineEpsilon // Im
(i2[3, 0.1, 1 + 0.5 Sqrt@$$MachineEpsilon, 0] - i2[3, 0.1, 1 - 0.5 Sqrt@$$MachineEpsilon, 0])/Sqrt@$MachineEpsilon (* 77.8076 77.8076 77.8076 *)  • The integration now makes a lot more sense, thank you. Could you also please try running the 3d plots? they are very slow to generate. Jul 1, 2020 at 16:29 • The y differentiation you could actually get an analytical expression for. I was more stuck with the x and xp differentiation, since it doesn't seem that Mathematica handles the differentiation of EllipticTheta very well. Jul 1, 2020 at 16:30 • @dedekindCuttage I have no problems with changing the y derivative to either x or xp. For the second derivative, I'd have to think about whether that's possible with the complex-step, but putting differentiating under the integral sign with respect to both x and xp and the central difference formula (with step size Sqrt@Sqrt@$MachineEpsilon) work for me. Jul 1, 2020 at 17:20

The integral is zero for Element[{x,y},Reals] (Thanks to the answer Michael E2)

Integrate[integrand[3, x, y, 0, x0, T], {T, 0, \[Infinity]}, {x0, 0, 1}]
(*ConditionalExpression[0, Re[y^2] > 0]*)


The integral depending on x,y,xp is zero for Element[{y},Reals]

 Integrate[integrand[3, x, y, xp, x0, T], {T, 0, \[Infinity]}, {x0, 0, 1}]

(*ConditionalExpression[0, Re[y^2] > 0]*)

• curious...could you specify what version of Mathematica you are using? Mine is 11.3 and it takes very long to run this, without an outcome at the end. Also, would it be possible to tweak this to be able to use Manipulate for xp? Jul 1, 2020 at 9:05
• @dedekindCuttage My Mathematica version is v12 (Windows 10). Have a look at my edited answer. Jul 1, 2020 at 9:06
• OK. I will update mine and try again. Jul 1, 2020 at 9:12
• And thanks for the addition. It is actually unexpected to me that this is always identical to 0. So what does Integrate exactly do here, and is it more reliable than NIntegrate? Is there a way to confirm this result numerically as well? Jul 1, 2020 at 9:14
• @dedekindCuttage Integrate is surely more reliable, the symbolical capabilities of Mathematica are the reason why I use it. You might verify the results pointwise for x,y,xp but Nintegrate gives warning messages. Jul 1, 2020 at 9:29

I've found similar issues when doing high-dimension integrals. A reliable method is QuasiMonteCarlo, since the set of sampling points it uses are more evenly distributed than in MonteCarlo, and therefore it will converge faster. However, if your integral receives the most contribution from a single point, e.g. a spike/singularity, then an adaptive method would work better, since it will preferentially sample the singularity (as long as your initial grid refinement is fine enough to see it in the first place), and therefore will converge faster.

In your case, identify any singularities and then do some integrals focused around them to see whether they will make a large contribution to the integral. If they don't make a large contribution, then QuasiMonteCarlo should be fine. If they do contribute a lot, then I recommend breaking your integral down into several domains, so that you can integrate the singularities separately from the rest of the domain.

• Thank you for the suggestion. I have actually tried breaking down the interval to separate the regions that include the singularities (which I think is around T=0 and x0=0,1`) - without much success. The main contribution is, I think, coming from around these points and I don't know how to choose the best method to numerically integrate them. Jul 1, 2020 at 8:08
• Try using AdaptiveQuasiMonteCarlo near singularities; it will automatically refine the grid in the vicinity of the singularity, so it should converge quickly. To check your answer, see if you can Taylor expand your integrand at the singularity (if you know its exact location).
– Guy
Jul 2, 2020 at 5:06