Roots of BesselJ

Question

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!

Answer 1

$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 *)

Answer 2

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)]]

Answer 3

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.