Commit 30bc76c2 authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

New bisection implementation

For now it's in the cc file, because I didn't come up
with a nice name. The new version also exactly documents
what it does.
parent 2187fbb8
Pipeline #8677 passed with stage
in 5 minutes and 27 seconds
......@@ -130,6 +130,200 @@ public:
};
class NewBisection
{
template <class K>
static inline K dist(const Solvers::Interval<K>& I, K x)
{
if (x < I[0])
return I[0]-x;
if (x > I[1])
return x-I[1];
return 0;
}
public:
/** \brief Constructor a bisection solver
*
* This computes an approximation of a root of the given function
* using bisection.
*
* It is assumed that the given initial interval I_0 contains
* at least one root and that the initial interval
* provides a sign change, i.e. I_0[0]*I_0[1] <= 0.
*
* Termination can be done using absolute and relative errors.
* Assume that x^* is the solution, x some iterate, and x_ref some
* reference value. Then the absolute and relative error are defined
* by |x-x^*| and |x-x^*|/|x_ref-x^*|.
*
* Notice that, if x_ref is 'on the I_0[0] side'
*
* Assuming that x is from the interval I that brackets the root,
* computable upper bounds for the absolute and relative error are
* given by |x_k-x^*| <= absErr := |I[1] - I[0]|
* and |x_k-x^*|/|x_ref-x^*| <= relErr := absErr/dist(x_ref,I).
* Notice that both bounds only depend on I and are valid for
* all points in I.
*
* The algorithm constructs a sequence of bracketing intervals I_k.
* It will terminate if either of the following is true:
* 1) The residual satisfies |f(I_k[1])| <= residualTol.
* 2) The residual satisfies |f(I_k[0])| <= residualTol.
* 3) The absolute error bound satisfies absErr <= absTol.
* 4) The relative error bound satisfies relErr <= relTol.
* 5) I_k[0] and I_k[1] are consecutive floating point numbers.
*
* In the first case the alrorithm will return I_k[1] in all other
* cases I_k[0]. Assuming that I_0 contains a unique root x^*, this
* will always be the point in I_k that is closest to I_0[0].
*
* Remark: The algorithm will overflow, if the initial interval
* produces an overflow for I_0[0]+I_0[1]. Avoiding this would
* require a more expenisive algorithm for the midpoint computation.
*
* \param absTol Absolute error tolerance
* \param relTol Relative error tolerance
* \param residualTol Residual tolerance
*/
NewBisection(double absTol = 0.0, double relTol = 0.0, double residualTol = 0.0) :
absTol_(absTol),
relTol_(relTol),
residualTol_(residualTol)
{};
/** \brief Find root of function
*
* \param f The function, e.g. the derivative of a convex functional
* \param I Initial interval
* \param fI Values at initial interval
* \param x_ref Reference value. This is used to determine the relative error |x-x^*|/|x_ref-x^*|
* \return Approximate solution
*/
template <class F>
double solve(const F& f, Solvers::Interval<double> I, Solvers::Interval<double> fI, double x_ref, int verbosity=0) const
{
double mid;
double fmid;
double absErr;
double d;
std::size_t k = 0;
if (verbosity>0)
{
std::cout << " # | left | mid | right | f(left) | f(mid) | f(right) |" << std::endl;
std::cout << "----|--------------|--------------|--------------|--------------|--------------|--------------|" << std::endl;
auto p = std::cout.precision();
std::cout << std::setw(4) << k << "|"
<< std::setprecision(8)
<< std::setw(14) << I[0] << "|"
<< " |"
<< std::setw(14) << I[1] << "|"
<< std::setw(14) << fI[0] << "|"
<< " |"
<< std::setw(14) << fI[1] << "|" << std::endl;
std::cout << std::setprecision(p);
}
// Terminate if residual bound is satisfied
if (std::fabs(fI[1]) <= residualTol_)
return I[1];
if (std::fabs(fI[0]) <= residualTol_)
return I[0];
// Check if we have a sign change. Notice that signbit
// identifies negative values. While the criterion
// would not be satisfied if we have a non-negative
// value and zero, that's not a problem, because we
// excluded zero by the previous termination check.
if (std::signbit(fI[0]) == std::signbit(fI[1]))
DUNE_THROW(Exception, "Initial interval provided to bisection has no sign change");
// Iterate while midpoint does not coincide with boundary numerically
while (std::nextafter(I[0],I[1]) != I[1])
{
++k;
// Here we would like to check for the residual criterion as follows:
//
// if (std::fabs(fI[1]) <= residualTol_)
// return I[1];
// if (std::fabs(fI[0]) <= residualTol_)
// return I[0];
//
// This would result in always double checking for the unchanged
// interval bound. To avoid this we can equivalently only check
// for the midpoint at the end of the previous loop (see below).
// Terminate if absolute error bound is satistied.
absErr = std::fabs(I[1]-I[0]);
if (absErr <= absTol_)
return I[0];
// Terminate if relative error bound is satistied.
d = dist(I, x_ref);
if (absErr <= relTol_*d)
return I[0];
// Compute midpoint
double mid = (I[0] + I[1]) / 2;
double fmid = f(mid);
if (verbosity>0)
{
auto p = std::cout.precision();
std::cout << std::setw(4) << k << "|"
<< std::setprecision(8)
<< std::setw(14) << I[0] << "|"
<< std::setw(14) << mid << "|"
<< std::setw(14) << I[1] << "|"
<< std::setw(14) << fI[0] << "|"
<< std::setw(14) << fmid << "|"
<< std::setw(14) << fI[1] << "|" << std::endl;
std::cout << std::setprecision(p);
}
// Terminate if residual bound is satisfied. By checking for
// the midpoint here, we don't need to check for |f(I[0])|
// and |f(I[1])| in the next loop (see above).
if (std::fabs(fmid) <= residualTol_)
return mid;
// Replace suitable bound by midpoint. Since the value zero
// is excluded by the previous residual checks, signbit
// is here really an indicator for strict negative or positive
// value.
if (std::signbit(fmid) == std::signbit(fI[0]))
{
I[0] = mid;
fI[0] = fmid;
}
else
{
I[1] = mid;
fI[1] = fmid;
}
}
return I[0];
}
private:
double absTol_;
double relTol_;
double residualTol_;
};
/* Class to project a correction onto the admissible set.
* All we have is bound constraints. Therefore, in an ideal world, we would
* simply use the ObstacleDefectProjection class instead of this one. However,
......@@ -635,11 +829,60 @@ int main (int argc, char *argv[]) try
ScalarSumFunctional<decltype(std::get<0>(restriction.functions())),
decltype(std::get<1>(restriction.functions()))> sumFunctional(std::get<0>(restriction.functions()), std::get<1>(restriction.functions()));
#if 0
Bisection bisection;
xi = 1;
int bisectionSteps;
// Old version why do we use xi here? Especially the second one
// does not make any sense.
xi = bisection.minimize(sumFunctional, xi, xi, bisectionSteps);
// It should probably by one of the following:
// xi = bisection.minimize(sumFunctional, 0, 0, bisectionSteps);
// xi = bisection.minimize(sumFunctional, 1, 0, bisectionSteps);
#else
auto f = [&](auto z) {
return sumFunctional.subDifferential(z).projectIn(0);
};
auto I = sumFunctional.domain();
I[0] = 0;
auto fI = I;
fI[0] = f(I[0]);
if (f(0) < 0)
{
fI[1] = f(I[1]);
// Absolute error tolerance 0
// Relative error tolerance 0.1
// If possible, start from interval [1,3].
NewBisection bisection(0, 0.2);
if (I[1] >= 1)
{
auto f1 = f(1);
if (f1 <= 0)
{
I[0] = 1;
fI[0] = f1;
}
}
if (I[1] >= 3)
{
auto f3 = f(3);
if (f3 >= 0)
{
I[1] = 3;
fI[1] = f3;
}
}
xi = bisection.solve(f, I, fI, 0);
}
else
xi = 0;
#endif
};
using QuadraticPartLinearization = TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<QuadraticPartFunctional, BitVector>;
......
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