Can you suggest a software to solve this equation?

Question

Can you please suggest a free software or website that would allow me to approximate, numerically, the first $n$ roots of the equation $ 2J_{0}(2\alpha)+2\alpha J'_{0}(2\alpha)=0$ ? I'm trying to find a Bessel-Fourier expansion for a function in the interval $[0,2]$ and that's the boundary condition given. As you know, I need the roots of that equation to find the eigenvalues. I've checked Maxima and Octave but I think they don't offer this-or at least it wasn't clear to me they had this feature. I need this ASAP, Thanks in advance!

Note: By $J_{0}$ we mean the Bessel function of order $0$

Answer 1

Realizing my approach in the comments, using numerical integration of $xy''+y'+xy=0$, $y(0)=1$, $y'(0)=0$ with a counting/registering event function for the roots of $2y(x)+xy'(x)$ can be done via the code

a = 1e-5
ya = 1-a*a/4
dya = -a/2

n=15

def rootfun(x,y): return 2*y[0]+x*y[1]
def ode(x,y): return y[1], -y[1]/x-y[0]
res = solve_ivp(ode, [a,n*3.2], [ya,dya], events=rootfun, atol=1e-10, rtol=1e-12, method="DOP853")
print(f"{res.message}, {len(res.t)} nodes, {len(res.t_events[0])} roots")
print(res.t_events[0]/2)

which prints

The solver successfully reached the end of the integration interval., 125 nodes, 16 roots
[ 0.7997246   2.14547923  3.64419446  5.18291555  6.73594101  8.29551651
  9.8586033  11.42378277 12.99030737 14.55774679 16.1258362  17.6944038
 19.26333335 20.83254342 22.40197526 23.97158545]

or as suggested by Claude Leibovici in terms of multiples of $\frac\pi2$ and an offset

print(', '.join(f"{r:.10f}+{n:.0f}*pi/2" for n,r in zip((res.t_events[0]/2)//(np.pi/2), (res.t_events[0]/2)%(np.pi/2))))

resulting in

0.7997246032+0*pi/2, 0.5746829035+1*pi/2, 0.5026018019+2*pi/2, 
0.4705265694+3*pi/2, 0.4527557017+4*pi/2, 0.4415348768+5*pi/2, 
0.4338253348+6*pi/2, 0.4282084863+7*pi/2, 0.4239367597+8*pi/2, 
0.4205798533+9*pi/2, 0.4178729303+10*pi/2, 0.4156442033+11*pi/2, 
0.4137774296+12*pi/2, 0.4121911701+13*pi/2, 0.4108266857+14*pi/2, 
0.4096405508+15*pi/2

Answer 2

Javier Segura has written a Fortran implementation of his 2012 work with Amparo Gil concerning the roots of $x \mathcal{C}'_{\nu}(x)+\gamma \mathcal{C}_{\nu}(x)$, $\gamma\in\mathbb{R}$, $\mathcal{C}$ a linear combination of Bessel functions.

Gil, Amparo, and Javier Segura. 2012. “Computing the Real Zeros of Cylinder Functions and the Roots of the Equation $x\mathcal{C}_{\nu}′(x)+\gamma \mathcal{C}_{\nu}(x)=0$.” Computers & Mathematics with Applications 64 (1): 11–21.

Answer 3

The $n^{\text{th}}$ zero of function $$f(x)=J_0(2 x)-2 x J_1(2 x)$$ is quite close to $$x^{(n)}_0= -1+ \frac \pi 2 n$$

Performing one iteration of Halley method $$x^{(n)}_1=\frac t2\,\,\frac{2 t \left(t^2+2\right) J_0(t){}^2-\left(t^2-1\right) J_1(t) \,J_0(t)+t \left(t^2+1\right) J_1(t){}^2 } {2 t\left(t^2+1\right) J_0(t){}^2+\left(t^2-1\right) J_1(t)\, J_0(t)+t \left(t^2+3\right) J_1(t){}^2 }$$ where $t=(\pi\,n-2)$.

For example $x^{(123)}_1=192.0347$ while the solution is $192.0307$.

For sure, better could be done using instead Householder method (the formula is too long to be typed here). For the above case, the result would be $192.0304$.

Tabulating once the values of $J_0(t)$ and $J_1(t)$ could make life quite easy except if you are looking for very accurate results.