52
$\begingroup$

I have some questions for multiroot search for transcendental equations. Is there any clever solution to find all the roots for a transcendental equation in a specific range?

Perhaps FindRoot is the most efficient way to solve transcendental equations, but it only gives one root around a specific value. For example,

FindRoot[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]], {x, 10}]

Of course, one can first Plot the equation and then choose several start values around each root and then use FindRoot to get the exact value.

  1. Is there any elegant way to find all the roots at once?

    Actually, I come up with this question when I solve the eigenequation for optical waveguides and I want to get the dispersion relation. I find ContourPlot is very useful to get the curve of the dispersion relation. For example,

    ContourPlot[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin@Sin[a*x] == 0, {x, 0, 10}, {a, 0, 4}]

    You can get

    enter image description here

  2. Is there any elegant way to get all the values in the ContourPlot for x when a==0 ?

  3. Is it possible to know how the ContourPlot gets all the points shown in the figure? Perhaps we can harness it to get all the roots for the transcendental equation.

Improve this question
$\endgroup$
13
  • 2
    $\begingroup$ You might be interested in this question. Anyway, for finding the roots of a function of a single variable, you can directly use the output of Plot[] to find initial approximations for FindRoot[]. If you're interested in that approach, I can write up an answer. $\endgroup$ Commented May 17, 2012 at 8:09
  • $\begingroup$ @J.M. I'm expecting for your answer ~~ $\endgroup$ Commented May 17, 2012 at 8:53
  • 2
    $\begingroup$ A tiny comment on Reduce[]-based methods: all of these hinge on the fact that Reduce[] apparently knows quite a fair bit about the transcendental functions built within Mathematica. If, however, you are dealing with a black-box function that can only evaluate at numerical values (you can simulate this behavior with something like f[x_?NumericQ] := Haversine[Pi x]), then Reduce[] won't be able to do much. $\endgroup$ Commented May 17, 2012 at 9:15
  • $\begingroup$ @J.M. Yes, I have tried Reduce[], and I find it is very slow, especially for a much more complex equation. $\endgroup$ Commented May 17, 2012 at 9:50
  • 2
    $\begingroup$ @yulinlinyu It seems you wanted to accept a few answers. You can accept ONLY ONE of the answers (clicking the green tick beside an answer), which you find the most helpful. You can upvote EVERY answer you find helpful or downvote if you find an answer misleading. $\endgroup$ Commented May 18, 2012 at 1:57

9 Answers 9

Reset to default
47
$\begingroup$

Borrowing almost verbatim from a recent response about finding extrema, here is a method that is useful when your function is differentiable and hence can be "tracked" by NDSolve.

