# Using Python to Find Angles Where Graphs of Bessel Functions Cross

### We explore the mathematical and scientific power of Python and its libraries (like SciPy) by using them to calculate and visualize Bessel functions.

Join the DZone community and get the full member experience.

Join For FreeA previous post of mine looked at the angles that graphs make when they cross. For example, sin() and cos(*x*) always cross with the same angle. The same holds for sin() and cos(*kx*) since the *k* simply rescales the *x*-axis.

The post ended with wondering about functions analogous to sine and cosine, such as Bessel functions. This post will look at that question in more detail. Specifically we'll look at the functions *J*_{ν} and *Y*_{ν}.

Because these two Bessel functions satisfy the same second order linear homogeneous differential equation, the Strum separation theorem says that their zeros are interlaced: between each pair of consecutive zeros of *J*_{ν} is exactly one zero of *Y*_{ν}, and between each pair of consecutive zeros of *Y*_{ν} there is exactly one zero of

In the following Python code, we find zeros of *J*_{ν}, then look in between for places where *J*_{ν} and *Y*_{ν} cross. Next we find the angle the two curves make at each intersection and plot the angles.

```
from scipy.special import jn_zeros, jv, yv
from scipy.optimize import bisect
from numpy import empty, linspace, arccos
import matplotlib.pyplot as plt
n = 3 # bessel function order
N = 100 # number of zeros
z = jn_zeros(n, N) # Zeros of J_n
crossings = empty(N-1)
f = lambda x: jv(n,x) - yv(n,x)
for i in range(N-1):
crossings[i] = bisect(f, z[i], z[i+1])
def angle(n, x):
# Derivatives of J_nu and Y_nu
dj = 0.5*(jv(n-1,x) - jv(n+1,x))
dy = 0.5*(yv(n-1,x) - yv(n+1,x))
top = 1 + dj*dy
bottom = ((1 + dj**2)*(1 + dy**2))**0.5
return arccos(top/bottom)
y = angle(n, crossings)
plt.plot(y)
plt.xlabel("Crossing number")
plt.ylabel("Angle in radians")
plt.show()
```

This shows that the angles steadily decrease, apparently quadratically.

This quadratic behavior is what we should expect from the asymptotics of *J*_{ν} and *Y*_{ν}: For large arguments they act like shifted and rescaled versions of sin(*x*)/√*x*. So if we looked at √*xJ*_{ν} and √*xY*_{ν} rather than *J*_{ν} and *Y*_{ν} we'd expect the angles to reach some positive asymptote, and they do, as shown below.

Published at DZone with permission of John Cook, DZone MVB. See the original article here.

Opinions expressed by DZone contributors are their own.

Comments