diff --git a/MathLib/Nonlinear/Root1D.h b/MathLib/Nonlinear/Root1D.h
new file mode 100644
index 0000000000000000000000000000000000000000..f015d66417ca1d8addfc54f3152e7f0254fbeda5
--- /dev/null
+++ b/MathLib/Nonlinear/Root1D.h
@@ -0,0 +1,165 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_NONLINEAR_ROOT1D_H
+#define MATHLIB_NONLINEAR_ROOT1D_H
+
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <type_traits>
+#include <logog/include/logog.hpp>
+
+namespace MathLib
+{
+namespace Nonlinear
+{
+
+namespace detail
+{
+
+//! Tells if \c a and \c b have the same sign.
+inline bool same_sign(double a, double b)
+{
+ // note: signbit() casts integers to double and distinguishes between +0.0 and
+ // -0.0. However, the latter is not a problem, because if f(x0) == 0, we've
+ // already found the root x0 that we are searching for.
+ return std::signbit(a) == std::signbit(b);
+}
+
+inline bool almost_zero(double a)
+{
+ return std::abs(a) <= std::numeric_limits<double>::epsilon();
+}
+
+}
+
+
+//! Use the regula falsi method to find the root of some scalar function of one
+//! variable.
+template<typename SubType, typename Function>
+class RegulaFalsi
+{
+public:
+ //! Initializes finding a root of the \c Function \c f in the interval
+ //! [\c a, \c b].
+ RegulaFalsi(Function const& f, double a, double b)
+ : _f(f), _a(a), _b(b), _fa(f(a)), _fb(f(b))
+ {
+ static_assert(
+ std::is_same<double, decltype(f(0.0))>::value,
+ "Using this class for functions that do not return double"
+ " involves a lot of casts. Hence it is disabled.");
+
+ if (detail::almost_zero(_fa)) {
+ _b = _a;
+ } else if (detail::almost_zero(_fb)) {
+ _a = _b;
+ } else if (detail::same_sign(_fa, _fb)) {
+ ERR("Regula falsi cannot be done, because the function values"
+ " at the interval ends have the same sign.");
+ std::abort();
+ }
+ }
+
+ //! Do \c num_steps iteration of regula falsi.
+ void step(const unsigned num_steps)
+ {
+ for (unsigned i=0; i<num_steps; ++i)
+ {
+ if (_a == _b) return;
+
+ const double s = (_fb - _fa)/(_b - _a);
+ const double c = _a - _fa/s;
+ const double fc = _f(c);
+
+ if (detail::almost_zero(fc)) {
+ _a = _b = c;
+ return;
+ } else if (!detail::same_sign(fc, _fb)) {
+ _a = _b; _fa = _fb;
+ _b = c; _fb = fc;
+ } else {
+ const double m = SubType::get_m(_fa, _fb, fc);
+ _fa *= m;
+ _b = c;
+ _fb = fc;
+ }
+ }
+ }
+
+ //! Returns the current estimate of the root.
+ double getResult() const
+ {
+ if (_a == _b) return _a;
+
+ const double s = (_fb - _fa)/(_b - _a);
+ const double c = _a - _fa/s;
+
+ return c;
+ }
+
+ //! Returns the size of the current search interval.
+ double getRange() const { return std::fabs(_a - _b); }
+
+private:
+ Function const& _f;
+ double _a, _b, _fa, _fb;
+};
+
+
+/*! Creates a new instance of \c RegulaFalsi that uses the given modified iteration
+ * method \c SubType.
+ *
+ * \see https://en.wikipedia.org/wiki/False_position_method#Improvements_in_regula_falsi
+ */
+template<typename SubType, typename Function>
+RegulaFalsi<SubType, Function>
+makeRegulaFalsi(Function const& f, double const a, double const b)
+{
+ return RegulaFalsi<SubType, Function>(f, a, b);
+}
+
+
+//! Used by \c RegulaFalsi in the original regula falsi algorithm.
+struct Unmodified
+{
+ static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
+ { return 1.0; }
+};
+
+//! Used by \c RegulaFalsi in a modified version of the regula falsi algorithm.
+struct Illinois
+{
+ static double get_m(const double /*fa*/, const double /*fb*/, const double /*fc*/)
+ { return 0.5; }
+};
+
+//! Used by \c RegulaFalsi in a modified version of the regula falsi algorithm.
+struct Pegasus
+{
+ static double get_m(const double /*fa*/, const double fb, const double fc)
+ { return fb / (fb+fc); }
+};
+
+//! Used by \c RegulaFalsi in a modified version of the regula falsi algorithm.
+struct AndersonBjorck
+{
+ static double get_m(const double /*fa*/, const double fb, const double fc)
+ {
+ const double f = 1.0 - fc / fb;
+ return (f >= 0.0) ? f : 0.5;
+ }
+};
+
+}
+
+} // namespace MathLib
+
+#endif // MATHLIB_NONLINEAR_ROOT1D_H
diff --git a/Tests/MathLib/TestNonlinear1D.cpp b/Tests/MathLib/TestNonlinear1D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7968acb414f070af66339c91c4b57abf8c8e8284
--- /dev/null
+++ b/Tests/MathLib/TestNonlinear1D.cpp
@@ -0,0 +1,59 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <limits>
+#include <type_traits>
+#include <gtest/gtest.h>
+#include <logog/include/logog.hpp>
+#include "MathLib/Nonlinear/Root1D.h"
+
+double f(double x)
+{
+ return x*x-1;
+}
+
+template <typename T>
+class MathLibRegulaFalsi : public ::testing::Test
+{
+};
+
+namespace NL = MathLib::Nonlinear;
+using RegulaFalsiTypes = ::testing::Types<NL::Unmodified, NL::Illinois, NL::Pegasus, NL::AndersonBjorck>;
+
+TYPED_TEST_CASE(MathLibRegulaFalsi, RegulaFalsiTypes);
+
+TYPED_TEST(MathLibRegulaFalsi, QuadraticFunction)
+{
+ auto rf = NL::makeRegulaFalsi<TypeParam>(f, -0.1, 1.1);
+ double old_range = rf.getRange();
+
+ DBUG(" 0 -- x ~ %23.16g, range = %23.16g", rf.getResult(), old_range);
+
+ for (unsigned n=0; n<10; ++n)
+ {
+ rf.step(1);
+ double range = rf.getRange();
+ // expect that the interval of the root search shrinks
+ EXPECT_GT(old_range, range);
+ old_range = range;
+ DBUG("%2i -- x ~ %23.16g, range = %23.16g", n+1, rf.getResult(), range);
+
+ if (range < std::numeric_limits<double>::epsilon()) break;
+ }
+
+ auto const error = std::abs(f(rf.getResult()));
+
+ if (!std::is_same<NL::Unmodified, TypeParam>::value) {
+ EXPECT_GT(std::numeric_limits<double>::epsilon(), old_range);
+ EXPECT_GT(std::numeric_limits<double>::epsilon(), error);
+ } else {
+ // The unmodified regula falsi method converges very slowly.
+ EXPECT_GT(100.0*std::numeric_limits<double>::epsilon(), error);
+ }
+}