OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add combined regula falsi + bisection method, and test.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2019-11-05 10:33:02 +01:00
parent 53619980c5
commit 2951df2bfb
3 changed files with 185 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -258,6 +258,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_nonuniformtablelinear.cpp
tests/test_OpmLog.cpp
tests/test_param.cpp
tests/test_RootFinders.cpp
tests/test_SimulationDataContainer.cpp
tests/test_sparsevector.cpp
tests/test_uniformtablelinear.cpp

139
opm/common/utility/numeric/RootFinders.hpp
View File

@ -1,21 +1,6 @@
//===========================================================================
//
// File: RootFinders.hpp
//
// Created: Thu May 6 19:59:42 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Jostein R Natvig <jostein.r.natvig@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010 Statoil ASA.
Copyright 2010, 2019 SINTEF Digital
Copyright 2010, 2019 Equinor ASA
This file is part of the Open Porous Media project (OPM).
@ -295,6 +280,126 @@ namespace Opm
template <class ErrorPolicy = ThrowOnError>
class RegulaFalsiBisection
{
public:
inline static std::string name()
{
return "RegulaFalsiBisection";
}
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford. Current variant is the 'Pegasus' method.
/// Combines this method with the bisection method, inspired by
/// http://phillipmfeldman.org/Python/roots/find_roots.html
template <class Functor>
inline static double solve(const Functor& f,
const double a,
const double b,
const int max_iter,
const double tolerance,
int& iterations_used)
{
using namespace std;
const double sqrt_half = std::sqrt(0.5);
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps * max(max(fabs(a), fabs(b)), 1.0);
double x0 = a;
double x1 = b;
double f0 = f(x0);
const double epsF = tolerance + macheps * max(fabs(f0), 1.0);
if (fabs(f0) < epsF) {
return x0;
}
double f1 = f(x1);
if (fabs(f1) < epsF) {
return x1;
}
if (f0 * f1 > 0.0) {
return ErrorPolicy::handleBracketingFailure(a, b, f0, f1);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
double width = fabs(x1 - x0);
double contraction = 1.0;
while (width >= 1e-9 * eps) {
// If we are contracting sufficiently to at least halve
// the interval width in two iterations we use regula
// falsi. Otherwise, we take a bisection step to avoid
// slow convergence.
const double xnew = (contraction < sqrt_half)
? regulaFalsiStep(x0, x1, f0, f1)
: bisectionStep(x0, x1, f0, f1);
const double fnew = f(xnew);
++iterations_used;
if (iterations_used > max_iter) {
return ErrorPolicy::handleTooManyIterations(x0, x1, max_iter);
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi1/(1.0 - phi0);
// 3. The method called 'Method 4' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi0 - phi1;
// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
// " gamma = " << gamma << endl;
// 4. The 'Pegasus' method
const double gamma = f1 / (f1 + fnew);
// 5. Straight unmodified Regula Falsi
// const double gamma = 1.0;
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
const double widthnew = fabs(x1 - x0);
contraction = widthnew/width;
width = widthnew;
}
return 0.5 * (x0 + x1);
}
private:
inline static double regulaFalsiStep(const double a, const double b, const double fa, const double fb)
{
assert(fa * fb < 0.0);
return (b * fa - a * fb) / (fa - fb);
}
inline static double bisectionStep(const double a, const double b, const double fa, const double fb)
{
static_cast<void>(fa);
static_cast<void>(fb);
assert(fa * fb < 0.0);
return 0.5*(a + b);
}
};
/// Attempts to find an interval bracketing a zero by successive
/// enlargement of search interval.
template <class Functor>

62
tests/test_RootFinders.cpp Normal file
View File

@ -0,0 +1,62 @@
/*
Copyright 2019 Equinor ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE RootFindersTest
#include <boost/test/unit_test.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
using namespace Opm;
template<class Method>
struct Test
{
template <class Functor>
static int run(const Functor& f,
const double a,
const double b,
const int max_iter,
const double tolerance)
{
int iter = 0;
Method::solve(f, a, b, max_iter, tolerance, iter);
return iter;
}
};
BOOST_AUTO_TEST_CASE(SimpleFunction)
{
auto f = [](double x) { return (x - 1.0)*x; };
BOOST_CHECK_EQUAL(Test<RegulaFalsi<>>::run(f, -0.5, 0.99, 50, 1e-12), 10);
BOOST_CHECK_EQUAL(Test<RegulaFalsiBisection<>>::run(f, -0.5, 0.99, 50, 1e-12), 15);
}
BOOST_AUTO_TEST_CASE(ToughFunction)
{
auto f = [](double x) { return x < 50.0 ? -0.0001 : x - 50.0001; };
BOOST_CHECK_THROW(Test<RegulaFalsi<>>::run(f, 0.0, 100.0, 50, 1e-12), std::exception);
BOOST_CHECK_EQUAL(Test<RegulaFalsi<>>::run(f, 0.0, 100.0, 1000000, 1e-12), 90);
BOOST_CHECK_EQUAL(Test<RegulaFalsiBisection<>>::run(f, 0.0, 100.0, 1000000, 1e-12), 2);
}