0
$\begingroup$

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][[2]]
dom1[gh_, ts_, ms_] := 
 ImplicitRegion[
  cosa[ts, th, ph] - cosacrittemp[gh, ms][[2]] > 
   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)]
$\endgroup$
  • $\begingroup$ 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]}}[[2]]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. $\endgroup$ – Alex Trounev May 11 at 12:41
  • $\begingroup$ @AlexTrounev : I apologize. I have already corrected the code. $\endgroup$ – John Taylor May 11 at 17:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.