# Surprising error message from NIntegrate

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ϕ][];
θhSol[θS_, a_] =
θh /. Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] == 1, θh] /. {C -> 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]][], θhSol[θS, A[ES, mS]][]], 0],
Min[Max[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 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?

• The limits of integration depend also on \[Theta]h]! – Ulrich Neumann May 28 at 10:14
• @UlrichNeumann That kind of dependency is not a problem. – Anton Antonov May 28 at 10:31
• As per my posted answer/comment voting to close as a "simple mistake". But I do not feel strongly about it, ~65%. – Anton Antonov May 28 at 10:43
• @AntonAntonov Thanks, I got it! – Ulrich Neumann May 28 at 10:47
• @AntonAntonov The doc page for NIntegrate states NIntegrate[f[x, y], {x, a, b}, {y, c, d}] $\int_a^b \int _c^df(x,y)\,dy\,dx$, and in normal math notation, $c$ and $d$ may depend on $x$. Aside from that, I don't know of a dedicated discussion. It comes up in some of the singularities handling and preprocessors, for instance "UnitCubeRescaling". – Michael E2 May 28 at 13:19

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ϕ][];

θhSol[θS_,
a_] := θh /.
Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] ==
1, θh] /. {C -> 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]][], θhSol[θS, A[ES, mS]][]], 0],
Min[Max[θhSol[θS,
A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 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:= 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:= 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:= 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:= General::stop: Further output of Solve::ifun will be suppressed during this calculation.

During evaluation of In:= 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]][]], 0],
Min[Max[θhSol[0.01, A[1500, 40]][[
2]], θhSol[0.01, A[1500, 40]][]], π]}, {ϕ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]][], θhSol[θS, A[ES, mS]][]], 0],
Min[Max[θhSol[θS,
A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 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:= 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:= 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:= 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:= General::stop: Further output of Solve::ifun will be suppressed during this calculation. *)

(* 0.0981353 *)

• Thank you! Your suggestion about Re solved a problem. – John Taylor May 28 at 10:41
• @JohnTaylor Good to hear. Good luck! – Anton Antonov May 28 at 10:44

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]][], θhSol[θS, A[ES, mS]][]], 0],
Min[Max[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 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[], θh0[]},
{ϕ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  *)


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[], θh0[]},
MinRecursion -> 5, MaxRecursion -> 20,
Method -> "GaussKronrodRule"]
];

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

• In numerical problems, the branch points of arc cosine and arc sine are trouble spots. Spurious complex results are just one symptom of the numerical instability there. It is usually better to compute both the sine and cosine of the unknown angle, and then use the two argument form of ArcTan to get the angle. – John Doty May 29 at 12:55
• @JohnDoty I think I understand what you're saying, but I don't see how it applies here. The input to ArcCos[]` goes well outside the domain $[-1,1]$, and by my reading (i.e. guess), the intention is to exclude regions outside the domain from integration. – Michael E2 May 29 at 19:45