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

enter image description here

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

enter image description here

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

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

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

Increasing the WorkingPrecision improves the performance of NSolve

 erg = NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals,
   WorkingPrecision -> wp];
   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}]

enter image description here

  • $\begingroup$ Thanks for your effort. I'm still wondering that such a smooth function needs such an amount of WorkingPrecision to calculate the roots... $\endgroup$ – Ulrich Neumann Jun 29 '18 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}]

Mathematica graphics

rr=FindRoot[BesselJ[0, x], {x, # - 1, # + 1}] & /@ tt

Mathematica graphics

Show[Plot[BesselJ[0, x], {x, 0, 100}], 
 Graphics[ Point /@ ({x, 0} /. rr)]]

enter image description here


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

enter image description here

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

 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"

enter image description here

  NSolve[{BesselJ[0, x], 0 < x < 100}, x, Reals],
  HoldPattern[Select[_, Function[_Equal]]],
  TraceInternal -> True] // Flatten

Mathematica graphics

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.

  • 1
    $\begingroup$ Thank you for your detailed explanation, which shows the imprecision of numerics. I'm still wondering how to check additionaly numerical calculations... $\endgroup$ – Ulrich Neumann Dec 28 '18 at 7:30
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Dec 28 '18 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.