3
$\begingroup$

Is there an elegant way to nest a function in a similar fashion to tree graphs; i.e., create and follow branches until some criterion is met (until we reach a leaf)? An example:

We want to numerically find all roots of some function within a specific interval, for instance all roots (between -4 and 4) of the polynomial

poly[x_] := (x - 3) (x - 3.2) (x - 2) (x - 1) (x) (x + 1) (x + 2) (x + 3) 

We search for a root via

FindRoot[poly@x, {x ,0 ,-4, 4}]
(* {x -> 0.} *)

We now split our search interval into two intervals [-4,0] and [0,-4] and search again for roots

FindRoots[poly@x, {x,-2, -4, 0}]
FindRoots[poly@x, {x, 2, 0, 4}]
(* {x->-2}*)
(* {x->2}*)

Repeating the search within [-4, -2], [-2, 0], [0, 2] and [2, 4] yields the additional roots: -3, -1 , 1, 3

Yet another round reveals the additional root at 3.2 via FindRoot[poly@x, {x, 3.5, 3, 4 }].

If we now check the new intervals [2, 3.2] and [3.2, 4] we don't get any new roots we didn't already know. To look for roots that we might have missed due to picking the wrong initial search values we could now vary all search intervals for instance by excluding already found roots, e.g searching for roots between [3.21, 3.99] which yields a warning from FindRoot indicating that there are no more roots inside the interval. (In this case there are no additional roots inside the other intervals as well)

FindRoot[poly@x, {x, 3.5, 3.21, 3.99}]
(* The point {3.21} is at the edge of the search region {3.21,3.99} in coordinate 1 and the computed search direction points outside the region. >>*)

How does one implement such an algorithm in an elegant way?

$\endgroup$
  • $\begingroup$ Do you know that you can use Roots[poly[x] == 0, x]? $\endgroup$ – C. E. Apr 4 '15 at 21:42
  • $\begingroup$ Finding roots was actually just the first example that came to my mind to illustrate the question (and Roots would only work for polynomials and not for any arbitrary functions ) $\endgroup$ – Sascha Apr 4 '15 at 21:55
1
$\begingroup$

Normal recursion is about as elegant as it's gonna get, I believe:

poly[x_] := (x - 3) (x - 3.2) (x - 2) (x - 1) (x) (x + 1) (x + 2) (x + 3)
sol[solutions_, {min_, max_}] := Module[{root, eps = 0.1},
  root = Quiet@Check[x /. FindRoot[poly[x], {x, (max + min)/2, min, max}], $Failed];
      Flatten@If[root === $Failed,
    solutions,
    Union[solutions, {
      root, sol[solutions, {root + eps, max}], 
      sol[solutions, {min, root - eps}]
      }]
    ]
  ]
sol[{}, {-4, 4}]

{0., -2., -3., -1., 2., 1., 3., 3.2}

$\endgroup$
  • $\begingroup$ Any idea why this fails when the interval is larger e.g. sol[{}, {-100, 100}] (*{0., 3.2, 3., -3., -1., -2.}*) ? $\endgroup$ – Sascha Apr 6 '15 at 19:50
  • $\begingroup$ @Sascha Probably because (max + min)/2 is not always a good starting point for FindRoot. $\endgroup$ – C. E. Apr 6 '15 at 20:03

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.