Message Boards Message Boards

Extending FindRoot

GROUPS:

The function FindRoot can only find one root; there is no way of finding multiple roots. User J. M. extended the FindRoot by the function FindAllCrossings by changing Stan Wagon's book Mathematica in Action slightly.

The implementation is pretty slick. You plot the function on a certain range and tell the Plot function to mesh the zeros, then you extract the zeros and use FindRoot with them as an initial condition.

But this implementation has a major drawback; it fails in cases where the function doesn't cross zero, as x^2 for example.

To extend J. M. implementation I added some more lines of code, and I've made the code more verbose. The new lines of code calculate the zeros of the derivative of the function and then test it is close to zero and later refine it with FindRoot, this works even with functions as Abs[x], which has discontinuous derivative and never crosses zero, but thanks to the interpolated result, it takes care of that.

Options@FindRoots = Sort@Join[Options@FindRoot, {MaxRecursion -> Automatic, PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic, Debug -> False, ZeroTolerance -> 10^-2}];

FindRoots[fun_, {var_, min_, max_}, opts:OptionsPattern[]] := Module[{PlotRules, RootRules, g, g2, pts, pts2, lpts, F, sol},
    (* Extract the Options *)
    PlotRules = Sequence @@ FilterRules[Join[{opts}, Options@FindRoots], Options@Plot];
    RootRules = Sequence @@ FilterRules[Join[{opts}, Options@FindRoots], Options@FindRoot];

    (* Plot the function and "mesh" the point with y-coordinate 0 *)
    g = Normal@Plot[fun, {var, min, max}, MeshFunctions -> (#2 &), Mesh -> {{0}}, Method -> Automatic, Evaluate@PlotRules];

    (* Get the meshes zeros *)
    pts = Cases[g, Point[p_] :> SetPrecision[p[[1]], OptionValue@WorkingPrecision], Infinity];
    (* Get all plot points *)
    lpts = Join@@Cases[g, Line[p_] :> SetPrecision[p, OptionValue@WorkingPrecision], Infinity];

    (* Derive the interpolated data to find other zeros *)
    F = Interpolation[lpts, InterpolationOrder->2];
    g2 = Normal@Plot[Evaluate@D[F@var, var], {var, min, max}, MeshFunctions -> (#2 &), Mesh -> {{0}}, Method -> Automatic, Evaluate@PlotRules];

    (* Get the meshes zeros and retain only small ones *)
    pts2 = Cases[g2, Point[p_] :> SetPrecision[p[[1]], OptionValue@WorkingPrecision], Infinity];
    pts2 = Select[pts2, Abs[F@#] < OptionValue@ZeroTolerance &];
    pts = Join[pts, pts2]; (* Join all zeros *)

    (* Refine zeros by passing each point through FindRoot *)
    If[Length@pts > 0,
       pts = Map[FindRoot[fun, {var, #}, Evaluate@RootRules]&, pts];
       sol = Union@Select[pts, min <= Last@Last@# <= max &];
       (* For debug purposes *)
       If[OptionValue@Debug, Print@Show[g, Graphics@{PointSize@0.02, Red, Point[{var, fun} /. sol]}]];
       sol
    ,
       If[OptionValue@Debug, Print@g];
       {}
    ]
]

Example:

enter image description here

My primary use of this function is to find the eigenvalues of dielectric waveguides. They are found as the zeros of the boundary conditions.

POSTED BY: Thales Fernandes
Answer
5 months ago

Another use, finding zeroes of Zeta-function:

y /. FindRoots[Abs@Zeta[1/2 + I y], {y, 0, 100}, PlotPoints -> 500] // Quiet
Im@*ZetaZero /@ Range@Length@% // N;
(* Output *)
{14.1347, 21.022, 25.0109, 30.4249, 32.9351, 37.5862, 40.9187, \
43.3271, 48.0052, 49.7738, 52.9703, 56.4462, 59.347, 60.8318, \
65.1125, 67.0798, 69.5464, 72.0672, 75.7047, 77.1448, 79.3374, \
82.9104, 84.7355, 87.4253, 88.8091, 92.4919, 94.6513, 95.8706, \
98.8312}
(* Calculate error *)
%% - % // Abs // Mean;
(* Output *)
1.58635*10^-8
POSTED BY: Thales Fernandes
Answer
5 months ago

There is a 2006 package by Ted Ersek, RootSearch, which has a similar purpose. I haven't tried it lately to see if it still works. There may be more recent versions somewhere.

POSTED BY: Gianluca Gorni
Answer
4 months ago

Since Thales wanted an extension of the improved FindAllCrossings[] to roots that are also extrema, here is an implementation based on PeakDetect[] to find the non-crossing roots:

Options[FindAllRoots] = Sort[Join[Options[FindRoot], {MaxRecursion -> Automatic, PerformanceGoal :> $PerformanceGoal,
                                  PlotPoints -> Automatic, Tolerance -> Automatic}]];

FindAllRoots[f_, {t_, a_, b_}, opts : OptionsPattern[]] := 
    Module[{op, plt, prec, r, r1, r2, tol, v, xa, ya},
           prec = OptionValue[WorkingPrecision];
           plt = Normal[Plot[f, {t, a, b}, MeshFunctions -> (#2 &), Mesh -> {{0}}, Method -> Automatic, 
                             Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllRoots]], Options[Plot]]]]];
           r1 = Cases[plt, Point[p_] :> SetPrecision[p[[1]], prec], Infinity];
           {xa, ya} = Transpose[First[Cases[plt, Line[l_] :> l, Infinity]]];
           op = Join[Pick[xa, PeakDetect[-ya], 1], Pick[xa, PeakDetect[ya], 1]];
           v = Function[t, f] /@ op;
           tol = OptionValue[Tolerance]; If[tol === Automatic, tol = Max[v] $MachineEpsilon^(1/3)];
           r2 = SetPrecision[Pick[op, Thread[Abs[v] <= tol]], prec];
           r = Join[r1, r2];
           If[r =!= {}, 
              Union[Select[t /. Map[FindRoot[f, {t, #}, Evaluate[Sequence @@ FilterRules[Join[{opts},
                                             Options[FindAllRoots]], Options[FindRoot]]]] &,
                                    r], a <= # <= b &]], {}]]

This should work on Thales's original example.

POSTED BY: J. M.
Answer
4 months ago

enter image description here - Congratulations! This post is now a Staff Pick as distinguished on your profile! Thank you, keep it coming!

POSTED BY: Moderation Team
Answer
4 months ago

Nice. To me it seems more like an extension of NSolve:

NSolve[Abs[x] Sin[Abs[x]]^2 == 0 && -5 < x < 5, x]
(*  {{x -> 0}, {x -> -3.14159}, {x -> 3.14159}}  *)

For some reason Solve/NSolve won't strip the trivial Abs from Abs[Zeta[...]] == 0, so human thought is needed:

Solve[Zeta[1/2 + I y] == 0 && 0 < Re@y < 100 && Im[y] == 0, y, GeneratedParameters -> c]

Table[
 Chop[N@First@%, $MinMachineNumber],
 Evaluate@First@Cases[%, _[a_, ___, x_c, ___, b_] :> {x, a, b}, Infinity]]
(*
{{y -> ConditionalExpression[-(1/2) I (-1 + 2 ZetaZero[c[1]]), 
    c[1] \[Element] Integers && 1 <= c[1] <= 29]}}

{{y -> 14.1347}, {y -> 21.022}, {y -> 25.0109}, {y -> 30.4249}, {y -> 32.9351},
 {y -> 37.5862}, {y -> 40.9187}, {y -> 43.3271}, {y -> 48.0052}, {y -> 49.7738},
 {y -> 52.9703}, {y -> 56.4462}, {y -> 59.347}, {y -> 60.8318}, {y -> 65.1125},
 {y -> 67.0798}, {y -> 69.5464}, {y -> 72.0672}, {y -> 75.7047}, {y -> 77.1448},
 {y -> 79.3374}, {y -> 82.9104}, {y -> 84.7355}, {y -> 87.4253}, {y -> 88.8091},
 {y -> 92.4919}, {y -> 94.6513}, {y -> 95.8706}, {y -> 98.8312}}
*)
POSTED BY: Michael Rogers
Answer
4 months ago

For your first code, your NSolve only returns the zero solution. And the second code is kind of a cheat since Mathematica has implemented ZetaZero.

NSolve isn't the best solution, try to run the following code:

Solve[Zeta[1/2 + I y] Cos[y] == 0 && 0 < Re@y < 100 && Im[y] == 0, y, GeneratedParameters -> c] // AbsoluteTiming

The solution is pretty straightforward. I don't know if Mathematica is able to solve it, but my lack of patience had my abort the code.

NSolve and Solve are too general, hence they are inherently slow.

POSTED BY: Thales Fernandes
Answer
4 months ago

My point, fwiw, was that FindRoots seems more an extension of NSolve than of FindRoot, based on the output. I think it's clear from my example (i.e., your original Zeta example) that FindRoots outperforms Solve/NSolve on the zeta function, but my intention was not to compare them. (Yes, it's unfortunate that it can rely on ZetaZero, since there are other methods it can use for analytic functions, which might work.)

(For my first code, I get all three solutions in V9.0.1, so you must be using a different version.)

POSTED BY: Michael Rogers
Answer
4 months ago

Hi Fernandes,

Thanks very nice and useful code.

Can I ask you a further question: can you also specified the the dots as example: y-5 or y-10 line with x -10,10, Currently y=0 in your code. My intention is also to count the dots with the intersection of the plot?

I have tried but without success !

f[x_?NumericQ] := Abs[x] Sin[Abs@x]^2
count = FindRoots[f[x*2], {x, -10, 10}
   , Debug -> True
   , ImageSize -> 300] // Quiet
Length[count]

{{x -> -9.42478}, {x -> -7.85398}, {x -> -6.28319}, {x -> -4.71239}, \ {x -> -3.14159}, {x -> -1.5708}, {x -> 2.72823*10^-8}, {x -> 1.5708}, {x -> 3.14159}, {x -> 4.71239}, {x -> 6.28319}, {x -> 7.85398}, {x -> 9.42478}}

And the length is13 at y=0

I hope to understand my question. Thanks.

enter image description here

POSTED BY: Jos Klaps
Answer
4 months ago

Just shift your function. y - y0.

POSTED BY: Thales Fernandes
Answer
4 months ago

Hi Fernandes,

Thanks for your quick reply. But can you give me some more details where to shift the function in the code?

Best Regards,...Jos

POSTED BY: Jos Klaps
Answer
4 months ago
FindRoots[f[x*2] - 5, {x, -10, 10}, Debug -> True, ImageSize -> 300] // Quiet

Just shift the y-axis, it has the same desirable result. $y=y_0$ can be modified to $y-y_0=0$.

POSTED BY: Thales Fernandes
Answer
4 months ago

Yes, Ted Ersek's RootSearch package is GREAT if you're looking for all possible roots along a line (whether they cross the line or not) - a common application that I don't think that WRI addresses. I have a 2010 version and tried out all the sample problems and they all worked.

POSTED BY: David J M Park Jr
Answer
4 months ago

To All,

Thanks for your excellent support !

Best Regards,....Jos

POSTED BY: Jos Klaps
Answer
4 months ago

Good demo. I should note, however, that in practice (i.e. part of what ZetaZero[] does under the hood) one uses the function RiemannSiegelZ[] for finding zeroes of Riemann $\zeta$; this has the advantage of the function (supposedly) actually crossing the axis instead of just touching it.

POSTED BY: J. M.
Answer
4 months ago

Group Abstract Group Abstract