Root Finding Without Derivatives

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

概要

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

namespace boost{ namespace math{ namespace tools{

template <class F, class T, class Tol>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

template <class F, class T, class Tol>
std::pair<T, T> 
   bracket_and_solve_root(
      F f, 
      const T& guess, 
      const T& factor, 
      bool rising, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   bracket_and_solve_root(
      F f, 
      const T& guess, 
      const T& factor, 
      bool rising, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

template <class F, class T, class Tol>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

template <class F, class T, class Tol>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      const T& fa, 
      const T& fb, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      const T& fa, 
      const T& fb, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

// 結束條件:
template <class T>
struct eps_tolerance;

struct equal_floor;
struct equal_ceil;
struct equal_nearest_integer;

}}} // namespaces

描述

這些函數求解某些函數 f(x) 的根而不需要計算函數 f(x) 的導數。這些函數在這裡使用TOMS 算法 748，這個算法是已知的最有效的漸近算法，並且對於一些光滑的函數(smooth functions)進行了優化。

另一方面，在一些情況下，有一個平分函數是很有用的，在調用一個更高級的算法之前，用於將包含根的範圍進行縮小。

在這一部分中的所有算法使用同樣的漸近效率(asymptotic efficiency)來減小包含根的區間的直徑。這與基於導數的算法是相反的，使用導數的算法可能永遠都不會明顯地減小包含根的區間的大小，即使算法快速地趨近於根。這與其它的無導數算法(derivative-free methods)(例如Brent 或 Dekker 算法)只在計算的最後一步減小包含根的區間的大小也是相反的。因此，這些算法返回一個包含找到的區間的std::pair ，並接收一個函數對像作為結束條件。對於已生成(ready-made)的結束條件已提供了3個函數對像：eps_tolerance 會在相對誤差小於一個特定值的時候結束計算過程，而equal_floor 和 equal_ceil 對於一些計算結果已知是整數的統計分佈是很有用的。其它用戶自定義的結束條件也類似地使用，這種情況很少，但是在一些特定的情況下可能很有用。

template <class F, class T, class Tol>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   bisect(
      F f, 
      T min, 
      T max, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

這些函數使用二分法來定位函數的根，函數參數為：

f: 一個一元的函數對象，為這個函數查找根。
min: 包含這個根的左邊界值。
max: 包含這個根的右邊界值。前提條件為min < max 且 f(min)*f(max) <= 0，如果違反了前提條件，那麼這個函數會產生一個計算錯誤。這種行為由計算錯誤處理策略來控制。返回一個最佳的猜測值，可能明顯是錯誤的。
tol: 一個用於指定結束條件的二元函數對像：當tol(min,max) 變為真的時候，這個函數返回包含根的區間。
max_iter: 當查找根的時候調用函數f(x) 的最大次數。

最後一個策略參數是可選的，並可以用來控制函數的行為：如何處理錯誤，使用哪種層次的精度等等，參考策略文檔瞭解詳細信息。

返回包含根的數對r使得：

f(r.first) * f(r.second) <= 0

或者

tol(r.first, r.second) == true

或

max_iter >= m

其中 m 是傳遞給函數的參數 max_iter 的初始值。

換句話說，由調用者來確保是否由超過max_iter 次函數調用而導致計算結束(通過在函數返回時檢查max_iter 的值)，而不是因為結束條件tol 得到滿足。

template <class F, class T, class Tol>
std::pair<T, T> 
   bracket_and_solve_root(
      F f, 
      const T& guess, 
      const T& factor, 
      bool rising, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   bracket_and_solve_root(
      F f, 
      const T& guess, 
      const T& factor, 
      bool rising, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

這是一個方便易用的函數，在庫的內部稱之為toms748_solve ，用於查找函數f(x)的根。僅當函數f(x)是一個單調函數，並且根的位置是近似知道的時候才可以使用函數toms748_solve，尤其是知道根為正數還是為負數。參數為：

f: 一個一元的函數對象，為這個函數查找根。 f(x)必須在x上一致遞增或一致遞減。
guess: 根的初始逼近值。
factor: 一個用於包含根的尺度因子：guess 的值乘以 (或者如果恰當則除以)factor直到找到包括(bracket)根的兩個值。類似於2的值通常是因子factor的一個選擇。
rising: 如果函數f(x) 關於x遞增則設置為true，如果函數f(x)關於x遞減則設置為false。這個值與f(guess) 的值一起使用來判斷guess 是大於還是小於根。
tol: 一個用於判定結束條件的二元的函數對像：tol 在每一步都傳遞給當前的區間，直到返回真的時候就將當前的區間作為返回值。
max_iter: 在查找根的過程中所使用的最大的函數調用次數。

最後一個策略參數是可選的，並可以用來控制函數的行為：如何處理錯誤，使用哪種層次的精度等等，參考策略文檔瞭解詳細信息。

返回包含根的數對r使得：

f(r.first) * f(r.second) <= 0

或者

tol(r.first, r.second) == true

或

max_iter >= m

其中 m 是傳遞給函數的參數 max_iter 的初始值。

換句話說，由調用者來確保是否由超過max_iter 次函數調用而導致計算結束(通過在函數返回時檢查max_iter 的值)，而不是因為結束條件tol 得到滿足。

template <class F, class T, class Tol>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      Tol tol, 
      boost::uintmax_t& max_iter);

template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

template <class F, class T, class Tol>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      const T& fa, 
      const T& fb, 
      Tol tol, 
      boost::uintmax_t& max_iter);
      
template <class F, class T, class Tol, class Policy>
std::pair<T, T> 
   toms748_solve(
      F f, 
      const T& a, 
      const T& b, 
      const T& fa, 
      const T& fb, 
      Tol tol, 
      boost::uintmax_t& max_iter,
      const Policy&);

這兩個函數實現 TOMS Algorithm 748：它是三次方程，二次方程以及線性插值( cubic, quadratic and linear interpolation)的混合使用，用於定位函數f(x)的根。這兩個函數的區別在於f(a) 和 f(b) 的值是否是已經可用的：

f: 一個被查找根的二元函數。 f(x) 不需要在x 上一致遞增或一致遞減，並且可以有多個根。
a: 初始的包含根的區間的下邊界值。
b: 初始的包含根的區間的上邊界值。前提條件為a < b 且 a 和b 包括(bracket)要查找的根，使得f(a)*f(b) < 0。
fa: 可選的： f(a)的值。
fb: 可選的： f(b)的值。
tol: 一個用於判定結束條件的二元的函數對像：tol 在每一步都傳遞給當前的區間，直到返回真的時候就將當前的區間作為返回值。
max_iter: 在查找根的過程中所使用的最大的函數調用次數。在退出的時候max_iter 的值設置為實際調用的函數的次數。

最後一個策略參數是可選的，並可以用來控制函數的行為：如何處理錯誤，使用哪種層次的精度等等，參考策略文檔瞭解詳細信息。

返回包含根的數對r使得：

f(r.first) * f(r.second) <= 0

或者

tol(r.first, r.second) == true

或

max_iter >= m

其中 m 是傳遞給函數的參數 max_iter 的初始值。

換句話說，由調用者來確保是否由超過max_iter 次函數調用而導致計算結束(通過在函數返回時檢查max_iter 的值)，而不是因為結束條件tol 得到滿足。

template <class T>
struct eps_tolerance
{
   eps_tolerance(int bits);
   bool operator()(const T& a, const T& b)const;
};

這是這些根查找算法通常所使用的結束條件：當a 和b 之間的相對距離小於兩倍的T類型的機器最小精度值( machine epsilon for T)，或者 2^1-bits，無論這兩個值中哪一個更大。換句話說，你將bits 設置為結果中你所想要的bit數目。大小為T類型的最小機器精度的2倍的最小容許度( minimal tolerance)確保我們得到一個括起來的區間(bracketing interval)：因為這個區間至少為1個機器精度大小。

struct equal_floor
{
   equal_floor();
   template <class T> bool operator()(const T& a, const T& b)const;
};

當你想要查找一個真實根向下捨入的整數值的結果時，你使用這個結束條件。當區間的兩個邊界值向下捨入值相同時這個函數就立即返回。

struct equal_ceil
{
   equal_ceil();
   template <class T> bool operator()(const T& a, const T& b)const;
};

當你想要查找一個真實根向上捨入的整數值的結果時，你使用這個結束條件。當區間的兩個邊界值向上捨入值相同時這個函數就立即返回。

struct equal_nearest_integer
{
   equal_nearest_integer();
   template <class T> bool operator()(const T& a, const T& b)const;
};

當你想要查找一個真實根捨入到最接近的整數值的結果時，你使用這個結束條件。當區間的兩個邊界值捨入到最接近的整數值相同時這個函數就立即返回。

實現

二分算法的實現相當直接並且沒有在這裡詳細說明。TOMS algorithm 748 在下面的文獻中進行了詳細的描述：

Algorithm 748: Enclosing Zeros of Continuous Functions, G. E. Alefeld, F. A. Potra and Yixun Shi, ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. Pages 327-344.

這裡的實現是將這篇論文的算法如實地用C++實現。