Bessel functions

From Wikipedia

Membrane vibration

Iterative implementation

Bessel Jn source code

I came up with the above implementation, with which I plotted these.

J0..4 in [0..2pi]

Bessel Jn plot

2*abs(J1..4(x)/x) in [0..2pi]

Bessel Jn plot (abs)

Iterative methods are a terrible choice for many use cases because of error accumulation, jitter, the very presence of a loop, etc… But in this case this method also happens to be pretty weak and way too easy to break:

J0..4 in [0..8pi]

Bessel error accumulation

Deus ex machina

After some unsuccessful googling I was about to work on my own Taylor Series approximation or at least study the proposed implementations in Numerical Recipes. But I learnt to my surprise :o) that <math.h> provides j0/j1/jn/y0/y1/yn and even CUDA provides j0f/j1f/jnf/y0f/y1f/ynf.

I checked whether the C/CUDA implementations match my iterative implementation, and they do.

For reasons I may discuss in a future post, I am after j1 only, so… my job here is done.