Surprising error message from NIntegrate

Question

I have the following integral:

cosϕSol[θS_, sinθh_, cosθh_, a_] =
  cosϕ /. Solve[cosϕ*Sin[θS] sinθh + Cos[θS]*cosθh == a, cosϕ][[1]];
θhSol[θS_, a_] =
  θh /. Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] == 1, θh] /. {C[1] -> 0};
A[ES_, mS_] = Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];
integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh,
   Max[Min[θhSol[θS, A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0], 
   Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0, 
   Min[Pi, ArcCos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]]}]

When evaluating it at some point, say, integral1[1500, 0.01, 40], I got an error:

NIntegrate::inumr: "The integrand 1 has evaluated to non-numerical values for all \ sampling points in the region with boundaries {{0,0.25},{0,1}}"

What can be a reason for this?

Answer 1

Not an answer, extended comment...

First some redefinitions:

Clear[cosϕSol, θhSol, A, integral1]

cosϕSol[θS_, sinθh_, cosθh_, a_] := 
  cosϕ /. 
   Solve[cosϕ*Sin[θS] sinθh + 
       Cos[θS]*cosθh == a, cosϕ][[1]];

θhSol[θS_, 
  a_] := θh /. 
   Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] ==
      1, θh] /. {C[1] -> 0}

A[ES_, mS_] := Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];

Because of the complicated symbolic computations for the integration ranges, I first switched off the symbolic preprocessing. I got the following output:

integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh, 
   Max[Min[θhSol[θS, 
       A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0], 
   Min[Max[θhSol[θS, 
       A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0, 
   Min[Pi, ArcCos[
     cosϕSol[θS, Sin[θh], Cos[θh], 
      A[ES, mS]]]]}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]

integral1[1500, 0.01, 40]

(* During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[132]:= General::stop: Further output of Solve::ifun will be suppressed during this calculation.

During evaluation of In[132]:= NIntegrate::nlim: ϕh = Min[3.14159,3.14159 -0.545349 I] is not a valid limit of integration. *)

(* NIntegrate[1, {θh, 
  Max[Min[θhSol[0.01, A[1500, 40]][[
     1]], θhSol[0.01, A[1500, 40]][[2]]], 0], 
  Min[Max[θhSol[0.01, A[1500, 40]][[
     2]], θhSol[0.01, A[1500, 40]][[1]]], π]}, {ϕh, 0,
   Min[π, 
   ArcCos[cosϕSol[0.01, Sin[θh], Cos[θh], 
     A[1500, 40]]]]}, 
 Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}] *)

Using Re in the second integration range produces a numerical result:

integral1[ES_, θS_, mS_] :=
 NIntegrate[1,
  {θh, 
   Max[Min[θhSol[θS, 
       A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0], 
   Min[Max[θhSol[θS, 
       A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
  {ϕh, 0, 
   Min[Pi, Re@
     ArcCos[cosϕSol[θS, Sin[θh], Cos[θh], 
       A[ES, mS]]]]}, 
  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]

integral1[1500, 0.01, 40]

(* During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.

During evaluation of In[134]:= General::stop: Further output of Solve::ifun will be suppressed during this calculation. *)

(* 0.0981353 *)

Answer 2

Another extended comment:

I wonder if the following might be another "fix." The question I have is that when code depends on an erroneous comparison of real and complex numbers, all means that make the error go away, be it by applying Re[], Abs[] or whatever, cannot all be correct. Since the illegal comparison arises because of inputs to ArcCos[] that lie outside its domain as a real-valued function, it seems better to me to fix the problem by restricting the input to the domain and defining the integral to be zero outside this domain. This can be done in this case by making the upper ArcCos[] limit of ϕh be equal to the lower limit 0. There is, of course, not enough information in the question to decide which of the possible fixes is correct. That is for the OP to figure out.

integral1[ES_, θS_, mS_] :=
  NIntegrate[1,
   {θh,
    Max[Min[θhSol[θS, A[ES, mS]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
    Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], Pi]},
   {ϕh, 0, 
    Min[Pi, ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
       cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]
     ]}];

integral1[1500, 0.01, 40]

NIntegrate::slwcon: Numerical integration converging too slowly....

(*  0.0289182  *)

The code can be simplified, too, and the NIntegrate::slwcon warning disappears by including MinRecursion -> 2:

integral1[ES_, θS_, mS_] :=
  With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
   NIntegrate[1,
    {θh, θh0[[1]], θh0[[2]]},
    {ϕh,0,
     ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
       cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]
      ]},
    MinRecursion -> 2]
   ];

integral1[1500, 0.01, 40]
(*  0.0289182  *)

Addendum

There's another issue that I neglected to consider. The integral is equivalent to an area-under-the-curve single integral:

ClearAll[realACos];  (* another way to code a real arc cosine *)
realACos[x_?NumericQ /; -1 <= x <= 1] := ArcCos[x];
realACos[_?NumericQ] := 0;

integral2[ES_, θS_, mS_] :=
  With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
   NIntegrate[
    realACos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]],
    {θh, θh0[[1]], θh0[[2]]},
    MinRecursion -> 5, MaxRecursion -> 20,
    Method -> "GaussKronrodRule"]
   ];

integral2[1500, 0.01, 40]
(*  0.0289182  *)