# Difficult Numerical Integral with special functions

## 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.

• Instead of Plot, do a ListLinePlot on a Table e.g 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. – flinty Jul 17 '20 at 22:57
• @flinty thanks, I was thinking of just a normal Table, but probably a parallel one is faster. And yes I tried MC ones, but the results are too unreliable and I don't know if I can make sense of them. – SaMaSo Jul 17 '20 at 23:01
• @SaMaSo It looks like integral does not depend on x and diverges at T->0. – Alex Trounev Jul 21 '20 at 20:55
• @AlexTrounev Ah my bad! the mathematica code is the correct formula, but I had made a mistake in the formula I gave in for context. Will correct it, thanks – SaMaSo Jul 22 '20 at 1:50
• @SaMaSo Unfortunately integral diverges at T->0. So we can regularized it by cutting integral over ypp like {ypp, eps, Infinity} with eps>0 or over T with eps=10^-2. With the last assumption I have a solution. – Alex Trounev Jul 22 '20 at 15:12

We can integrate in 3 steps:

Integrate[(yp/(b + yp^2)^(3/2)), {yp, 0, Infinity},
Assumptions -> b > 0]*(x - xp) /. {b ->
g (x - xp)^2} //Simplify

Out[]: (x - xp)/Sqrt[g (x - xp)^2]


So we have intyp=1/Sqrt[g] as results and it means that Q[g,x] not depends on x. Next step:

Integrate[(Exp[-ypp^2/(8 T)])*(-1 +
ypp^2/(8 T)) (ypp/(g*(xp - xpp)^2 + ypp^2)^(3/2)), {ypp, 0, Infinity},  Assumptions ->{...}]


I made substitutions s->ypp/Sqrt[8 T], a->g*(xp - xpp)^2/(8 T), it turns into

Integrate[
Exp[-s^2] (-1 + s^2) s/(a + s^2)^(3/2), {s, 0, Infinity},
Assumptions -> {a > 0}]

Out[]= -((1 + a)/Sqrt[a]) +
1/2 (3 + 2 a) E^a Sqrt[\[Pi]] Erfc[Sqrt[a]]


Restoring all coefficient coming from ypp normalization on Sqrt[8 T] we have

intypp=
With[{a = g*(xp - xpp)^2/(8 T)},
Sqrt[8 T]/(8 T)^(3/2) Sqrt[
8 T] (-((1 + a)/Sqrt[a]) + 1/2 (3 + 2 a) E^a Sqrt[\[Pi]] Erfc[Sqrt[a]]) //
Simplify]

Out[]=
(-((8*(1 + (g*(xp - xpp)^2)/(8*T)))/Sqrt[(g*(xp - xpp)^2)/T]) +
E^((g*(xp - xpp)^2)/(8*T))*Sqrt[2*Pi]*(3 + (g*(xp - xpp)^2)/(4*T))*
Erfc[Sqrt[(g*(xp - xpp)^2)/T]/(2*Sqrt[2])])/(8*Sqrt[T])


Therefore we get integrand

intp intpp (xp - xpp)/T^2 (EllipticTheta[3, 1/2 Pi (xp + x0), Exp[-Pi^2 T]] +
EllipticTheta[3, 1/2 Pi (xp - x0), Exp[-Pi^2 T]])*(EllipticTheta[3,
1/2 Pi (xpp + x0), Exp[-Pi^2 T]] +
EllipticTheta[3, 1/2 Pi (xpp - x0), Exp[-Pi^2 T]])


And finally we have

int2[g_, x0_, T_, xp_,
xpp_] := (EllipticTheta[3, 1/2 Pi (xp + x0), Exp[-Pi^2 T]] +
EllipticTheta[3, 1/2 Pi (xp - x0),
Exp[-Pi^2 T]])*(EllipticTheta[3, 1/2 Pi (xpp + x0),
Exp[-Pi^2 T]] +
EllipticTheta[3, 1/2 Pi (xpp - x0),
Exp[-Pi^2 T]])/(8 T^3) (-2 Sqrt[
2 T] (1 + (g (xp - xpp)^2)/(8 T))/Sqrt[g ] +
1/2 E^((g (xp - xpp)^2)/(8 T))
Sqrt[\[Pi]] (3 + (g (xp - xpp)^2)/(4 T)) Erfc[Sqrt[(
g (xp - xpp)^2)/T]/(2 Sqrt[2])]*(xp - xpp))/Sqrt[g] Sqrt[8 T]


This is what we can work with. But it diverges at T->0. We can perform numerical integration cutting of temperature limits as follows

hardintegral[g_?NumericQ] :=
NIntegrate[
int2[g, x0, T, xp, xpp], {x0, 0, 1}, {xp, 0, 1}, {xpp, 0, 1}, {T,
10^-2, 10}, AccuracyGoal -> 2, PrecisionGoal -> 2]


