# Convergence of a multiple numerical integration

I would like to define the following function $F(A,\nu)$, which is the result of a numerical integration, and needs the following function definitions.

χ[x_] := (1 - x^2) Exp[-x^2/2.]
v1[x_, T_, A_, ν_] :=
Sqrt[2 ν A] (-Tanh[Sqrt[A/ 2. ν] (x - Sqrt[-Sqrt[A] T])] +
Tanh[Sqrt[A/ 2. ν] (x + Sqrt[-Sqrt[A] T])])
u1[x_, T_, A_, ν_] := v1[x, T, A, ν]/v1[0., T, A, ν]
P[x_, T_] := -x Abs[x]/(8. T^2)

F[A_?NumericQ, ν_?NumericQ] :=
NIntegrate[
P[x, T] u1[x, T, A, ν] P[y, T] u1[y, T, A, ν] χ[
x - y], {y, - Infinity, Infinity}, {x, - Infinity,
Infinity}, {T, -Sqrt[A], -1/Sqrt[A]}, WorkingPrecision -> 6]


$F(A,\nu)$ is the integral of the product of the function $P(x,T)\, u_1(x,T,A,\nu)$ (displayed in the image) at two different points, x and y, and a kernel $(1-x^2) \exp(-x^2/2)$. The image was generated for $T=-10.$, $A=10.$ and $\nu=100.$ (these are typical values for these variables, $T$ is always negative and $A$ is always larger than 1).

The evaluation of this function takes a long time and returns some numerical warnings. For instance,

F[1000., 100.]


returns these messages:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

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 288.7638201609691666670328291224718571098961230923382242956. and 5.401118408723842142699384811776130162688582048448406609656. for the integral and error estimates.

What options can I use for better convergence of this integral, and how can I decide which method is the best, from looking at the properties of the integrand?

• Your working precision is too small. Experiment using higher working precisions and small precision goals. – Anton Antonov Jul 17 '16 at 20:37
• what are you trying to accomplish setting WorkingPrecision less than machine precision? In any case if you want to go to high precision you,ll need to specify exact values for your literals ( make 8. 8, etc) – george2079 Jul 17 '16 at 21:04
• @MichaelE2 Thanks! I thought something like that might be in place... – Anton Antonov Jul 17 '16 at 21:56
• @AntonAntonov Oops, I misread the limits. 0 is not in the T interval, so the integrand is ok (or if not, I haven't discovered the pathological part). – Michael E2 Jul 17 '16 at 22:04

This is not a complete answer, the code provides results that demonstrate the need for more detailed investigation.

### Re-definition

The re-defintion uses exact numbers and adds options argument to F.

Clear["Global*"];

χ[x_] := (1 - x^2) Exp[-x^2/2]
v1[x_, T_, A_, ν_] :=
Sqrt[2 ν A] (-Tanh[Sqrt[A/2 ν] (x - Sqrt[-Sqrt[A] T])] +
Tanh[Sqrt[A/2 ν] (x + Sqrt[-Sqrt[A] T])])
u1[x_, T_, A_, ν_] := v1[x, T, A, ν]/v1[0, T, A, ν]
P[x_, T_] := -x Abs[x]/(8 T^2)

F[A_?NumericQ, ν_?NumericQ, opts : OptionsPattern[]] :=
NIntegrate[
P[x, T] u1[x, T, A, ν] P[y, T] u1[y, T, A, ν] χ[x - y],
{y, -∞, ∞}, {x, -∞, ∞}, {T, -Sqrt[A], -1/Sqrt[A]}, opts]


### Experiments

Using the default working precision produces result that seems to be wrong:

In[63]:= AbsoluteTiming[F[1000, 100]]

Out[63]= {5.68169, 288.231}


Here is a result with higher working precision and denser sampling points:

In[64]:= AbsoluteTiming[
F[1000, 100, WorkingPrecision -> 30, PrecisionGoal -> 6,
MinRecursion -> 2,
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]
]

Out[64]= {72.3681, 292.119919516539460106135415685}


And one more without singularity handling:

In[65]:= AbsoluteTiming[
F[1000, 100, WorkingPrecision -> 60, PrecisionGoal -> 6,
MinRecursion -> 2,
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0,
"SingularityHandler" -> None}]
]

Out[65]= {65.9541, 293.123327189195721528823498643979792005706149710632238819309}


Having these kind of different results means that the integral has some pathological parts (narrow peaks, oscillations, singular points/curves...) This have to be investigated further and NIntegrate`'s advanced documentation provides a good start.