{{ message }}

Instantly share code, notes, and snippets.

Ben1980/bisection.h

Last active Jun 10, 2019
 class Bisection : public Iteration { public: Bisection(double epsilon, const std::function &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 &mf; };
 class Brent : public Iteration { public: Brent(double epsilon, const std::function &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::max(); double fs = std::numeric_limits::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::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 &mf; };
 class Dekker : public Iteration { public: Dekker(double epsilon, const std::function &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 &mf; };
 class Newton : public Iteration { public: Newton(double epsilon, const std::function &f, const std::function &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::min()); return x - fx/fxPrime; } const std::function &mf; const std::function &mfPrime; };
 class Secant : public Iteration { public: Secant(double epsilon, const std::function &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::min()); return x - fx*(x-lastX)/functionDifference; } const std::function &mf; };