I define the following function $$ f(x,y,z) = \frac{\sin x}{\pi \sqrt{\Big[ \cos x - \cos (y-z) \Big] \Big[ \cos (y+z) - \cos (x) \Big]}} $$ where $y, z$ are all between $0, \pi/2$ and $ |y-z| < x < \text{min}(y+z, \pi/2)$. Additionally, if $y+z > \pi/2$ we reflect the function in the $x$ argument around from $\pi$, $x \to \pi - x$. Thus: $$ g(x,y,z) = \begin{cases} f (x,y,z); \; &y + z < \pi/2 \\ f (x,y,z) + f (\pi-x,y,z); \; &y + z > \pi/2 \end{cases} $$
The problem with $g$ is that it has some singularities (all integrable, though).
In Mathematica, I define the function like this
test[x_, y_, z_] := Sin[x] ((HeavisideTheta[x - Abs[y - z]] HeavisideTheta[y + z - x])/(Pi Sqrt[(Cos[x] - Cos[y - z]) (Cos[y + z] - Cos[x])]) + (HeavisideTheta[y + z - Pi/2] HeavisideTheta[x + y + z - Pi])/(Pi Sqrt[-(Cos[x] + Cos[y - z]) (Cos[y + z] + Cos[x])]));
All the thetas are to ensure we don't get any complex numbers from the function. My goal is to do a parametric integral of test(x,y,z)
multiplied by some simple function of y
and z
(no singularities in this additional function) over a square {y, 0, Pi/2}
, {z, 0, Pi/2}
with x
as a parameter.
With[{x = Pi/50}, NIntegrate[ Sin[y] Sin[z] test[x, y, z], {y, 0, Pi/2}, {z, 0, Pi/2}]] // Timing
With[{x = Pi/16}, NIntegrate[ Sin[y] Sin[z] test[x, y, z], {y, 0, Pi/2}, {z, 0, Pi/2}]] // Timing
With[{x = Pi/4}, NIntegrate[ Sin[y] Sin[z] test[x, y, z], {y, 0, Pi/2}, {z, 0, Pi/2}]] // Timing
With[{x = 7 Pi/16}, NIntegrate[ Sin[y] Sin[z] test[x, y, z], {y, 0, Pi/2}, {z, 0, Pi/2}]] // Timing
Here I took Sin[y] Sin[z]
as the additional function, as an example.
The timings are absolutely atrocious,
{8.26563, 0.0617321}
{7.25, 0.193324}
{5.875, 0.700855}
{7.35938, 0.972718}
And I also keep getting the same two errors: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
and 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 0.06173208991557969 and 0.0006314495223153685 for the integral and error estimates.
I don't know why I'm getting these, the integral isn't 0, the integrand isn't highly oscillatory and even when I increase the WorkingPrecision
(tried 16, 50, 100) I still get the same error (but now it takes much more time).
If I try to limit the integration domain to points that already satisfy $ |y-z| < x < \text{min}(y+z, \pi/2)$
reg[x_] := ImplicitRegion[(x > Abs[y - z]) && (x < y + z) && (y > 0) && (y < \[Pi]/2) && (z > 0) && (z < \[Pi]/2), {y, z}];
NIntegrate[Sin[y] Sin[z] test[x, y, z], {y, z} \[Element] reg[x]]
Now the first integral won't even evaluate (the integrand has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0,1},{0.,5.84851*10^-15}}
) and the other three integrals give slightly different results, and, despite limiting the region to even smaller area, the timings are worse
{8.54688, 0.19472}
{6.92188, 0.704321}
{8.625, 0.97823}
The error on the first integral is especially weird, as MMA seems to have make up its own weird domain (the region I specified is not $[0,1] \times [0, 5.84851 \times 10^{-15}]$, where is this number even coming from?)
There are so many weird things happening here, can someone shed some light on it?
u == y + z
,v == y - z
. $\endgroup$u+v = 2y>0
,u+v < Pi / 4
,u-v = 2z > 0
,u - v < Pi/4
. This seems to be a tilted in the u-v world. Should I do it via integral {u, 0, pi/4}, {v, Abs[Pi/8 - u] - Pi/8, Pi/8 - Abs[Pi/8 - u]} (not sure if MMA supports this kind of functional bound), or via region function? Should I also remove the theta functions from the definition, since now it's easier to catch? $\endgroup$HeavisideTheta
stands for the certain distribution,UnitStep
is the right command here. $\endgroup$