I use Mathematica 9 and am having trouble understanding 1-dimensional FindRoot, which ought to be using Brent’s interval method. Problems arise when my real functions have involved Undefined as an auxiliary.
For example,
step := If[x < 1, -1, 1] + x/7;
stepa := If[NumberQ[If[x < 1, 0, Undefined]], -1, 1] + x/7;
stepb := If[NumberQ[If[x > 1, 0, Undefined]], 1, -1] + x/7;
all define the same real function and give the same Plots, but
{FindRoot[step, {x, 0, 2}],
FindRoot[stepa, {x, 0, 2}],
FindRoot[stepb, {x, 0, 2}]}
yields 3 different answers, the last two not being in the initial interval [0,2]:
{{x → 1.}, {x → -7.}, {x → 7.}}
Can anyone help? Howand why does Mathematica use the innards of the function definitions?
Mathematica uses the innards of function definitions to decide which method to use, whether or not to do symbolic or numerical differentiation, etc. In the section “Some Notes on Internal Implementation” the description of FindRoot gives 3 methods, damped Newton’s method, secant method and Brent’s method.
Frank,
Thanks for your quick response. As I understood things, Mathematica should be using Brent’s method because I gave FindRoot an interval of starting values. It ought to be pretty infallible if the values at the ends of the interval have opposite signs and it certainly should produce a result in the interval. What goes wrong here?
Is it possible that under some conditions Mathematica might not use Brent’s method even when two starting values are given? I found the same problem even when I specified: Method->“Brent”. Can you point me to documentation on this? What I found under “Some Notes on Internal Implementation” didn’t suggest such a possibility.
Michael,
Yes, an interval of starting methods should force Brent’s method. For further exploration, you might want to try defining your function as f[x_?NumberQ], rather than putting NumberQ inside the function.
Frank,
The functions I am really using don’t have NumberQ explicitly. I was just trying to produce the simplest functions to illustrate the problem. All these functions evaluate to reals without apparent problems and plot nicely. What I can’t understand is why Mathematica, which is supposed to be using Brent’s method, uses anything other than function evaluations and produces an answer outside the initial interval. Maybe we need a Mathematica insider to explain??
this works
In[87]:= Clear[step, stepa, stepb]
step[x_] := If[x < 1, -1, 1] + x/7;
stepa[x_?NumericQ] := If[NumberQ[If[x < 1, 0, Undefined]], -1, 1] + x/7;
stepb[x_?NumericQ] := If[NumberQ[If[x > 1, 0, Undefined]], 1, -1] + x/7;
In[78]:= FindRoot[step[x], {x, 0, 2}]
During evaluation of In[78]:= FindRoot::brmp: The root has been bracketed as closely as possible with machine precision but the function value exceeds the absolute tolerance 1.0536712127723497`*^-8. >>
Out[78]= {x -> 1.}
In[93]:= FindRoot[stepa[x], {x, 0, 2}]
During evaluation of In[93]:= FindRoot::brmp: The root has been bracketed as closely as possible with machine precision but the function value exceeds the absolute tolerance 1.0536712127723497`*^-8. >>
Out[93]= {x -> 1.}
In[94]:= FindRoot[stepb[x], {x, 0, 2}]
During evaluation of In[94]:= FindRoot::brmp: The root has been bracketed as closely as possible with machine precision but the function value exceeds the absolute tolerance 1.0536712127723497`*^-8. >>
Out[94]= {x -> 1.}
You could use EvaluationMonitor to see where FindRoot is going for hte different cases.
In[10]:= step[x_] := If[x < 1, -1, 1] + x/7;
In[11]:= FindRoot[step[x], {x, 0, 2}, EvaluationMonitor :> Print[x]]
During evaluation of In[11]:= 0.
During evaluation of In[11]:= 2.
During evaluation of In[11]:= 0.875
During evaluation of In[11]:= 1.4375
During evaluation of In[11]:= 1.11159
During evaluation of In[11]:= 0.976787
During evaluation of In[11]:= 1.04419
During evaluation of In[11]:= 1.00565
During evaluation of In[11]:= 0.989177
During evaluation of In[11]:= 0.997412
During evaluation of In[11]:= 1.00153
During evaluation of In[11]:= 0.999176
During evaluation of In[11]:= 1.00035
During evaluation of In[11]:= 0.99968
During evaluation of In[11]:= 1.00002
During evaluation of In[11]:= 0.999824
During evaluation of In[11]:= 0.999921
During evaluation of In[11]:= 0.999969
During evaluation of In[11]:= 0.999993
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 0.999998
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 0.999999
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= 1.
During evaluation of In[11]:= FindRoot::brmp: The root has been bracketed as closely as possible with machine precision but the function value exceeds the absolute tolerance 1.0536712127723497`*^-8. >>
Out[11]= {x -> 1.}
Udo, Thanks for ?NumericQ magic. I’ll try to understand what it’s doing; I haven’t used that before.
Frank, A good tip to use EvaluationMonitor. Thanks.
I’d like to understand why FindRoot doesn’t seem to do “what it says on the tin”, i.e., do a proper Brent method for an interval. Maybe you have given me the tools to do it.
x_?NumberQ and x_?NumericQ seem to deal with my actual problems nicely. Thank you guys for helping me with this.