Rational minimax approximation for J₀(x)

How I found a small error in the Digital Library of Mathematical Functions ¹.

[DLMF, §3.11(iii)] describes a rational minimax approximation to the Bessel function J₀(x) (emphasis mine)

With w(x)=1 and 14-digit computation, we obtain the following rational approximation of type [3,3] to the Bessel function J₀(x) (§10.2(ii)) on the interval 0≤x≤j₀,₁, where j₀,₁ is the first positive zero of J₀(x):

The legend of [DLMF, Figure 3.11.1] states j₀,₁ = 0.89357…, but I seem to recall that the first positive zero of J₀ is closer to 2.405. [AS, Table 9.5, p. 409]² gives the value j₀,₁ = 2.40485 5577 for the zero and [DMLF, Figure 10.3.1] shows the zero close to x=2.5. Since both [AS] and [DLMF, §10.21(i)] use j₀,₁ to denote the first positive zero of J₀, it appears to have a different meaning here. The symbol b will be used below for the value 0.89357…. I shall look at the approximation to resolve this contradiction.

Plotting R₃,₃(x), J₀(x) and the approximation error e(x) = R₃,₃(x) - J₀(x) (using the logarithmic y-axis on the right) shows that the approximation works well for 0 ≤ x ≤ b, but that the error increases rapidly outside this interval.

# rational minimax approximation to J0(x), from
# https://dlmf.nist.gov/3.11#iii
# coefficients
# https://dlmf.nist.gov/3.11.T1
p0 =  0.99999998917854
p1 = -0.34038938209347	
q1 = -0.34039052338838
p2 = -0.18915483763222
q2 =  0.06086501629812
p3 =  0.06658319420166	
q3 = -0.01864476809090
# https://dlmf.nist.gov/3.11.E19
R33(x) = (p0 + x*(p1 + x*(p2 + x*p3))) / (1 + x*(q1 + x*(q2 + x*q3)))
b = 0.89357
set xrange [-0.2:2.5]
set yrange [-0.2:1.1]
set logscale y2
set y2range [1e-9:1]
set y2tics
set samples 500
set key outside top center

e(x) = R33(x) - besj0(x)
print 0, e(0)
print b, e(b)

set term pngcairo transparent
set output 'plot/j0rat-01.png'
plot R33(x) t 'R₃,₃(x)' w lines, \
     besj0(x) t 'J₀(x)' w lines, \
     abs(R33(x) - besj0(x)) axes x1y2 \
     t '|R₃,₃(x) - J₀(ax)|', \
     '-' using 1:2 axes x1y2 w points t 'Boundary points'
     0        1.08214599503498e-08
     0.89357 +1.08214395222461e-08
 e

set output
set term svg
set output 'plot/j0rat-01.svg'
replot
set output

set xrange [-0.05:0.95]
set yrange [-2.0e-8:2.0e-8]
unset y2tics
#set samples 10000
set output 'plot/j0rat-02.svg'

plot R33(x) - besj0(x) axes x1y1 \
     t 'R₃,₃(x) - J₀(ax)', \
     '-' using 1:2 axes x1y1 w points t 'Extrema and boundary points'
     0       -1.08214599503498e-08
     0.046   +1.08216265948258e-08
     0.175   -1.08208797477971e-08
     0.355   +1.08209655680369e-08
     0.555   -1.08227087292079e-08
     0.73    +1.08232458551072e-08
     0.851   -1.08223704442523e-08
     0.89357 +1.08214395222461e-08
 e
set output

the gnuplot script produces the plots shown above and below and prints the approximation error e(x) = R₃,₃(x) - J₀(x) for x = 0 and x = b.

gnuplot> load 'j0rat.gp'
0 -1.08214599503498e-08
0.89357 1.08214395222461e-08
gnuplot>   

The plot below reproduces [DLMF, Figure 3.11.1], the rightmost point has the x-coordinate b=0.89357…, so this value is the upper limit if the interval where the absolute value of the approximation error does not exceed ϵ.

The value b = 0.89357… is the upper limit of the interval where the rational approximation maintains its error below ϵ ≈ 1.08e-8 and is unrelated to the first zero j₀,₁ of J₀.