//
//
//
//
// Lab. Calcolo II - Esempio di codice
//
//
//
//
//
//
// <![CDATA[
// Example code: Function zero search.
// (francesco.prelz@mi.infn.it 20060814)
#include <iostream>
#include <iomanip>
#include <vector>
#include <functional>
#include <cmath>
template <typename DT>
class FunctionZeroSearch
{
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;
// Creator
FunctionZeroSearch(functor_type &f);
FunctionZeroSearch(function_ptr_type f);
// Destructor
~FunctionZeroSearch();
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;
bracket_container_type m_brackets;
};
template <typename DT>
FunctionZeroSearch<DT>::FunctionZeroSearch(functor_type &f) : m_f(f) {}
template <typename DT>
FunctionZeroSearch<DT>::FunctionZeroSearch(function_ptr_type f) : m_f(f) {}
template <typename DT>
FunctionZeroSearch<DT>::~FunctionZeroSearch() {}
template <typename DT>
DT FunctionZeroSearch<DT>::operator() (DT argument)
{
return m_f(argument);
}
template <typename DT>
int FunctionZeroSearch<DT>::find_bracket_expand(bracket &bra)
{
m_brackets.clear();
const DT expansion_factor=1.6;
const int 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 FunctionZeroSearch<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 FunctionZeroSearch<DT>::clear_brackets()
{
m_brackets.clear();
}
template <typename DT>
const typename FunctionZeroSearch<DT>::bracket_container_type &
FunctionZeroSearch<DT>::get_brackets()
{
return m_brackets;
}
template <typename DT>
std::vector<DT> FunctionZeroSearch<DT>::get_zeros(double tolerance)
{
const int maximum_bisections=100;
std::vector<DT> results;
for (typename FunctionZeroSearch<DT>::bracket_container_type::const_iterator bra = m_brackets.begin();
bra != m_brackets.end(); ++bra)
{
double xmiddle;
double fleft, fmiddle, fright;
// Computing f might be expensive
fleft = m_f(bra->x1);
fright = m_f(bra->x2);
if (fleft*fright >= 0) continue; // Not a valid bracket.
bracket cbra(*bra); // Modifiable copy.
for (int i=0; i<maximum_bisections; i++)
{
xmiddle = (cbra.x1 + cbra.x2) / 2;
fmiddle = m_f(xmiddle);
if (fmiddle == 0)
{
results.push_back(xmiddle); // We cannot expect
break; // to be this lucky too often.
}
if (std::fabs(cbra.x2 - cbra.x1) < std::fabs(tolerance))
{
results.push_back(xmiddle);
break;
}
if (fleft * fmiddle < 0)
{
cbra.x2 = xmiddle;
fright = fmiddle;
}
else
{
cbra.x1 = xmiddle;
fleft = fmiddle;
}
}
}
return results;
}
int
main(int argc, char *argv[])
{
FunctionZeroSearch<double> fzcos(std::cos);
FunctionZeroSearch<double>::bracket test;
test.x1 = 0;
test.x2 = M_PI/4;
int result = fzcos.find_bracket_expand(test);
std::cout << "find_bracket_expand cos 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 = fzcos.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 = 0;
test.x2 = 4*M_PI;
result = fzcos.find_bracket_partition(test, 20);
std::cout << "get_bracket_partition cos with initial interval (0,4*PI) returns "
<< result << std::endl;
const FunctionZeroSearch<double>::bracket_container_type &brackets =
fzcos.get_brackets();
FunctionZeroSearch<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 = fzcos.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>