Zeros of Bessel functions
Dominic Berry - 21 October 2019
The following code is for the asymptotic expansions of the zeros of Bessel functions, and the zeros of the derivatives of Bessel functions. We determine the numbers α and α’ as in 10.21.25 and 10.21.30 of the NIST DLMF, except to an arbitrary number of terms.
https://dlmf.nist.gov/10.21
The method is from the paper “A further method for the evaluation of zeros of Bessel functions and some new asymptotic expansions for zeros of functions of large order” by FWJ Olver in 1951.
https://doi.org/10.1017/S0305004100027158
The first block of code is for the zeros of the Bessel function. The value of mmax is the range of values of m that we find ϕ_m for, and smax is the maximum value of s and therefore the number of α that we attempt to solve for. On the LHS of (8.7) we have orders up to 1/u^{2*smax}, and so we need to expand the RHS to at least that order.  Given m, the minimum order on the RHS is 1/u^{6m}, which is with s=0 on the RHS.  That means that if we have 6*mmax=2*smax, then the range of m is sufficient to give orders up to 1/u^{2*smax} correctly.  For the sum over s on the RHS, the worst case is m=0, in which case the sum over s needs to go up to smax to correctly give orders up to 1/u^{2*smax}.

Clear[ϵ];
$Assumptions={ϵ≥0,ϵ<1,u∈Reals,u>0,x∈Reals};
mmax=8; (*This sets the number of terms solved for. We get α up to α[3*mmax].*)
smax=3*mmax;
ϕ=ConstantArray[0,mmax+1];
For[m=0,m≤ mmax,m++,
(*This corresponds to Eq. (2.5), where we have used the residue theorem instead of taking the integral. Then rs is the same as a_m(λ*Sech[u]) with λ=1/(1+ϵ). Here we use x=1+Sech[u].*)
rs=Residue[(t*(x-1)/(1+ϵ)+Sinh[t])^(-(2*m+1)),{t,0}];
(*We now take the only λ-dependent parts in (2.7) and replace λ with 1/(1+ϵ). We perform an expansion in ϵ to make the integral easier. Multipliying by x^(m*3+smax+1) and expanding gives an expression with no negative powers of x=1+Sech[u], so Expand greatly simplies the expression.*)
s2=Expand[Normal[Series[rs/(2*(1+ϵ))^(2*m),{ϵ,0,smax}]]*x^(m*3+smax+1)];
(*Next we perform the integral in (2.7). Now this ϕ corresponds to that ϕ_m(1+ϵ) in (8.4), and ϕ_m(λ)=f_m(1/λ). That is why we take λ=1/(1+ϵ). The expansion in (8.4) is up to order smax in ϵ.*)
ϕ[[m+1]]=(-1)^(m)*Factorial[2*m]*Expand[Total[Integrate[MonomialList[s2/.x-> 1+Sech[u],Sech[u]]*Sech[u]^(2*m+1)/(1+Sech[u])^(m*3+smax+1),{u,0,∞}]]]]
(*We now replace x with 1+Sech[u] in s2, and use MonomialList which separates out the powers of Sech[u]. Then each term has a power of Sech[u] divided by (1+Sech[u])^(m*3+smax+1), which Mathematica can integrate much more easily.*)
αv=Array[α,smax];
ϵ=Sum[αv[[s]]/u^(2*s),{s,1,smax}]; (*We use this substitution for ϵ because that is what is used in (8.7).*)
lhs=1-Sum[(2*s-3)/3*αv[[s]]/u^(2*s),{s,1,smax}];(*This is the left-hand side of (8.7).*)
rhs=Sum[(-1)^m/u^(6*m)*ϕ[[m+1]],{m,0,mmax}];(*This is the right-hand side of (8.7).*)
s3=Series[rhs-lhs,{u,∞,2*smax}];(*We take the series of the difference of the sides so we can match the left and right sides.*)
eqs=ConstantArray[0,smax-1];
For[m=1,m<smax,m++,eqs[[m]]=Simplify[SeriesCoefficient[s3,2*m+2]]==0](*We set the equations to be those solving (8.7) for α.*)
Simplify[Solve[eqs,Part[αv,2;;smax]]]

