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