Numerical Integral of a "Complicated" Function

Question

I am trying to calculate the numerical integral of a 4D function, but NIntegrate keeps saying that is not able to perform a good estimate.

My function comes from an exponential of four derivatives of something else so it is a "big" chunck (around 3 lines), but it is not oscillatory, it does not have divergences and the final value cannot be zero because it is always positive. It actually looks more or less like a 4D Gaussian.

The function is:

L[a_,c_,w_,u_]=    Exp[
  -(Y[1] - f[X[1], a, c, w, u])^2/(2*0.005^2) - (Y[2] - f[X[2], a, c, w, u])^2/(2*0.005^2)
  -(Y[3] - f[X[3], a, c, w, u])^2/(2*0.005^2) - (Y[4] - f[X[4], a, c, w, u])^2/(2*0.005^2)
    ]

where:

X={0.266667, 0.5, 1, 1.16667}
Y={0.867, 0.596, 0.0689, -0.00554}

and f[q,a,c,w,u] is given by:

f[q_,a_,c_,w_,u_]=(1/(c (c^2 + a^2 π^2) q)
   3 a π Csch[a π q] (c Cos[
      c q] (-1 - c^4 u - 5 a^4 π^4 u + 
       c^2 (10 a^2 π^2 u - w) + 3 a^2 π^2 w + 
       2 a^2 π^2 Csch[a π q]^2 (10 (c^2 - 5 a^2 π^2) u + 3 w - 
          60 a^2 π^2 u Csch[a π q]^2)) + 
    a π Coth[a π q] (1 + 5 c^4 u + a^4 π^4 u - a^2 π^2 w + 
       c^2 (-10 a^2 π^2 u + 3 w) + 
       6 a^2 π^2 Csch[a π q]^2 (-10 (c - a π) (c + a π) u - w + 
          20 a^2 π^2 u Csch[a π q]^2)) Sin[c q]))/1.653

In theory I am interested in the integral from -infinity to infinity of the four parameters a,c,w,u. But I know that most of the time the function is zero, so it will be enough to integrate {a, 0.5,0.7}, {c, 3.05, 3.25}, {w, 0.02, 0.2}, {u, -0.003,-0.0001}.

I have tried with the automatic configurations for NIntegrate and it says:

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 6.0484948162078716`*^-9 and 1.3350642633041837`*^-12 for the integral and error estimates.

Do you know a better NIntegrate strategy for my problem? I was thinking of replacing my function by an approximate interpolating function, but I don't know if this is the best I can do. I have tried MonteCarlo methods but after increasing the Iterations until it finally says it converges, the results changes by 10% if you try the command again.

Answer 1

If the only thing you need is the value of this integral, I would try to use Monte Carlo integration by using Mathematica for generation of random input, evaluation of your function and averaging results. If you are not familiar with this technique, please see simple explanation on the Internet, for example, a good review with illustrations in https://www.cs.dartmouth.edu/~wjarosz/publications/dissertation/appendixA.pdf. You may also look at presentation in http://www.cs.utah.edu/~edwards/research/mcIntegration.pdf. Monte Carlo method is one of the most popular when it comes to computation of multi-dimensional integral. Nick L.

Answer 2

The general rule with multivariate integrals should be, to do as much integration as you can analytically.

Here the u und w integration can be done analytically and then the remaining two-dimensional numerical integration is better to do.

For analytical integration it is better to rationalize the functions

X = Rationalize[{0.266667, 0.5, 1, 1.16667}, 0]
Y = Rationalize[{0.867, 0.596, 0.0689, -0.00554}, 0]

f[q_, a_, c_, w_, u_] = (1/(c (c^2 + a^2 \[Pi]^2) q) 3 a \[Pi] Csch[
  a \[Pi] q] (c Cos[
     c q] (-1 - c^4 u - 5 a^4 \[Pi]^4 u + 
      c^2 (10 a^2 \[Pi]^2 u - w) + 3 a^2 \[Pi]^2 w + 
      2 a^2 \[Pi]^2 Csch[
         a \[Pi] q]^2 (10 (c^2 - 5 a^2 \[Pi]^2) u + 3 w - 
         60 a^2 \[Pi]^2 u Csch[a \[Pi] q]^2)) + 
   a \[Pi] Coth[
     a \[Pi] q] (1 + 5 c^4 u + a^4 \[Pi]^4 u - a^2 \[Pi]^2 w + 
      c^2 (-10 a^2 \[Pi]^2 u + 3 w) + 
      6 a^2 \[Pi]^2 Csch[
         a \[Pi] q]^2 (-10 (c - a \[Pi]) (c + a \[Pi]) u - w + 
         20 a^2 \[Pi]^2 u Csch[a \[Pi] q]^2)) Sin[c q]))/(1653/
 1000) // Simplify

