# How to handle the instability of the integral evaluation in the following case?

I want to integrate a function \[Kappa][gh, ms, ts, th, ph] (for all definitions, see a code below) with the distribution distr[gh,th] over variables gh,th,ph in a region $$ph\in (0,2\pi), \quad th \in (0,\pi) \quad \text{if} \quad gh \in (1,\text{ghcrit}[ms]),$$ $$\{ph, th\} \in \text{dom1} \quad \text{if} \quad gh \in (\text{ghcrit}[ms],25)$$ The integral, Average[ms,ts], is defined in a domain $$ms \in (0,62.5), \quad ts \in (0,\pi)$$ I compute it using MonteCarlo method.

For small values of ms (say, ms = 19, ts = Pi/8) the integral evaluates "stable" (i.e. always returns a result), while for large ms (say, ms = 60 and ts = Pi/8 as an example) it behaves "unstable": sometimes it gives the result, while sometimes it stuck (i.e. evaluates infinitely long). I do not know how to handle this unstable behavior. My goal is to compute

Table[{ts, Average[60, ts]},{ts,0,Pi,0.01}]


and the stability of the evaluation is crucial for me.

Could you please tell me how to fix this integral?

The code:

sol = Simplify[Solve[Sqrt[1 - x^2]/x == Sqrt[1 - y^2]/(A (y + B)), y]]
y1 = (A^2 B (1 - x^2) -
x Sqrt[A^2 (B^2 - 1) (x^2 - 1) + x^2])/(A^2 (x^2 - 1) - x^2);
y2 = (A^2 B (1 - x^2) +
x Sqrt[A^2 (B^2 - 1) (x^2 - 1) + x^2])/(A^2 (x^2 - 1) - x^2);
b[gh_, ms_] = Sqrt[gh^2 - 1]/(gh*Sqrt[1 - 4*ms^2/125^2]);
cosacrittemp[gh_, ms_] =
x /. Solve[A^2*(B^2 - 1)*(x^2 - 1) + x^2 == 0, x] /. {A -> gh,
B -> b[gh, ms]}
cosa[ts_, th_, ph_] = Cos[ph]*Sin[th]*Sin[ts] + Cos[th]*Cos[ts];
\[Kappa][gh_, ms_, ts_, th_, ph_] =
Evaluate[Simplify[D[y1, x] /. {x -> Cos\[Alpha]}]]*
UnitStep[ghcrit[ms] - gh] +
UnitStep[gh - ghcrit[ms]]*
Evaluate[
Simplify[
D[y1 - y2, x] /. {x -> Cos\[Alpha]}]] /. {Cos\[Alpha] ->
cosa[ts, th, ph]} /. {A -> gh, B -> b[gh, ms]};
ghcrit[ms_] = gh /. Solve[b[gh, ms] == 1, gh][]
dom1[gh_, ts_, ms_] :=
ImplicitRegion[
cosa[ts, th, ph] - cosacrittemp[gh, ms][] >
0, {{th, 0, Pi}, {ph, 0, 2*Pi}}]
distr[gh_, th_] = Exp[-Sqrt[gh]*(1 + Cos[th]^4)]*Cos[th]^20;
Average[ms_, ts_] :=
If[ghcrit[ms] > 25,
NIntegrate[
distr[gh, th]*\[Kappa][gh, ms, 125, cosa[ts, th, ph]], {gh, 1,
25}, {th, 0, Pi}, {ph, 0, 2*Pi},
Method -> "MonteCarlo"]/(2*Pi*Pi),
NIntegrate[
distr[gh, th]*\[Kappa][gh, ms, ts, th, ph], {gh, 1,
ghcrit[ms]}, {th, 0, Pi}, {ph, 0, 2*Pi},
Method -> "MonteCarlo"]/(2*Pi*Pi) +
2*NIntegrate[
distr[gh, th]*\[Kappa][gh, ms, ts, th, ph], {gh, ghcrit[ms],
25}, {th, ph} \[Element] dom1[gh, ts, ms],
Method -> "MonteCarlo"]/(2*Pi*Pi)]

• I have three messages Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. Part::partw: Part 2 of {{gh->InverseFunction[Bb,1,2][1,ms]}} does not exist. ReplaceAll::reps: {{{gh->InverseFunction[Bb,1,2][1,ms]}}[]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. – Alex Trounev May 11 '19 at 12:41
• @AlexTrounev : I apologize. I have already corrected the code. – John Taylor May 11 '19 at 17:12