Context
I am trying to calculate some transport coefficients for a heat equation in confinement. The Boundaries are in the $x$ direction, and $y$ represents the parallel directions. This function essentially boils down to the following
$Q ( g , x ) = \int_0^1 dx' \int_0^\infty \, dy' \int_0^1 dx'' \int_0^\infty \, dy'' \, \int_0^1 dx_0 \int_0^\infty dT \, f(x' , x'') \, \times\frac{ y'( x - x')}{ ( g \,(x-x')^2 + {y'}^2 )^{3/2} }\frac{ y''(x' - x'')}{ ( g \,(x'-x'')^2 + {y''}^2 )^{3/2} },$
where we have
$f(x',x'') = \frac{\partial^2}{\partial x' \partial x''} \frac{ e^{-\frac{{y''}^2}{8T}}}{T^2} \left( -1 + \frac{{y''}^2}{8 T} \right)[ \theta_3 ( \frac{\pi ( x' + x_0 )}{2},e^{-\pi^2 T}) +\theta_3 ( \frac{\pi ( x' - x_0 )}{2},e^{-\pi^2 T}) ] \times[ \theta_3 ( \frac{\pi ( x'' + x_0 )}{2},e^{-\pi^2 T}) +\theta_3 ( \frac{\pi ( x'' - x_0 )}{2},e^{-\pi^2 T})]$
and $\theta_3$ represents the Jacobi Theta function which solves the heat equation in confinement.
I want to plot the behavior of $Q(g,x=0)$ and $Q(g,x=1)$ for $ 0 < g < 2$
Mathematica code
As a continuation of a previous question, I am now trying to numerically calculate the following integral:
hardintegral [ g_?NumericQ , x_?NumericQ ] :=
NIntegrate[
(Exp[-ypp^2/(8T)] / T^2) * ( -1 + ypp^2/(8T) ) *
( EllipticThetaPrime[3, 1/2 Pi (xp + x0), Exp[-Pi^2 T] ] +
EllipticThetaPrime[3, 1/2 Pi (xp - x0), Exp[-Pi^2 T] ] ) *
( EllipticThetaPrime[3, 1/2 Pi (xpp + x0), Exp[-Pi^2 T] ] +
EllipticThetaPrime[3, 1/2 Pi (xpp - x0), Exp[-Pi^2 T] ] ) *
( yp*(x-xp) / ( g*(x-xp)^2 + yp^2 )^(3/2) ) * ( ypp*(xp-xpp) / ( g*(xp-xpp)^2 + ypp^2 )^(3/2) ),
{x0, 0, 1} , {T, 0, ∞}, {xp, 0, 1} , {xpp, 0, 1} , {yp, 0, ∞}, {ypp, 0, ∞} ]
I want to get the following plots: Plot[ hardintegral [g,0] , {g,0,2} ] and
Plot[ hardintegral [g,1] , {g,0,2} ]. However, even obtaining a single result, for say g=1.1 is taking a very long time on my computer. Using Method->"GlobalAdaptive" I get 2.83493*10^6 with the following error
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times.
The global error is expected to decrease monotonically after a number of integrand evaluations.
Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0.
Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration.
NIntegrate obtained 2.8349279022111776`*^6 and 7.683067946598636`*^7 for the integral and error estimates.
Also, With Method->"GaussKronrodRule the computation goes on forever with no outcome.
Is there a way to speed up these integrations? I guess a possible solution for the plot then will be to use ListPlot.
PS
The yp and ypp integrations can be done using Integrate. For example
Integrate[ Exp[-z^2/8T] * ( z / (a + z^2)^(3/2) ) , {z, 0, ∞}, Assumptions-> a>0 && T>0 ]
gives
( Gamma[1/2 (-1 + d)] HypergeometricU[ 1/2 (-1 + d), 1/2, a/(8 T) ] ) / (2 Sqrt[a])
Also for
Integrate[ Exp[-z^2/8T] * ( z^3 / (a + z^2)^(3/2) ) , {z, 0, ∞}, Assumptions-> a>0 && T>0 ]
the result is
1/2 Sqrt[a] * ( Gamma[1/2 (1 + d)] HypergeometricU[ 1/2 (1 + d), 3/2, a/(8 T) ] )
I tried plugging these back into the NIntegrate but it doesn't seem to do much in terms of the speed.
ParallelTable[{g,hardintegral [g,0]} , {g,1,3,0.1}]. The advantage is you control the step size and you can do it in parallel. Also have you tried the MonteCarlo methods for NIntegrate? They may be faster. $\endgroup$xand diverges atT->0. $\endgroup$T->0. So we can regularized it by cutting integral overypplike{ypp, eps, Infinity}witheps>0or overTwitheps=10^-2. With the last assumption I have a solution. $\endgroup$