Handle case x = 0 in the modified gradient method

src/bisection-example-flexible.cc

@@ -49,6 +49,8 @@ class SampleFunctional {
 
   SmallVector minimise(const SmallVector x, unsigned int iterations) const {
     SmallVector descDir = ModifiedGradient(x);
+    if (descDir == SmallVector(0.0))
+      return SmallVector(0.0);
 
     Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
     Dune::dverb << "Minimizing in direction w with dJ(x,w) = "
@@ -129,9 +131,16 @@ class SampleFunctional {
   }
 
   SmallVector ModifiedGradient(const SmallVector x) const {
-    if (x == SmallVector(0.0))
-      // TODO
-      DUNE_THROW(Dune::Exception, "The case x = 0 is not yet handled.");
+    if (x == SmallVector(0.0)) {
+      SmallVector d = SmoothGrad(x);
+      // Decline of the smooth part in the negative gradient direction
+      double smoothDecline = -(d * d);
+      double nonlinearDecline =
+          func_.rightDifferential(0.0) * d.two_norm(); // TODO: is this correct?
+      double combinedDecline = smoothDecline + nonlinearDecline;
+
+      return (combinedDecline < 0) ? d : SmallVector(0.0);
+    }
 
     SmallVector const pg = PlusGrad(x);
     SmallVector const mg = MinusGrad(x);