Elliptic integrals arise occasionally in contexts relevant to us, particularly in geometric calculations involving ellipses or ellipsoids.

Here we provide functions for calculating the complete elliptic integrals of the first and second kinds, which arise in Hertz contact theory for point contacts where the contact area is elliptical. We use the following definitions for these two integrals:

    K(m) = integ[0,Pi/2] {1 / sqrt(1 - m sin^2(t)) dt}     1st kind
    E(m) = integ[0,Pi/2] {    sqrt(1 - m sin^2(t)) dt}     2nd kind
    0 <= m <= 1

Elliptic integrals are defined only for arguments in range [0,1] inclusive.

We provide a function that calculates K(m) and E(m) to machine precision (float or double) with a fast-converging iterative method adapted from ref. 1, which was in turn adapted from ref. 2. A much faster approximate version is also available, using the higher-precision approximation of the two provided in ref. 2. The approximate version provides a smooth function that gives at least 7 digits of accuracy (in either float or double precision) across the full range at about 1/4 the cost of the machine precision version. For many applications, including engineering- or scientific-quality contact, 7 digits is more than adequate and in float precision that's all you can expect anyway.

Warning: In the literature there are two different definitions for elliptic integrals. The other definition (call them K'(k) and E'(k)) uses the argument k=sqrt(m) so that K'(k)=K(m^2), E'(k)=E(m^2). Our definition is used by Matlab's ellipke() function and much engineering literature while the K',E' definitions seem to be preferred by mathematicians.

For an ellipse with semimajor axis a and semiminor axis b (so a >= b), the eccentricity e=sqrt(1-(b/a)^2). Our argument to the elliptic integrals in that case is m = e^2 = 1-(b/a)^2. In constrast, K.L. Johnson uses the mathematicians' definition in Chapter 4, pg. 95 of his book (ref. 4) where he writes K(e) and E(e), where e is eccentricity as defined above, that is e=sqrt(m), so we would call his K'(e)=K(e^2) and E'=E(e^2).

Author: Michael Sherman

References

(1) Dyson, Evans, Snidle. "A simple accurate method for calculation of stresses and deformations in elliptical Hertzian contacts", J. Mech. Eng. Sci. Proc. IMechE part C 206:139-141, 1992.

(2) Abramovitz, Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, NY, 1972.

(3) Antoine, Visa, Sauvey, Abba. "Approximate analytical model for Hertzian Elliptical Contact Problems", ASME J. Tribology 128:660, 2006.

(4) Johnson, K.L. Contact Mechanics. Cambridge University Press 1987 (corrected edition).