# Roots of BesselJ

Bug introduced in 11.0 or earlier and fixed in 14.0.0 or earlier

BesselJ[0,x] is a smooth quasiperiodic function. NSolve is able to find all roots in the given range 0<x<50:

erg = NSolve[{BesselJ[0, x], 0 < x < 50}, x, Reals]
Show[{Plot[BesselJ[0, x], {x, 0, 50}] , Graphics[Point[{x, 0} /. erg]]}]


But if I increase the range 0<x<100 NSolve misses many of the expected 32 roots.

erg = NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals]
Show[{Plot[BesselJ[0, x], {x, 0, 100}] ,Graphics[Point[{x, 0} /. erg]]}]


Any idea why NSolve fails? Is there some kind of critical argumentsize x in the definition of BesselJ[0,x]? Thanks!

• FYI it doesn't fail on mathematica 11.3 macOS Commented Jun 29, 2018 at 10:50
• @ chris:Thanks. My version 11.0.1.0 (Windows 64) only detects 18 (of the expected 31?) roots in NSolve[{BesselJ[0, x], 0 < x < 100 }, x, Reals] Commented Jun 29, 2018 at 10:59
• It fails on Version 11.3, Windows 64, only 18 roots.
– rmw
Commented Jun 29, 2018 at 13:52
• Version 8.0 on Windows 32 has no problem to find all 32 roots. Commented Jun 29, 2018 at 15:50
• Linux Version does not fail, even for 200 where it finds 63 roots ($Version 10.0 for Linux x86 (64-bit) (June 29, 2014)) Commented Dec 27, 2018 at 18:17 ## 3 Answers $Version

(* "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)" *)


Increasing the WorkingPrecision improves the performance of NSolve

Manipulate[
erg = NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals,
WorkingPrecision -> wp];
Show[{
Plot[BesselJ[0, x], {x, 0, 100}],
Graphics[{Red, AbsolutePointSize[4], Point[{x, 0} /. erg]}]}],
{{wp, Automatic, "WorkingPrecision"},
{Automatic, Range[20, 30, 2]} // Flatten,
ControlType -> SetterBar}]


EDIT: NSolve in the current version of Mathematica performs as expected for much larger values of xmax

$Version (* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *) Clear["Global*"]; xmax = 2000; (erg = Solve[{BesselJ[0, x] == 0, 0 < x < xmax}, x, Reals])[[1 ;; 4]] (* {{x -> BesselJZero[0, 1]}, {x -> BesselJZero[0, 2]}, {x -> BesselJZero[0, 3]}, {x -> BesselJZero[0, 4]}} *) Length@erg (* 636 *) ergN = NSolve[{BesselJ[0, x] == 0, 0 < x < xmax}, x, Reals]; (x /. erg) == (x /. ergN) (* True *)  • Thanks for your effort. I'm still wondering that such a smooth function needs such an amount of WorkingPrecision to calculate the roots... Commented Jun 29, 2018 at 14:53 I am not sure if this would qualify as a workaround for what you want to do but this works. tt = Table[BesselJZero[0, i] // N, {i, 25}]  rr=FindRoot[BesselJ[0, x], {x, # - 1, # + 1}] & /@ tt  Show[Plot[BesselJ[0, x], {x, 0, 100}], Graphics[ Point /@ ({x, 0} /. rr)]]  I would like to argue that this is a bug. The roots are roughly located and polished with FindRoot at a WorkingPrecision of 10 plus the WorkingPrecision of NSolve. In this case that is ten more than machine precision, or MachinePrecision + 10, which is about 26.95 digits. We can see that all 32 roots are found. fn = BesselJ[0, #] &; (* for convenience *) rts0 = Trace[ NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals], HoldPattern[f : FindRoot[_, {v_, _}, ___]] :> (v /. f), TraceInternal -> True] // Flatten // ReleaseHold // Sort; Plot[fn[x], {x, 0, 100}, PlotLabel -> Row[{Length@rts, "roots"}, " "], Epilog -> {Red, Point@Thread[{rts0, 0}]}]  NSolve then set the roots to the arbitrary-precision equivalent of the working precision, which in this case is $MachinePrecision and deletes the value of x for which fn[x] is nonzero:

rts = SetPrecision[rts0, $MachinePrecision]; ListPlot[ fn[rts] /. r_ /; r == 0 :> 0 // RealExponent, GridLines -> {Flatten@Position[fn[rts], r_ /; r != 0], None}, PlotLabel -> "Indices of final roots, rejected with nonzero residual" ]  Trace[ NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals], HoldPattern[Select[_, Function[_Equal]]], TraceInternal -> True] // Flatten  The real problem is that at an arbitrary-precision of $MachinePrecision there is no number near each discarded root that can make the Bessel function evaluate to an arbitrary-precision zero. (It's the nature of a finite machine that not every real number can be represented.) Usually what happens at the closest arbitrary-precision number to a root is that the round-off error bound calculated by the arbitrary-precision code is greater than the computed value of fn[x] at a zero x, and an arbitrary-precision zero is returned. I suspect the assumption was that this always happens, but apparently it does not.

• Thank you for your detailed explanation, which shows the imprecision of numerics. I'm still wondering how to check additionaly numerical calculations... Commented Dec 28, 2018 at 7:30
• @UlrichNeumann You're welcome. I don't know your last comment is supposed to be a question for me or not. I can add that BesselJ[] seems to exhibit interesting precision-loss issues as the working precision changes. It suggests to me that it switches algorithms depending on the argument and precision. Nothing very dramatic. I don't understand it at all and could not figure out how to incorporate into my answer. It may have nothing to do with the problem at hand, or it possibly might explain why NSolve` fails here. Commented Dec 28, 2018 at 23:50