f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] In[191]:= zeros = Reap[soln = y[x] /. First[ NDSolve[{y'[x] == Evaluate[D[f[x], x]], y[10] == (f[10])}, y[x], {x, 10, 0}, Method -> {"EventLocator", "Event" -> y[x], "EventAction" :> Sow[{x, y[x]}]}]]][[2, 1]] During evaluation of In[191]:= NDSolve::mxst: Maximum number of 10000 steps reached at the point x == 1.5232626281716416`*^-124. >> Out[191]= {{9.39114, 8.98587*10^-16}, {6.32397, -3.53884*10^-16}, {3.03297, -8.45169*10^-13}, {0.886605, -4.02456*10^-15}} Plot[f[x], {x, 0, 10}, Epilog -> {PointSize[Medium], Red, Point[zeros]}]

function and roots

If it were a trickier function, one might use Method -> {"Projection", ...} to enforce the condition that y[x] is really the same as f[x]. This method may be useful in situations (if you can find them) where we have one function in one variable, and Reduce either cannot handle it or takes a long time to do so.

Addendum by J. M.

WhenEvent is now the documented way to include event detection in NDSolve, so using it along with the trick of specifying an empty list where the function should be, here's how to get a pile of zeroes:

f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] zeros = Reap[NDSolve[{y'[x] == D[f[x], x], WhenEvent[y[x] == 0, Sow[{x, y[x]}]], y[10] == f[10]}, {}, {x, 10, 0}]][[-1, 1]]; Plot[f[x], {x, 0, 10}, Epilog -> {PointSize[Medium], Red, Point[zeros]}]

function and roots

Improve this answer
$\endgroup$
6
  • $\begingroup$ @ Daniel Lichtblau, it's a very clever, efficient and general way to get all the roots and can be harnessed in many other cases. Thank you a lot. $\endgroup$ Commented May 18, 2012 at 1:39
  • 2
    $\begingroup$ Excellent! However, this method doesn't work in v10, possibly a bug. $\endgroup$ Commented Nov 25, 2014 at 4:00
  • 1
    $\begingroup$ @luyuwuli Thanks for pointing that out. I tracked this to an option setting change wherein the default became MaxSteps -> Infinity. You can restore the behavior shown above with an explicit MaxSteps -> 10000. I will ask in house about whether the change is for the better given that it can cause an example like this to hang. $\endgroup$ Commented Nov 25, 2014 at 20:44
  • 2
    $\begingroup$ @DanielLichtblau, this is a bug that I think is not related to the MaxSteps -> Automatic change - at least not directly. You can get this to work with NDSolve[{y'[x] == Evaluate[D[f[x], x]], y[10] == (f[10])}, y[x], {x, 10, 10^-16}, Method -> {"EventLocator", "Event" -> y[x], "EventAction" :> Sow[{x, y[x]}]}] where the 0 is replaced with a small number. $\endgroup$ Commented Nov 26, 2014 at 8:58
  • 1
    $\begingroup$ @Nigel I think for complex-valued functions one will need to track real and imaginary parts both as a function of a real parameter. $\endgroup$ Commented Jun 29, 2017 at 17:05
35
$\begingroup$

I might as well elaborate on my comment. Here is a modification of Stan Wagon's FindAllCrossings[] function (from his book Mathematica in Action, second edition) that uses Plot[] to generate the initial approximations to be subsequently polished by FindRoot[]:

Options[FindAllCrossings] = Sort[Join[Options[FindRoot], {MaxRecursion -> Automatic, PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic}]]; FindAllCrossings[f_, {t_, a_, b_}, opts___] := Module[{r, s, s1, ya}, {r, ya} = Transpose[First[Cases[Normal[ Plot[f, {t, a, b}, Method -> Automatic, Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings]], Options[Plot]]]]], Line[l_] :> l, Infinity]]]; s1 = Sign[ya]; If[ ! MemberQ[Abs[s1], 1], Return[{}]]; s = Times @@@ Partition[s1, 2, 1]; If[MemberQ[s, -1] || MemberQ[Take[s, {2, -2}], 0], Union[Join[Pick[r, s1, 0], Select[t /. Map[FindRoot[f, {t, r[[#]], r[[# + 1]]}, Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings]], Options[FindRoot]]]] &, Flatten[Position[s, -1]]], a <= # <= b &]]], {}]]

Try it out:

FindAllCrossings[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]], {x, 0, 100}, WorkingPrecision -> 20] {0.88660463531346207679, 3.0329660890136683539, 6.3239665137114782212, 9.3911434075850854017, 12.589067252797192964, 15.687789316501627036, 18.865248000326751595, 21.976728589463954937, 25.144727536576135544, 28.263111694495812775, 31.425621611972587076, 34.548333253230934213, 37.707250582575859710, 40.832929091244028281, 43.989309899299199529, 47.117149292753968158, 50.271642808648651326, 53.401126310508732375, 56.554160423425108364, 59.684936863656120488, 62.836808589969946989, 65.968628466404205710, 69.119552435670107405, 72.252232119042699356, 75.402368490602128759, 78.535768909563822046, 81.685240379776511855, 84.819253682565487033, 87.968156330842562535, 91.102697190752433240, 94.251107663305908556, 97.386107414041584404}

A different implementation involves the use of the MeshFunctions option of Plot[] to generate the seeds. Here's how it looks:

FindAllCrossings[f_, {t_, a_, b_}, opts : OptionsPattern[]] := Module[{r}, r = Cases[Normal[Plot[f, {t, a, b}, MeshFunctions -> (#2 &), Mesh -> {{0}}, Method -> Automatic, Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings]], Options[Plot]]]]], Point[p_] :> SetPrecision[p[[1]], OptionValue[WorkingPrecision]], Infinity]; If[r =!= {}, Union[Select[t /. Map[FindRoot[f, {t, #}, Evaluate[Sequence @@ FilterRules[Join[{opts}, Options[FindAllCrossings]], Options[FindRoot]]]] &, r], a <= # <= b &]], {}]]

This version might be a bit faster in some cases, but it no longer has the safety feature of root bracketing in the previous version.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Of course, this method won't find even-order multiple roots (i.e. tangencies to the horizontal axis), but almost any numerical method will have difficulty with them anyway. If you know for certain that there are tangencies within the range of interest, you might consider using FindMinimum[]/FindMaximum[] instead to find them. $\endgroup$ Commented May 17, 2012 at 8:49
  • $\begingroup$ Usage example here mathematica.stackexchange.com/a/94881/193 $\endgroup$ Commented Sep 17, 2015 at 16:09
  • $\begingroup$ @J. M. Nice answer. It seems this method with FindAllCrossings[] does not work for zeroes of complex functions. For example I get the error "First::first: "{} has a length of zero and no first element." with pure imaginary Sin[x]/Sin[I*x]. Can it be used for complex functions? $\endgroup$ Commented Jun 29, 2017 at 16:44
  • $\begingroup$ @Nigel, the function here would not work directly since Plot[] does not work for complex functions. One might try using FindAllCrossings2D[] on the real and imaginary parts, tho. $\endgroup$ Commented Jul 26, 2017 at 4:50
  • $\begingroup$ @J.M.'spersistentexhaustion maybe in the case of even-order multiple roots one could perturb the original function find the splitted roots, decrease the perturbation parameter by two and find them again (maybe repeat that one more time) then extrapolate (or just use the initial splitted roots and take the middle as an approximation). $\endgroup$ Commented Aug 16, 2022 at 16:38
23
$\begingroup$

One can use Solve as well, e.g. 

s = Solve[ BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[ Sin[x]] == 0 && 0 < x < 10, x] 
Solve::incs: Warning: Solve was unable to prove that the solution set found is complete. >> {{x -> Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 0.886604635313462076794393681674}]}, {x -> Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 3.03296608901366835385376172847}]}, {x -> Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 6.32396651371147786252003752922}]}, {x -> Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 9.39114340758508579766919382120}]}}

Solve similarly as Reduce cannot prove that the set of solutions is complete. It returns the result in the form of rules but we can show that they return the same set solutions e.g. : 

r = Reduce[ BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[ Sin[x]] == 0 && 0 < x < 10, x]; s[[All, 1, 2]] == List @@ r[[All, 2]] 
Reduce::incs: Warning: Reduce was unable to prove that the solution set found is complete. >> True

Edit

Defining 

f[x_, a_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[a x]]

we can visualise the graphs of functions f[x,a] making use of ParametricPlot3D (look at answers to this question), e.g.

Show[ ParametricPlot3D[ Evaluate[ Table[{x, a, f[x, a]}, {a, 0, 5}]], {x, 0, 10}, BoxRatios -> {10, 8, 4}], ParametricPlot3D[{x, 1, f[x, 1]}, {x, 0, 10}, PlotStyle->{Thick, Darker@Red}, BoxRatios -> {10, 8, 4}]

enter image description here

Red thick curve is the graph of f[x,1],

or we can make use of Plot3D as well in the following way :

Plot3D[ f[x, a], {x, 0, 10}, {a, 0, 4}, MeshFunctions -> {#1 &, #2 &, #3 &}, Mesh -> {9, 3, 5}, Filling -> 0, PlotPoints -> 200, MaxRecursion -> 5]

]

enter image description here

The option Filling ->0 makes an impression that the level f[x,a] == 0 is like the surface of water.

Improve this answer
$\endgroup$
19
$\begingroup$

For finding all roots in a given interval I use the following function:

Clear[findRoots] Options[findRoots] = Options[Reduce]; findRoots[gl_Equal, {x_, von_, bis_}, prec : (_Integer?Positive | MachinePrecision | Infinity) : MachinePrecision, wrap_: Identity, opts : OptionsPattern[]] := Module[{work, glp, vonp, bisp}, {glp, vonp, bisp} = {gl, von, bis} /. r_Real :> SetPrecision[r, prec]; work = wrap@Reduce[{glp, vonp <= x <= bisp}, opts]; work = {ToRules[work]}; If[prec === Infinity, work, N[work, prec]]];

Example 

bes[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] findRoots[bes[x] == 0, {x, 0, 100}] During evaluation of In[35]:= Reduce::incs: Warning: Reduce was unable to prove that the solution set found is complete. >> Out[39]= {{x -> 0.886605}, {x -> 3.03297}, {x -> 6.32397}, {x -> 9.39114}, {x -> 12.5891}, {x -> 15.6878}, {x -> 18.8652}, {x -> 21.9767}, {x -> 25.1447}, {x -> 28.2631}, {x -> 31.4256}, {x -> 34.5483}, {x -> 37.7073}, {x -> 40.8329}, {x -> 43.9893}, {x -> 47.1171}, {x -> 50.2716}, {x -> 53.4011}, {x -> 56.5542}, {x -> 59.6849}, {x -> 62.8368}, {x -> 65.9686}, {x -> 69.1196}, {x -> 72.2522}, {x -> 75.4024}, {x -> 78.5358}, {x -> 81.6852}, {x -> 84.8193}, {x -> 87.9682}, {x -> 91.1027}, {x -> 94.2511}, {x -> 97.3861}}

The parameter prec gives the precision for the calculation. prec=Infinity gives exact solutions.
wrap (default: Identity) is a function you can apply to the solutions.

Improve this answer
$\endgroup$
1
  • $\begingroup$ With this example don't work: findRoots[Im[Zeta[x]^(1/Log[x])] == 0, {x, 0, 1}]? $\endgroup$ Commented May 16, 2018 at 19:29
17
$\begingroup$

Regarding your first question:

For a certain set of equations, Reduce is able to find all roots and prove that no more roots exist in a given interval. I am not sure how this works, but there's an interesting blog post about it here.

When it is not able to guarantee (using symbolic methods) that there are no more solutions, it may still return a set, but it will give a warning. This is the case for your equation.

Plot[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]], {x, 0, 10}]

Mathematica graphics

Reduce[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] == 0 && 0 < x < 10, x] (* Reduce::incs: Warning: Reduce was unable to prove that the solution set found is complete. >> *) (* ==> x == Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 0.886604635313462076794393681674}] || x == Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 3.03296608901366835385376172847}] || x == Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 6.32396651371147786252003752922}] || x == Root[{BesselJ[1, #1]^2 + BesselK[1, #1]^2 - Sin[Sin[#1]] &, 9.39114340758508579766919382120}] *) ToRules@N[%] (* ==> Sequence[{x -> 0.886605}, {x -> 3.03297}, {x -> 6.32397}, {x -> 9.39114}] *)

Note that it is important to restrict the search to an interval. Otherwise Reduce will tell you that it can't find any solutions.

Improve this answer
$\endgroup$
0
16
$\begingroup$

Very recently, I learned a useful procedure due to J. P. Boyd (see also this and this) that involves expanding a function as a Chebyshev polynomial series, forming the so-called "colleague matrix", and then finding the eigenvalues of this matrix, which are often good approximations to the roots of the original function. I shall present how to do a barebones version of this strategy in Mathematica.

The first step is to form the Chebyshev series coefficients of the function concerned. One could either start from the Taylor series and use any number of methods to convert these coefficients to Chebyshev coefficients. Another route, which I shall present here, is to evaluate the function at the so-called "Chebyshev nodes", suitably transformed to the domain of the original function, and then using the (type-1) discrete cosine transform to produce the Chebyshev coefficients. It goes like this:

f[x_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]]; {xmin, xmax} = {1/2, 100}; n = 128; (* arbitrarily chosen integer *) prec = 20; (* precision *) cnodes = Rescale[N[Cos[Pi Range[0, n]/n], prec], {-1, 1}, {xmin, xmax}]; cc = Sqrt[2/n] FourierDCT[f[cnodes], 1]; cc[[{1, -1}]] /= 2;

n here should be chosen to be greater than or equal to the number of roots known to be in the given interval. That is usually determined through a preliminary Plot[]. (A more sophisticated version, as used in the Chebfun package for MATLAB, does an adaptive increase in the order of the underlying Chebyshev approximation using a technique similar to Clenshaw-Curtis quadrature. I have chosen not to implement this for now.)

With the Chebyshev coefficients now available, one can then form the colleague matrix like so:

colleague = SparseArray[{{i_, j_} /; i + 1 == j :> 1/2, {i_, j_} /; i == j + 1 :> 1/(2 - Boole[j == 1])}, {n, n}] - SparseArray[{{i_, n} :> cc[[i]]/(2 cc[[n + 1]])}, {n, n}];

We can then find the eigenvalues of this matrix, filter out irrelevant eigenvalues, and then use Rescale[] to transform the remaining eigenvalues back to the domain of the original function:

rts = Sort[Select[DeleteCases[ Rescale[Eigenvalues[colleague], {-1, 1}, {xmin, xmax}], _Complex | _DirectedInfinity], xmin <= # <= xmax &]] {0.8866466, 3.0327390, 6.3243343, 9.3928065, 12.5892495, 15.6908929, 18.8634014, 21.9762069, 25.1562125, 28.3003285, 31.4743983, 34.5197235, 37.6746911, 40.8716527, 43.9572481, 47.1394116, 50.2776398, 53.3723848, 56.5916636, 59.6440741, 62.8638369, 66.0001767, 69.0778338, 72.2121395, 75.3874360, 78.5335813, 81.6858635, 84.8171664, 87.9688503, 91.1015749, 94.2504610, 97.3866141}

The eigenvalues obtained look a bit rough here; this is due to the ill-behavior of the original function near the origin. For more well-behaved functions, the eigenvalues are often quite accurate roots of the original function. Of course, one can always use FindRoot[] to polish off these approximations.

Improve this answer
$\endgroup$
12
$\begingroup$

For your first question one can use Ted Ersek's package here :

RootSearch[BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]] == 0, {x, 0, 10}] (* {{x -> 0.886605}, {x -> 3.03297}, {x -> 6.32397}, {x -> 9.39114}} *)
Improve this answer
$\endgroup$
0
9
$\begingroup$

For a finite range of interest, NSolve works well

f[x_] = BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin[Sin[x]]; Manipulate[ Module[{sol}, Column[{ sol = NSolve[{f[x] == 0, 0 <= x <= xmax}, x], Plot[f[x], {x, 0, xmax}, Epilog -> {Red, AbsolutePointSize[6], Point[{x, f[x]} /. sol]}, ImageSize -> 360]}]], {{xmax, 16}, 1, 31, 3, Appearance -> "Labeled"}]

enter image description here

EDIT: With later Mathematica versions (e.g., 10.4.1), some of the roots can be missed unless you use higher precision. The //N after NSolve is used to reduce the displayed precision in the results.

Manipulate[ Module[{sol}, Column[ {sol = NSolve[{f[x] == 0, 0 <= x <= xmax}, x, WorkingPrecision -> 30] // N, Plot[f[x], {x, 0, xmax}, Epilog -> {Red, AbsolutePointSize[6], Point[{x, f[x]} /. sol]}, ImageSize -> 360]}]], {{xmax, 16}, 1, 31, 3, Appearance -> "Labeled"}]
Improve this answer
$\endgroup$
3
$\begingroup$

using my myRootSearch and myRootSearchShow

ClearAll[f]; f[x_, a_] := BesselJ[1, x]^2 + BesselK[1, x]^2 - Sin@Sin[a*x] Manipulate[ Column[{Quiet@myRootSearchShow[f[#, a] &, 0, 50, ImageSize -> 250], Quiet@myRootSearch[f[#, a] &, 0, 50]}], {a, 0, 4}]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.