Длина дуги эллипса с радиусамиaиbдана:

L(a, b) = 4aE(k)

с:

k = √(1 - b2/a2)

гдеE(k)— полный эллиптический интеграл второго рода — см.ellint_2.

Предположим, что мы знаем длину дуги и один радиус, а затем можем вычислить другой радиус, перевернув формулу выше. Начнем с кодирования приведенной выше формулы в функтор, который могут вызвать наши алгоритмы поиска корней.

Обратите внимание, что, хотя это не совсем очевидно из приведенной выше формулы, функция полностью симметрична в двух радиусах, которые могут быть заменены по желанию, в этом случае нам нужно убедиться, что<a >=b>, чтобы мы случайно не взяли квадратный корень из отрицательного числа:

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

Нам также понадобится приличная оценка, чтобы начать поиск, приближение:

L(a, b) ≈ 4√(a2 + b2)

Легко перевернуть, чтобы дать нам то, что нам нужно, что с помощью поиска корней без производных приводит к алгоритму:

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;        // Initally 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)

Эта функция обычно находит корень в пределах 8-10 итераций, поэтому, учитывая, что время выполнения полностью зависит от стоимости вызова эллиптического интеграла, было бы неплохо несколько уменьшить это количество. Мы постараемся сделать это, используя метод, основанный на производных; производные этой функции довольно трудно выработать вручную, но, к счастьюWolfram Alphaможет выполнить работу, которую мы можем дать:

d/da L(a, b) = 4(a2E(k) - b2K(k)) / (a2 - b2)

Обратите внимание, что теперь у нас естьдваэллиптических интегральных вызова, чтобы получить производную, так что наш функтор будет по крайней мере в два раза дороже звонить, чем производный-свободный один выше: нам придется значительно уменьшить количество итераций, чтобы сделать разницу!

Вот пересмотренный функтор:

template <class T = double>
struct elliptic_root_functor_1deriv
{ // Functor also returning 1st derviative.
   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

Код поиска корней теперь почти такой же, как и раньше, но мы будем использовать Ньютон-итерацию, чтобы получить результат:

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

Количество итераций, требуемых для точности<double>, теперь обычно составляет около 4 - так что мы чуть более чем вдвое сократили количество итераций, но сделали функтор вдвое дороже!

Интересно, что вторая производная не требует более дорогих эллиптических интегральных вызовов, чем первая, другими словами, она поставляется по существу «бесплатно», и в этом случае мы могли бы также использовать ее и использовать Галлей-итерацию. Это довольно типичная ситуация при инвертировании специальных функций. Вот пересмотренный функтор:

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
};

Фактический код поиска корней почти такой же, как и раньше, за исключением того, что мы можем использовать итерацию Галлея, а не Ньютона:

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

В то время как эта функция использует только немного меньше итераций (обычно около 3), чтобы найти корень, по сравнению с исходным методом без производных, мы перешли от 8-10 эллиптических интегральных вызовов к 6.

Полный код этого примера находится по адресуroot_elliptic_finding.cpp.