# Numerical Integral of a "Complicated" Function

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.

• What error message do you get? If it's a convergence error, then probably the function has been evaluated hundreds or thousands of times, and it's probably not because of simplification. How bad is the error estimate in the message? -- "SymbolicProcessing" -> 0 would prevent symbolic manipulation, but that is usually effective only when you know what Method strategy/rule to use. -- 2D integrals are already numerically hard, and >2D are really hard. It may simply be hard. -- You could increase MaxRecursion. Sep 25, 2016 at 16:35
• It will be easier to help you if you provide the complete code. If the function is too big and you can't reduce it, you can post the code on a website, such as pastebin.com. Sep 25, 2016 at 19:22
• Thanks, I have added more detail. I tried the "SymbolicProcessing" -> 0 and didn't work, so you were right that it was not that. Sep 26, 2016 at 19:09
• The largest value of the function is around 10^(-94703) within the range of values that you give and that maximum occurs at the maximum values of a and d and minimum values of w and u. So the range of values is not at all centered near the maximum value of the function. You'll likely need to rationalize all of the coefficients, center the range of values better, and scale the resulting value of the function.
– JimB
Oct 26, 2016 at 20:38
• Is the issue that in the definition of L you mean, for example, X[[1]] with double brackets and not X[1] with single brackets? Oct 27, 2016 at 1:39

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]

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.

• I have used my own code of Metropolis MonteCarlo in python to do the Integral and it was not stable (the result changed for every trial even for 100,000 points. I read the reviews you suggested to see if I can improve it. Thank you very much Sep 26, 2016 at 19:11