291 lines
5.2 KiB
C++
291 lines
5.2 KiB
C++
|
|
|
|
#ifndef TZEROFINDER_H
|
|
#define TZEROFINDER_H
|
|
|
|
#include <cmath>
|
|
#include <assert.h>
|
|
|
|
enum {
|
|
FZ_OK = 0,
|
|
FZ_FTOL = 1,
|
|
FZ_XTOL = 2,
|
|
FZ_MAX_ITER = 3,
|
|
FZ_NOT_ZERO = 4
|
|
};
|
|
|
|
/*
|
|
|
|
//Daniele
|
|
// Use the one at the end instead. This may overshoot - plus, it does not
|
|
// bisecate the interval... (I suppose it has never been used much :) )
|
|
|
|
template< class unaryOp ,
|
|
class unaryOpDerivate>
|
|
bool findZero_Newton ( double x0 ,
|
|
double x1 ,
|
|
unaryOp f ,
|
|
unaryOpDerivate fd1 ,
|
|
double xtol ,
|
|
double ftol ,
|
|
int maxIter ,
|
|
double &x ,
|
|
int &err )
|
|
{
|
|
double f_x_n_plus_1 = 0.0 ;
|
|
double f_x_n = f(x1) ;
|
|
double f_x_n_minus_1 = f(x0) ;
|
|
|
|
if( f_x_n*f_x_n_minus_1 > 0 )
|
|
{
|
|
err = FZ_NOT_ZERO;
|
|
return false;
|
|
}
|
|
|
|
int i;
|
|
err=FZ_OK;
|
|
|
|
if( f_x_n == 0 ||
|
|
f_x_n_minus_1 == 0 )
|
|
{
|
|
x = (f_x_n == 0.0) ? x1 : x0;
|
|
err = FZ_OK;
|
|
return true;
|
|
}
|
|
|
|
double x_n_plus_1 ;
|
|
double x_n ;
|
|
double x_n_minus_1 ;
|
|
|
|
x_n = x1;
|
|
x_n_minus_1 = x0;
|
|
|
|
for(i=0;i<maxIter;++i)
|
|
{
|
|
// the firts time x_n = x1
|
|
double fd1_x_n = fd1(x_n);
|
|
assert( fd1_x_n != 0.0 );
|
|
|
|
x_n_plus_1 = x_n - f_x_n / fd1_x_n;
|
|
if( fabs(x_n - x_n_plus_1) < xtol )
|
|
{
|
|
err=FZ_XTOL;
|
|
break;
|
|
}
|
|
|
|
f_x_n_plus_1 = f(x_n_plus_1);
|
|
if( fabs(f_x_n_plus_1) < ftol )
|
|
{
|
|
err=FZ_FTOL;
|
|
break;
|
|
}
|
|
|
|
x_n_minus_1 = x_n;
|
|
x_n = x_n_plus_1;
|
|
f_x_n = f_x_n_plus_1;
|
|
}
|
|
|
|
x=x_n_plus_1;
|
|
err = (i == maxIter ? FZ_MAX_ITER : err);
|
|
return true;
|
|
}
|
|
*/
|
|
|
|
//-----------------------------------------------------------------------------
|
|
|
|
/*!
|
|
|
|
Bisection algorithm.
|
|
Get an object and a method of object class.
|
|
x1 and x2 are extremes of test and xtol
|
|
is the accuracy.
|
|
|
|
Return:
|
|
the minimum or,
|
|
-2 if min it's unreacheable.
|
|
-1 if min does not exist
|
|
Remark:
|
|
Max iterations numbers is fixed to 100 in in maxIter.
|
|
This value fix the accuracy to 2e-100.
|
|
|
|
From Numerical Recipes.
|
|
*/
|
|
|
|
template <class unaryOp>
|
|
bool findZero_bisection(double x0,
|
|
double x1,
|
|
unaryOp f,
|
|
double xtol,
|
|
int maxIter,
|
|
double &x,
|
|
int &err)
|
|
{
|
|
int i;
|
|
err = FZ_OK;
|
|
double dx, fx0, fmid, xmid, rtb;
|
|
|
|
fx0 = f(x0);
|
|
fmid = f(x1);
|
|
|
|
if (fx0 * fmid > 0.0) {
|
|
err = FZ_NOT_ZERO;
|
|
return false;
|
|
}
|
|
|
|
if (fx0 == 0 ||
|
|
fmid == 0) {
|
|
x = (fmid == 0.0) ? x1 : x0;
|
|
err = FZ_OK;
|
|
return true;
|
|
}
|
|
|
|
rtb = fx0 < 0.0 ? (dx = x1 - x0, x0) : (dx = x0 - x1, x1);
|
|
|
|
for (i = 1; i <= maxIter; i++) {
|
|
fmid = f(xmid = rtb + (dx *= 0.5));
|
|
if (fmid <= 0.0)
|
|
rtb = xmid;
|
|
if (fabs(dx) < xtol || fmid == 0.0) {
|
|
x = rtb;
|
|
err = FZ_XTOL;
|
|
break;
|
|
}
|
|
}
|
|
err = (i == maxIter ? FZ_MAX_ITER : err);
|
|
return true;
|
|
}
|
|
|
|
template <class unaryOp>
|
|
bool findZero_secant(double x0,
|
|
double x1,
|
|
unaryOp f,
|
|
double xtol,
|
|
double ftol,
|
|
int maxIter,
|
|
double &x,
|
|
int &err)
|
|
{
|
|
double f_x_n_plus_1 = 0.0,
|
|
f_x_n = f(x1),
|
|
f_x_n_minus_1 = f(x0);
|
|
|
|
if (f_x_n * f_x_n_minus_1 > 0) {
|
|
err = FZ_NOT_ZERO;
|
|
return false;
|
|
}
|
|
|
|
int i;
|
|
err = FZ_OK;
|
|
|
|
if (f_x_n == 0 ||
|
|
f_x_n_minus_1 == 0) {
|
|
x = (f_x_n == 0.0) ? x1 : x0;
|
|
return true;
|
|
}
|
|
|
|
double x_n_plus_1 = x1,
|
|
x_n = x1,
|
|
x_n_minus_1 = x0;
|
|
|
|
for (i = 0; i < maxIter; ++i) {
|
|
double den = f_x_n - f_x_n_minus_1;
|
|
assert(den != 0.0);
|
|
if (den == 0.0) {
|
|
x_n = x_n_plus_1;
|
|
err = FZ_FTOL;
|
|
break;
|
|
}
|
|
|
|
x_n_plus_1 = x_n - (x_n - x_n_minus_1) * f_x_n / den;
|
|
if (fabs(x_n - x_n_plus_1) < xtol) {
|
|
err = FZ_XTOL;
|
|
break;
|
|
}
|
|
|
|
f_x_n_plus_1 = f(x_n_plus_1);
|
|
if (fabs(f_x_n_plus_1) < ftol) {
|
|
err = FZ_FTOL;
|
|
break;
|
|
}
|
|
|
|
x_n_minus_1 = x_n;
|
|
x_n = x_n_plus_1;
|
|
f_x_n_minus_1 = f_x_n;
|
|
f_x_n = f_x_n_plus_1;
|
|
}
|
|
|
|
x = x_n_plus_1;
|
|
err = (i == maxIter ? FZ_MAX_ITER : err);
|
|
return true;
|
|
}
|
|
|
|
template <class unaryOp1, class unaryOp2>
|
|
bool findZero_Newton(double x0,
|
|
double x1,
|
|
unaryOp1 f,
|
|
unaryOp2 f1,
|
|
double xtol,
|
|
double ftol,
|
|
int maxIter,
|
|
double &x,
|
|
int &err)
|
|
{
|
|
double f_0 = f(x0);
|
|
double f_1 = f(x1);
|
|
double f_x_n = f_0, f_x_n_plus_1 = 0.0;
|
|
|
|
bool s_0 = (f_0 > 0);
|
|
bool s_1 = (f_1 > 0);
|
|
bool s_n;
|
|
|
|
if (s_0 && s_1) {
|
|
err = FZ_NOT_ZERO;
|
|
return false;
|
|
}
|
|
|
|
int i;
|
|
err = FZ_OK;
|
|
|
|
if (f_0 == 0 ||
|
|
f_1 == 0) {
|
|
x = (f_1 == 0.0) ? x1 : x0;
|
|
return true;
|
|
}
|
|
|
|
double x_n = x1, x_n_plus_1;
|
|
|
|
for (i = 0; i < maxIter; ++i) {
|
|
double den = f1(x_n);
|
|
if (den > 1e-2)
|
|
x_n_plus_1 = x_n - f_x_n / den;
|
|
|
|
if (den <= 1e-2 || x0 > x_n_plus_1 || x1 < x_n_plus_1)
|
|
x_n_plus_1 = 0.5 * (x0 + x1);
|
|
|
|
if (fabs(x_n - x_n_plus_1) < xtol) {
|
|
err = FZ_XTOL;
|
|
break;
|
|
}
|
|
|
|
f_x_n_plus_1 = f(x_n_plus_1);
|
|
if (fabs(f_x_n_plus_1) < ftol) {
|
|
err = FZ_FTOL;
|
|
break;
|
|
}
|
|
|
|
x_n = x_n_plus_1;
|
|
f_x_n = f_x_n_plus_1;
|
|
s_n = (f_x_n > 0);
|
|
|
|
if (s_n == s_0)
|
|
x0 = x_n;
|
|
else
|
|
x1 = x_n;
|
|
}
|
|
|
|
x = x_n_plus_1;
|
|
err = (i == maxIter ? FZ_MAX_ITER : err);
|
|
return true;
|
|
}
|
|
|
|
#endif
|