|
Simbody
3.5
|
|
Simbody
3.5
|
Elliptic integrals arise occasionally in contexts relevant to us, particularly in geometric calculations involving ellipses or ellipsoids. More...
Functions | |
| std::pair< double, double > | SimTK::approxCompleteEllipticIntegralsKE (double m) |
| Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m), approximated but with a maximum error of 2e-8 so at least 7 digits are correct (same in float or double precision). See Elliptic integrals for a discussion. More... | |
| std::pair< float, float > | SimTK::approxCompleteEllipticIntegralsKE (float m) |
| This is the single precision (float) version of the approximate calculation of elliptic integrals, still yielding about 7 digits of accuracy even though all calculations are done in float precision. More... | |
| std::pair< double, double > | SimTK::approxCompleteEllipticIntegralsKE (int m) |
| This integer overload is present to prevent ambiguity; it converts its argument to double precision and then calls approxCompleteEllipticIntegralsKE(double). More... | |
| std::pair< double, double > | SimTK::completeEllipticIntegralsKE (double m) |
| Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m), calculated to (roughly) machine precision (float or double). See Elliptic integrals for a discussion. More... | |
| std::pair< float, float > | SimTK::completeEllipticIntegralsKE (float m) |
| This is the single precision (float) version of the machine-precision calculation of elliptic integrals, providing accuracy to float precision (about 7 digits) which is no better than you'll get with the much faster approximate version, so use that instead! More... | |
| std::pair< double, double > | SimTK::completeEllipticIntegralsKE (int m) |
| This integer overload is present to prevent ambiguity; it converts its argument to double precision and then calls completeEllipticIntegralsKE(double). More... | |
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.
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).
(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).
|
inline |
Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m), approximated but with a maximum error of 2e-8 so at least 7 digits are correct (same in float or double precision). See Elliptic integrals for a discussion.
Note that a full-precision computation of these integrals is also available; see completeEllipticIntegralsKE(). However, if you are using float precision there is no point in using the exact computation since this approximation is just as good and much faster.
| [in] | m | The argument to the elliptic integrals. Must be in range [0,1] although we allow for a very small amount of numerical error (see SignificantReal) that might put m outside that range. |
You can find the approximating formula we're using in ref. 2, sections 17.3.34 and 17.3.36, pp. 591-2, or you can look at the code for this function (in the header file). The formulas are accurate to 2e-8 over the full [0-1] argument range, according to the authors. I checked our implementation against Matlab's ellipke() function and the results are very good, providing at least 7 correct digits for a range of sample values.
The cost is about 90 flops.
(1) Antoine, Visa, Sauvey, Abba. "Approximate analytical model for Hertzian Elliptical Contact Problems", ASME J. Tribology 128:660, 2006.
(2) Abramovitz, Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, NY, 1972.
|
inline |
This is the single precision (float) version of the approximate calculation of elliptic integrals, still yielding about 7 digits of accuracy even though all calculations are done in float precision.
|
inline |
This integer overload is present to prevent ambiguity; it converts its argument to double precision and then calls approxCompleteEllipticIntegralsKE(double).
Note that the only legal values here are 0 and 1.
|
inline |
Given 0<=m<=1, return complete elliptic integrals of the first and second kinds, K(m) and E(m), calculated to (roughly) machine precision (float or double). See Elliptic integrals for a discussion.
Note that we also provide a faster approximate method for calculating these functions (see approxCompleteEllipticIntegralsKE()). The approximate method is good enough for many scientific and engineering applications and is always preferred if you are using float precision.
| [in] | m | The argument to the elliptic integrals. Must be in range [0,1] although we allow for a very small amount of numerical error (see SignificantReal) that might put m outside that range. |
Cost here is about 50 + 50*n flops, where n is the number of iterations required. The number of iterations n you can expect to get a double-precision result is 4 for a:b < 2, 5 for a:b < 20, 6 for a:b < 1000, 7 after that; for single-precision it will take one fewer iterations. In flops that's 250, 300, 350, 400 – 300 will be typical for double, 250 for float.
|
inline |
This is the single precision (float) version of the machine-precision calculation of elliptic integrals, providing accuracy to float precision (about 7 digits) which is no better than you'll get with the much faster approximate version, so use that instead!
|
inline |
This integer overload is present to prevent ambiguity; it converts its argument to double precision and then calls completeEllipticIntegralsKE(double).
Note that the only legal values here are 0 and 1.
1.8.8