5
$\begingroup$

I have an equation of the form f(x) = 0 which I need to solve numerically because f is absolutely horrible. I'm using FindRoot to achieve my goal. Here's my code:

FindRoot[
(43.76760935027122*(ArcCos[(3.6958544099301736-0.09097576027495978*x^2-1.*x^4)/(3.695854409930174-0.30306895456563354*x^2-1.*x^4)]-1.*ArcTan[(1223.7099426464301*x^4-725.4797933949869*x^6-69.76907700248954*x^8+105.90945888977912*Sqrt[x^4*(22.772281766782058-51.683884026721294*x^2+43.76287298560037*x^4-16.38570295772005*x^6+2.2889328165770233*x^8)]+x^2*(507.6064341807468-45.96118379610641*Sqrt[x^4*(22.772281766782058-51.683884026721294*x^2+43.76287298560037*x^4-16.38570295772005*x^6+2.2889328165770233*x^8)]))/(x*(16.765159281029014-4851.4404675526785*Sqrt[x^4*(3.157334518321-3.553778*x^2+x^4)*(3.151032666129048-3.60488687607268*x^2+1.*x^4)]+x*(1.8755408291095235*^-15+x*(33.141011104612154+x*(6.514741496419986*^-15+x*(-41.50048368818499+x*(3.2465458701698195*^-15+x*(8.093976+x*(2.2903016461860367*^-16+1.*x)))))))))]))/x,
{x, 1.4, 1, 2}]

(I told you it was horrible.)

When I try to run this, I get the following error:

FindRoot::reged: The point {1.4} is at the edge of the search region {1.,2.} in coordinate 1 and the computed search direction points outside the region. >>

Obviously, 1.4 is not outside the region [1,2]. Furthermore, I know that the solution is 1.553, which is not outside the region [1,2] either. I know this because if I play around with the starting value I get {x -> 1.553} as the output if the starting value is close enough to that value and an error if it is not; I've also plotted f in the interval [1.5,1.6], and it looks like the root is indeed around 1.553. What's going on?

$\endgroup$
  • $\begingroup$ It works for me. Maybe it's your f? Hard to say without the code. Please post a complete minimal working example (MWE). $\endgroup$ – Michael E2 Jun 17 '17 at 21:18
  • $\begingroup$ Thanks for the reply! I've edited my question to include the function. $\endgroup$ – Rain Jun 17 '17 at 22:40
  • 1
    $\begingroup$ Your expression is complex at x->1.4. $\endgroup$ – John Doty Jun 17 '17 at 23:13
  • $\begingroup$ That's true. Actually, it's complex all the way up to some number between 1.454 and 1.455, at which point it becomes real. Still, though, why does Mathematica say 1.4 is at the edge of [1,2] instead of saying f is complex at 1.4? Sounds like the error message is poorly chosen. $\endgroup$ – Rain Jun 17 '17 at 23:26
  • 1
    $\begingroup$ The error message is rather inscrutable. The part about the "computed search direction" is a hint. When the expression is complex, the search wanders off into the complex plane. That, of course, puts it outside your interval. $\endgroup$ – John Doty Jun 17 '17 at 23:34
5
$\begingroup$

The Newton update step at x == 1.4 is

-fn/D[fn, x] /. x -> 1.4
(*  0.0424827 + 0.0445837 I  *)

This points outside the search region which is the line segment, or interval, from 1 to 2.

One way to fix it is to give FindRoot a rectangle in the complex plane for the search region:

FindRoot[fn = fn, {x, 1.4, 1 - I, 2 + I}]
(*  {x -> 1.55322 - 2.69371*10^-17 I}  *)

The negligible imaginary part can be deleted with Chop@FindRoot[...].

That's assuming you do not want to use the information you have at hand, namely that the root is approximately x -> 1.553.

