libcruft-util/roots/bisection.hpp

59 lines
1.7 KiB
C++
Raw Normal View History

/*
2018-08-04 15:14:06 +10:00
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Copyright 2016 Danny Robson <danny@nerdcruft.net>
*/
#pragma once
#include "../maths.hpp"
#include <cmath>
namespace cruft::roots {
/// find a root of a function using the bisection method
///
/// the user is responsible for ensuring there is in fact a root. there
/// is no iteration limit currently, so poorly considered arguments could
/// result in an infinite loop.
///
/// a: lower bound
/// b: upper bound
/// f: real function
/// tolerance: minimum range termination condition
template <typename T>
T
bisection (T a, T b, T (*f)(T), T tolerance)
{
// to simplify the exit conditions we assume that a..b is an
// increasing range.
if (a > b)
std::swap (a, b);
while (b - a > tolerance) {
// calculate the midpoint
auto c = (a + b) / T{2};
auto f_c = f (c);
// return early if we happened across the root
if (exactly_zero (f_c))
return c;
// check which interval we fell into, and divide for the next
// iteration
auto f_a = f (a);
if (samesign (f_a, f_c))
a = c;
else
b = c;
}
// we didn't find it exactly, but it's somewhere in this range so take
// a punt and return the middle of the range. not sure if this is
// valid mathematically.
return (a + b) / T{2};
}
2017-01-05 15:06:49 +11:00
}