# 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 1 year 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 1 year 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. Answer 1 year ago  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. Answer 1 year ago  - Congratulations! This post is now a Staff Pick as distinguished on your profile! Thank you, keep it coming! Answer 1 year 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}} *) 
1 year 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.
1 year 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.)
1 year 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.
1 year ago
 Just shift your function. y - y0.
1 year 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
1 year 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$.
 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.