L[a_, c_, w_, u_] = 
   Simplify[Rationalize[
    Exp[-(Y[[1]] - f[X[[1]], a, c, w, u])^2/(2*0.005^2) - (Y[[2]] - 
     f[X[[2]], a, c, w, u])^2/(2*0.005^2) - (Y[[3]] - 
     f[X[[3]], a, c, w, u])^2/(2*0.005^2) - (Y[[4]] - 
     f[X[[4]], a, c, w, u])^2/(2*0.005^2)], 0], 
   Assumptions -> {a, c, w, u} \[Element] Reals]

L is of the form Exp[L[a,c,w,u][[2]]] and can be writen as Exp[r + s u + t u^2] Calculate the coefficients r,s,t

cl = CoefficientList[Collect[L[a, c, w, u][[2]], u], u]

(*  A very large output was generated ...   *)

cl // Length

(*  3   *)

rule1 = Thread[{r, s, t} -> 
  Simplify[cl, Assumptions -> {a, c, w, u} \[Element] Reals]];

The integration over u can be done easily. According to graphical tests t is assumed to be t<0, although no deep proof is done here.

 int11 = Integrate[Exp[r + s u + t u^2], {u, -\[Infinity], \[Infinity]}, 
    Assumptions -> t < 0]

 (*  (E^(r - s^2/(4 t)) Sqrt[\[Pi]])/Sqrt[-t]    *)

Nearly the same can be done with the w-integration

 exp2 = int11 /. rule1;

 exp2 // Length

 (*  4   *)

 exp2 == exp2[[1]]*Exp[exp2[[2, 2]]]*exp2[[3]]*exp2[[4]] // Simplify

 (*   True   *)

 fac = exp2[[1]]*exp2[[3]]*exp2[[4]];

Fortunately fac does not depend on w

 clw = Simplify[CoefficientList[exp2[[2, 2]], w], 
    Assumptions -> {a, c, w} \[Element] Reals]

 (*  A very large output was generated ...   *)

 clw // Length
 (*  3  *)

 rule2 = Thread[{o, p, q} -> clw];

 int12 = Integrate[
   fact*Exp[o + p*w + q*w^2], {w, -\[Infinity], \[Infinity]}, 
     Assumptions -> q < 0]

 (*    (E^(o - p^2/(4 q)) fact Sqrt[\[Pi]])/Sqrt[-q]   *)

 fact = Simplify[fac, Assumptions -> {a, c} \[Element] Reals];

Finaly you get the integrand for the twodimensional numerical integration

 integrand[a_, c_] = int12 /. rule2

Since this integrand has a few thousands of terms, it takes quite a long time to calculate. And it is very sensitive to input. MachinePrecision is not enough here. Examples:

 N[integrand[4/10, 1/2], 3]

 (*   0.0000400 + 0.*10^-8 I    *)

 N[integrand[3/10, 1/2], 3] // Timing

 (*    {2.688, 7.78*10^-38 + 0.*10^-41 I}    *)

Calculate an array of the integrand in order to get a rough overview where to integrate

 (tab4 = Table[{a, c, N[integrand[a, c], 3]}, {a, 1/40, 4, 1/20}, {c, 
  1/10, 20, 2/10}];) // Timing

  ListPlot3D[Flatten[Re[tab4], 1], PlotRange -> All, 
    AxesLabel -> {a, c, integrand}]

But attention, this and the following integration takes several hours (I aborted therefore)

 int1 = NIntegrate[integrand[a, c], {a, 0, 1}, {c, 0, 5}, 
    Method -> {"MultidimensionalRule", "Generators" -> 9}, 
      MaxRecursion -> 100, WorkingPrecision -> 45]