3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
findscalarroot.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
25#define DUMUX_COMMON_SCALAR_ROOT_FINDING_HH
26
27#include <cmath>
28#include <limits>
29#include <type_traits>
30
34
35namespace Dumux {
36
46template<class Scalar, class ResFunc, class DerivFunc,
47 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar> &&
48 std::is_invocable_r_v<Scalar, DerivFunc, Scalar>>...>
49Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual, const DerivFunc& derivative,
50 const Scalar tol = 1e-13, const int maxIter = 200)
51{
52 Scalar xNew = xOld;
53 Scalar r = residual(xNew);
54
55 int n = maxIter;
56 Scalar relativeShift = std::numeric_limits<Scalar>::max();
57 while (relativeShift > tol && n > 0)
58 {
59 xNew = xOld - r/derivative(xOld);
60 r = residual(xNew);
61
62 using std::abs; using std::max;
63 relativeShift = abs(xOld-xNew)/max(abs(xOld), abs(xNew));
64 xOld = xNew;
65 n--;
66 }
67
68 using std::isfinite;
69 if (!isfinite(r))
70 DUNE_THROW(NumericalProblem, "Residual is not finite: " << r << " after " << maxIter - n << " iterations!");
71
72 if (relativeShift > tol)
73 DUNE_THROW(NumericalProblem, "Scalar newton solver did not converge after " << maxIter << " iterations!");
74
75 return xNew;
76}
77
87template<class Scalar, class ResFunc,
88 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar>>...>
89Scalar findScalarRootNewton(Scalar xOld, const ResFunc& residual,
90 const Scalar tol = 1e-13, const int maxIter = 200)
91{
92 const auto eps = NumericDifferentiation::epsilon(xOld);
93 auto derivative = [&](const auto x){ return (residual(x + eps)-residual(x))/eps; };
94 return findScalarRootNewton(xOld, residual, derivative, tol, maxIter);
95}
96
109template<class Scalar, class ResFunc,
110 typename std::enable_if_t<std::is_invocable_r_v<Scalar, ResFunc, Scalar>>...>
111Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc& residual,
112 const Scalar tol = 1e-13, const int maxIter = 200)
113{
114 // precompute the residuals (minimize function evaluations)
115 Scalar fa = residual(a);
116 Scalar fb = residual(b);
117 Scalar fs = 0.0;
118
119 // check if the root is inside the given interval
120 using std::signbit;
121 if (!signbit(fa*fb))
122 DUNE_THROW(NumericalProblem, "Brent's algorithm failed: [a,b] does not contain any, or no uniquely defined root!");
123
124 // sort bounds
125 using std::abs; using std::swap;
126 if (abs(fa) < abs(fb))
127 {
128 swap(a, b);
129 swap(fa, fb);
130 }
131
132 Scalar c = a;
133 Scalar fc = fa;
134 Scalar d = 0.0;
135 Scalar s = 0.0;
136 bool flag = true;
137
138 for (int i = 0; i < maxIter; ++i)
139 {
140 // stopping criterion
141 using std::max;
142 if (abs(b-a) < tol*max(abs(a), abs(b)))
143 return b;
144
145 // inverse quadratic interpolation
146 if (fa != fc && fb != fc)
147 {
148 const auto fab = fa-fb;
149 const auto fac = fa-fc;
150 const auto fbc = fb-fc;
151 s = a*fb*fc/(fab*fac) - b*fa*fc/(fab*fbc) + c*fa*fb/(fac*fbc);
152 }
153
154 // secant method
155 else
156 {
157 s = b - fb*(b-a)/(fb-fa);
158 }
159
160 // bisection method
161 if ( (s < (3*a + b)*0.25 || s > b)
162 || (flag && abs(s-b) >= abs(b-c)*0.5)
163 || (!flag && abs(s-b) >= abs(c-d)*0.5)
164 || (flag && abs(b-c) < tol*max(abs(b), abs(c)))
165 || (!flag && abs(c-d) < tol*max(abs(c), abs(d))) )
166 {
167 s = (a+b)*0.5;
168 flag = true;
169 }
170 else
171 flag = false;
172
173 // compute residual at new guess
174 fs = residual(s);
175 d = c;
176 c = b;
177 fc = fb;
178
179 // check on which side of the root s falls
180 if (fa*fs < 0.0)
181 {
182 b = s;
183 fb = fs;
184 }
185 else
186 {
187 a = s;
188 fa = fs;
189 }
190
191 // sort if necessary
192 if (abs(fa) < abs(fb))
193 {
194 swap(a, b);
195 swap(fa, fb);
196 }
197 }
198
199 DUNE_THROW(NumericalProblem, "Brent's algorithm didn't converge after " << maxIter << " iterations!");
200}
201
202} // end namespace Dumux
203
204#endif
Some exceptions thrown in DuMux
A class for numeric differentiation.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Scalar findScalarRootNewton(Scalar xOld, const ResFunc &residual, const DerivFunc &derivative, const Scalar tol=1e-13, const int maxIter=200)
Newton's root finding algorithm for scalar functions (secant method)
Definition: findscalarroot.hh:49
Scalar findScalarRootBrent(Scalar a, Scalar b, const ResFunc &residual, const Scalar tol=1e-13, const int maxIter=200)
Brent's root finding algorithm for scalar functions.
Definition: findscalarroot.hh:111
Definition: adapt.hh:29
Exception thrown if a fixable numerical problem occurs.
Definition: exceptions.hh:39
static Scalar epsilon(const Scalar value, const Scalar baseEps=1e-10)
Computes the epsilon used for numeric differentiation.
Definition: numericdifferentiation.hh:47