Different Methods in NIntegrate

Question

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, 
  Method -> "LocalAdaptive"]

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.

Answer 1

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
*)

Answer 2

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]*)

addition

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]*)

Answer 3

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.