The upper limits T does not matter since integrand very fast vanished at T>1 , but T=10^-2 is essential for fast calculations. So we make a table and plot

 lst = Table[{g, hardintegral[g]}, {g, .1, 2, .1}]
ListLinePlot[lst, PlotRange -> All, FrameLabel -> {"g", "Q"},
Frame -> True]


I can recommend to use function Q[g,T] for future research. We can define function

Q[g_?NumericQ, T_?NumericQ] :=
NIntegrate[
int2[g, x0, T, xp, xpp], {x0, 0, 1}, {xp, 0, 1}, {xpp, 0, 1},
AccuracyGoal -> 2, PrecisionGoal -> 2]


Now we plot it to check singularity at g->0 and T->0:

Plot3D[Q[g, T], {g, .1, 2}, {T, .1, 2}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic, PlotRange -> All]


• +1. Nice breakdown, However, I think the cut-off T=10^-2 is arbitrary and unwarranted. Nonetheless, it does bring focus on the problem, and that is helpful. Not convinced that it diverges at T->0, but there are problems around 10^-6 < T < 10^-2. – Michael E2 Jul 22 '20 at 18:18
• @MichaelE2 I agree with you that my approach is only 3 steps of 5 to solve this problem. Next 2 steps could be around properties of function Q[g_, T_]:=NIntegrate[ int2[g, x0, T, xp, xpp], {x0, 0, 1}, {xp, 0, 1}, {xpp, 0, 1}]. – Alex Trounev Jul 22 '20 at 19:24
• @AlexTrounev Thank you very much. This looks on the right track. I have a question: as it is, you are performing the T integration after ypp (which I like) - but doesn't this mean there shouldn't be a singularity in T? – SaMaSo Jul 23 '20 at 22:14
• @SaMaSo There is a singularity 1/T^2on a hypersurface ypp=0. And we can't remove it, but we can cut of T as I did. It is a reason why your hardintegral has no stable result. – Alex Trounev Jul 24 '20 at 0:21
• @AlexTrounev I am not fully sure again, because ypp=0 is also multiplied in the numerator, so it is ambiguous. – SaMaSo Jul 24 '20 at 7:35

To reduce the number of integration in NIntegrate seems reasonable. The effects are somehow dependent on the choices of options for NIntegrate.

Choices are

Values for lower integral bounds larger than zero. Values replacing the infinite upper bound of the integral to a meaningful numerical integration value.

The default method is GlobalAdaptive. This can be changed.

GlobalAdaptive has the method option MaxErrorIncreases that governs the time needed or spend by NIntegrate in this question very much. MaxErrorIncreases takes time and is always used to full extend.

WorkingPrecision should be set high according to persisting error messages.

Leave most of the work to NIntegrate is in general very good advice. A best practice recommendation from Wolfram Inc and its competitors.

This works moderately crude:

Nn = 10^6; eps = 10^-8; Table[
ListPlot[Last[
Reap[NIntegrate[(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]])*(0.2727575560073645 -
0.6266570686577505 E^(1.6801824043209879 T) Sqrt[T] +
0.6266570686577502 E^(1.6801824043209879 T) Sqrt[T]
Erf[1.2962185017661907 Sqrt[T]])*(Integrate[
Exp[-z^2/8 T]*(z^3/(g + z^2)^(3/2)), {z, 0, \[Infinity]},
Assumptions -> g > 0 && T > 0]), {x0, eps, 1}, {T, eps,
Nn}, {xp, eps, 1}, {xpp, eps, 1}, {g, eps, 2},
Method -> {str, "MaxErrorIncreases" -> 15},
WorkingPrecision -> 50, EvaluationMonitor :> Sow[g]]]],
"Trapezoidal", "DoubleExponential"}}]


This output is a bunch of error messages and the desired plot.

Higher upper bound and smaller lower bound offer convergence independent of the value of MaxErrorIncrease.

This might perform even better with compilations and parallel processing.

I soon will continue this answer.

• Thank you for this. There was a mistake in the formula I had provided for context (the one for Mathematica code was correct though), and I have edited this. – SaMaSo Jul 22 '20 at 1:59
• I don't seem to understand to understand your approach, especially what those numbers in the exponentials etc are, and why there is sqrt of T. Would be good if you could add explanations later – SaMaSo Jul 22 '20 at 2:03
• for Methods, based on the previous post that I linked in the question, it could be that the GaussKronrod is the best. – SaMaSo Jul 22 '20 at 2:04
• The output of you code is i.stack.imgur.com/vrmjJ.png, which seems useless. Why do you think it's worth posting? – Michael E2 Jul 22 '20 at 13:31
• @SteffenJaeschke Option WorkingPrecision -> 50 is not applicable to the integrand with numbers like 0.2727575560073645, – Alex Trounev Jul 22 '20 at 14:45