+ View the NASA Portal

Search JPL

Search CLARAty

- Algorithms- Infrastructure - Doxygen (public)- Doxygen (private)- Tutorials

- Repository- Nightly Build- Nightly Regression- Feature/Bug Tracking- Build Targets- Procedure- Conventions- Tools

1d_function_solver.cc

Go to the documentation of this file.

// -*-c++-*-
//---------------------------< /-/ CLARAty /-/ >------------------------------
/**
 * @file  1d_function_solver.cc
 *
 * Provides algorithms to (a) find the roots for a single variable
 * math function (i.e. zero crossings), and (b) the local
 * minimum/maximum of a single variable function over a specified
 * interval. (a) roots is implemented using a Newton Raphson
 * minimization iterative technique.
 *
 * <br>@b Designer(s):  Hari Nayar
 * <br>@b Author(s):    Issa A.D. Nesnas
 * <br>@b Date:         May 15, 2001
 *
 * <b>Software License:</b><br>
 * <code>  http://claraty.jpl.nasa.gov/license/open_src/  or
 *         file: license/open_src.txt  </code>
 *
 * &copy; 2006, Jet Propulsion Laboratory, California Institute of Technology<br>
 *
 * $Revision: 1.4 $
 */
//-----------------------------------------------------------------------------

#include "claraty/1d_function_solver.h"
#include <iostream>
#include <math.h>

using namespace std;

namespace claraty {

//-----------------------------------------------------------------------------
/**
 * 1-D function solver.
 *
 * @param[in]   function  Function which accepts double arguments
 * @param[in]   guess
 * @param[in]   max_iterations
 * @param[in]   tolerance
 * @param[in]   perturbation
 * @param[in]   step_size
 */

N_1D_Function_Solver::
N_1D_Function_Solver(double (*function)(double x),
                      double guess,
                      int    max_iterations,
                      double tolerance,
                      double perturbation,
                      double step_size)
  :
  _function1(NULL)
{ 
  set_function(function);
  set_initial_guess(guess);
  set_maximum_iterations(max_iterations);
  set_solution_tolerance(tolerance);
  set_perturbation(perturbation);
  set_step_size(step_size);
}

//-----------------------------------------------------------------------------
/**
 * 1-D generic function solver.
 *
 * @param[in]   function  Function which accepts <> arguments.
 * @param[in]   guess
 * @param[in]   max_iterations
 * @param[in]   tolerance
 * @param[in]   perturbation
 * @param[in]   step_size
 */
N_1D_Function_Solver::
N_1D_Function_Solver(Function<> * function,
                     double       guess,
                     int          max_iterations,
                     double       tolerance,
                     double       perturbation,
                     double       step_size)
  :
  _function2(NULL)
{
  set_function(function);
  set_initial_guess(guess);
  set_maximum_iterations(max_iterations);
  set_solution_tolerance(tolerance);
  set_perturbation(perturbation);
  set_step_size(step_size);
}

//-----------------------------------------------------------------------------
/**
 * Find root of a function.
 *
 * @param [out] result
 * @return      Return code.
 */
int N_1D_Function_Solver::
find_root(double& result)
{
  int  iterations = 0;
  bool terminate  = false;
  int return_code = CONTINUE_ITERATION;
  double current_value, next_value=0, delta_value;
  double derivative;
  
  while (!terminate)
    {
      // Compute the function value at the current guess
      current_value = _function(_guess);
      
      // If the current value is below the threshold,
      // set termination condition and quit the loop
      
      if (fabs(current_value) < _error_tolerance)
        {
          terminate = true;
          return_code = SOLVED;
        }
      else
        {
          // Compute the function value at a delta perturbation from
          // the current guess to determine the derivative of the function
          next_value = _function(_guess+_perturbation);
          delta_value = next_value - current_value;
            
          // Compute the derivative 
          derivative = delta_value/_perturbation;
          
          // Use the derivative to determine the next guess
          _guess -= current_value / derivative;
          
          // Increment iteration counter
          iterations++;
          
          // Check for termination
          if (iterations > _max_iterations)
            {
              terminate = true;
              return_code = MAX_ITERATIONS_EXCEEDED;
            }
        }
#ifdef DEBUG
      cout << "Iteration #: "      << iterations
           << "  Current guess = " << _guess 
           << "  Current value = " << current_value
           << "  Next value = "    << next_value;
           << endl;
#endif // DEBUG
    }
    result = _guess;
    return return_code;
}

//-----------------------------------------------------------------------------
/**
 * Bisection.
 *
 * @param[out]  x0
 * @param[out]  x2
 * @param[out]  x4
 * @param[out]  y0
 * @param[out]  y2
 * @param[out]  y4
 */
void N_1D_Function_Solver::
_bisect(double &x0, double &y0,
        double &x2, double &y2,
        double &x4, double &y4)
{
  // Bisect the first interval
  double x1 = ( x0 + x2 )/2;
  double y1 = _function(x1);
  
  // if new point is lower, use as new middle
  if(y1 < y2)
    {
      x4 = x2; y4 = y2; x2 = x1; y2 = y1;
      return;
    }
    
  // if it's higher, move in left limit
  if(y1 > y2)
    {
      x0 = x1; y0 = y1;
    }
    
  // Bisect the second interval
  double x3 = ( x2 + x4 )/2;
  double y3 = _function(x3);
    
  // if new point is lower, use as new middle
  if(y3 < y2)
    {
      x0 = x2; y0 = y2; x2 = x3; y2 = y3;
      return;
    }
  
  // if it's higher, move in right limit
  if(y3 > y2)
    {
      x4 = x3; y4 = y3;
    }
  
  // if y1 was equal to y2, replace y2
  if(y1 == y2)
    {
      x2 = x1; y2 = y1;
    }
}

/*
 * Solve for minimum of function, using n successive bisections.
 * Initial points _MUST_ represent legal  configuration, i.e.,
 * y0 > y1, y2 > y1.
 * The points need not be evenly spaced.
 *
 * Example of usage of the minimization function:
 *  double x0 = 0;
 *  double y0 = f(x0);
 *  double x1 = 0.5;
 *  double y1 = f(x1);
 *  double x2 = 1;
 *  double y2 = f(x2);
 *  cout << minimize(x0, y0, x1, y1, x2, y2, 20);
 */

int N_1D_Function_Solver::
minimize(double &x0, double &y0,
         double &x1, double &y1,
         double &x2, double &y2,
         long n, double &result)
{
  for(long i=0; i<n; i++) _bisect(x0, y0, x1, y1, x2, y2);
  
  result = x1; 
  return SOLVED;
}

} // namespace claraty

Last modified: June 05 2007 10:19:53

JPL PRIVACY STATEMENT GLOSSARY SITEMAP FEEDBACK CONTACT US

Site Manager: Issa A.D. Nesnas
Site Maintainer: Kelly Breed
Document Clearance #: CL 07-0635