home *** CD-ROM | disk | FTP | other *** search
- //===================================================================
- // rootfind.hpp
- //
- // Version 1.0
- //
- // Written by:
- // Brent Worden
- // WordenWare
- // email: Brent@Worden.org
- //
- // Copyright (c) 1999 WordenWare
- //
- // Created: May 09, 1999
- // Revised:
- //===================================================================
-
- #ifndef _ROOTFIND_HPP_
- #define _ROOTFIND_HPP_
-
- #include "numerics.h"
-
- NUM_BEGIN
-
- template<class Function, class T>
- T bisection(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
- //-------------------------------------------------------------------
- // Returns the root of func bracketed between x0 and x1 using the
- // bisection method. The returned root is within tol of the true
- // value of the root.
- //-------------------------------------------------------------------
- {
- T f0 = func(x0);
- T f1 = func(x1);
- T zero(0);
-
- if(f0 * f1 > zero){
- throw Exception("bisection", "func(x0) and func(x1) must have oppisite signs");
- }
-
- T fmid;
- T xmid;
- while(!areEqual(x0, x1, tol)){
- xmid = average(x0, x1);
- fmid = func(xmid);
- if(fmid * f0 > zero){
- f0 = fmid;
- x0 = xmid;
- } else {
- f1 = fmid;
- x1 = xmid;
- }
- }
- return average(x0, x1);
- };
-
- template<class Function, class T>
- T falsePosition(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
- //-------------------------------------------------------------------
- // Returns the root of func using the false position method with
- // initial approximations x0 and x1. The returned root is within tol
- // of the true value of the root.
- //-------------------------------------------------------------------
- {
- T f0 = func(x0);
- T f1 = func(x1);
- T zero(0);
-
- if(f0 * f1 > zero){
- throw Exception("falsePosition", "func(x0) and func(x1) must have oppisite signs");
- }
-
- T f;
- T x;
- int n = 0;
-
- do {
- x = x1 - f1 * (x1 - x0) / (f1 - f0);
- f = func(x);
- if(f * f1 < zero){
- x0 = x1;
- f0 = f1;
- }
- x1 = x;
- f1 = f;
- ++n;
- } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
- if(n >= NUMERICS_ITMAX){
- throw Exception("falsePosition", "Failed to converge to root.");
- }
-
- return x1;
- };
-
- template<class Function, class T>
- T newton(Function func, Function der, T x0, T tol = NUMERICS_MAX_ERROR)
- //-------------------------------------------------------------------
- // Returns the root of func with first derivative der using Newton's
- // method with initial approximation x0. The returned root is within
- // tol of the true value of the root.
- //-------------------------------------------------------------------
- {
- int n = 0;
- T x1 = x0;
-
- do {
- x0 = x1;
- x1 = x0 - func(x0) / der(x0);
- ++n;
- } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
- if(n >= NUMERICS_ITMAX){
- throw Exception("newton", "Failed to converge to root.");
- }
-
- return x1;
- };
-
- template<class Function, class T>
- T secant(Function func, T x0, T x1, T tol = NUMERICS_MAX_ERROR)
- //-------------------------------------------------------------------
- // Returns the root of func using the secant method with initial
- // approximations x0 and x1. The returned root is within tol of the
- // true value of the root.
- //-------------------------------------------------------------------
- {
- T f0 = func(x0);
- T f1 = func(x1);
- int n = 0;
- T x;
-
- do {
- x = x1 - f1 * (x1 - x0) / (f1 - f0);
- x0 = x1;
- f0 = f1;
- x1 = x;
- f1 = func(x);
- ++n;
- } while(n < NUMERICS_ITMAX && !areEqual(x0, x1, tol));
- if(n >= NUMERICS_ITMAX){
- throw Exception("secant", "Failed to converge to root.");
- }
-
- return x1;
- };
-
- NUM_END
-
- #endif
-
- //===================================================================
- // Revision History
- //
- // Version 1.0 - 05/09/1999 - New.
- //===================================================================
-