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?

Improve this question
$\endgroup$
6
  • $\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$ Commented Jun 17, 2017 at 21:18
  • $\begingroup$ Thanks for the reply! I've edited my question to include the function. $\endgroup$ Commented Jun 17, 2017 at 22:40
  • 1
    $\begingroup$ Your expression is complex at x->1.4. $\endgroup$ Commented Jun 17, 2017 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$ Commented Jun 17, 2017 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$ Commented Jun 17, 2017 at 23:34

2 Answers 2

Reset to default
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

Improve this answer
$\endgroup$
1
  • $\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$ Commented Jun 18, 2017 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}}

Improve this answer
$\endgroup$
1
  • $\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$ Commented Jun 18, 2017 at 0:29

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.