1
$\begingroup$

It's almost the same question with "multiple root finding", but the order of roots is important in my case.

Basically, I'm trying to find $n$-th crossing point of the following graph:

Plot[{x BesselJ[1, x], 2 BesselJ[0, x]}, {x, 0, 15}, PlotRange -> All]

For example, I need to find out the 1005th (positive and the order count from the smallest) root of the equation below:

x BesselJ[1, x]-2 BesselJ[0, x]==0

Until now, FindRoot seems to be the fastest and clearest method, but I couldn't find a way to specify the $n$-th root. 

FindRoot[x BesselJ[1, x] == A BesselJ[0, x], {x, 8, 12}, Method -> "Brent"]
Improve this question
$\endgroup$

2 Answers 2

Reset to default
4
$\begingroup$

x BesselJ[1, x] very soon dominates the behavior and you can use BesselJZero[n, k] to constrain FindRoot[]:

f[x_] = x BesselJ[1, x] - 2 BesselJ[0, x] zero[n_] := BesselJZero[1, n - 1] // N FindRoot[f[x] == 0, {x, zero[1005] - 1, zero[1005] + 1}]

Of course you need to be a little careful about the interval you search for the root in. A more thorough approach where we aren't guaranteed to have nicely spaced roots would be to step forward one root at a time.

Another approach if you know that roots are at least dx apart would be to observe the sign of your function as you step forward on the x-axis with dx/2 steps and simply count the number of sign changes, then again limit the search with two values - FindRoot[f[x], {x, x1, x2}] - to find the root in the identified interval.

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

Another possibility when trying to find the $k$-th root of a function would be to adapt the NDSolve[]-based event detection strategy by Daniel given in this answer; that is, keep a running count as each zero is crossed, and then bail when you've hit the one you wanted. For instance:

f[x_] := x BesselJ[1, x] - 2 BesselJ[0, x] k = 0; j = 1005; (* index of zero *) sol = Reap[NDSolve[{y'[x] == f'[x], y[0] == f[0], WhenEvent[y[x] == 0, k++; If[k == j, Sow[x]; "StopIntegration"], "DetectionMethod" -> "Interpolation", "IntegrateEvent" -> True, "LocationMethod" -> "Brent"]}, {}, {x, 0, ∞}, AccuracyGoal -> ∞, WorkingPrecision -> 20];][[-1, 1, 1]]

which yields

3154.9449374317543982

You can polish further with FindRoot[] if wanted:

sol = \[FormalX] /. First[FindRoot[f[\[FormalX]], {\[FormalX], sol}, WorkingPrecision -> 20]] 3154.9449374319666351

If you wish for an entire range (say, the 26th to the 30th zeros), the event inside WhenEvent[] is easily modified:

k = 0; {kmin, kmax} = {26, 30} Reap[NDSolve[{y'[x] == D[f[x], x], y[0] == f[0], WhenEvent[y[x] == 0, k++; If[kmin <= k <= kmax, Sow[x], If[30 < k, StopIntegration"]], "DetectionMethod" -> "Interpolation", "IntegrateEvent" -> True, "LocationMethod" -> "Brent"]}, {}, {x, 0, ∞}, AccuracyGoal -> ∞, WorkingPrecision -> 20];][[-1, 1]] {79.345692016011214192, 82.486505133453918880, 85.627375391257917317, 88.768296740115704530, 91.909263963793310359}

and again one can post-process with FindRoot[] if need be.

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.