TransWikia.com

How can I reduce the error estimate of numerical integration in Mathematica?

Mathematica Asked on December 31, 2021

Before introducing the integral I want to go through some definitions.

Define triangles $T_1$ and $T_2$ by the points ${a,b,c}$ and ${d,e,f}$ respectively.

Define $$D(u,v,x,y) = sign(det (y-u,y-v,y-x)),$$

and define
$$n_cap (T_1, x, y) = frac{1}{8} (D(a, b, x, y) + D(b, c, x, y) + D(c, a, x, y)
+ \D(a, b, x, y)D(b, c, x, y)D(c, a, x, y))
(1 – D(a, b, c, d)D(a, b, c, e)).$$

Say $lk(T_1,T_2)^2$ is defined as

$$lk(T_1, T_2)^2 = [n_cap (T_1, d, e) + n_cap(T_1, e, f) + n_cap (T_1, f, d)]^2.$$

How can I reduce the error of numerical integration of the following integral?
$$int_{Omega^3} int_{Omega^3} lk(T_1, T_2)^2,$$

where $Omega = [0,1]^3$ is the space in which we take the points $a,b,c,d,e,f$ to generate the triangles $T_1$ and $T_2$. We integrate over the points ${a, b, c}$ and ${d, e, f}$ (hence the $Omega^3$).
Note that the points $a,b,c,d,e,f$ are points in the unit box.

I have tried to calculate the value of this integral numerically using Mathematica and the
result is approximately 0.15… the error estimate is around 0.0016042 and I haven’t been able
to reduce this error.

I am not used to work with Mathematica but I have tried to use Global and Local adaptive strategies and haven’t been succesful. I have
also tried setting an Acuracy Goal, Max and Min recursions but this has not worked at all.
I even tried to change the integration method but this has not worked either.

Any advice of how to reduce the error of the numerical integration or how to calculate the integral symbolically would be appreciated.

The code for the integral in Mathematica is

a = {a1, a2, a3};
b = {b1, b2, b3};
c = {c1, c2, c3};
d = {d1, d2, d3};
e = {e1, e2, e3};
f = {f1, f2, f3};
x = {x1, x2, x3};
y = {y1, y2, y3};

lk2 := ((1/
       8 (Sign[Det[{e - a, e - b, e - d}]] + 
        Sign[Det[{e - b, e - c, e - d}]] + 
        Sign[Det[{e - c, e - a, e - d}]] + (Sign[
           Det[{e - a, e - b, e - d}]]*
          Sign[Det[{e - b, e - c, e - d}]]*
          Sign[Det[{e - c, e - a, e - d}]])) (1 - (Sign[
           Det[{d - a, d - b, d - c}]]*
          Sign[Det[{e - a, e - b, e - c}]]))) + (1/
       8 (Sign[Det[{f - a, f - b, f - e}]] + 
        Sign[Det[{f - b, f - c, f - e}]] + 
        Sign[Det[{f - c, f - a, f - e}]] + (Sign[
           Det[{f - a, f - b, f - e}]]*
          Sign[Det[{f - b, f - c, f - e}]]*
          Sign[Det[{f - c, f - a, f - e}]])) (1 - (Sign[
           Det[{e - a, e - b, e - c}]]*
          Sign[Det[{f - a, f - b, f - c}]]))) + (1/
       8 (Sign[Det[{d - a, d - b, d - f}]] + 
        Sign[Det[{d - b, d - c, d - f}]] + 
        Sign[Det[{d - c, d - a, d - f}]] + (Sign[
           Det[{d - a, d - b, d - f}]]*
          Sign[Det[{d - b, d - c, d - f}]]*
          Sign[Det[{d - c, d - a, d - f}]])) (1 - (Sign[
           Det[{f - a, f - b, f - c}]]*
          Sign[Det[{d - a, d - b, d - c}]]))))^2

NIntegrate[lk2, {a1, 0, 1}, {a2, 0, 1}, {a3, 0, 1}, {b1, 0, 1}, {b2, 
  0, 1}, {b3, 0, 1}, {c1, 0, 1}, {c2, 0, 1}, {c3, 0, 1}, {d1, 0, 
  1}, {d2, 0, 1}, {d3, 0, 1}, {e1, 0, 1}, {e2, 0, 1}, {e3, 0, 1}, {f1,
   0, 1}, {f2, 0, 1}, {f3, 0, 1}]
```

One Answer

On such a high dimensional integral, the default rule is the Monte Carlo rule. You can increase the number of points. I also increased the PrecisionGoal, so that the error estimate will be reported.

NIntegrate[..., 
 Method -> {"MonteCarloRule", "Points" -> 10^6}, PrecisionGoal -> 6]

NIntegrate::maxp: The integral failed to converge after 3000000 integrand evaluations. NIntegrate obtained 0.15226900000000002 and 0.0002590579838772007 for the integral and error estimates.

(*  0.152269  *)

The Monte Carlo Rule error is proportional to $1/sqrt{N}$ where $N$ is the number of "Points" (under certain assumptions). It converges slowly. To get another digit of precision, use about 100 times as many points and wait about 100 times as along.

Answered by Michael E2 on December 31, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP