Introduction and summary
The main point is that feeding inexact input to an exact solver may produce incorrect results. This should be more generally known. It is discussed in the last paragraph of this answer, Understand the difference between exact and approximate (Real) numbers.
A secondary point is that there is a bug in Limit[]
(imo). The problem is extremely localized, but maybe it is worth exposing. The exact reasons for it and the extent of its generality are unclear. It can be avoided by using exact or arbitrary-precision input, so the main point above undermines its importance as a bug. In connection with the main point, it is worth pointing out that the Solve
/Reduce
suite of solvers take a different approach. They rationalize the input and numericize the output, emitting a warning message along the way. The same approach would work on this example. It may be worth considering whether Limit
, Integrate
and others might benefit from such an approach.
My original guess was that rounding error led to an unpredictable choice, probably choosing the wrong branch. I'm pretty sure now that it is not a wrong branch being chosen. There is a discontinuity at Pi/3
, but the result does not correspond to choosing the wrong branch, AFAICT. Instead, I believe I have localized the problem to Limit
failing to find the limit at a single machine floating-point number near Pi/3
. This may or may not be due to a rounding error; for instance, if arbitrary-precision reals are used, Limit
succeeds, and one of the principal effects of arbitrary-precision is to reduce rounding error, especially with respect to quantities canceling or testing equal. If it is rounding error, it's a hidden, internal error, which might be dismissed by referring to the main point.
On the other hand, since there is not such a problem at other floating point numbers near Pi/3
, including Pi/3.
itself, the user point of view should be to label this a bug and report it to WRI, the main point notwithstanding.
Analyzing the problem
A trace will reveal several limits, some of which test convergence, some which are used to apply the Fundamental Theorem of Calculus, and some which seem to test for branch cuts or other discontinuities.
Trace[
Integrate[
Sqrt[1/(2*98/10*(Sin[Pi/3] - Sin[θ]))], {θ,
ArcSin[2/3 Sin[Pi/3]], Pi/3.}],
_Limit,
TraceInternal -> True
] // Flatten
Below we show the last four and rewrite them in terms of the antiderivative used by Integrate
internally (which has been stripped of a constant factor to be an antiderivative for the given integral), which I call ff[]
. The first two limits evaluate ff[]
at (or near) the end points of the interval of integration. The last two seem to be checking the antiderivative at the limit of integration for a branch discontinuity.
ff[θ_] := -2 EllipticF[1/4 (π - 2 θ), 4 (2 + Sqrt[3])] *
Sqrt[1/(Sqrt[3] - 2 Sin[θ])] *
Sqrt[(Sqrt[3] - 2 Sin[θ])/(-2 + Sqrt[3])];
{Limit[ff[θ], θ -> 1.0471975511965974, Analytic -> False,
Assumptions -> True, Direction -> 1, Method -> "InternalClassic"],
Limit[ff[θ], θ -> 0.6154797086703874, Analytic -> False,
Assumptions -> True, Direction -> -1, Method -> "InternalClassic"],
Limit[ff[θ], θ -> 1.0471975511965976, Analytic -> False,
Assumptions -> True, Direction -> -1, Method -> "InternalClassic"],
Limit[ff[θ], θ -> 1.0471975511965976, Analytic -> False,
Assumptions -> True, Direction -> 1, Method -> "InternalClassic"]}
(*
{0.,
-1.24682 - 1.59814 I,
-2.39469*10^-8 - 1.59814 I,
-2.39469*10^-8 - 1.59814 I}
*)
If we look at the limits of ff[]
at consecutive floating-point numbers near the limit of integration, we see there is some misbehavior at θ -> 1.0471975511965974
, which is equal to Pi/3. - $MachineEpsilon
. This incorrect limit is what leads to the incorrect integral. There is also trouble with the limit at the number just below it also in one of the one-sided limits, but otherwise there is no problem with limits at other nearby numbers.
Limit[ff[θ], θ -> 1.0471975511965974 - 2 $MachineEpsilon]
{Limit[ff[θ], θ -> 1.0471975511965974 - $MachineEpsilon, Direction -> 1],
Limit[ff[θ], θ -> 1.0471975511965974 - $MachineEpsilon, Direction -> -1]} ->
Limit[ff[θ], θ -> 1.0471975511965974 - $MachineEpsilon]
Limit[ff[θ], θ -> 1.0471975511965974]
Limit[ff[θ], θ -> 1.0471975511965974 + $MachineEpsilon]
Limit[ff[θ], θ -> 1.0471975511965974 + 2 $MachineEpsilon]
(*
-5.86663*10^-8 - 1.59814 I
{0., -4.65708*10^-8 - 1.59814 I} -> Indeterminate
0.
-2.39469*10^-8 - 1.59814 I
0. - 1.59814 I
*)
One might wonder if there is an evaluation problem near Pi/3
. There does not seem to be one:
ff[θ] /. θ -> 1.0471975511965974 - 2 $MachineEpsilon
ff[θ] /. θ -> 1.0471975511965974 - $MachineEpsilon
ff[θ] /. θ -> 1.0471975511965974
ff[θ] /. θ -> 1.0471975511965974 + $MachineEpsilon (* == Pi/3. *)
ff[θ] /. θ -> 1.0471975511965974 + 2 $MachineEpsilon
(*
-5.56222*10^-8 - 1.59814 I
-4.62805*10^-8 - 1.59814 I
-3.44954*10^-8 - 1.59814 I
Indeterminate <-- θ -> Pi/3.
-5.94877*10^-16 - 1.59814 I
*)
The Indeterminate
result comes from Sqrt[3]-2 Sin[θ]
evaluating to zero at Pi/3.
.
An irrelevant point: The round-off error in Sqrt[3]-2 Sin[θ]
for θ
around Pi/3.
is nearly catastrophic in itself, and arbitrary-precision numbers lose about 16 digits of precision. However, since we multiply and divide by the same thing, (the square-root of) this quantity, the rounding error cancels out, and the results are much more accurate than they appear.
ff[θ] /. θ -> 1.0471975511965974`24 - 2*2^-52
ff[θ] /. θ -> 1.0471975511965974`24 - 2^-52
ff[θ] /. θ -> 1.0471975511965974`24
ff[θ] /. θ -> 1.0471975511965974`24 + 2^-52
ff[θ] /. θ -> 1.0471975511965974`24 + 2*2^-52
(*
-5.622254*10^-8 - 1.59814200 I
-4.767384*10^-8 - 1.59814200 I
-3.72104*10^-8 - 1.59814200 I
-2.22809*10^-8 - 1.5981420 I
0.*10^-14 - 1.5981420 I
*)
Discontinuity
The discontinuity alluded to above occurs in the complex plane. It's conceivable that the problem relates to it, but it is unclear how exactly. Here is a plot of the real and imaginary parts of ff[]
along a small path encircling the point where Limit
fails:
Plot[
ReIm[ff[θ]] /. θ ->
Rationalize[1.0471975511965974, 0] + 0.001 Exp[I t] // Evaluate,
{t, 0, 2 Pi}]
Workarounds
For the limit:
N@Limit[ff[θ], θ -> Rationalize[1.0471975511965974`, 0]]
(* -3.44954*10^-8 - 1.59814 I *)
N@Limit[ff[θ], θ -> 1.0471975511965974`16]
(* -3.54651*10^-8 - 1.59814 I *)
Limit[Simplify[ff[θ], Sqrt[3] - 2 Sin[θ] > 0],
θ -> 1.0471975511965974]
(* -3.44954*10^-8 - 1.59814 I *)
For the integral:
N@Integrate[Sqrt[1/(2*98/10*(Sin[Pi/3] - Sin[θ]))],
{θ, ArcSin[2/3 Sin[Pi/3]], Pi/3}]
(* 0.398283 - 1.41859*10^-16 I *)
N@Integrate[Sqrt[1/(2*98/10*(Sin[Pi/3] - Sin[θ]))],
{θ, ArcSin[2/3 Sin[Pi/3]], Pi/3.`16}]
(* 0.398283 + 0. I *)
Integrate[Sqrt[1/(2*98/10*(Sin[Pi/3] - Sin[θ]))],
{θ, ArcSin[2/3 Sin[Pi/3]], Pi/3.}, Assumptions -> Sin[Pi/3] - Sin[θ] > 0]
(* 0.398283 - 2.83718*10^-16 I *)
Integrate[Sqrt[1/(2*98/10*(Sin[Pi/3] - Sin[θ]))],
{θ, ArcSin[2/3 Sin[Pi/3]], 1., Pi/3.}]
(* 0.398283 + 0. I *)
Integrate[integrand, {\[Theta], ArcSin[2/3 Sin[Pi/3.]], Pi/3}]
gives you0.398283 + 0.510508 I
, which is a significant difference from the other answers, but wouldn't be expected to differ. The upper bound of the integral is singular, and causes problems. $\endgroup${θ, ArcSin[2./3 Sin[Pi/3]], Pi/3}
with an exact upper bound but numerical lower, that I get a significantly different answer from exact upper and lower bounds of{θ, ArcSin[2/3 Sin[Pi/3]], Pi/3}
, and the imaginary part is not round-off error. It's just the wrong answer. $\endgroup$