# NIntegrate gives $0$ for a positive integrand

I have a positive function called "integrand" and I want to perform numerical integration. The following code creates the function "integrand" and performs the numerical integration:

NN = 300; n = 30; M = 5;
k = a*(M + 1)^2/NN;
ei[l_] := 1 + 2*k*Cos[Pi*l/(M + 1)] - 2*k;
prod = 1;
For[l = 1, l <= M, l++,
prod = prod*ei[l];
];
dA = Simplify[prod];
R = ConstantArray[0, {M, M}];
For[jj = 1, jj <= M, jj++,
For[l = 1, l <= M, l++,
R[[jj, l]] = Sin[l*Pi*jj/(M + 1)];
];
];
Q = Simplify[Sqrt[2/(M + 1)]*R];
Dinvn = ConstantArray[0, {M, M}];
For[l = 1, l <= M, l++,
Dinvn[[l, l]] = 1/ei[l]^n;
];
Ainvn = Q.Dinvn.Transpose[Q];
dAi = 1/dA^n;
f[x_] := PDF[
ProductDistribution[NormalDistribution[0, 1],
NormalDistribution[0, 1], NormalDistribution[0, 1],
NormalDistribution[0, 1], NormalDistribution[0, 1]], x];
w = {x, x, v, x, x};
integrand = f[Ainvn.w]*dAi;
NIntegrate[integrand /. v -> 1,
{x, -∞, ∞},{x,-∞,∞}, {x, -∞, ∞}, {x,-∞, ∞}, {a, 1, 2}]


The output:

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option.

0.


We can analyze "integrand" by doing some examples:

dAi /. a -> 1.4

2.06591*10^30

Ainvn.w /. {x -> 1, x -> 1, x -> 1, x -> 1, a -> 1.5, v -> -0.5}

{-6.71576*10^13, 1.1632*10^14, -1.34315*10^14, 1.1632*10^14, -6.71576*10^13}


Notice that $f$ is the joint density of five independent $\text{Normal}(0,1)$ distributions. Anyway, the thing is that $f(\{10^{13},10^{13},10^{13},10^{13},10^{13}\})\approx 0$, but not exactly zero. Mathematica gives underflow. The result of the evaluation of $f$ should be a number very close to $0$, that with dAi should compensate. Or maybe it is a problem of the integration performed by NIntegrate. How can I modify my code so that the result obtained by NIntegrate is correct?

• @bills yes, it has. I have increased the AccuracyGoal, WorkingPrecission and it returns a very small number -3.49697474357660*10^-221213430006. I also work with double precision NN=300 and n=30 and with Method -> "MultidimensionalRule", AccuracyGoal -> 10, \ WorkingPrecision -> 15. I do not know how reliable (and useful) is this result. – José Antonio Díaz Navas Feb 13 '18 at 17:57
• @JoséAntonioDíazNavas Do you mean to put Method -> "MultidimensionalRule", AccuracyGoal -> 10, WorkingPrecision -> 15 in NIntegrate? What do you mean by double precision in NN=300 and n=30? – user39756 Feb 13 '18 at 18:05
• You should evaluate integrand, not the pieces. That's the way to find out the problem. – Michael E2 Feb 13 '18 at 18:05
• @MichaelE2 If I evaluate integrand, I obtain underflow. I do not know how to solve that. – user39756 Feb 13 '18 at 18:07
• It means your problem in not well scaled for machine precision. The values computed, including intermediate ones, that are to be treated as nonzero should fit comfortably between the magnitudes 10^±308. The culprit is the argument to Exp in integrand. Haven't thought of a fix yet. – Michael E2 Feb 13 '18 at 18:27

This is probably more accurate, since the integrand underflows outside the region:

NIntegrate[
integrand /. v -> 1,
{x, -2, -1, 0}, {x, -1, 0, 1}, {x, -1, 0, 1}, {x, -2, -1, 0}, {a, 1.8, 1.9, 2},
MinRecursion -> 2, PrecisionGoal -> 4]
(*  1.4796*10^41  *)


This gives a positive result:

NIntegrate[
integrand /. v -> 1,
{x, -∞, -1, ∞}, {x, -∞, 0, ∞}, {x, -∞, 0, ∞}, {x, -∞, -1, ∞}, {a, 1, 1.9, 2},
MinRecursion -> 1, PrecisionGoal -> 4]
(*  5.22114*10^43  *)


THIS IS AN EXTENDED COMMENT RATHER THAN AN ANSWER.

For efficiency and brevity of code, recommend that you avoid use of For constructs, e.g.,

NN = 300; n = 30; M = 5;
k = a*(M + 1)^2/NN;
ei[l_] := 1 + 2*k*Cos[Pi*l/(M + 1)] - 2*k;

dA = Times @@ ei@Range[M] // Simplify;

R = Array[Sin[#1*#2*Pi/(M + 1)] &, {M, M}];

Q = Sqrt[2/(M + 1)]*R;

Dinvn = DiagonalMatrix[1/(ei@Range[M])^n];

• Thank you. I do not know much about Mathematica. Your comment has been very interesting. – user39756 Feb 13 '18 at 19:56