|
Home | Libraries | People | FAQ | More |
The arc length of an ellipse with radii a and b is given by:
L(a, b) = 4aE(k)
with:
k = √(1 - b2/a2)
where E(k) is the complete elliptic integral of the second kind - see ellint_2.
Let's suppose we know the arc length and one radii, we can then calculate the other radius by inverting the formula above. We'll begin by encoding the above formula into a functor that our root-finding algorithms can call.
Note that while not completely obvious from the formula above, the function
is completely symmetrical in the two radii - which can be interchanged at
will - in this case we need to make sure that a
>= b
so that we don't accidentally take the square root of a negative number:
template <typename T = double> struct elliptic_root_functor_noderiv { // Nth root of x using only function - no derivatives. elliptic_root_functor_noderiv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius) { // Constructor just stores value a to find root of. } T operator()(T const& x) { using std::sqrt; // return the difference between required arc-length, and the calculated arc-length for an // ellipse with radii m_radius and x: T a = (std::max)(m_radius, x); T b = (std::min)(m_radius, x); T k = sqrt(1 - b * b / (a * a)); return 4 * a * boost::math::ellint_2(k) - m_arc; } private: T m_arc; // length of arc. T m_radius; // one of the two radii of the ellipse }; // template <class T> struct elliptic_root_functor_noderiv
We'll also need a decent estimate to start searching from, the approximation:
L(a, b) ≈ 4√(a2 + b2)
Is easily inverted to give us what we need, which using derivative-free root finding leads to the algorithm:
template <class T = double> T elliptic_root_noderiv(T radius, T arc) { // return the other radius of an ellipse, given one radii and the arc-length using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For bracket_and_solve_root. T guess = sqrt(arc * arc / 16 - radius * radius); T factor = 1.2; // How big steps to take when searching. const boost::uintmax_t maxit = 50; // Limit to maximum iterations. boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. bool is_rising = true; // arc-length increases if one radii increases, so function is rising // Define a termination condition, stop when nearly all digits are correct, but allow for // the fact that we are returning a range, and must have some inaccuracy in the elliptic integral: eps_tolerance<T> tol(std::numeric_limits<T>::digits - 2); // Call bracket_and_solve_root to find the solution, note that this is a rising function: std::pair<T, T> r = bracket_and_solve_root(elliptic_root_functor_noderiv<T>(arc, radius), guess, factor, is_rising, tol, it); // Result is midway between the endpoints of the range: return r.first + (r.second - r.first) / 2; } // template <class T> T elliptic_root_noderiv(T x)
This function generally finds the root within 8-10 iterations, so given that the runtime is completely dominated by the cost of calling the elliptic integral it would be nice to reduce that count somewhat. We'll try to do that by using a derivative-based method; the derivatives of this function are rather hard to work out by hand, but fortunately Wolfram Alpha can do the grunt work for us to give:
d/da L(a, b) = 4(a2E(k) - b2K(k)) / (a2 - b2)
Note that now we have two elliptic integral calls to get the derivative, so our functor will be at least twice as expensive to call as the derivative-free one above: we'll have to reduce the iteration count quite substantially to make a difference!
Here's the revised functor:
template <class T = double> struct elliptic_root_functor_1deriv { // Functor also returning 1st derivative. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); elliptic_root_functor_1deriv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius) { // Constructor just stores value a to find root of. } std::pair<T, T> operator()(T const& x) { using std::sqrt; // Return the difference between required arc-length, and the calculated arc-length for an // ellipse with radii m_radius and x, plus it's derivative. // See http://www.wolframalpha.com/input/?i=d%2Fda+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29] // We require two elliptic integral calls, but from these we can calculate both // the function and it's derivative: T a = (std::max)(m_radius, x); T b = (std::min)(m_radius, x); T a2 = a * a; T b2 = b * b; T k = sqrt(1 - b2 / a2); T Ek = boost::math::ellint_2(k); T Kk = boost::math::ellint_1(k); T fx = 4 * a * Ek - m_arc; T dfx = 4 * (a2 * Ek - b2 * Kk) / (a2 - b2); return std::make_pair(fx, dfx); } private: T m_arc; // length of arc. T m_radius; // one of the two radii of the ellipse }; // struct elliptic_root__functor_1deriv
The root-finding code is now almost the same as before, but we'll make use of Newton-iteration to get the result:
template <class T = double> T elliptic_root_1deriv(T radius, T arc) { using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For newton_raphson_iterate. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); T guess = sqrt(arc * arc / 16 - radius * radius); T min = 0; // Minimum possible value is zero. T max = arc; // Maximum possible value is the arc length. // Accuracy doubles at each step, so stop when just over half of the digits are // correct, and rely on that step to polish off the remainder: int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.6); const boost::uintmax_t maxit = 20; boost::uintmax_t it = maxit; T result = newton_raphson_iterate(elliptic_root_functor_1deriv<T>(arc, radius), guess, min, max, get_digits, it); return result; } // T elliptic_root_1_deriv Newton-Raphson
The number of iterations required for double
precision is now usually around 4 - so we've slightly more than halved the
number of iterations, but made the functor twice as expensive to call!
Interestingly though, the second derivative requires no more expensive elliptic integral calls than the first does, in other words it comes essentially "for free", in which case we might as well make use of it and use Halley-iteration. This is quite a typical situation when inverting special-functions. Here's the revised functor:
template <class T = double> struct elliptic_root_functor_2deriv { // Functor returning both 1st and 2nd derivatives. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); elliptic_root_functor_2deriv(T const& arc, T const& radius) : m_arc(arc), m_radius(radius) {} std::tuple<T, T, T> operator()(T const& x) { using std::sqrt; // Return the difference between required arc-length, and the calculated arc-length for an // ellipse with radii m_radius and x, plus it's derivative. // See http://www.wolframalpha.com/input/?i=d^2%2Fda^2+[4+*+a+*+EllipticE%281+-+b^2%2Fa^2%29] // for the second derivative. T a = (std::max)(m_radius, x); T b = (std::min)(m_radius, x); T a2 = a * a; T b2 = b * b; T k = sqrt(1 - b2 / a2); T Ek = boost::math::ellint_2(k); T Kk = boost::math::ellint_1(k); T fx = 4 * a * Ek - m_arc; T dfx = 4 * (a2 * Ek - b2 * Kk) / (a2 - b2); T dfx2 = 4 * b2 * ((a2 + b2) * Kk - 2 * a2 * Ek) / (a * (a2 - b2) * (a2 - b2)); return std::make_tuple(fx, dfx, dfx2); } private: T m_arc; // length of arc. T m_radius; // one of the two radii of the ellipse };
The actual root-finding code is almost the same as before, except we can use Halley, rather than Newton iteration:
template <class T = double> T elliptic_root_2deriv(T radius, T arc) { using namespace std; // Help ADL of std functions. using namespace boost::math::tools; // For halley_iterate. BOOST_STATIC_ASSERT_MSG(boost::is_integral<T>::value == false, "Only floating-point type types can be used!"); T guess = sqrt(arc * arc / 16 - radius * radius); T min = 0; // Minimum possible value is zero. T max = arc; // radius can't be larger than the arc length. // Accuracy triples at each step, so stop when just over one-third of the digits // are correct, and the last iteration will polish off the remaining digits: int get_digits = static_cast<int>(std::numeric_limits<T>::digits * 0.4); const boost::uintmax_t maxit = 20; boost::uintmax_t it = maxit; T result = halley_iterate(elliptic_root_functor_2deriv<T>(arc, radius), guess, min, max, get_digits, it); return result; } // nth_2deriv Halley
While this function uses only slightly fewer iterations (typically around 3) to find the root, compared to the original derivative-free method, we've moved from 8-10 elliptic integral calls to 6.
Full code of this example is at root_elliptic_finding.cpp.