Skip to content

Instantly share code, notes, and snippets.

@Ben1980

Ben1980/bisection.h

Last active Jun 10, 2019
Embed
What would you like to do?
class Bisection : public Iteration {
public:
Bisection(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}
double solve(double a, double b) override {
resetNumberOfIterations();
checkAndFixAlgorithmCriteria(a, b);
fmt::print("Bisection -> [{:}, {:}]\n", a, b);
fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "x", "f(x)");
fmt::print("------------------------------------------------------------------------------------ \n");
double x = 0.5*(a+b);
double fx = mf(x);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, x, fx);
while(fabs(fx) >= epsilon()) {
x = calculateX(x, a, b, fx);
fx = mf(x);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, x, fx);
}
fmt::print("\n");
return x;
}
private:
void checkAndFixAlgorithmCriteria(double &a, double &b) const {
//Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
assert(mf(a)*mf(b) < 0);
if(mf(a) < mf(b)) {
std::swap(a,b);
}
}
static double calculateX(double x, double &a, double &b, double fx) {
if(fx >= 0) a = x;
if(fx < 0) b = x;
return 0.5*(a+b);
}
const std::function<double (double)> &mf;
};
class Brent : public Iteration {
public:
Brent(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}
double solve(double a, double b) override {
resetNumberOfIterations();
double fa = mf(a);
double fb = mf(b);
checkAndFixAlgorithmCriteria(a, b, fa, fb);
fmt::print("Brent -> [{:}, {:}]\n", a, b);
fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)", "f(s)");
fmt::print("---------------------------------------------------------------------------------------------------------- \n");
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb, "");
double lastB = a; // b_{k-1}
double lastFb = fa;
double s = std::numeric_limits<double>::max();
double fs = std::numeric_limits<double>::max();
double penultimateB = a; // b_{k-2}
bool bisection = true;
while(fabs(fb) > epsilon() && fabs(fs) > epsilon() && fabs(b-a) > epsilon()) {
if(useInverseQuadraticInterpolation(fa, fb, lastFb)) {
s = calculateInverseQuadraticInterpolation(a, b, lastB, fa, fb, lastFb);
}
else {
s = calculateSecant(a, b, fa, fb);
}
if(useBisection(bisection, a, b, lastB, penultimateB, s)) {
s = calculateBisection(a, b);
bisection = true;
}
else {
bisection = false;
}
fs = mf(s);
penultimateB = lastB;
lastB = b;
if(fa*fs < 0) {
b = s;
}
else {
a = s;
}
fa = mf(a);
lastFb = fb;
fb = mf(b);
checkAndFixAlgorithmCriteria(a, b, fa, fb);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb, fs);
}
fmt::print("\n");
return fb < fs ? b : s;
}
private:
static double calculateBisection(double a, double b) {
return 0.5*(a+b);
}
static double calculateSecant(double a, double b, double fa, double fb) {
//No need to check division by 0, in this case the method returns NAN which is taken care by useSecantMethod method
return b-fb*(b-a)/(fb-fa);
}
static double calculateInverseQuadraticInterpolation(double a, double b, double lastB, double fa, double fb, double lastFb) {
return a*fb*lastFb/((fa-fb)*(fa-lastFb)) +
b*fa*lastFb/((fb-fa)*(fb-lastFb)) +
lastB*fa*fb/((lastFb-fa)*(lastFb-fb));
}
static bool useInverseQuadraticInterpolation(double fa, double fb, double lastFb) {
return fa != lastFb && fb != lastFb;
}
static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) {
//Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
assert(fa*fb < 0);
if (fabs(fa) < fabs(fb)) {
std::swap(a, b);
std::swap(fa, fb);
}
}
bool useBisection(bool bisection, double a, double b, double lastB, double penultimateB, double s) const {
static const double DELTA = epsilon() + std::numeric_limits<double>::min();
return (bisection && fabs(s-b) >= 0.5*fabs(b-lastB)) || //Bisection was used in last step but |s-b|>=|b-lastB|/2 <- Interpolation step would be to rough, so still use bisection
(!bisection && fabs(s-b) >= 0.5*fabs(lastB-penultimateB)) || //Interpolation was used in last step but |s-b|>=|lastB-penultimateB|/2 <- Interpolation step would be to small
(bisection && fabs(b-lastB) < DELTA) || //If last iteration was using bisection and difference between b and lastB is < delta use bisection for next iteration
(!bisection && fabs(lastB-penultimateB) < DELTA); //If last iteration was using interpolation but difference between lastB ond penultimateB is < delta use biscetion for next iteration
}
const std::function<double (double)> &mf;
};
class Dekker : public Iteration {
public:
Dekker(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}
double solve(double a, double b) override {
resetNumberOfIterations();
double fa = mf(a);
double fb = mf(b);
checkAndFixAlgorithmCriteria(a, b, fa, fb);
fmt::print("Dekker -> [{:}, {:}]\n", a, b);
fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)");
fmt::print("------------------------------------------------------------------------------------ \n");
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb);
double lastB = a;
double lastFb = fa;
while(fabs(fb) > epsilon() && fabs(b-a) > epsilon()) {
const double s = calculateSecant(b, fb, lastB, lastFb);
const double m = calculateBisection(a, b);
lastB = b;
b = useSecantMethod(b, s, m) ? s : m;
lastFb = fb;
fb = mf(b);
if (fa * fb > 0 && fb * lastFb < 0) {
a = lastB;
}
fa = mf(a);
checkAndFixAlgorithmCriteria(a, b, fa, fb);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb);
}
fmt::print("\n");
return b;
}
private:
static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) {
//Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
assert(fa*fb < 0);
if (fabs(fa) < fabs(fb)) {
std::swap(a, b);
std::swap(fa, fb);
}
}
static double calculateSecant(double b, double fb, double lastB, double lastFb) {
//No need to check division by 0, in this case the method returns NAN which is taken care by useSecantMethod method
return b-fb*(b-lastB)/(fb-lastFb);
}
static double calculateBisection(double a, double b) {
return 0.5*(a+b);
}
static bool useSecantMethod(double b, double s, double m) {
//Value s calculated by secant method has to be between m and b
return (b > m && s > m && s < b) ||
(b < m && s > b && s < m);
}
const std::function<double (double)> &mf;
};
class Newton : public Iteration {
public:
Newton(double epsilon, const std::function<double (double)> &f, const std::function<double (double)> &fPrime) : Iteration(epsilon), mf(f), mfPrime(fPrime) {}
double solve(double x) override {
resetNumberOfIterations();
fmt::print("Newton -> [x0 = {:}]\n", x);
fmt::print("{:<5}|{:<20}|{:<20}|{:<20}\n", "K", "x", "f(x)", "f'(x)");
fmt::print("------------------------------------------------------------------ \n");
double fx = mf(x);
double fxPrime = mfPrime(x);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, fx, fxPrime);
while(fabs(fx) >= epsilon()) {
x = calculateX(x, fx, fxPrime);
fx = mf(x);
fxPrime = mfPrime(x);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, fx, fxPrime);
}
fmt::print("\n");
return x;
}
private:
static double calculateX(double x, double fx, double fxPrime) {
assert(fabs(fxPrime) >= std::numeric_limits<double>::min());
return x - fx/fxPrime;
}
const std::function<double (double)> &mf;
const std::function<double (double)> &mfPrime;
};
class Secant : public Iteration {
public:
Secant(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}
double solve(double a, double b) override {
resetNumberOfIterations();
fmt::print("Secant -> [{:}, {:}]\n", a, b);
fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)");
fmt::print("--------------------------------------------------------------------------------------\n");
if(mf(a) > mf(b)) {
std::swap(a,b);
}
double x = b;
double lastX = a;
double fx = mf(b);
double lastFx = mf(a);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, lastX, fx, lastFx);
while(fabs(fx) >= epsilon()) {
const double x_tmp = calculateX(x, lastX, fx, lastFx);
lastFx = fx;
lastX = x;
x = x_tmp;
fx = mf(x);
fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, lastX, fx, lastFx);
}
fmt::print("\n");
return x;
}
private:
static double calculateX(double x, double lastX, double fx, double lastFx) {
const double functionDifference = fx - lastFx;
assert(fabs(functionDifference) >= std::numeric_limits<double>::min());
return x - fx*(x-lastX)/functionDifference;
}
const std::function<double (double)> &mf;
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.