Convergence of numerical integration with parameter in the region of integration

Question

I would like to solve the integral:

$ f(a)=\iiint\limits_{\mathcal{R}(a)}\frac{\sqrt{12\left(\left(\frac{x-y}{x+y}\right)^2+\left(\frac{y-z}{y+z}\right)^2+\left(\frac{x-z}{x+z}\right)^2\right)}}{\exp(x+y+z)\sqrt{\left(\frac{4y(x-y)}{(x+y)^3}+\frac{4z(x-z)}{(x+z)^3}\right)^2+\left(\frac{-4x(x-y)}{(x+y)^3}+\frac{4z(y-z)}{(y+z)^3}\right)^2+\left(\frac{-4y(y-z)}{(y+z)^3}-\frac{4x(x-z)}{(x+z)^3}\right)^2}}\rm{d} x \rm{d} y \rm{d}z $

The region of integration $\mathcal{R}(a)$ is the surface:

$ \mathcal{R}(a)=\left\{(x,y,z)\in \mathbb{R}^{3}, x\geq 0, y\geq 0, z\geq 0,z\leq x,z\leq y, \sqrt{\frac{1}{3}\left(\left(\frac{x-y}{x+y}\right)^2+\left(\frac{y-z}{y+z}\right)^2+\left(\frac{x-z}{x+z}\right)^2\right)}=a\right\} $

With parameter $a\in[0,1]$.

All my efforts to solve this analytically have failed. I'm now struggling to get an accurate numerical solution.

I implemented it in Mathematica (v. 10.3) in this manner:

region[a_?NumericQ] =  ImplicitRegion[Sqrt[1/3 ((x - y)^2/(x + y)^2 + (x - z)^2/(x + z)^2 + (y - z)^2/(y + z)^2)] == a && x > 0 && y > 0 &&  z > 0 && z <= x && z <= y, {x, y, z}]

result = Table[{a,Re[NIntegrate[Sqrt[12 ((x - y)^2/(x + y)^2 + (x - z)^2/(x + z)^2 + (y - z)^2/(y + z)^2)]/(E^(x + y + z) Sqrt[(-((4 x (x - z))/(x + z)^3) - (4 y (y - z))/(y + z)^3)^2 + ((4 z (y - z))/(y + z)^3 - (4 x (x - y))/(x + y)^3)^2 + ((4 y (x - y))/(x + y)^3 + (4 z (x - z))/(x + z)^3)^2]), {x, y, z} [\Element] region[a]]]}, {a, 0, 1, .01}]

But I'm having all sort of convergence problems, 1/0 errors, accuracy, instabilities, etc you name it...

Do you know what I might be doing wrong? Do you know any work-around? Any advice is welcome.

Answer 1

I have rewritten your definition of the region as

region2[a_?NumericQ]:=ImplicitRegion[2 x^3 (y-z)^2 (y+z)+2 x y (y-z)^2 z (y+z)+x^4 (3 y^2+2 y z+3 z^2)+y^2 z^2 (3 y^2+2 y z+3 z^2)+x^2 (3 y^4-2 y^3 z-18 y^2 z^2-2 y z^3+3 z^4)==3 a^2 (x+y)^2 (x+z)^2 (y+z)^2&&x>0&&y>0&&z>0&&z<=x&&z<=y,{x,y,z}]

That makes at least ContourPlot3D for these regions considerably faster.

Then Mathematica 10.3 can do the integration. Of course, it is time consuming.

NIntegrate[ Sqrt[12 ((x-y)^2/(x+y)^2+(x-z)^2/(x+z)^2+(y-z)^2/(y+z)^2)]/(E^(x+y+z) Sqrt[(-((4 x (x-z))/(x+z)^3)-(4 y (y-z))/(y+z)^3)^2+((4 z (y-z))/(y+z)^3-(4 x (x-y))/(x+y)^3)^2+((4 y (x-y))/(x+y)^3+(4 z (x-z))/(x+z)^3)^2]), {x,y,z} \[Element] region2[0.3]] // Timing

(* {166.765, 0.321009 -1.75152*10^-8 I} *)