Numerical (in)stability of recurrence relations
The previous post gave several examples of three-term recurrence relations for special functions. These relations can be computationally useful, but they have to be applied carefully.
Several years ago I wrote a post on stable and unstable recurrences. In that post I show that the stability of the recurrence relation for Bessel functions produces depends on which kind of Bessel function and which direction the recurrence is applied.
In the forward direction, computing higher order values from lower order values, works well for Bessel functions of the second kind _Y_ _n_ but not for Bessel functions of the first kind _J_ _n_. In the reverse direction, the recurrence is stable for _J_ _n_ but not for _Y_ _n_.
I didn’t explain in that post why this is. In this post I will.
Second order linear difference equations have two independent solutions, just like second order linear differential equations. For both kinds of equations, all solutions are linear combinations of the two solutions. Suppose one solution grows with _n_ and the other decays. You may want to compute the decaying solution, but in doing so you might pick up a small component of the growing solution due to rounding error. This post illustrates this phenomena for differential equations, and this post illustrates it for difference equations.
When you look at a plot of Bessel functions in a text book, you’ll probably see a few plots of _J_ _n_(_x_) and _Y_ _n_(_x_) for a few small values of _n_. The functions seem to behave roughly the same way, like sine and cosine. And that’s true,**as functions of _x_**.
But it’s not true for _J_ _n_(_x_) and _Y_ _n_(_x_) as functions of _n_ for fixed _x_. As _n_ increases, _J_ _n_(_x_) decays to zero and _Y_ _n_(_x_) goes off to −∞.
That’s the source of numerical instability. And there will be similar instability problems for other recurrences where the ratios of the two independent solutions goes to zero or infinity as a function of _n_.
There are techniques for computing the solution that does not diverse, the so-called minimal solution, such as Miller’s algorithm mentioned here.