LegendreEllipticIntegral.bigE(x, m)

[axkr]

Can you add a link to the Wikipedia definition to LegendreEllipticIntegral.bigE(x, m) and describe the relaton of your LegendreEllipticIntegral.bigE(x, m) to the definition in Mathematica?

• https://reference.wolfram.com/language/ref/EllipticE.html

Can you implement both "Konventionen" from this table:

• https://de.wikipedia.org/wiki/Elliptisches_Integral#Definition_der_vollst%C3%A4ndigen_elliptischen_Integrale

[maisonobe]

I can surely add the link. Do you think it is urgent enough to cancel the ongoing vote for releasing 2.0 or can it wait until 2.1?

About the convention, I am reluctant to implement both, as there is really a simple conversion between the two (as explained in the class javadoc, m = k²) and since there is a single parameter (either k or m), we would have to rely on different method names, like bigEk/bigEm and this for all integrals (E, K, K', D, F, Π), and for all signatures (double, field, complex, field complex), and for both complete and incomplete integrals, so this would double the combinatorics which is already huge in this class (there are 36 different integrals in the class currently). I have tried to follow the most widespread convention, i.e. to use elliptic modulus k for the integrals, parameter m for the Jacobi functions and nome q for the Jacobi theta functions.

[axkr]

Can you verify if I m doing something wrong for the example from the MMA page: EllipticE[3 + 2.5*I, 2.3 - 1.5*I]

If I use LegendreEllipticIntegral.bigE(x, m.sqrt()) it should evaluate to 3.05969 + I* 11.1650 but gives -0.912102+I*(-9.13895)

BTW is the name bigE suitable? Should it be capitalE?

• https://en.wiktionary.org/wiki/capital_letter

[axkr]

Regarding the question which convention (kor m) will be better to implement; if you want to implement the m convention with the current implementation you would have a sqrt() call for your second parameter which will be "squared" inside your bigE method. If the methods in hipparchus would implement the m convention only the "square" operation for the second parameter would be necessary.

[maisonobe]

If I use LegendreEllipticIntegral.bigE(x, m.sqrt()) it should evaluate to 3.05969 + I* 11.1650 but gives -0.912102+I*(-9.13895)

I will look at this. There are also different conventions for the first argument in the case of incomplete integrals. In mathematical text, the chosen convention for both arguments is shown by using different separators between the arguments. For example when a comma is used as in E(φ, k), the first argument φ is an angle and the second argument k is the elliptic modulus: this is the trigonometric form of the integral. When a semicolon is used as in E(φ; m) the first argument φ is an angle and the second argument m=k² is the parameter: this is also a trigonometric form of the integral. When a vertical bar is used as in E(x|m) the first argument x=sin(φ) is not an angle anymore and the second argument m=k² is the parameter: this is the Legendre form. When a reverse solidus is used as in E(φ\α) the first argument φ is an angle and the second argument α is the modular angle.

I think Wolfram uses the Legendre form, and we use the trigonometric form. However, even using x.asin() to go back to our convention or changing the convention internally to Legendre form did not give the expected result, so I have to investigate further. I checked the complete integral (i.e. only one parameter), and found the same result as Wolfram Alpha.

I will change the convention from k to m as you suggest, and also change the second parameter to x (as internally we indeed always compute the sine of the first parameter. So in terms of mathematical separators, I will change from trigonometric form E(φ, k) to Legendre form E(x|m).

I canceled the ongoing vote on the forum for 2.0, until this issue is fixed.

Concerning the name bigE instead of capitalE, I think the former is consistent with another classical naming scheme in math, the big-O notation.

[maisonobe]

From the documentation, it seems Mathematica is using the trigonometric form of the integral with parameter m (hence it computes E(φ; m)). For real values of the Jacobi amplitude φ between 0 and π/2, we get similar results for the incomplete integrals as them. For the complete integral (φ=π/2), we also get consistent results.

However, our implementation fails for Jacobi amplitude outside of the [0 ; π/2] interval. I have to look at argument reduction. Here is a plot of out implementation of E(z, 2.3-1.5i). It clearly shows discontinuities at kπ/2. Do you have access to a reference program that could produce a comparable display?

[maisonobe]

The issue seems to be fixed by now. I have pushed a number of changes in the release-2.0 branch. Beware this branch is still in development, the vote for release has not occurred yet. I still have some problems with incomplete elliptic integral of the third kind Π(φ, α², m) when I compare the to Wolfram alpha results. Two tests fail (so I have put an @ignore annotation temporarily to avoid breaking the continuous integration). The tests compare three different computations of the same integral: one from wolfram alpha, one using the new numerical complex integration feature and one using the Carlson-based implementation. Two cases out of four give consistent results. One case give consistent results between Wolfram alpha and numerical integration and inconsistent with Carlson. One case give consistent results between numerical integration and Carlson and inconsistent with Wolfram alpha. So I don't understand these tests for now.

I am totally exhausted and demoralized by this and am seriously considering either shipping the library as is with a warning about Π(φ, α², m) being unreliable or even remove this integral from 2.0.

[axkr]

Does it make sense to push this also to 2.0-SNAPSHOT ? Or what is the policy here?

[maisonobe]

I have also pushed it on master now, reverting the version id to 2.0-SNAPSHOT so the continuous integration takes it into account. It should be available on packages.orekit.org.

[axkr]

Does it make sense to compare with ellipticPi

• https://github.com/axkr/symja_android_library/blob/e4b2bb1822c67f729f127619a619950bcb125478/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/functions/EllipticIntegralsJS.java#L465

which was ported from this Javascript lib:

• https://github.com/paulmasson/math/blob/40fdc802f9891c6b0e0849f9b11cbaf17c9f134e/src/functions/elliptic-integrals.js#L326

[maisonobe]

Yes, but beware to the order an types of arguments.

• Abramowitz and Stegun (A&S) uses Π(n; φ\α)
• Carlson 1995 uses Π(φ, k, n) but the n as the opposite sign as the one in A&S
• DLMF uses Π(φ, α², k) but the α² from DLFM is equal to the n of A&S and the k from DLMF is equal to sin²(α) in A&S (so it is not the same α)
• Wolfram uses π(n; φ, k) in Mathworld but implements EllipticPI(n, φ, m) with m=k² in the Wolfram language
• at first, Hipparchus signature was bigPI(φ, α², k), which was chosen for consistency with DLMF. However, since following your suggestion we already switched from k to m=k2 convention, our current notation became inconsistent. I therefore changed it again to bigPI(n, φ, m) so it now matches Wolfram language (i.e. I renamed our α² to n and switched the first two arguments); this change has just been pushed to master and continuous integration already pushed the compiled snapshot artifacts.

I looked at ellipticPi and it seems to handle the quasi-periodicity differently from what I did in Hipparchus (there is a sign switch every quadrant in my implementation and none in yours). My implementation was made for consistency with F, E and D according to DLMF equation 19.2.10. This equation does not mention Π (only F, E, and D) but it seemed to be just an oversight. On the other hand your implementation seems consistent with Wolfram reference functions (EllipticPI periodicity or (EllipticF periodicity. When I don't put the sign change, I get discontinuities (it was part of solving this issue). I think it is due to one method handling only the [0; π/2] interval and the other one able to handle [-π;π] due to function parity, but am not sure (I am not sure of anything by now).

There are two leads that could explain the problem with incomplete Π in the complex case. Both are compatible with the fact that for some input values the results are really good (to about 15 significant digits) and for some other they are completely out.

The first lead is the proble described in Carlson[2000] paper, which fixes some issues in the algorithm from his previous paper Carlson[1995]. I am looking at this right now.

The second lead is roots selection hidden in some computation. As an example, when developing the Carlson integrals, I have noticed spurious roots switch in the middle of nowhere. I tracked it down to the DLMF equation 19.21.10, where the last term \sqrt{\frac{xy}{z}} should really be computed as \frac{\sqrt{x}\sqrt{y}}{\sqrt{z}} (i.e compute three roots independently and compute the fraction afterwards. Computing the fraction fist and a single root afterwards leads to the wrong root to be selected in some complex cases, despite none of x, y, and z changes quadrant (but the fraction does). The same problem occurs in equations 19.2.7 and 19.2.8. I have notified DLMF editors about this and one of the mathematical editors from the University of Edinburgh replied me and confirmed this observation was correct and these equations would be corrected in a future version. So perhaps once again there are some root changes for some input values that throw the result far away.

If anybody could help me understand what happens, it would be greatly appreciated.

[maisonobe]

Making some progress. The issue seems to lie in Carlson integral R_C, which is used to compute Carlson integral R_J using the algorithm in Carlson[2000] (not the same as Carlson[1995]), and which combined with Carlson integral R_F allows to compute Π(n; φ, k).

[maisonobe]

More progress.

I found that in Hipparchus current implementation, Π (n = 3.4-1.3i, φ=1.2, m=0.2+0.6i) was correct but Π (n = 3.4-1.3i, φ=1.2+0.75i, m=0.2+0.6i) was very wrong, so I made a simple loop and found a switch point at φ=1.2+0.2505i with result jumping far away. What changed at this point was integral R_C, which is used inside the computation of integral R_J.

R_J is computed using a duplication formula (i.e. an iterative algorithm which converges in just a few iterations) plus an underlying computation of the R_C integral. In Carlson[1995], the R_C integral is computed at each iteration (this is what is done in elliptic_integrals.js). In Carlson[2000] an improvement was made and only one R_C integral is computed at the end thus enhancing performances (this is what is cone in Hipparchus).

I first thought that this was the culprit and changed the implementation to compute R_C at each iteration. It did not change the jump, which occurred with both the slower 1995 algorithm and the faster 2000 algorithm.

Further investigation showed that the change in R_C was due to changes early in the iterations in R_J(x, y, z, p). Regardless of R_C being computed once or several times, the arguments passed to R_C depend on intermediate variables in the iterations. The switch is caused by intermediate variable dₘ=(√pₘ+√yₘ)(√pₘ+√yₘ)(√pₘ+√zₘ), where for all square roots the one with nonnegative real part must be selected. At the switch point, pₘ is almost real and negative and its imaginary part changes sign, so the square root selection switches from a point along the imaginary axis from negative to positive. The iterations change a lot and the result on both sides of the switch are very different.

Reading again and more accurately the 1995 paper, I noticed an hypothesis I missed:

"Let x, y, z have nonnegative real part and at most one of them be 0, while Rep > 0.
 Alternatively, if p ≠ 0 and Iph pl < π, either let x, y, z be real and nonnegative and
 at most one of them be 0, or else let two of the variables x, y, z be nonzero and
 conjugate complex with phase less in magnitude than π and the third variable be real
 and nonnegative."

So as we have everything complex, the fourth parameter in R_f(x, y, z, p) must have a nonnegative real part!

In the case that causes problem, the p computed to convert Legendre integrals to Carlson integrals is p≈-2.82711+0.99869i. most Carlson integrals are symmetric, but R_F is not symmetric on its fourth parameter, so we cannot change parameters order at will. I will try to find what can be done tomorrow.

[maisonobe]

The part related to the integral of the third kind has its own issue number: #152.