Boost C++ Libraries Home Libraries People FAQ More

PrevUpHomeNext

使用導數查找根(Root Finding With Derivatives)

概要

#include <boost/math/tools/roots.hpp>

namespace boost{ namespace math{ namespace tools{

template <class F, class T>
T newton_raphson_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

template <class F, class T>
T halley_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

template <class F, class T>
T schroeder_iterate(F f, T guess, T min, T max, int digits);

template <class F, class T>
T schroeder_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter);

}}} // namespaces
說明

這些函數都進行迭代根查找:newton_raphson_iterate 進行二階Newton Raphson iteration,而halley_iterateschroeder_iterate 分別進行三階HalleySchroeder 迭代。

這些函數都帶有同樣的參數:

根查找函數的參數

F f

類型F必須是一個可調用的函數對象,接收一個參數並返回一個元組(tuple):

對於二階迭代算法 (Newton Raphson) 這個元組必須包含兩個元素:函數的值以及它的一階導數。

對於三階迭代算法 (Halley and Schroeder)這個元組應當有三個元素:包含函數的值以及它的一階和二階導數。

T guess

起始值。

T min

結果的最小可能值,這個值作為區間的下邊界值。

T max

結果的最大可能值,這個值作為區間的上邊界值。

int digits

預期的二進制數字個數。

uintmax_t max_iter

一個可選的進行迭代的最大次數。

當使用這些函數時,你應當注意:

Newton Raphson 算法

給定一個初始的估測值 x0 ,後續的值使用下面的等式計算:

超出範圍之後 回退到當前範圍的二分逼近。

在理想情況下,在每次迭代中,正確的數字個數都在成倍增加。

Halley's 算法

給定一個初始的估測值 x0 ,後續的值使用下面的等式計算:

由二階導數的過度補償(Over-compensation)(會在錯誤的方向進行下去)引起這個算法退回到 Newton-Raphson 。

超出範圍之後 回退到當前範圍的二分逼近。

在理想情況下,在每次迭代中,正確的數字個數都在成三倍地增加。

Schroeder's 算法

給定一個初始的估測值 x0 ,後續的值使用下面的等式計算:

由二階導數的過度補償(Over-compensation)(會在錯誤的方向進行下去)引起這個算法退回到 Newton-Raphson 。類似地,無論何時如果Newton步驟可以將下一個值改變超過10%,那麼將會使用 Newton 步驟。

超出範圍之後 回退到當前範圍的二分逼近。

在理想情況下,在每次迭代中,正確的數字個數都在成三倍地增加。

例子

假設我們想要計算一個數的立方根,我們將求解的等式以及它的導數是:

讓我們使用 Newton Raphson 迭代開始,我們將定義一個返回求解的函數的計算方法的函數對象,以及它的一階導數:

template <class T>
struct cbrt_functor
{
   cbrt_functor(T const& target) : a(target){}
   std::tr1::tuple<T, T> operator()(T const& z)
   {
      T sqr = z * z;
      return std::tr1::make_tuple(sqr * z - a, 3 * sqr);
   }
private:
   T a;
};

實現立方根現在非常簡單了,最難的部分是查找一個好的逼近值:在這種情況下,我們將僅僅把指數除以3:

template <class T>
T cbrt(T z)
{
   using namespace std;
   int exp;
   frexp(z, &exp);
   T min = ldexp(0.5, exp/3);
   T max = ldexp(2.0, exp/3);
   T guess = ldexp(1.0, exp/3);
   int digits = std::numeric_limits<T>::digits;
   return tools::newton_raphson_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits);
}

使用 libs/math/test/cbrt_test.cpp 中的測試數據,在每種情況下,發現這種計算精確到最後一個數字,並且在double精確度時不超過6次迭代。然而,你會發現在這個例子中使用了一個更高的精度,精度是在早一些的時候在這個文檔中提醒的那樣!在這個特殊的例子中,可以精確地計算函數f(x)而不會有過度的消去錯誤(cancellation error ),所以,一個高的限制並不是一個太大的問題。然而將限制降低到std::numeric_limits<T>::digits * 2 / 3 在所有的測試情況下(除了一個)給出了完全的精度。最大的迭代次數仍為6,但在大多數的情況下減少了一次迭代。

注意,上面的代碼忽略了錯誤處理,並且沒有正確地處理z的負值。這將留作讀者的一個練習!

現在讓我們稍微修改一下這個函數對像來返回二階導數:

template <class T>
struct cbrt_functor
{
   cbrt_functor(T const& target) : a(target){}
   std::tr1::tuple<T, T, T> operator()(T const& z)
   {
      T sqr = z * z;
      return std::tr1::make_tuple(sqr * z - a, 3 * sqr, 6 * z);
   }
private:
   T a;
};

並且將函數cbrt 修改為使用 Halley 迭代:

template <class T>
T cbrt(T z)
{
   using namespace std;
   int exp;
   frexp(z, &exp);
   T min = ldexp(0.5, exp/3);
   T max = ldexp(2.0, exp/3);
   T guess = ldexp(1.0, exp/3);
   int digits = std::numeric_limits<T>::digits / 2;
   return tools::halley_iterate(detail::cbrt_functor<T>(z), guess, min, max, digits);
}

注意,在精度僅為完全精度的一半時,迭代就被設置為終止,即使是這樣,沒有一個測試情況有一個bit錯誤。更重要的是,最大的迭代次數現在僅為4。

讓我進行進一步的完善,在最後一個例子中,我們調用schroeder_iterate :而實際上,在這個特殊的例子中,在精度和迭代次數上並沒有產生差別。然而,這兩種方法的相關的性能依賴於函數 f(x)的性質,以及初始的估測值可以計算到的精度。除了「嘗試並觀察」,似乎沒有通用的方法。

最後,我們已經使用精確度為1000 bit 的 NTL::RR 庫調用函數cbrt,只需要7次迭代就可以獲得完全的精度。將精度增加一個數量為20的因子,迭代次數的增加少於2倍。這僅僅是為了強調所有的迭代次數都在獲取最開始的正確數字時用完了:在這之後,這些方法可以以很高的效率計算出更多的正確數字。或者用另一種方式來陳述:沒有什麼可以打敗一個好的初始估測值!


PrevUpHomeNext