(It's always confusing that Mathematica tends to work in the complex plane, but the error message is accurate. It assumes the user realizes that the search is in the complex plane when the function value is complex.)

Visualizations:

Objective function:

Plot[fn // ReIm // Evaluate, {x, 1, 2}]

Mathematica graphics

The discontinuities are at

excl = Simplify`FunctionSingularities[fn, {x}, {"ALL"}];

DeleteDuplicates@
 Flatten[Cases[excl, eq_Equal :> NSolve[eq && 1 < x < 2, x], 
   Infinity], 1]
(*
  {{x -> 1.333}, {x -> 1.333}, {x -> 1.22055}, {x -> 1.22056}, {x -> 1.333},
   {x -> 1.45435}, {x -> 1.22056}, {x -> 1.333}, {x -> 1.45435}, {x -> 1.35147}}
*)

This suggests you should not start your root-search below x == 1.455 or so.

The FindRoot search:

{steps} = ReIm@Last@Reap[
    Sow[1.4];
    FindRoot[fn, {x, 1.4, 1 - I, 2 + I}, StepMonitor :> Sow[x]]
    ]
(*
  {{{1.4, 0}, {1.44248, 0.0445837}, {1.49397, 0.0689971}, {1.55744, 0.0558957},
   {1.57813, -0.000239207}, {1.54782, 0.000107546}, {1.55298, 9.27012*10^-6},
   {1.55322, 3.52565*10^-8}, {1.55322, 2.54977*10^-13}, {1.55322, -2.69371*10^-17}}}
*)

Graphics[
 {Line[steps], Point[steps], 
  MapIndexed[Text[#2, #1, {0, -1.5}] &, steps]},
 Axes -> True, AxesLabel -> {HoldForm[x], HoldForm[I y]}]

Mathematica graphics

The small steps at the end are hard to see without some sort of log-scale around the root. A 3D graphics shows that last few steps seems to line up as the reduction in the error of the root is unseeable.

plot = With[{steps3D = MapIndexed[Join, steps]},
  Graphics3D[
   {Thickness[Medium], Darker@Blue, Line[steps3D], Point[steps3D], 
    MapIndexed[Text[#2, #1, {1.3, -1.1}] &, steps3D]},
   Axes -> True, AxesLabel -> {HoldForm[x], HoldForm[I y], "step"},
   BoxRatios -> {1.6, 1, 0.6}]
  ]

GraphicsRow[Show[plot, ViewPoint -> #] & /@ {Top, Front}]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Thanks! This takes care of the problem. I'll mark this as the answer because it effectively gets rid of the problem and also explains why it was there in the first place. $\endgroup$ – Rain Jun 18 '17 at 0:50
0
$\begingroup$

Instead of using FindRoot with a 4-arg second argument, you could use NSolve and add a bracketing condition. Let f[x] be your function:

f[x_] := (43.76760935027122*(ArcCos[(3.6958544099301736-0.09097576027495978*x^2-1.*x^4)/(3.695854409930174-0.30306895456563354*x^2-1.*x^4)]-1.*ArcTan[(1223.7099426464301*x^4-725.4797933949869*x^6-69.76907700248954*x^8+105.90945888977912*Sqrt[x^4*(22.772281766782058-51.683884026721294*x^2+43.76287298560037*x^4-16.38570295772005*x^6+2.2889328165770233*x^8)]+x^2*(507.6064341807468-45.96118379610641*Sqrt[x^4*(22.772281766782058-51.683884026721294*x^2+43.76287298560037*x^4-16.38570295772005*x^6+2.2889328165770233*x^8)]))/(x*(16.765159281029014-4851.4404675526785*Sqrt[x^4*(3.157334518321-3.553778*x^2+x^4)*(3.151032666129048-3.60488687607268*x^2+1.*x^4)]+x*(1.8755408291095235*^-15+x*(33.141011104612154+x*(6.514741496419986*^-15+x*(-41.50048368818499+x*(3.2465458701698195*^-15+x*(8.093976+x*(2.2903016461860367*^-16+1.*x)))))))))]))/x

Then, you can use:

NSolve[f[x]==0 && 1 < x < 2, x]

{{x -> 1.55322}}

$\endgroup$
  • $\begingroup$ Thanks! I gave this a try. Apparently, NSolve takes more than 10 times longer than FindRoot to compute, so when I have 30,000 repetitions it's going to take a lot of time to complete. Still, it seems to produce the right answer, which is nice. $\endgroup$ – Rain Jun 18 '17 at 0:29

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.