16
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ When I used NDSolve on the question you cited in your first line, I too encountered a problem with RegionPlot, although I am uncertain whether it was the same problem. For this reason, I also used ContourPlot. In any case, the problem went away, when I used ParameterNDSolve. $\endgroup$ – bbgodfrey Dec 25 '14 at 21:23
  • 2
    $\begingroup$ No (obvious) problem on version 9.0.1, with no messages produced after about 10 minutes of run-time. So, I think it looks like a possible bug in version 10. I don't have that installed on this computer, I'm not feeling patient enough to let version 9 run to completion, and also not much in the mood for troubleshooting, so this is of course not conclusive, but clearly we should not be seeing FindMinimum's throw come back to the top level. I'd report it, if I were you. $\endgroup$ – Oleksandr R. Dec 25 '14 at 23:10
  • 2
    $\begingroup$ This is most certainly a bug. I can only emphasize what Olek mentioned at the end: A Throw inside a built-in function should never reach the user. This would be very bad style. Btw, I get the same output as you did on Linux 64 with Mathematica 10.0.2. $\endgroup$ – halirutan Dec 25 '14 at 23:20
  • $\begingroup$ As suggested in my previous comment, I had hoped that using ParametricNDSolve with RegionPlot might provide better results. It did not, instead producing much the same errors described above. $\endgroup$ – bbgodfrey Dec 26 '14 at 0:39
  • 1
    $\begingroup$ @OleksandrR. Sorry about that - my original lhmin only rarely called FindMinimum. I didn't think about it actually working and running for a good long time. As far as I can tell, on V10, it fails in the first few calls. Thanks for checking V9. P.S. I have reported it and will update when I hear back. $\endgroup$ – Michael E2 Dec 26 '14 at 1:06
6
$\begingroup$

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).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.