Bug in plotting function using FindMinimum?

Question

Bug introduced in 10.0 and persisting through 11.0.1 or later

---

While working on an answer to RegionPlot from NDSolve, I came across this strange behavior of RegionPlot and ParametricPlot involving a function that sometimes called FindMinimum. Below is a simplified version of the code, and in it the function lhmin always calls FindMinimum.

I'm thinking it is probably a bug. However, I wonder if I'm breaking a rule somewhere, and something I'm doing is causing this apparent bug. The following is the complete code, run on a new kernel. I'm using V10.0.1 (on Mac OSX 10.9.5).

ode3[y_?NumericQ, M_?NumericQ] := Block[{t},
  First@NDSolve[{
     Dt[lh[t], t] == -ls[t]/2 - (lh[t] - 1)^3/4,
     Dt[ls[t], t] == lh[t]/2 - ls[t]^3,
     lh[0] == y, ls[0] == M}, {lh, ls}, {t, 0, 40}]]

lhmin[y0_?NumericQ, M0_?NumericQ] := With[{sol = ode3[y0, M0]},
   With[{values = lh["ValuesOnGrid"] /. sol},
    Print[{y0, M0}];
    With[{minPos = First@Ordering[values, 1]},
     First@FindMinimum[{lh[t] /. sol, {0 <= t <= 40}},
       {t, (lh["Coordinates"] /. sol)[[1, minPos]]}]
     ]]
   ];

Note that lhmin prints its input before calling FindMinimum. You will note below that it does not produce a slew of output, because the whole thing pretty quickly comes to a halt because of an uncaught Throw.

RegionPlot. In this plot, lhmin is used to define the region. There are messages about Break[] being abused and some delinquent function acting up and throwing things.

RegionPlot[lhmin[y, M] > 0, {y, 0.1, 1.}, {M, 0, 10}]

{0.10004741578947368`, 0.0005268421052631579`}  (* output of Print *)
{0.10000004736842105`, 5.263157894736842`*^-7}
{0.10000004736842105`, 5.263157894736842`*^-7}

Break::nofwd: No enclosing For, While, or Do found for Break[]. >>

Break::nofwd: No enclosing For, While, or Do found for Break[]. >>

Break::nofwd: No enclosing For, While, or Do found for Break[]. >>

General::stop: Further output of Break::nofwd will be suppressed during this calculation. >>

FindMinimum::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {0.0423358,2.58505*10^-10,0.0211679}, is returned. >>

Throw::nocatch: Uncaught Throw[{{8.34096},{0.0635038,{0.0423358,2.58505*10^-10,0.0211679}},500,NOTCONVERGED},FindMinimum`InteriorPoint] returned to top level. >>

Throw::nocatch: Uncaught Throw[Null,FindMinimum`InteriorPoint] returned to top level. >>

(* Hold[Throw[{{8.34096}, {0.0635038, {0.0423358, 2.58505*10^-10, 
       0.0211679}}, 500, "NOTCONVERGED"}, FindMinimum`InteriorPoint], 
    Throw[Null, FindMinimum`InteriorPoint]] *)

However, note that evaluating the function lhmin on the point indicated works fine and without generating messages.

lhmin @@ {0.10000004736842105`, 5.263157896957288`*^-7}

{0.1,5.26316*10^-7}   (* output of Print *)

(*  0.0913811  *)

ParametricPlot. I used lhmin in the RegionFunction to ParametricPlot, and it generates the same messages as above except no FindMinimum::eit, perhaps because FindMinimum did not have a convergence problem on the particular input {0.1`, 0.`}. The output is different, but that's probably not significant, since it's a different function.

ParametricPlot[{y, M}, {y, 0.1, 1.}, {M, 0, 10}, 
  RegionFunction -> (lhmin[#1, #2] > 0 &)]

{0.1`, 0.`}   (* output of Print *)

Break::nofwd: No enclosing For, While, or Do found for Break[]. >>
...
Throw::nocatch: Uncaught Throw[{{8.34095},{0.0635036,{0.0423357,3.05593*10^-10,0.0211679}},500,NOTCONVERGED},FindMinimum`InteriorPoint] returned to top level. >>

Throw::nocatch: Uncaught Throw[Null,FindMinimum`InteriorPoint] returned to top level. >>

Hold[Break[], Break[], <<7499>>, Throw[Null, FindMinimum`InteriorPoint]]

Again evaluating the function lhmin on the point indicated works fine and without generating messages.

lhmin @@ {0.1`, 0.`}

{0.1,0.}

(*  0.0913811  *)

Summary. On the one hand lhmin and FindMinimum work without a problem when called by themselves. On the other hand, when lhmin is called from within RegionPlot or ParametricPlot, FindMinimum seems to generate fatal errors (on the same input values) that kill the program. Is this a bug, or am I doing something wrong?

---

Confirmed by WRI.

Answer 1

As mentioned in the question, this is a bug.

For a possible workaround, try

RegionPlot[lhmin[y, M] > 0, {y, 0.1, 1.}, {M, 0, 10}, "NumericalFunction" -> False]

Some related questions: (1), (2), (3).