# Compute unstable integral with high precision

I want to calculate $$\int_{\mathbb{R}^n} e^{-t_1^4-t_2^4- ...-t_n^4-(1-t_1-t_2-...-t_n)^4} \text{d}t_1 \text{d}t_2 ... \text{d}t_n$$ with the highest possible n.

I tried with

Module[{n = 3, intVars, poly},
intVars = Table[{Subscript[t, i], -∞, ∞}, {i, 1, n - 1}];
poly = Sum[Subscript[t, i]^4, {i, 1, n - 1}] + (1 - Sum[Subscript[t, i], {i, 1, n - 1}])^4;

Print[E^-poly];

NIntegrate[E^-poly, ##, Method -> {Automatic, "SymbolicProcessing" -> 0}] & @@ intVars
]


which seems to work well only for n<6. Is there any trick I can use?

• Play with NIntegrate[E^-poly, ##, Method -> {"QuasiMonteCarlo", "SymbolicProcessing" -> 0}, AccuracyGoal -> 3, PrecisionGoal -> 3, WorkingPrecision -> 30] & @@ intVars]. It works for $n=19$. Commented Jun 23, 2019 at 17:35
• You cannot calculate a high-dimensional integral with high precision. Read the tutorial: "For low-dimensional integrals, the default setting for PrecisionGoal is related to WorkingPrecision. For high-dimensional integrals, it is typically taken to be a fixed value, usually 2. " Commented Jun 23, 2019 at 17:38
• @Alex Trounev: You are right. My idea is unsufficiently considered. Commented Jun 23, 2019 at 17:42

One idea I like to use for these kinds of integrals is to add an auxiliary variable and a Dirac delta function, convert the Dirac delta function to it's integral formulation, and then do a bunch of simple 1D integrals. For your case, this would proceed as follows, starting from:

$$\underset{t\in \mathbb{R}^n}{\int }e^{-t_1^4-t_2^4-\ldots \ -t_n^4-\left(1-t_1-t_2-\ldots -t_n\right){}^4}$$

Introduce the auxiliary variable $$s$$ and add a Dirac delta function:

$$\underset{t\in \mathbb{R}^n}{\int }\underset{s\in \ \mathbb{R}}{\int }e^{-t_1^4-t_2^4-\ldots -t_n^4-s^4} \delta(1-t_1-t_2-\ldots -t_n-s)$$

Next, introduce the integral formulation of the Dirac delta function:

$$\frac{1}{2 \pi }\underset{t\in \ \mathbb{R}^n}{\int }\underset{s\in \mathbb{R}}{\int }\underset{u\in \ \mathbb{R}}{\int }e^{-t_1^4-t_2^4-\ldots -t_n^4-s^4} e^{i u-i t_1 u-i t_2 \ u-\ldots -i t_n u-i s u}$$

Finally, we can do all of the $$t$$ and $$s$$ integrals to obtain:

$$\frac{1}{2 \pi } \int_{-\infty \ }^{\infty } e^{i u} g(u)^{n+1} \, du$$

where:

$$g(u)=\int_{-\infty }^{\infty } e^{-i t u-t^4} \, dt$$

Now, let's have Mathematica do these integrals:

g[u_] = Sqrt[2 Pi] FourierTransform[Exp[-t^4], t, u]


2 Gamma[5/4] HypergeometricPFQ[{}, {1/2, 3/4}, u^4/256] - 1/4 u^2 Gamma[3/4] HypergeometricPFQ[{}, {5/4, 3/2}, u^4/256]

So, the desired integral has become:

int[n_, opts:OptionsPattern[NIntegrate]] := (1/(2 Pi)) NIntegrate[
Cos[u] g[u]^(n+1),
{u, -Infinity, Infinity},
opts,
]


Now, g[u] is real:

Refine[g[u] ∈ Reals, u ∈ Reals]


True

so I use Cos[u] instead of Exp[I u] since the integral is real. I also customize the integration method. Let's check:

int[2, WorkingPrecision->20]
int[2, WorkingPrecision->40]
int[2, WorkingPrecision->60]


1.4733172914977911077

1.473317291497785926905017339845596712841

1.47331729149778592690501733984559670949096610342311667206502

The result seems correct, and improves with higher working precision. Now, for higher orders:

int[4, WorkingPrecision -> 20]
int[5, WorkingPrecision -> 20]
int[6, WorkingPrecision -> 20]
int[40, WorkingPrecision -> 20]


4.4732305211180348293

7.7543594355221995796

13.461688085347942892

4.0351905913672630176*10^9

• You change the order of the integration in the above. This is not a simple matter for improper integrals. Next, in the integral presentation of DiracDelta the principal value is taken. Commented Jun 23, 2019 at 19:54
• @user64494 What makes you say g[u] is complex valued? I don't think it is. As for order of integration, I haven't been rigorous in proving that it is possible to do so, but the results seem correct, and that is sufficient for me. Commented Jun 23, 2019 at 19:58
• You are right concerning g[u]. However, other remarks of me remain open. The statement "the results seem correct" is not based. Commented Jun 23, 2019 at 20:02
• You are right: I don't execute g[u_]. It works well. Sorry for the trouble. Commented Jun 23, 2019 at 20:37
• @CarlWoll your solution is very elegant and it seems to do the job excellently. Thank you very much for the effort invested! Just for curiosity, where is used the method of inserting delta and then Fourier transform? To decouple integration variables? Do you have an intuitive interpretation for this?
– m137
Commented Jun 23, 2019 at 21:28

We can increase n, reducing the accuracy of calculations, for example

f[m_, p_] :=
Block[{n = m, $$MinPrecision = p,$$MaxPrecision = p},
intVars =
Table[{Subscript[t, i], -\[Infinity], \[Infinity]}, {i, 1, n}];
poly = Sum[
Subscript[t, i]^4, {i, 1,
n}] + (1 - Sum[Subscript[t, i], {i, 1, n}])^4;
NIntegrate[E^-poly, ##, PrecisionGoal -> p/2,
AccuracyGoal -> p/2] & @@ intVars]
f[4, 6]

(* 4.47148*)

f[5, 4]

(* 7.75615*)

f[6, 2]

(* 15.8289*)