# Extending FindRoot

GROUPS:
 Thales Fernandes 5 Votes 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:My primary use of this function is to find the eigenvalues of dielectric waveguides. They are found as the zeros of the boundary conditions. Answer 9 months ago 14 Replies  Thales Fernandes 2 Votes 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  Answer 9 months ago  - Congratulations! This post is now a Staff Pick as distinguished on your profile! Thank you, keep it coming! Answer 8 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}} *) 
8 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.
8 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.)
8 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=0I hope to understand my question. Thanks.
8 months ago
 Just shift your function. y - y0.
8 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
8 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$.
8 months ago
 Gianluca Gorni 1 Vote 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.
8 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.
 J. M. 1 Vote 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.
 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.