Commit 34b8b965 authored by Timo Koch's avatar Timo Koch
Browse files

[test] Add test for scalar root finding algorithms

parent 6fa5f714
......@@ -36,14 +36,13 @@
namespace Dumux {
/*!
* \file
* \ingroup Common
* \brief Newton's root finding algorithm for scalar functions (secant method)
* \param xOld initial guess
* \param residual Residual function
* \param derivative Derivative of the residual
* \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
* \param maxIter Maximum number of iterations
* \todo replace std::result_of by std::invoke_result for C++17
*/
template<class Scalar, class ResFunc, class DerivFunc,
......@@ -68,7 +67,8 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFun
n--;
}
if (!std::isfinite(r))
using std::isfinite;
if (!isfinite(r))
DUNE_THROW(NumericalProblem, "Residual is not finite: " << r << " after " << maxIter - n << " iterations!");
if (relativeShift > tol)
......@@ -78,14 +78,13 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFun
}
/*!
* \file
* \ingroup Common
* \brief Newton's root finding algorithm for scalar functions (secant method)
* \note The derivative is numerically computed. If the derivative is know use signature with derivative function.
* \param xOld initial guess
* \param residual Residual function
* \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
* \param maxIter Maximum number of iterations
*/
template<class Scalar, class ResFunc,
typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
......@@ -98,7 +97,6 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
}
/*!
* \file
* \ingroup Common
* \brief Brent's root finding algorithm for scalar functions
* \note Modified from pseudo-code on wikipedia: https://en.wikipedia.org/wiki/Brent%27s_method
......@@ -107,8 +105,8 @@ Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
* \param a Lower bound
* \param b Upper bound
* \param residual Residual function
* \param maxIter Maximum number of iterations
* \param tol Relative shift tolerance
* \param maxIter Maximum number of iterations
*/
template<class Scalar, class ResFunc,
typename std::enable_if_t<Dune::Std::is_invocable<ResFunc, Scalar>::value>...>
......@@ -123,7 +121,7 @@ Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc& residual,
// check if the root is inside the given interval
using std::signbit;
if (!signbit(fa*fb))
DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain the root!");
DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain any, or no uniquely defined root!");
// sort bounds
using std::abs; using std::swap;
......
......@@ -4,5 +4,6 @@ add_subdirectory(freeflow)
add_subdirectory(io)
add_subdirectory(material)
add_subdirectory(multidomain)
add_subdirectory(nonlinear)
add_subdirectory(porousmediumflow)
add_subdirectory(discretization)
add_subdirectory(findscalarroot)
dumux_add_test(SOURCES test_findscalarroot.cc
LABELS unit nonlinear)
#include <config.h>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/float_cmp.hh>
#include <dumux/nonlinear/findscalarroot.hh>
int main(int argc, char* argv[]) try
{
using namespace Dumux;
const auto func = [](const double x){ return x*x - 5.0; };
const auto deriv = [](const double x){ return 2*x; };
const auto absTol = 1e-15;
const auto relTol = 1e-15;
const auto rootNewton = findScalarRootNewton(5.0, func, deriv, absTol);
const auto rootNewtonNumDiff = findScalarRootNewton(5.0, func, absTol);
const auto rootBrent = findScalarRootBrent(0.0, 5.0, func, relTol);
const auto rootExact = std::sqrt(5.0);
if (Dune::FloatCmp::ne(rootExact, rootNewton, absTol))
DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (exact derivative)! Relative error: " << std::abs(rootExact-rootNewton)/rootExact);
if (Dune::FloatCmp::ne(rootExact, rootNewtonNumDiff, absTol))
DUNE_THROW(Dune::Exception, "Wrong root for Newton algorithm (numeric differentiation)! Relative error: " << std::abs(rootExact-rootNewtonNumDiff)/rootExact);
if (Dune::FloatCmp::ne(rootExact, rootBrent, absTol))
DUNE_THROW(Dune::Exception, "Wrong root for Brent algorithm! Relative error: " << std::abs(rootExact-rootBrent)/rootExact);
std::cout << "Positive root of f(x) = (x*x - 5) is\n" << std::setprecision(15)
<< "-- Exact: " << rootExact << "\n"
<< "-- Newton (exact derivative): " << rootNewton << "\n"
<< "-- Newton (numeric derivative): " << rootNewtonNumDiff << "\n"
<< "-- Brent: " << rootBrent << std::endl;
return 0;
}
catch (const Dune::Exception& e)
{
std::cout << e << std::endl;
return 1;
}
catch (...)
{
std::cout << "Unknown exception thrown!" << std::endl;
return 1;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment