Commit 54d820b7 authored by Gräser, Carsten's avatar Gräser, Carsten
Browse files

Improve algorithm for initial bisection interval

The old version was to pessimistic and led to numerical errors
if the domain was very large (e.g., because all constrained dofs
have been zero).
parent 923098cb
Pipeline #8762 passed with stage
in 3 minutes and 49 seconds
......@@ -840,7 +840,7 @@ int main (int argc, char *argv[]) try
// It should probably by one of the following:
// xi = bisection.minimize(sumFunctional, 0, 0, bisectionSteps);
// xi = bisection.minimize(sumFunctional, 1, 0, bisectionSteps);
#else
#elif 0
auto f = [&](auto z) {
return sumFunctional.subDifferential(z).projectIn(0);
......@@ -882,6 +882,71 @@ int main (int argc, char *argv[]) try
}
else
xi = 0;
#else
auto f = [&](auto z) {
return sumFunctional.subDifferential(z).projectIn(0);
// auto result = sumFunctional.subDifferential(z).projectIn(0);
// std::cout << "f(" << z << ") = " << result << std::endl;
// return result;
};
auto contains = [](const Solvers::Interval<double>& J, double t) {
return (J[0] <= t) and (t <= J[1]);
};
auto dom = sumFunctional.domain();
auto I = dom;
auto fI = I;
I[0] = 0;
fI[0] = f(I[0]);
// Ensure that we have a descent direction. Otherwise use step length = 0
if (fI[0] >= 0)
{
xi = 0;
return;
}
// If the domain contains 1, try to find a better initial interval guess.
// Otherwise use domain as initial guess.
if (contains(I, 1))
{
// Start from z=1 and successively double the value of z
// until we have found f(z)>=0 for using z as upper bound.
// While doing so, replace the lower bound by z whenever
// f(z)<0. In the special case f(z)=0, we're done and
// can instantly return the root z.
double z = 1;
double fz = f(z);
while(fz < 0)
{
I[0] = z;
fI[0] = fz;
z = dom.projectIn(z*2);
fz = f(z);
}
if (fz == 0)
{
xi = z;
return;
}
I[1] = z;
fI[1] = fz;
}
else
fI[1] = f(I[1]);
// Start bisection with:
// * absolute error tolerance 0
// * relative error tolerance 0.2
// * residual tolerance 0
// * initial interval guess I
// * x_ref=0 as refenence value to determine the relative error
NewBisection bisection(0, 0.2);
xi = bisection.solve(f, I, fI, 0);
#endif
};
......
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