maisonobe
commented
3 years ago There are probably two different problems here. One is discussed in the forum (see Incomplete elliptic integral of the third kind Π(n, Φ, m) sometimes selects the wrong sheet of the Riemann surface in the complex plane . The forum discussion is about an integral of the third kind whereas EllipticF is an integral of the first kind, but the problem is similar: when using Carlson transforms, at some points the algorithm changes the sheet of the Riemann surface it selects, and jumps instantly from one value to the other.
The second problem is argument reduction which seems wrong.
Here are several plots I made using either the small program in the hipparchus-sample module or Wolfram cloud.
The first plot was computed using numerical integration, using a Gauss-Legendre integrator with 24 points and convergence settings set to 1.0e-6 (both absolute and relative tolerance), and letting the integrator call the integrand function up to 100000 times before giving up. When it gives up, it corresponds to a white point in the diagram:
As you see, numerical integration is not an ideal solution, here the computation is very unstable, perhaps due to m being real and very close to π/2. This computation took more than one hour to complete.
The second plot was computed using numerical integration too, but I relaxed tolerances to a much cruder value of 1.0e-4. Then the computation completed for more points, but the accuracy is probably not very good.
I don't know how much time it took to complete but it was certainly faster than the previous one because the integrator converged for more point, so it did not waste time trying again and again until it reaches 100000 calls to the integrand. We see however that the integral should at least be almost smooth left to right (i.e. along the real direction), despite there are still some sharp corners at ±π/2 in the real part, that looks more like computation artifacts than reality to me. There is only one location of realistic sharp hue transitions which is exactly along the real axis (i.e. there is a discontinuity when crossing the real axis from negative imaginary values to positive imaginary values). Our aim would be to have this kind of result, without the white dots, and perhaps without the sharp corners that appear here and there, and fast.
The third plot was computed using the current implementation of LegendreEllipticIntegral.bigE for complex numbers.
It computes vary fast (a few seconds), and exhibits obvious continuity problems when the real part of the numbers are ±kπ/2. The center band is good, though.
The last plot was in fact my reference, it was computed using Wolfram cloud, you see the command in Wolfram language used just on top of the image.
As you see, Wolfram as the same problem as us! At first, I would have thought their result is correct and the integral has discontinuity, but there are two arguments that make me think this is a problem that affect both our and their implementations. The first argument is that for an integral to be discontinuous, the integrand has to cross a singularity, and the integrand in this case has only isolated singularities (at arcsin(1/sqrt(m) and their repetitions due to arcsin), so we could have a singular point and a ray extending radially from that point to infinity as the integral starts from origin 0 + 0i (see the first plot in the forum discussion for an example in the case of an integral of the third kind). The second argument is that despite it is highly numerically unstable and is inaccurate as depicted by the first two plots, numerical integration shows the integral is almost smooth when crossing the ±kπ/2 lines (except for a few sharp corners).
Hipparchus implementation and Wolfram implementation therefore seem to have the same kind of problems related to the argument reduction used to bring the computation back inside the [-π/2;+π/2] interval (for the real part of the argument) where everything is fine.
Also note that the test point in your issue report is at a very specific and unstable point. It is both:
- just at the boundary of the [-π/2;+π/2] interval
- exactly opposite to the ellipse parameter m
- just in the only discontinuity we think is correct: the one corresponding to the real axis itself
So I think it is really a difficult point, and due to the discontinuity along real axis, we could expect two different values. We should in theory need to set up a branch cut differentiating +0 and -0 imaginary part, but I think we already have too many problems.
As explained in the forum discussion, I don't think I will be able to solve this fast (Wolfram did not either). I still want to try, but think it should wait for next version.
What do you think?
axkr
Please Verify
EllipticF(-1.5708,1.5708)See these discussions: