//
//
//
//
// Lab. Calcolo II - Esempio di codice
//
//
//
//
//
//
// <![CDATA[
// Example code: Function zero search with Newton-Raphson method.
// (francesco.prelz@mi.infn.it 20060814)
#include <iostream>
#include <iomanip>
#include <vector>
#include <functional>
#include <cmath>
template <typename DT>
class FunctionZeroNRSearch
{
public:
struct bracket
{
DT x1;
DT x2;
};
typedef DT value_type;
typedef DT (*function_ptr_type)(DT); // Normal - nonmember function pointer
typedef std::function<DT (DT arg)> functor_type;
typedef std::vector<bracket> bracket_container_type;
// Creators
FunctionZeroNRSearch(functor_type &f, functor_type &df);
FunctionZeroNRSearch(function_ptr_type f, function_ptr_type df);
FunctionZeroNRSearch(functor_type &f, function_ptr_type df);
FunctionZeroNRSearch(function_ptr_type f, functor_type &df);
// Destructor
~FunctionZeroNRSearch();
DT operator() (DT argument);
int find_bracket_expand(bracket &bra);
int find_bracket_partition(const bracket &bra, int n_partitions);
void clear_brackets();
const bracket_container_type &get_brackets();
std::vector<DT> get_zeros(double tolerance=1e-6);
private:
functor_type m_f;
functor_type m_df;
bracket_container_type m_brackets;
};
template <typename DT>
FunctionZeroNRSearch<DT>::FunctionZeroNRSearch(functor_type &f, functor_type &df) : m_f(f), m_df(df) {}
template <typename DT>
FunctionZeroNRSearch<DT>::FunctionZeroNRSearch(function_ptr_type f, function_ptr_type df) : m_f(f), m_df(df) {}
template <typename DT>
FunctionZeroNRSearch<DT>::FunctionZeroNRSearch(functor_type &f, function_ptr_type df) : m_f(f), m_df(df) {}
template <typename DT>
FunctionZeroNRSearch<DT>::FunctionZeroNRSearch(function_ptr_type f, functor_type &df) : m_f(f), m_df(df) {}
template <typename DT>
FunctionZeroNRSearch<DT>::~FunctionZeroNRSearch() {}
template <typename DT>
DT FunctionZeroNRSearch<DT>::operator() (DT argument)
{
return m_f(argument);
}
template <typename DT>
int FunctionZeroNRSearch<DT>::find_bracket_expand(bracket &bra)
{
m_brackets.clear();
const DT expansion_factor=1.6;
const DT expansion_max_retries=50;
for (int i=0;i<expansion_max_retries;i++)
{
DT fx1 = m_f(bra.x1);
DT fx2 = m_f(bra.x2);
if (fx1*fx2 < 0)
{
m_brackets.push_back(bra);
return 1;
}
if (std::fabs(fx1) < std::fabs(fx2))
{
// Try expanding the interval at the x1 side
bra.x1 = bra.x2 + (bra.x1-bra.x2) * expansion_factor;
}
else
{
// Try expanding the interval at the x2 side
bra.x2 = bra.x1 + (bra.x2-bra.x1) * expansion_factor;
}
}
return 0; // 0 means no good brackets were found
}
template <typename DT>
int FunctionZeroNRSearch<DT>::find_bracket_partition(const bracket &bra, int n_partitions)
{
m_brackets.clear();
int n_results = 0;
DT cur_x1,cur_x2;
DT segment_size = (bra.x2 - bra.x1) / static_cast<DT>(n_partitions);
cur_x1 = bra.x1;
cur_x2 = cur_x1 + segment_size;
for (int i=0; i<n_partitions; i++)
{
// Last segment
if (i == (n_partitions-1)) cur_x2 = bra.x2;
if (m_f(cur_x1)*m_f(cur_x2) < 0)
{
bracket new_result;
new_result.x1 = cur_x1;
new_result.x2 = cur_x2;
m_brackets.push_back(new_result);
n_results++;
}
cur_x1 += segment_size;
cur_x2 += segment_size;
}
return n_results;
}
template <typename DT>
void FunctionZeroNRSearch<DT>::clear_brackets()
{
m_brackets.clear();
}
template <typename DT>
const typename FunctionZeroNRSearch<DT>::bracket_container_type &
FunctionZeroNRSearch<DT>::get_brackets()
{
return m_brackets;
}
template <typename DT>
std::vector<DT> FunctionZeroNRSearch<DT>::get_zeros(double tolerance)
{
const int maximum_iterations=100;
std::vector<DT> results;
for (typename FunctionZeroNRSearch<DT>::bracket_container_type::const_iterator bra = m_brackets.begin();
bra != m_brackets.end(); ++bra)
{
DT xmiddle, newx;
DT fleft, fmiddle, fright;
DT increment,old_increment;
DT dfmiddle;
// Computing f might be expensive
fleft = m_f(bra->x1);
fright = m_f(bra->x2);
if (fleft*fright >= 0) continue; // Not a valid bracket.
xmiddle = (bra->x1 + bra->x2) / 2;
fmiddle = m_f(xmiddle);
old_increment = std::fabs(bra->x1 - bra->x2);
increment = old_increment;
bracket cbra(*bra); // Locally modifiable copy.
for (int i=0; i<maximum_iterations; i++)
{
if (fmiddle == 0)
{
results.push_back(xmiddle); // We cannot expect
break; // to be this lucky too often.
}
// Compute first derivative at xmiddle
dfmiddle = m_df(xmiddle);
if (dfmiddle != 0) newx = xmiddle - fmiddle/dfmiddle;
else newx = xmiddle;
// If we are trying to move out of the bracket boundaries
// (that is, if the tangent line has the same sign at the boundaries)
// or if the tangent line does not decrease at least by twice
// as the function value over the current interval, revert to
// normal bisection.
if ( ((xmiddle - cbra.x1)*dfmiddle - fmiddle) *
((xmiddle - cbra.x2)*dfmiddle - fmiddle) > 0 ||
(std::fabs(2*fmiddle) > std::fabs(old_increment*dfmiddle)) )
{
old_increment = increment;
if (fleft * fmiddle < 0)
{
cbra.x2 = xmiddle;
fright = fmiddle;
}
else
{
cbra.x1 = xmiddle;
fleft = fmiddle;
}
xmiddle = (cbra.x1 + cbra.x2) / 2;
increment = std::fabs(cbra.x1 - cbra.x2) / 2;
fmiddle = m_f(xmiddle);
}
else
{
old_increment = increment;
increment = std::fabs(newx - xmiddle);
xmiddle = newx;
fmiddle = m_f(xmiddle);
if (fleft * fmiddle < 0)
{
cbra.x2 = xmiddle;
fright = fmiddle;
}
else
{
cbra.x1 = xmiddle;
fleft = fmiddle;
}
}
if (std::fabs(cbra.x2 - cbra.x1) < std::fabs(tolerance))
{
results.push_back(xmiddle);
break;
}
}
}
// We get here only if we received a bad tolerance value
return results;
}
int
main(int argc, char *argv[])
{
FunctionZeroNRSearch<double> fzsin(std::sin,std::cos);
FunctionZeroNRSearch<double>::bracket test;
test.x1 = -M_PI/2;
test.x2 = M_PI/2;
int result = fzsin.find_bracket_expand(test);
std::cout << "find_bracket_expand sin with initial interval (0,PI/4) returns "
<< result << std::endl;
std::cout << " Bracket: [" << test.x1/M_PI << " PI," << test.x2/M_PI
<< " PI]" << std::endl;
std::vector<double> zeros = fzsin.get_zeros();
for (std::vector<double>::const_iterator it = zeros.begin();
it != zeros.end(); ++it)
{
std::cout << " zero = " << (*it)/M_PI << " PI (+/- 1E-6)" << std::endl;
}
test.x1 = M_PI/2;
test.x2 = 9*M_PI/2;
result = fzsin.find_bracket_partition(test, 20);
std::cout << "get_bracket_partition sin with initial interval (0,4*PI) returns "
<< result << std::endl;
const FunctionZeroNRSearch<double>::bracket_container_type &brackets =
fzsin.get_brackets();
FunctionZeroNRSearch<double>::bracket_container_type::const_iterator bra;
for (bra=brackets.begin();bra!=brackets.end();++bra)
{
std::cout << " Bracket: [" << bra->x1/M_PI << " PI," << bra->x2/M_PI
<< " PI]" << std::endl;
}
zeros = fzsin.get_zeros();
for (std::vector<double>::const_iterator it = zeros.begin();
it != zeros.end(); ++it)
{
std::cout << " zero = " << (*it)/M_PI << " PI (+/- 1E-6)" << std::endl;
}
return 0;
}
]]></plaintext></span>