regula_falsi.h

Go to the documentation of this file.

/*
 *  regula_falsi.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef REGULA_FALSI_H
#define REGULA_FALSI_H
 
 
namespace nest
{
 
template < typename CN >
double
regula_falsi( const CN& node, const double dt )
{
  double root;
  double threshold_dist_root;
 
  int last_threshold_sign = 0;
 
  double a_k = 0.0;
  double b_k = dt;
 
  double threshold_dist_a_k = node.threshold_distance( a_k );
  double threshold_dist_b_k = node.threshold_distance( b_k );
 
  if ( threshold_dist_a_k * threshold_dist_b_k > 0 )
  {
    throw NumericalInstability( "regula_falsi: time step too short to reach threshold." );
  }
 
  const int MAX_ITER = 500;
  const double TERMINATION_CRITERION = 1e-14;
 
  for ( int iter = 0; iter < MAX_ITER; ++iter )
  {
    assert( threshold_dist_b_k != threshold_dist_a_k );
 
    root = ( a_k * threshold_dist_b_k - b_k * threshold_dist_a_k ) / ( threshold_dist_b_k - threshold_dist_a_k );
    threshold_dist_root = node.threshold_distance( root );
 
    if ( std::abs( threshold_dist_root ) < TERMINATION_CRITERION )
    {
      return root;
    }
 
    if ( threshold_dist_a_k * threshold_dist_root > 0.0 )
    {
      // threshold_dist_a_k and threshold_dist_root have the same sign
      a_k = root;
      threshold_dist_a_k = threshold_dist_root;
 
      if ( last_threshold_sign == 1 )
      {
        // If threshold_dist_a_k and threshold_dist_root had the same sign in
        // the last time step, we halve the value of threshold_distance_(b_k) to
        // force the root in the next time step to occur on b_k's side. This is
        // done to improve the convergence rate.
        threshold_dist_b_k /= 2;
      }
      last_threshold_sign = 1;
    }
    else if ( threshold_dist_b_k * threshold_dist_root > 0.0 )
    {
      // threshold_dist_b_k and threshold_dist_root have the same sign
      b_k = root;
      threshold_dist_b_k = threshold_dist_root;
 
      if ( last_threshold_sign == -1 )
      {
        threshold_dist_a_k /= 2;
      }
      last_threshold_sign = -1;
    }
    else
    {
      throw NumericalInstability( "regula_falsi: Regula falsi method did not converge" );
    }
  }
  throw NumericalInstability( "regula_falsi: Regula falsi method did not converge during set number of iterations" );
}
 
}  // namespace nest
 
#endif  // REGULA_FALSI_H