home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / rootfind.hpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-17  |  3.8 KB  |  154 lines

  1. //===================================================================
  2. // rootfind.hpp
  3. //
  4. // Version 1.0
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1999 WordenWare
  12. //
  13. // Created:  May 09, 1999
  14. // Revised:  
  15. //===================================================================
  16.  
  17. #ifndef _ROOTFIND_HPP_
  18. #define _ROOTFIND_HPP_
  19.  
  20. #include "numerics.h"
  21.  
  22. NUM_BEGIN
  23.  
  24. template<class Function, class T>
  25. T bisection(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
  26. //-------------------------------------------------------------------
  27. // Returns the root of func bracketed between x0 and x1 using the
  28. // bisection method.  The returned root is within tol of the true
  29. // value of the root.
  30. //-------------------------------------------------------------------
  31. {
  32.     T f0 = func(x0);
  33.     T f1 = func(x1);
  34.     T zero(0);
  35.  
  36.     if(f0 * f1 > zero){
  37.         throw Exception("bisection", "func(x0) and func(x1) must have oppisite signs");
  38.     }
  39.  
  40.     T fmid;
  41.     T xmid;
  42.     while(!areEqual(x0, x1, tol)){
  43.         xmid = average(x0, x1);
  44.         fmid = func(xmid);
  45.         if(fmid * f0 > zero){
  46.             f0 = fmid;
  47.             x0 = xmid;
  48.         } else {
  49.             f1 = fmid;
  50.             x1 = xmid;
  51.         }
  52.     }
  53.     return average(x0, x1);
  54. };
  55.  
  56. template<class Function, class T>
  57. T falsePosition(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
  58. //-------------------------------------------------------------------
  59. // Returns the root of func using the false position method with
  60. // initial approximations x0 and x1.  The returned root is within tol
  61. // of the true value of the root.
  62. //-------------------------------------------------------------------
  63. {
  64.     T f0 = func(x0);
  65.     T f1 = func(x1);
  66.     T zero(0);
  67.  
  68.     if(f0 * f1 > zero){
  69.         throw Exception("falsePosition", "func(x0) and func(x1) must have oppisite signs");
  70.     }
  71.  
  72.     T f;
  73.     T x;
  74.     int n = 0;
  75.  
  76.     do {
  77.         x = x1 - f1 * (x1 - x0) / (f1 - f0);
  78.         f = func(x);
  79.         if(f * f1 < zero){
  80.             x0 = x1;
  81.             f0 = f1;
  82.         }
  83.         x1 = x;
  84.         f1 = f;
  85.         ++n;
  86.     } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
  87.     if(n >= NUMERICS_ITMAX){
  88.         throw Exception("falsePosition", "Failed to converge to root.");
  89.     }
  90.  
  91.     return x1;
  92. };
  93.  
  94. template<class Function, class T>
  95. T newton(Function func, Function der, T x0, T tol = NUMERICS_MAX_ERROR)
  96. //-------------------------------------------------------------------
  97. // Returns the root of func with first derivative der using Newton's
  98. // method with initial approximation x0.  The returned root is within
  99. // tol of the true value of the root.
  100. //-------------------------------------------------------------------
  101. {
  102.     int n = 0;
  103.     T x1 = x0;
  104.  
  105.     do {
  106.         x0 = x1;
  107.         x1 = x0 - func(x0) / der(x0);
  108.         ++n;
  109.     } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
  110.     if(n >= NUMERICS_ITMAX){
  111.         throw Exception("newton", "Failed to converge to root.");
  112.     }
  113.  
  114.     return x1;
  115. };
  116.  
  117. template<class Function, class T>
  118. T secant(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
  119. //-------------------------------------------------------------------
  120. // Returns the root of func using the secant method with initial
  121. // approximations x0 and x1.  The returned root is within tol of the
  122. // true value of the root.
  123. //-------------------------------------------------------------------
  124. {
  125.     T f0 = func(x0);
  126.     T f1 = func(x1);
  127.     int n = 0;
  128.     T x;
  129.  
  130.     do {
  131.         x = x1 - f1 * (x1 - x0) / (f1 - f0);
  132.         x0 = x1;
  133.         f0 = f1;
  134.         x1 = x;
  135.         f1 = func(x);
  136.         ++n;
  137.     } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
  138.     if(n >= NUMERICS_ITMAX){
  139.         throw Exception("secant", "Failed to converge to root.");
  140.     }
  141.  
  142.     return x1;
  143. };
  144.  
  145. NUM_END
  146.  
  147. #endif
  148.  
  149. //===================================================================
  150. // Revision History
  151. //
  152. // Version 1.0 - 05/09/1999 - New.
  153. //===================================================================
  154.