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