The integral Π(n, Φ, m) is defined as
φ is the amplitude, n is the elliptic characteristic, m is the ellipse parameter (some conventions use the elliptic modulus k such that m=k² instead of the parameter m, or they use the nome q). All variables that appear in this expression may be complex.
The integrand has poles corresponding to
and their periodic repetitions due to the inverse sines.
The integral is expected to be computed over the straight path from 0 to φ. With this assumption, the integral is a single-valued function. The following picture shows the value of the integral for n=3.4-1.3i and m=0.2+0.6i.
One can clearly see the pole near 0.537+0.110i as the start point of a ray separating a green and a purple zone. Domain coloring uses hue to represent phase plus periodic brightness to enhance both modulus and phase. The sharp green/purple transition on the right hand side shows there is a discontinuity when crossing this ray. There is another pole (less visible) near 0.785-0.909i and another ray. Drawing the integral on a larger range would show other rays dues to the repetitions of these poles.
If we take care to never cross the rays cast by poles, we could in fact compute the integral using other paths than the straight line from 0 to φ. We could for example start parallel to the real axis and then parallel to the imaginary axis, or any other path.
If on the other hand we compute the integral using a path that crosses these rays (i.e. a path that goes after the poles and then bends upwards or downward before reaching φ, then we don’t notice that we cross the rays because in fact there is no pole along the path, the poles are really isolated. So in this case, we compute a value of the integral that is finite and looks perfectly reasonable, but is different from the value that would be computed by the straight path. The following example shows one example of the results we could obtain.
We see that the purple region extended upwards in a wave-like shape. In fact, we have chosen a different sheet of the Riemann surface that represent the integral value.
The first image was really computed using numerical integration. The second image was in fact computed using the Carlson transformation that allows much faster results.
Numerical integration is not only very slow, it sometimes fails to converge. The tiny white points on the ray cast by the first pole correspond to failures: the integrator exceeded its maximum number of iterations (despite it was huge). The lower ray is also in fact probably false numerically, there should be a more pronounced color change, here we seem to just have some darker purple which I think is wrong.
Carlson transformations are extremely fast, but they obviously select the wrong value without notice. The changes occur when during internal iterations we compute a value λₘ = √xₘ √yₘ+ √xₘ √zₘ+ √yₘ √yₘ. According to Carlson, each root must be computed separately and in each case the root with nonnegative real part must be selected (this is important because √xₘ √yₘ and √(xₘyₘ) may lead to different roots selection in the complex plane). If during computation one the value crosses the real axis while being negative, then the root selection switches from one root to its conjugate: the imaginary part changes its sign instantly. After the switch, the iterative algorithm continues flawlessly and converges. It may even converge to the same median point that it had before the switch, but as the intermediate values are used in the computation of Π(n, Φ, m) (they are not used in the computations of the integrals of the first kind F(φ,m) and of the second kind E(φ, m)), the final value computed for the integral is not the correct one. The algorithm has switched from one sheet of the Riemann surface to another one, and we have no way to know it. It is even weird as when we are in the vicinity of the ray, the truth would be to see a discontinuity, but the result we get is that the algorithm somehow recreates a continuous surface, which when looking only at a few values and the evolution in the neighborhood seems correct.
I have checked the results against Wolfram Alpha as this is what I used to get reference values I trust. Unfortunately, I cannot draw the same picture using a free account at Wolfram cloud because this computation exceeds the allowed time. So I had to just use the forms on Wolfram Alpha site and get values one by one… My check was to compute the integrals with the end point φ
moving from 1.2-1.5i to 1.2+0.75i. This means that despite the initial point of my integral was always 0+0i, the end point of my integral was moving along a line parallel to the imaginary axis. So the first integrals were computed using an horizontal path and ended up in the purple region, whereas the last integral were computed using a slanted path that was above the pole and ended up in the green region. Of course, there were some specific value of φ that came very close to the pole and ended up on the singularity. The various columns of the table are:
- integral end point
- numerical integration using a straight path from 0 to φ
- numerical integration using a path with an intermediate point at pole+i (i.e. above the pole)
- numerical integration using a path with an intermediate point at pole-i (i.e. below first pole, and then above second pole)
- numerical integration using a path with two intermediate points below both poles
- computation using Carlson transforms
- reference values from Wolfram Alpha
| φ | straight integration | integration ⇗⇘ | integration ⇘⇒ | integration ⇘⇗ | Carlson-based | WolframAlpha |
|---|---|---|---|---|---|---|
| 1.2 -1.5000000000 | 0.067423 -0.689888 | -0.473719 0.953654 | 0.141937 -0.823344 | 0.034362 -0.584955 | 0.141937 -0.823344 | 0.033512 -0.575665 |
| 1.2 -1.4000000000 | 0.119416 -0.777162 | -0.470786 0.956006 | 0.144870 -0.820991 | 0.037511 -0.583373 | 0.144870 -0.820991 | 0.036445 -0.573313 |
| 1.2 -1.3891907650 | 0.151145 -0.812147 | -0.470405 0.956297 | 0.145251 -0.820700 | 0.037901 -0.583188 | 0.145251 -0.820700 | 0.036826 -0.573022 |
| 1.2 -1.3000000000 | 0.149002 -0.817999 | -0.466655 0.959000 | 0.149001 -0.817998 | 0.041726 -0.581407 | 0.149001 -0.817998 | 0.040576 -0.570319 |
| 1.2 -1.2000000000 | 0.154831 -0.814314 | -0.460825 0.962684 | 0.154831 -0.814314 | 0.047644 -0.578858 | 0.154831 -0.814314 | 0.046406 -0.566636 |
| 1.2 -1.0666819680 | 0.166406 -0.808529 | -0.449250 0.968468 | 0.166406 -0.808529 | 0.059532 -0.574296 | 0.166406 -0.808529 | 0.057981 -0.560851 |
| 1.2 -1.0666819670 | 0.166406 -0.808529 | -0.449250 0.968468 | 0.166406 -0.808529 | 0.059532 -0.574296 | 0.166406 -0.808529 | 0.166406 -0.808529 |
| 1.2 -1.0000000000 | 0.174332 -0.805526 | -0.441324 0.971472 | 0.174332 -0.805526 | 0.067533 -0.572163 | 0.174332 -0.805526 | 0.174332 -0.805526 |
| 1.2 -0.7500000000 | 0.219717 -0.798532 | -0.395938 0.978466 | 0.219717 -0.798532 | 0.113451 -0.568136 | 0.219717 -0.798532 | 0.219717 -0.798532 |
| 1.2 -0.5000000000 | 0.288633 -0.808119 | -0.327022 0.968878 | 0.288633 -0.808119 | 0.182994 -0.580982 | 0.288633 -0.808119 | 0.288633 -0.808119 |
| 1.2 0.0000000000 | 0.461999 -0.906689 | -0.153656 0.870309 | 0.461999 -0.906689 | 0.357959 -0.686343 | 0.461999 -0.906689 | 0.461999 -0.906689 |
| 1.2 0.0500000000 | 0.477681 -0.922628 | -0.137974 0.854370 | 0.477681 -0.922628 | 0.373667 -0.703568 | 0.477681 -0.922628 | 0.477681 -0.922628 |
| 1.2 0.0700000000 | 0.483659 -0.929206 | -0.131997 0.847792 | 0.483659 -0.929206 | 0.379664 -0.710612 | 0.483659 -0.929206 | 0.483659 -0.929206 |
| 1.2 0.0800000000 | 0.486579 -0.932531 | -0.129077 0.844467 | 0.486579 -0.932531 | 0.382609 -0.714177 | 0.486579 -0.932531 | 0.486579 -0.932531 |
| 1.2 0.0850000000 | 0.488021 -0.934202 | -0.127635 0.842796 | 0.488021 -0.934202 | 0.384064 -0.715968 | 0.488021 -0.934202 | 0.488021 -0.934202 |
| 1.2 0.0851810000 | 0.488073 -0.934263 | -0.127583 0.842735 | 0.488073 -0.934263 | 0.384116 -0.716033 | 0.488073 -0.934263 | 0.488073 -0.934263 |
| 1.2 0.0851820000 | 0.488073 -0.934263 | -0.127582 0.842735 | 0.488073 -0.934263 | 0.384117 -0.716034 | 0.488073 -0.934263 | 1.103729 -2.711260 |
| 1.2 0.0852449100 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | 1.103747 -2.711282 |
| 1.2 0.0852449200 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | -1.358876 4.396709 |
| 1.2 0.0852450100 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | -1.358876 4.396709 |
| 1.2 0.0852450200 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | -0.127564 0.842714 |
| 1.2 0.0852450500 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | -0.127564 0.842714 |
| 1.2 0.0852451000 | 0.488091 -0.934284 | -0.127564 0.842714 | 0.488091 -0.934284 | 0.384135 -0.716056 | 0.488091 -0.934284 | -0.127564 0.842714 |
| 1.2 0.0860000000 | 0.488308 -0.934537 | -0.127348 0.842461 | 0.488308 -0.934537 | 0.384353 -0.716327 | 0.488308 -0.934537 | -0.127348 0.842461 |
| 1.2 0.0870000000 | 0.488595 -0.934872 | -0.127061 0.842126 | 0.488595 -0.934872 | 0.384643 -0.716686 | 0.488595 -0.934872 | -0.127061 0.842126 |
| 1.2 0.0900000000 | 0.489451 -0.935878 | -0.126204 0.841119 | 0.489451 -0.935878 | 0.385571 -0.717296 | 0.489451 -0.935878 | -0.126204 0.841119 |
| 1.2 0.1000000000 | 0.492276 -0.939245 | -0.123380 0.837753 | 0.492276 -0.939245 | 0.388421 -0.720900 | 0.492276 -0.939245 | -0.123380 0.837753 |
| 1.2 0.2000000000 | 0.517589 -0.973584 | -0.097969 0.803432 | 0.517687 -0.973566 | 0.414396 -0.755743 | 0.517687 -0.973566 | -0.097969 0.803432 |
| 1.2 0.2049000000 | 0.518611 -0.975291 | -0.096861 0.801740 | 0.518795 -0.975258 | 0.415518 -0.757550 | 0.518795 -0.975258 | -0.096861 0.801740 |
| 1.2 0.2051631601 | 0.518664 -0.975383 | -0.096802 0.801649 | 0.518854 -0.975349 | 0.415578 -0.757647 | 0.518854 -0.975349 | -0.096802 0.801649 |
| 1.2 0.2051631602 | 0.518664 -0.975383 | -0.096802 0.801649 | 0.518854 -0.975349 | 0.415578 -0.757647 | -0.712458 2.578646 | -0.096802 0.801649 |
| 1.2 0.2400000000 | 0.526267 -0.987204 | -0.089309 0.789657 | 0.526347 -0.987341 | 0.423172 -0.770470 | -0.704965 2.566654 | -0.089309 0.789657 |
| 1.2 0.2462000000 | 0.527856 -0.990561 | -0.088045 0.787533 | 0.527611 -0.989464 | 0.424455 -0.772744 | -0.703700 2.564531 | -0.088045 0.787533 |
| 1.2 0.2467578160 | - - | -0.087932 0.787342 | 0.527724 -0.989655 | 0.424570 -0.772948 | -0.703588 2.564340 | -0.087932 0.787342 |
| 1.2 0.2475000000 | -0.087274 0.784663 | -0.087782 0.787088 | 0.527874 -0.989909 | 0.424721 -0.773220 | -0.703438 2.564086 | -0.087782 0.787088 |
| 1.2 0.2500000000 | -0.087271 0.786103 | -0.087280 0.786234 | 0.528376 -0.990764 | 0.425080 -0.774295 | -0.702936 2.563231 | -0.087280 0.786234 |
| 1.2 0.4500000000 | -0.057161 0.722395 | -0.057161 0.722395 | 0.558494 -1.054603 | 0.455863 -0.841550 | -0.672817 2.499392 | -0.057161 0.722395 |
| 1.2 0.7500000000 | -0.037762 0.653477 | -0.037762 0.653477 | 0.577893 -1.123520 | 0.476489 -0.915387 | -0.653418 2.430475 | -0.037762 0.653477 |
This table shows the discontinuities. We see that as φ goes from 1.2 -1.5i to 1.2+0.75i, all method have singularities. There are ranges in which several agree, and ranges in which some switch to the wrong sheet of the Riemann surface. Comparing column 2 (straight integration) and column 6 (Carlson transforms), corresponds to look in the first two pictures at the vertical line with real value 1.2.
We see that for φ=1.2+0.2467578160i, the straight integral failed to compute, and we see that below and above this value, the integral exhibits a discontinuity. This is correct and this is the result we want. For the low values of φ however, the values of the numerical integral seem strange to me, I think it failed to converge.
For low values, Carlson based computation is on fact quite good, but it fails after some time, and when it switches, it switches to some wrong sheet that corresponds to none of the three integration paths I used.
Surprisingly, Wolfram Alpha seems to be quite wrong in a number of places. Around φ=1.2+0.0852, it even switches twice in a very short interval, and then switch to the same value of the above integral. This means in Wolfram Alpha, the green area that should be above the ray extends much below.
A suggestion to fall back to numerical integration when Carlson fails has been attempted but the problem is that there is no sign when Carlson fails. It computes something, and what it computes is realistic. Carlson gives one sufficient (but not necessary) condition for success, with one of the intermediate variables that should have a positive real part. This condition however corresponds to a very small zone (roughly the red-yellow bubble near origin) and numerical integration sometimes fails dramatically (with an exception triggered by maximum iteration reached).
So here I am on this development:
- the elliptic integrals of the first, second and third kind are available, in incomplete and complete cases, with primitive double, fields, complex and field complex
- the integrals seem to work well on real case (or Field, as long as it is not the complex field)
- the incomplete integral of the third kind in complex (and probably field complex) has issues
- the incomplete integral of the first kind also has an issue (see issue 151), but it is rather related to argument reduction
- both issues are difficult to solve, the first one being extremely difficult
- even Wolfram Alpha, despite being a reference I trust, seems to have issues with these integrals and produces sometimes obviously wrong results
I don’t think I will be able to solve these issues.
I also really need to have 2.0 published.
I see four different ways to go forward
- find someone who could help me solve this (I am on the verge of a burnout with this stuff)
- drop support for Complex/FieldComplex in elliptical integrals for 2.0 and add it later
- drop support for all elliptical integrals for 2.0 and add it later
- publish the version as is, with big warnings that Complex/FieldComplex in elliptical integrals
are provided as experimental code with issues and should not be relied upon yet
What do you think?