0
$\begingroup$

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[1], x[2], v, x[4], x[5]};
integrand = f[Ainvn.w]*dAi;
NIntegrate[integrand /. v -> 1,
 {x[1], -∞, ∞},{x[2],-∞,∞}, {x[4], -∞, ∞}, {x[5],-∞, ∞}, {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] -> 1, x[2] -> 1, x[4] -> 1, x[5] -> 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?

$\endgroup$
  • $\begingroup$ @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. $\endgroup$ – José Antonio Díaz Navas Feb 13 '18 at 17:57
  • $\begingroup$ @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? $\endgroup$ – user39756 Feb 13 '18 at 18:05
  • $\begingroup$ You should evaluate integrand, not the pieces. That's the way to find out the problem. $\endgroup$ – Michael E2 Feb 13 '18 at 18:05
  • $\begingroup$ @MichaelE2 If I evaluate integrand, I obtain underflow. I do not know how to solve that. $\endgroup$ – user39756 Feb 13 '18 at 18:07
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Feb 13 '18 at 18:27
1
$\begingroup$

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

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

Original answer:

This gives a positive result:

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

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];
$\endgroup$
  • $\begingroup$ Thank you. I do not know much about Mathematica. Your comment has been very interesting. $\endgroup$ – user39756 Feb 13 '18 at 19:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.