Mathematica Asked by Pedro Morales on June 22, 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}]
```
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
obtained0.15226900000000002
and0.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.
Correct answer by Michael E2 on June 22, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP