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) (d1  y^2/(4T)) Exp(y^2/(8T))
( EllipticTheta(3, 1/2 Pi (xx0), Exp(Pi^2 T) ) + EllipticTheta(3, 1/2 Pi (x+x0), Exp(Pi^2 T) ) )
( EllipticTheta(3, 1/2 Pi (xpx0), 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,
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, (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

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
theEllipticTheta
function approaches to a sum of Dirac delta functions. Analytically, this does not seem to be a problem since theExp(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. 
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}

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