Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2012

[Date Index] [Thread Index] [Author Index]

Search the Archive

BesselJ[] with machine-precision input gives sub-machine-precision output

  • To: mathgroup at smc.vnet.net
  • Subject: [mg128412] BesselJ[] with machine-precision input gives sub-machine-precision output
  • From: Yen Lee Loh <yenleeloh1 at gmail.com>
  • Date: Sun, 14 Oct 2012 00:14:22 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • Delivered-to: l-mathgroup@wolfram.com
  • Delivered-to: mathgroup-newout@smc.vnet.net
  • Delivered-to: mathgroup-newsend@smc.vnet.net

Dear MathGroup readers:

Calling BesselJ[0,x] with machine-precision input (with absolute error on the order of 1E-15) gives an output whose absolute error may be as large as 1E-13.  This effect is most severe for x ~ 11, which is apparently where Wolfram's algorithm switches over from the Taylor series to an asymptotic expansion.  See diagnostics in Listing A (including a plot).

The Mathematica documentation for BesselJ does not warn the user about this.  The wording in guide/BesselRelatedFunctions ("Using original algorithms developed at Wolfram Research, Mathematica has full coverage of all standard Bessel-related functions - evaluating every function to arbitrary precision with optimized algorithms for arbitrary complex values of its parameters") may also give the user a false sense of security.  It is understandable that complicated functions like NSum[], NIntegrate[], FindMinimum[] may give sub-machine-precision.  However, special functions like BesselJ[] REALLY ought to be implemented with precision of a few times $MachineEpsilon.

I have constructed an interpolating function that reproduces the correct value of J_n(x) with errors as small as 4E-16.  See demo in Listing B.

The fact that a user-supplied interpolant works better than Wolfram's built-in implementation (at machine precision) suggests that Wolfram's implementation ought to be improved.  I presume that Numerical Recipes gives the Bessel J function accurate to at least 15 decimals over 0<x<20; I have not checked.

See also my post entitled "BesselK, BesselI, BesselJ, BesselY are unnecessarily slow."

-YLL

