Brash & Plucky

Bessel functions

Bessel functions

From Wikipedia

Iterative implementation

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

J0-4 in [0..2pi]

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

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]

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 will discuss in a future post, I am after j1 only, so… my job here is done.