roots_without_derivatives.qbk 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499
  1. [section:roots_noderiv Root Finding Without Derivatives]
  2. [h4 Synopsis]
  3. ``
  4. #include <boost/math/tools/roots.hpp>
  5. ``
  6. namespace boost { namespace math {
  7. namespace tools { // Note namespace boost::math::tools.
  8. // Bisection
  9. template <class F, class T, class Tol>
  10. std::pair<T, T>
  11. bisect(
  12. F f,
  13. T min,
  14. T max,
  15. Tol tol,
  16. boost::uintmax_t& max_iter);
  17. template <class F, class T, class Tol>
  18. std::pair<T, T>
  19. bisect(
  20. F f,
  21. T min,
  22. T max,
  23. Tol tol);
  24. template <class F, class T, class Tol, class ``__Policy``>
  25. std::pair<T, T>
  26. bisect(
  27. F f,
  28. T min,
  29. T max,
  30. Tol tol,
  31. boost::uintmax_t& max_iter,
  32. const ``__Policy``&);
  33. // Bracket and Solve Root
  34. template <class F, class T, class Tol>
  35. std::pair<T, T>
  36. bracket_and_solve_root(
  37. F f,
  38. const T& guess,
  39. const T& factor,
  40. bool rising,
  41. Tol tol,
  42. boost::uintmax_t& max_iter);
  43. template <class F, class T, class Tol, class ``__Policy``>
  44. std::pair<T, T>
  45. bracket_and_solve_root(
  46. F f,
  47. const T& guess,
  48. const T& factor,
  49. bool rising,
  50. Tol tol,
  51. boost::uintmax_t& max_iter,
  52. const ``__Policy``&);
  53. // TOMS 748 algorithm
  54. template <class F, class T, class Tol>
  55. std::pair<T, T>
  56. toms748_solve(
  57. F f,
  58. const T& a,
  59. const T& b,
  60. Tol tol,
  61. boost::uintmax_t& max_iter);
  62. template <class F, class T, class Tol, class ``__Policy``>
  63. std::pair<T, T>
  64. toms748_solve(
  65. F f,
  66. const T& a,
  67. const T& b,
  68. Tol tol,
  69. boost::uintmax_t& max_iter,
  70. const ``__Policy``&);
  71. template <class F, class T, class Tol>
  72. std::pair<T, T>
  73. toms748_solve(
  74. F f,
  75. const T& a,
  76. const T& b,
  77. const T& fa,
  78. const T& fb,
  79. Tol tol,
  80. boost::uintmax_t& max_iter);
  81. template <class F, class T, class Tol, class ``__Policy``>
  82. std::pair<T, T>
  83. toms748_solve(
  84. F f,
  85. const T& a,
  86. const T& b,
  87. const T& fa,
  88. const T& fb,
  89. Tol tol,
  90. boost::uintmax_t& max_iter,
  91. const ``__Policy``&);
  92. // Termination conditions:
  93. template <class T>
  94. struct eps_tolerance;
  95. struct equal_floor;
  96. struct equal_ceil;
  97. struct equal_nearest_integer;
  98. }}} // boost::math::tools namespaces
  99. [h4 Description]
  100. These functions solve the root of some function ['f(x)] -
  101. ['without the need for any derivatives of ['f(x)]].
  102. The `bracket_and_solve_root` functions use __root_finding_TOMS748
  103. by Alefeld, Potra and Shi that is asymptotically the most efficient known,
  104. and has been shown to be optimal for a certain classes of smooth functions.
  105. Variants with and without __policy_section are provided.
  106. Alternatively, __bisect is a simple __bisection_wikipedia routine which can be useful
  107. in its own right in some situations, or alternatively for narrowing
  108. down the range containing the root, prior to calling a more advanced
  109. algorithm.
  110. All the algorithms in this section reduce the diameter of the enclosing
  111. interval with the same asymptotic efficiency with which they locate the
  112. root. This is in contrast to the derivative based methods which may ['never]
  113. significantly reduce the enclosing interval, even though they rapidly approach
  114. the root. This is also in contrast to some other derivative-free methods
  115. (for example, Brent's method described at
  116. [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker)]
  117. which only reduces the enclosing interval on the final step.
  118. Therefore these methods return a `std::pair` containing the enclosing interval found,
  119. and accept a function object specifying the termination condition.
  120. Three function objects are provided for ready-made termination conditions:
  121. * ['eps_tolerance] causes termination when the relative error in the enclosing
  122. interval is below a certain threshold.
  123. * ['equal_floor] and ['equal_ceil] are useful for certain statistical applications
  124. where the result is known to be an integer.
  125. * Other user-defined termination conditions are likely to be used
  126. only rarely, but may be useful in some specific circumstances.
  127. [section:bisect Bisection]
  128. template <class F, class T, class Tol>
  129. std::pair<T, T>
  130. bisect( // Unlimited iterations.
  131. F f,
  132. T min,
  133. T max,
  134. Tol tol);
  135. template <class F, class T, class Tol>
  136. std::pair<T, T>
  137. bisect( // Limited iterations.
  138. F f,
  139. T min,
  140. T max,
  141. Tol tol,
  142. boost::uintmax_t& max_iter);
  143. template <class F, class T, class Tol, class ``__Policy``>
  144. std::pair<T, T>
  145. bisect( // Specified policy.
  146. F f,
  147. T min,
  148. T max,
  149. Tol tol,
  150. boost::uintmax_t& max_iter,
  151. const ``__Policy``&);
  152. These functions locate the root using __bisection_wikipedia.
  153. `bisect` function arguments are:
  154. [variablelist
  155. [[f] [A unary functor (or C++ lambda) which is the function ['f(x)] whose root is to be found.]]
  156. [[min] [The left bracket of the interval known to contain the root.]]
  157. [[max] [The right bracket of the interval known to contain the root.[br]
  158. It is a precondition that ['min < max] and ['f(min)*f(max) <= 0],
  159. the function raises an __evaluation_error if these preconditions are violated.
  160. The action taken on error is controlled by the __Policy template argument: the default behavior is to
  161. throw a ['boost::math::evaluation_error]. If the __Policy is changed to not throw
  162. then it returns ['std::pair<T>(min, min)].]]
  163. [[tol] [A binary functor (or C++ lambda) that specifies the termination condition: the function
  164. will return the current brackets enclosing the root when ['tol(min, max)] becomes true.
  165. See also __root_termination.]]
  166. [[max_iter][The maximum number of invocations of ['f(x)] to make while searching for the root. On exit, this is updated to the actual number of invocations performed.]]
  167. ]
  168. [optional_policy]
  169. [*Returns]: a pair of values ['r] that bracket the root so that:
  170. [:f(r.first) * f(r.second) <= 0]
  171. and either
  172. [:tol(r.first, r.second) == true]
  173. or
  174. [:max_iter >= m]
  175. where ['m] is the initial value of ['max_iter] passed to the function.
  176. In other words, it's up to the caller to verify whether termination occurred
  177. as a result of exceeding ['max_iter] function invocations (easily done by
  178. checking the updated value of ['max_iter] when the function returns), rather than
  179. because the termination condition ['tol] was satisfied.
  180. [endsect] [/section:bisect Bisection]
  181. [section:bracket_solve Bracket and Solve Root]
  182. template <class F, class T, class Tol>
  183. std::pair<T, T>
  184. bracket_and_solve_root(
  185. F f,
  186. const T& guess,
  187. const T& factor,
  188. bool rising,
  189. Tol tol,
  190. boost::uintmax_t& max_iter);
  191. template <class F, class T, class Tol, class ``__Policy``>
  192. std::pair<T, T>
  193. bracket_and_solve_root(
  194. F f,
  195. const T& guess,
  196. const T& factor,
  197. bool rising,
  198. Tol tol,
  199. boost::uintmax_t& max_iter,
  200. const ``__Policy``&);
  201. `bracket_and_solve_root` is a convenience function that calls __root_finding_TOMS748 internally
  202. to find the root of ['f(x)]. It is generally much easier to use this function rather than __root_finding_TOMS748, since it
  203. does the hard work of bracketing the root for you. It's bracketing routines are quite robust and will
  204. usually be more foolproof than home-grown routines, unless the function can be analysed to yield tight
  205. brackets.
  206. Note that this routine can only be used when:
  207. * ['f(x)] is monotonic in the half of the real axis containing ['guess].
  208. * The value of the inital guess must have the same sign as the root: the function
  209. will ['never cross the origin] when searching for the root.
  210. * The location of the root should be known at least approximately,
  211. if the location of the root differs by many orders of magnitude
  212. from ['guess] then many iterations will be needed to bracket the root in spite of
  213. the special heuristics used to guard against this very situation. A typical example would be
  214. setting the initial guess to 0.1, when the root is at 1e-300.
  215. The `bracket_and_solve_root` parameters are:
  216. [variablelist
  217. [[f][A unary functor (or C++ lambda) that is the function whose root is to be solved.
  218. ['f(x)] must be uniformly increasing or decreasing on ['x].]]
  219. [[guess][An initial approximation to the root.]]
  220. [[factor][A scaling factor that is used to bracket the root: the value
  221. /guess/ is multiplied (or divided as appropriate) by /factor/
  222. until two values are found that bracket the root. A value
  223. such as 2 is a typical choice for ['factor].
  224. In addition ['factor] will be multiplied by 2 every 32 iterations:
  225. this is to guard against a really very bad initial guess, typically these occur
  226. when it's known the result is very large or small, but not the exact order
  227. of magnitude.]]
  228. [[rising][Set to ['true] if ['f(x)] is rising on /x/ and /false/ if ['f(x)]
  229. is falling on /x/. This value is used along with the result
  230. of /f(guess)/ to determine if /guess/ is
  231. above or below the root.]]
  232. [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
  233. for the root. /tol/ is passed the current brackets at each step,
  234. when it returns true then the current brackets are returned as the pair result.
  235. See also __root_termination.]]
  236. [[max_iter] [The maximum number of function invocations to perform in the search
  237. for the root. On exit is set to the actual number of invocations performed.]]
  238. ]
  239. [optional_policy]
  240. [*Returns]: a pair of values ['r] that bracket the root so that:
  241. [:f(r.first) * f(r.second) <= 0]
  242. and either
  243. [:tol(r.first, r.second) == true]
  244. or
  245. [:max_iter >= m]
  246. where ['m] is the initial value of ['max_iter] passed to the function.
  247. In other words, it's up to the caller to verify whether termination occurred
  248. as a result of exceeding ['max_iter] function invocations (easily done by
  249. checking the value of ['max_iter] when the function returns), rather than
  250. because the termination condition ['tol] was satisfied.
  251. [endsect] [/section:bracket_solve Bracket and Solve Root]
  252. [section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
  253. template <class F, class T, class Tol>
  254. std::pair<T, T>
  255. toms748_solve(
  256. F f,
  257. const T& a,
  258. const T& b,
  259. Tol tol,
  260. boost::uintmax_t& max_iter);
  261. template <class F, class T, class Tol, class ``__Policy``>
  262. std::pair<T, T>
  263. toms748_solve(
  264. F f,
  265. const T& a,
  266. const T& b,
  267. Tol tol,
  268. boost::uintmax_t& max_iter,
  269. const ``__Policy``&);
  270. template <class F, class T, class Tol>
  271. std::pair<T, T>
  272. toms748_solve(
  273. F f,
  274. const T& a,
  275. const T& b,
  276. const T& fa,
  277. const T& fb,
  278. Tol tol,
  279. boost::uintmax_t& max_iter);
  280. template <class F, class T, class Tol, class ``__Policy``>
  281. std::pair<T, T>
  282. toms748_solve(
  283. F f,
  284. const T& a,
  285. const T& b,
  286. const T& fa,
  287. const T& fb,
  288. Tol tol,
  289. boost::uintmax_t& max_iter,
  290. const ``__Policy``&);
  291. These functions implement TOMS Algorithm 748: it uses a mixture of
  292. cubic, quadratic and linear (secant) interpolation to locate the root of
  293. ['f(x)]. The two pairs of functions differ only by whether values for ['f(a)] and
  294. ['f(b)] are already available.
  295. Generally speaking it is easier (and often more efficient) to use __bracket_solve
  296. rather than trying to bracket the root yourself as this function requires.
  297. This function is provided rather than [@http://en.wikipedia.org/wiki/Brent%27s_method Brent's method] as it is known to be more
  298. efficient in many cases (it is asymptotically the most efficient known,
  299. and has been shown to be optimal for a certain classes of smooth functions).
  300. It also has the useful property of decreasing the bracket size
  301. with each step, unlike Brent's method which only shrinks the enclosing interval in the
  302. final step. This makes it particularly useful when you need a result where the ends
  303. of the interval round to the same integer: as often happens in statistical applications,
  304. for example. In this situation the function is able to exit after a much smaller
  305. number of iterations than would otherwise be possible.
  306. The __root_finding_TOMS748 parameters are:
  307. [variablelist
  308. [[f] [A unary functor (or C++ lambda) that is the function whose root is to be solved.
  309. f(x) need not be uniformly increasing or decreasing on ['x] and
  310. may have multiple roots. However, the bounds given must bracket a single root.]]
  311. [[a] [The lower bound for the initial bracket of the root.]]
  312. [[b] [The upper bound for the initial bracket of the root.
  313. It is a precondition that ['a < b] and that ['a] and ['b]
  314. bracket the root to find so that ['f(a) * f(b) < 0].]]
  315. [[fa] [Optional: the value of ['f(a)].]]
  316. [[fb] [Optional: the value of ['f(b)].]]
  317. [[tol] [A binary functor (or C++ lambda) that determines the termination condition for the search
  318. for the root. ['tol] is passed the current brackets at each step,
  319. when it returns true, then the current brackets are returned as the result.
  320. See also __root_termination.]]
  321. [[max_iter] [The maximum number of function invocations to perform in the search
  322. for the root. On exit, ['max_iter] is set to actual number of function
  323. invocations used.]]
  324. ]
  325. [optional_policy]
  326. `toms748_solve` returns: a pair of values ['r] that bracket the root so that:
  327. [:['f(r.first) * f(r.second) <= 0]]
  328. and either
  329. [:['tol(r.first, r.second) == true]]
  330. or
  331. [:['max_iter >= m]]
  332. where ['m] is the initial value of ['max_iter] passed to the function.
  333. In other words, it's up to the caller to verify whether termination occurred
  334. as a result of exceeding ['max_iter] function invocations (easily done by
  335. checking the updated value of ['max_iter]
  336. against its previous value passed as parameter),
  337. rather than because the termination condition ['tol] was satisfied.
  338. [endsect] [/section:TOMS748 Algorithm TOMS 748: Alefeld, Potra and Shi: Enclosing zeros of continuous functions]
  339. [section:brent Brent-Decker Algorithm]
  340. The [@http://en.wikipedia.org/wiki/Brent%27s_method Brent-Dekker algorithm], although very well know,
  341. is not provided by this library as __root_finding_TOMS748 or
  342. its slightly easier to use variant __bracket_solve are superior and provide equivalent functionality.
  343. [endsect] [/section:brent Brent-Decker Algorithm]
  344. [section:root_termination Termination Condition Functors]
  345. template <class T>
  346. struct eps_tolerance
  347. {
  348. eps_tolerance();
  349. eps_tolerance(int bits);
  350. bool operator()(const T& a, const T& b)const;
  351. };
  352. `eps_tolerance` is the usual termination condition used with these root finding functions.
  353. Its `operator()` will return true when the relative distance between ['a] and ['b]
  354. is less than four times the machine epsilon for T, or 2[super 1-bits], whichever is
  355. the larger. In other words, you set ['bits] to the number of bits of precision you
  356. want in the result. The minimal tolerance of ['four times the machine epsilon of type T] is
  357. required to ensure that we get back a bracketing interval, since this must clearly
  358. be at greater than one epsilon in size. While in theory a maximum distance of twice
  359. machine epsilon is possible to achieve, in practice this results in a great deal of "thrashing"
  360. given that the function whose root is being found can only ever be accurate to 1 epsilon at best.
  361. struct equal_floor
  362. {
  363. equal_floor();
  364. template <class T> bool operator()(const T& a, const T& b)const;
  365. };
  366. This termination condition is used when you want to find an integer result
  367. that is the ['floor] of the true root. It will terminate as soon as both ends
  368. of the interval have the same ['floor].
  369. struct equal_ceil
  370. {
  371. equal_ceil();
  372. template <class T> bool operator()(const T& a, const T& b)const;
  373. };
  374. This termination condition is used when you want to find an integer result
  375. that is the ['ceil] of the true root. It will terminate as soon as both ends
  376. of the interval have the same ['ceil].
  377. struct equal_nearest_integer
  378. {
  379. equal_nearest_integer();
  380. template <class T> bool operator()(const T& a, const T& b)const;
  381. };
  382. This termination condition is used when you want to find an integer result
  383. that is the /closest/ to the true root. It will terminate as soon as both ends
  384. of the interval round to the same nearest integer.
  385. [endsect] [/section:root_termination Termination Condition Functors]
  386. [section:implementation Implementation]
  387. The implementation of the bisection algorithm is extremely straightforward
  388. and not detailed here.
  389. __TOMS748 is described in detail in:
  390. ['Algorithm 748: Enclosing Zeros of Continuous Functions,
  391. G. E. Alefeld, F. A. Potra and Yixun Shi,
  392. ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995.
  393. Pages 327-344.]
  394. The implementation here is a faithful translation of this paper into C++.
  395. [endsect] [/section:implementation Implementation]
  396. [endsect] [/section:roots_noderiv Root Finding Without Derivatives]
  397. [/
  398. Copyright 2006, 2010, 2015 John Maddock and Paul A. Bristow.
  399. Distributed under the Boost Software License, Version 1.0.
  400. (See accompanying file LICENSE_1_0.txt or copy at
  401. http://www.boost.org/LICENSE_1_0.txt).
  402. ]