-------- Listing A --------
x = 11 + 2/9;
10^(-Accuracy[N@x])  (* absolute error in floating-point rep of x *)
BesselJ[0, x] // N[#, 30] & // Print (* use extra-precision arithmetic *)
BesselJ[0, x] // N[#, 20] & // Print (* use extra-precision arithmetic *)
BesselJ[0, x] //  N[#, 17] & // Print (* use extra-precision arithmetic *)
BesselJ[0, x] // N // InputForm // Print  (* use machine-precision arithmetic *)
BesselJ[0, x // N] // InputForm // Print  (* use machine-precision arithmetic *)

xlist = Range[0, 20, 1/100];
errlist = Table[
   N[BesselJ[0, x]] - N[BesselJ[0, x], 24]
   , {x, xlist}
   ];
ListLogPlot[Transpose@{xlist, Abs@errlist}, PlotRange -> All,
 FrameLabel -> {"x", "absolute\nerror in\nBesselJ[0,x]"}]


-------- Listing B --------
myHankelH[x_] := Exp[I x]*
   InterpolatingFunction[{{0.`, 1.`}}, {4, 11, 0, {43}, {22}, 0, 0, 0,
       0, Automatic}, {{0.`, 0.02380952380952378`,
       0.047619047619047616`, 0.0714285714285714`,
       0.09523809523809523`, 0.11904761904761901`,
       0.14285714285714285`, 0.16666666666666663`,
       0.19047619047619047`, 0.21428571428571425`,
       0.23809523809523808`, 0.26190476190476186`,
       0.2857142857142857`, 0.3095238095238095`, 0.3333333333333333`,
       0.3571428571428571`, 0.38095238095238093`, 0.4047619047619047`,
        0.42857142857142855`, 0.45238095238095233`,
       0.47619047619047616`, 0.49999999999999994`,
       0.5238095238095237`, 0.5476190476190474`, 0.5714285714285714`,
       0.5952380952380951`, 0.6190476190476191`, 0.6428571428571428`,
       0.6666666666666666`, 0.6904761904761905`, 0.7142857142857142`,
       0.7380952380952379`, 0.7619047619047619`, 0.7857142857142856`,
       0.8095238095238095`, 0.8333333333333333`, 0.8571428571428571`,
       0.8809523809523809`, 0.9047619047619047`, 0.9285714285714284`,
       0.9523809523809523`, 0.976190476190476`,
       1.`}}, {{0.` + 0.` I}, {0.013432133128241612` -
        0.013434036909816599` I}, {0.026858545827246595` -
        0.026873776036828972` I}, {0.0402734815472207` -
        0.04032488287548035` I}, {0.05367111242358459` -
        0.053792948595401166` I}, {0.0670455059881269` -
        0.06728345074056401` I}, {0.0803905948620758` -
        0.08080171153596163` I}, {0.0937001506328062` -
        0.09435285331019386` I}, {0.10696776324101506` -
        0.10794175083120973` I}, {0.12018682729128204` -
        0.12157298073226136` I}, {0.13335053669856584` -
        0.13525076866561364` I}, {0.1464518889482846` -
        0.14897893533052886` I}, {0.1594836999425468` -
        0.16276084302258054` I}, {0.17243862992160156` -
        0.17659934476757114` I}, {0.18530922031759284` -
        0.1904967383559213` I}, {0.19808794068534233` -
        0.20445472762078912` I}, {0.21076724415559123` -
        0.2184743930804028` I}, {0.22333962926673886` -
        0.23255617361242553` I}, {0.23579770562906985` -
        0.24669986020623685` I}, {0.24813426070305328` -
        0.2609046021330251` I}, {0.26034232503306337` -
        0.2751689251732932` I}, {0.2724152335382705` -
        0.2894907609234868` I}, {0.2843466808699183` -
        0.3038674857206848` I}, {0.2961307693366367` -
        0.31829596740209437` I}, {0.3077620484178032` -
        0.3327726179554218` I}, {0.31923554538072274` -
        0.34729345009977786` I}, {0.33054678695612444` -
        0.3618541359365182` I}, {0.3416918123879788` -
        0.37645006599293473` I}, {0.35266717845007994` -
        0.39107640721745457` I}, {0.3634699572147776` -
        0.405728158745364` I}, {0.37409772747658643` -
        0.4204002045170751` I}, {0.38454856078657906` -
        0.43508736208049015` I}, {0.39482100305533335` -
        0.44978442713454164` I}, {0.4049140526454271` -
        0.4644862135665527` I}, {0.41482713581057723` -
        0.47918758889938623` I}, {0.4245600802573951` -
        0.4938835051957403` I}, {0.434113087515415` -
        0.5085690255685296` I}, {0.44348670470776197` -
        0.5232393465211681` I}, {0.4526817962231814` -
        0.5378898163934211` I}, {0.4616995157033277` -
        0.5525159502210358` I}, {0.47054127867928064` -
        0.5671134413342503` I}, {0.479208736119363` -
        0.5816781700247895` I}, {0.48770374908695635` -
        0.5962062096060041` I}}, {Automatic}][x^(-1/2)];
myBesselJ0[x_] := Re@myHankelH@x;

xlist = Range[2, 20, 1/100];
errlist = Table[
   myBesselJ0[x] - N[BesselJ[0, x], 24]
   , {x, xlist}
   ];
ListLogPlot[Transpose@{xlist, Abs@errlist}, PlotRange -> All,
 FrameLabel -> {"x", "absolute\nerror in\nmyBesselJ[0,x]"}]



  • Prev by Date: Re: plot of a function of two variables
  • Next by Date: BesselK, BesselI, BesselJ, BesselY are unnecessarily slow
  • Previous by thread: Re: plot function of two variables
  • Next by thread: BesselK, BesselI, BesselJ, BesselY are unnecessarily slow