# numerical integration – 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, (Infinity)}, {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, (Infinity)}, {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, (Infinity)}, {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. 