Next, this code is for determining the expansion for the zeros of the derivative of the Bessel function.  The key equation is now (9.4).  On the LHS we expand both sums to order smax to get maximum order 1/u^{2*smax}.  However, now the orders up to 1/u^4 cancel, and also the terms with a_smax cancel on both sides.  That means we can take the  terms for 1/u^6 up to 1/u^{2*smax}, and solve for a_2 up to a_{smax-1}.  For the RHS the expansion for m needs to be up to mmax with 6*mmax=2*smax.

Clear[ϵ,λ];
$Assumptions={m∈Integers,ϵ≥0,ϵ<1,u∈Reals,u>0,λ>0,λ∈Reals,x∈Reals};
mmax=8; (*This sets the number of terms solved for. We get α up to α[3*mmax]-1.*)
smax=3*mmax;
ψ=ConstantArray[0,mmax+1];
rs=ψ;
rs1=ψ;
rs2=ψ;
h=ψ;
For[m=0,m≤mmax,m++,
rs[[m+1]]=(-1)^(m)*Factorial[2*m]*(λ/2)^(2*m)*Simplify[Residue[(t*(x-1)*λ+Sinh[t])^(-(2*m+1)),{t,0}]]; (*This is giving f_m from (2.5) and (2.7), where we now include the other terms but do not integrate over u. Here we use x=1+Sech[u].*)
rs1[[m+1]]=Simplify[D[rs[[m+1]],λ]];(*We take the derivatives of f_m with respect to λ that we need.*)
rs2[[m+1]]=Simplify[D[rs1[[m+1]],λ]];
h[[m+1]]=Simplify[If[m==0,0,λ^2/(2*(1-λ^2))*(λ^2*rs2[[m]]+3*λ*rs1[[m]]+rs[[m]])]]] (*This expression is what we need to subtract from f_m in the first line of (5.6), except we still need to integrate over u.*)
λ=1/(1+ϵ);
For[m=0,m≤mmax,m++,
(*Multipliying by x^(m*3+smax+1) and expanding gives an expression with no negative powers of x=1+Sech[u], so Expand greatly simplies the expression.*)
s2=(-1)^m*Expand[Normal[Series[ϵ*(Sech[u]^2*rs[[m+1]]-h[[m+1]]),{ϵ,0,smax}]]*x^(m*3+smax)];
(*This is now the expression for ψ in (9.1), except we still need to integrate over u. We have introduced the factors of Sech that will be needed in the integral, though. A subtlety is that the powers of Sech need to be different for the two parts. The expansion over ϵ corresponds to the sum over s on the RHS of (9.4), which is why order smax is needed.*)
ψ[[m+1]]=Expand[Total[Integrate[MonomialList[s2/.x-> 1+Sech[u],Sech[u]]*Sech[u]^(2*m-1)/(1+Sech[u])^(m*3+smax),{u,0,∞}]]]](*Finally we perform the integral from (2.7). Again we replace x with 1+Sech[u] in s2 and split it up to make the integral easier.*)
av=Array[a,smax-1];
Clear[ϵ];
lhs=(1-Sum[(2*s-3)/3*av[[s]]/u^(2*s),{s,1,smax-1}])*Sum[av[[s]]/u^(2*s),{s,1,smax-1}]; (*The LHS of (9.4).*)
rhs=Sum[1/u^(6*m)*ψ[[m+1]],{m,0,mmax}]; (*The RHS of (9.4) where we will next redefine ϵ to get the expression in round brackets.*)
ϵ=Sum[av[[s]]/u^(2*s),{s,1,smax-1}]; (*This corresponds to the part in round brackets on the RHS of (9.4).*)
s3=Series[rhs-lhs,{u,∞,2*smax}];(*We take the series of the difference of the sides so we can match the left and right sides.*)
eqs=ConstantArray[0,smax-2];
For[m=1,m<smax-1,m++,eqs[[m]]=Simplify[SeriesCoefficient[s3,2*m+4]]==0];(*We set the equations to be those solving (9.4) for α.*)
Expand[Solve[eqs,Part[av,2;;smax-1]]]