Elliptic integrals

Before this page

Hill polarization tensors, whose Newton-potential integrals for a spheroid are expressed through the functions defined here.

The MeanFieldHomogenization.Elliptic submodule provides type-generic Legendre and Carlson integrals that flow end-to-end through automatic differentiation (ForwardDiff.Dual), arbitrary precision (BigFloat) and symbolic scalars (SymPyPythonCall.Sym, Symbolics.Num).

1. Legendre elliptic integrals

Throughout the submodule, m denotes the parameter (not the modulus), with the convention   and   .

Complete integrals

Incomplete integrals

Both reduce to the complete integrals at  .

2. Arithmetic–geometric mean (AGM)

For any Number subtype other than Float64, the complete integrals are evaluated using the classical AGM recursion (Abramowitz & Stegun 17.6, NIST DLMF 19.8). Starting from

the sequence

converges quadratically, and

Eight to twelve iterations are typically enough to reach Float64 precision; BigFloat needs only a few more.

3. Carlson symmetric forms

Incomplete integrals are delegated to Carlson's symmetric integrals (Carlson 1995, Numerical computation of real or complex elliptic integrals):

With these two primitives one recovers the Legendre integrals via

R_F and R_D are implemented by the duplication theorem (Carlson 1995, §2): each iteration divides the differences between by four, until their relative spread is small enough for a fifth-order Taylor series in the Carlson invariants (and for ). The recursion is arithmetic-only — no branch cuts, no transcendentals — so it extends to any Number subtype. A scalar whose comparisons do not return a Bool (a complex number) has no stopping test and runs a fixed number of duplications; a symbolic one never gets there, the SymPy and Symbolics extensions keeping the integrals as unexpanded calls.

4. Dispatch table

Scalar type	Backend
Float64	Elliptic.jl (GSL C binding, fastest)
ForwardDiff.Dual	pure-Julia AGM / Carlson (derivatives work)
BigFloat, generic	AGM / Carlson
SymPyPythonCall.Sym	sympy.elliptic_{k,e,f} via the SymPy weak ext
Symbolics.Num	unexpanded calls, via the Symbolics weak ext

Why a SymPy weak extension?

The AGM unrolls ~60 nested sqrt(a*b) operations. On a SymPyPythonCall.Sym input this builds a deeply nested symbolic tree that overflows SymPy's pretty-printer. The weak extension MeanFieldHomogenizationSymPyPythonCallExt (loaded automatically whenever SymPy is loaded alongside MeanFieldHomogenization) routes ell_K, ell_E, ell_F on Sym arguments directly to sympy.elliptic_{k,e,f}, returning the native closed form instead.

5. Special cases

•  :   ;   .

•  :     (logarithmic divergence);  . The Float64 fast path through Elliptic.jl throws at  .

•  : incomplete integrals coincide with complete ones.

6. References

• M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, §17.6, Dover 1972.

• NIST Digital Library of Mathematical Functions, §19.

• B.C. Carlson, Numerical computation of real or complex elliptic integrals, Numerical Algorithms 10 (1995) 13–26.

• B.C. Carlson and E.M. Notis, Algorithms for incomplete elliptic integrals, ACM Transactions on Mathematical Software 7 (1981) 398–403 — public-domain SLATEC routines DRF / DRD.

Where to go next

The submodule is exercised on numerical, automatic-differentiation, arbitrary precision and symbolic inputs in Elliptic integrals — examples, and its entry points are listed in Elliptic integrals — API. The theory chapter closes here; Getting started and the Manual put the whole chain to use.