# Multidimensional Numeric Integration

I can't manage to calculate numerically the following multidimensional integral: $I(X,Y,T)=\int_{0}^{T}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{-(T-\tau)(\lambda_1^2+\lambda_2^2)^2+I((X-\zeta_1)\lambda_1+(X-\zeta_2)\lambda_2)}(\frac{\partial^2}{\partial \zeta_1^2}+\frac{\partial^2}{\partial \zeta_2^2})e^{-\zeta_1^2-\zeta_2^2-\tau^2}d\lambda_1d\lambda_2d\zeta_1d\zeta_2d\tau$

The code i use is as follows:

In[12]:= T = 1/10;
X = 1/10;
Y = 0;

In[6]:= NIntegrate[
Exp[-(T - tau)*(lambda1^2 + lambda2^2)^3]*
Cos[(X - zeta1)*lambda1 + (Y - zeta2)*
lambda2]*(D[Exp[-zeta1^2 - zeta2^2 - tau^2], {zeta1, 2}] +
D[Exp[-zeta2^2 - zeta1^2 - tau^2], {zeta2, 2}])/(2*Pi)^(3/
5), {lambda1, -Infinity, Infinity}, {lambda2, -Infinity,
Infinity}, {zeta1, -Infinity, Infinity}, {zeta2, -Infinity,
Infinity}, {tau, 0, T}]

During evaluation of In[6]:= 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. >>

During evaluation of In[6]:= 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 -0.858336 and 0.014313811647397338 for the integral and error estimates. >>

Out[6]= -0.858336


I tryied modifying integration proporties such as:"WorkingPrecison","MaxErrorIncreases". Nothing has fixed the problem yet.

• I confess to not studying the integral closely, but something about it has me wondering if you could cast the problem geometrically then lay it out with something like a Cholesky Decomposition? Just a thought. Mar 23, 2015 at 13:18

I think it just does not converge. In general high dimensional (>4) integral always converge fastest with monte carlo method.

k = 0;
Dynamic[k]
NIntegrate[
Exp[-(T - tau)*(lambda1^2 + lambda2^2)^3]*
Cos[(X - zeta1)*lambda1 + (Y - zeta2)*
lambda2]*(D[Exp[-zeta1^2 - zeta2^2 - tau^2], {zeta1, 2}] +
D[Exp[-zeta2^2 - zeta1^2 - tau^2], {zeta2, 2}])/(2*Pi)^(3/
5), {lambda1, -Infinity, Infinity}, {lambda2, -Infinity,
Infinity}, {zeta1, -Infinity, Infinity}, {zeta2, -Infinity,
Infinity}, {tau, 0, T}, Method -> "AdaptiveMonteCarlo",
MaxPoints -> 100000, EvaluationMonitor :> k++]


Output:

100100
During evaluation of In[44]:= NIntegrate::maxp: The integral failed to converge after 100100 integrand evaluations. NIntegrate obtained -0.760721 and 0.019040561157440964 for the integral and error estimates. >>
-0.760721

• Thanks for the tip. Youre code doesnt work though
– user16247
Mar 23, 2015 at 17:02

To try to avoid some badness from numerical integration of the Cos, I tried to simplify things by first integrating out the tau

Integrate[
Exp[-(T - tau)*(lambda1^2 + lambda2^2)^3]*(D[
Exp[-zeta1^2 - zeta2^2 - tau^2], {zeta1, 2}] +
D[Exp[-zeta2^2 - zeta1^2 - tau^2], {zeta2, 2}])/(2*Pi)^(3/
5), {tau, 0, T}]


Which is:

1/\[Pi]^(1/10)2^(2/5) E^(
1/20 (lambda1^2 + lambda2^2)^3 (-2 + 5 (lambda1^2 + lambda2^2)^3) -
zeta1^2 -
zeta2^2) (-1 + zeta1^2 +
zeta2^2) (Erf[1/2 (lambda1^2 + lambda2^2)^3] +
Erf[1/10 - 1/2 (lambda1^2 + lambda2^2)^3])


Then plugging that back in,

NIntegrate[
Cos[(X - zeta1)*lambda1 + (Y - zeta2)*lambda2] 1/\[Pi]^(1/10) 2^(2/5)
E^(1/20 (lambda1^2 + lambda2^2)^3 (-2 +
5 (lambda1^2 + lambda2^2)^3) - zeta1^2 -
zeta2^2) (-1 + zeta1^2 +
zeta2^2) (Erf[1/2 (lambda1^2 + lambda2^2)^3] +
Erf[1/10 - 1/2 (lambda1^2 + lambda2^2)^3]), {lambda1, -Infinity,
Infinity}, {lambda2, -Infinity, Infinity}, {zeta1, -Infinity,
Infinity}, {zeta2, -Infinity, Infinity}]


Unfortunately, it appears that there is some problems, such as around

{lambda1,lambda2,zeta1,zeta2} = {-3.06813*10^6,0.,0.,0.}


Which appears to be due to the Erf going to 0 but the Exp going to infinity.

I suggest trying to simplify things other ways before throwing it into MMA.