# Three-Term Recurrence Relations and Bessel Functions

Posted by **Cleve Moler**,

Three-term recurrence relations are the basis for computing Bessel functions.

### Contents

#### A Familiar Three-Term Recurrence

Start with two large integers, say 514229 and 317811. For reasons that will soon become clear, I'll label them $f_{30}$ and $f_{29}$.

clear n = 30 f(30) = 514229; f(29) = 317811

n = 30 f = Columns 1 through 6 0 0 0 0 0 0 Columns 7 through 12 0 0 0 0 0 0 Columns 13 through 18 0 0 0 0 0 0 Columns 19 through 24 0 0 0 0 0 0 Columns 25 through 30 0 0 0 0 317811 514229

Now take repeated differences, and count backwards.

$$ f_{n-1} = f_{n+1} - f_n, \ n = 29, 28, ..., 2, 1 $$

for n = 29:-1:2 f(n-1) = f(n+1) - f(n); end n f

n = 2 f = Columns 1 through 6 0 1 1 2 3 5 Columns 7 through 12 8 13 21 34 55 89 Columns 13 through 18 144 233 377 610 987 1597 Columns 19 through 24 2584 4181 6765 10946 17711 28657 Columns 25 through 30 46368 75025 121393 196418 317811 514229

Recognize those integers? They're the venerable Fibonacci numbers. I've just reversed the usual procedure for computing them. And, I didn't pick the first two values at random. In fact, their choice is very delicate. Try picking any other six-digit integers, and see if you can take 30 differences before you hit a negative value. Of course, we usually start with $f_0 = 0$, $f_1 = 1$ and go forward with additions.

The relation

$$ f_n = f_{n+1} - f_{n-1} $$

is an example of a *three-term recurrence*. I've written it in such a way that I can go either forward or backward.

#### Friedrich Bessel

Friedrich Bessel was a German astronomer, mathematician and physicist who lived from 1784 until 1846. He was a contemporary of Gauss, but the two had some sort of falling out. Although he did not have a university education, he made major contributions to precise measurements of planetary orbits and distances to stars. He has a crater on the Moon and an asteroid named after him.

In mathematics, the special functions that bear his name were actually discovered by someone else, Daniel Bernoulli, and christened after Bessel's death. They play the same role in polar coordinates that the trig functions play in Cartesian coordinates. They are applicable in a wide range of fields, from physics to finance.

#### Bessel Functions

Bessel functions are defined by an ordinary differential equation that generalizes the harmonic oscillator to polar coordinates. The equation contains a parameter $n$, the *order*.

$$ x^2 y'' + x y' + (x^2 - n^2) y = 0 $$

Solutions to this equation are also sometimes called *cylindrical harmonics*.

For my purposes here, their most important property is the three-term recurrence relation:

$$ \frac{2n}{x} J_n(x) = J_{n-1}(x) + J_{n+1}(x) $$

This has a variable coefficient that depends upon the two parameters, the order $n$ and the argument $x$.

To complete the specification, it would be necessary to supply two initial conditions that depend upon $x$. But, as we shall see, using the recurrence in the forward direction is disastrous numerically. Alternatively, if we could somehow supply two end conditions, we could use the recurrence in the reverse direction. This the idea behind a classical method known as Miller's algorithm.

#### Miller's Algorithm

Miller's algorithm also employs this interesting identity, which I'll call the *normalizer*.

$$ J_0(x) + 2 J_2(x) + 2 J_4(x) + 2 J_6(x) + \ ... = 1 $$

It's not obvious where this identity comes from and I won't try to derive it here.

We don't know the end conditions, but we can pick arbitrary values, run the recurrence backwards, and then use the normalizer to rescale.

Step 1. Pick a large $N$ and temporarily let

$$ \tilde{J}_N = 0 $$

$$ \tilde{J}_{N-1} = 1 $$

Then for $n = N-1, ..., 1$ compute

$$ \tilde{J}_{n-1} = 2n/x \tilde{J}_n \ - \tilde{J}_{n+1} $$

Step 2. Normalize:

$$ \sigma = \tilde{J}_0 + 2 \tilde{J}_2 + 2 \tilde{J}_4 + ... $$

$$ J_n = \tilde{J}_n/\sigma, \ n = 0, ..., N $$

I'll show some results in a moment, after I introduce an alternative formulation.

#### Matrix Formulation

You know that I like to describe algorithms in matrix terms. The three-term recurrence generates a tridiagonal matrix. But which three diagonals? I can put the three below the diagonal, above the diagonal, or centered on the diagonal. This is accomplished with the second argument of the sparse matrix generator, `spdiags`. Here is the code.

```
type blt
```

function B = blt(n,x,loc) % BesseL Three term recurence. % B = blt(n,x,loc) is an n-by-n sparse tridiagonal coefficent matrix % that generates the Bessel functions besselj((0:n-1)',x). % loc specifies the location of the diagonals. % loc = 'center', the default, centers the three diagonals. % loc = 'lower' puts the three diagonals in the lower triangle. % loc = 'upper' puts the three diagonals in the upper triangle. if nargin == 2 loc = 'center'; end switch loc case 'center', locd = -1:1; case 'lower', locd = -2:0; case 'upper', locd = 0:2; end % Three term recurrence: 2*n/x * J_n(x) = J_(n-1)(x) + J_n+1(x). d = -2*(0:n-1)'/x; e = ones(n,1); B = spdiags([e d e],locd,n,n); end

#### Lower

If there were no roundoff error, and if we knew the values of the first two Bessel functions, $J_0(x)$ and $J_1(x)$, we could use backslash, the linear equation solver, with the diagonals in the lower triangle to compute $J_n(x)$ for $n > 1$. Here, for example, is $x = 1$.

format long J = @(n,x) besselj(n,x); n = 8 x = 1.0 B = blt(n,x,'lower'); rhs = zeros(n,1); rhs(1:2) = J((0:1)',x); v = B\rhs

n = 8 x = 1 v = 0.765197686557967 0.440050585744934 0.114903484931901 0.019563353982668 0.002476638964110 0.000249757730213 0.000020938338022 0.000001502326053

These values are accurate to almost all the digits shown. Except for the last one. It's beginning to show the effects of severe numerical instability. Let's try a larger value of `n`.

n = 18 B = blt(n,x,'lower'); rhs = zeros(n,1); rhs(1:2) = J((0:1)',x); v = B\rhs; semilogy(0:n-1,v,'o-') title('J_n(1)') xlabel('n')

n = 18

The values beyond about `n = 9` are suspect and, in fact, totally wrong. The forward elimination that backslash uses to solve this lower triangular system is just the three-term recurrence for increasing `n`. And that recurrence has two solutions -- the one we're trying to compute and a second one corresponding to $Y_n(x)$, the *Bessel function of second kind*. $Y_n(x)$ is stimulated by roundoff error and completely takes over for $n > 9$.

So the lower triangular option -- the forward recurrence -- is unsatisfactory for two reasons: it requires two Bessel values to get started and it is unstable.

#### Upper

We need to make one modification to the upper triangular coefficient matrix.

```
n = 30;
B = blt(n,x,'upper');
B(n-1,n) = 0;
```

Now the last two rows and columns are a 2-by-2 identity.

full(B(n-1:n,n-1:n))

ans = 1 0 0 1

The back substitution process just runs the three-term recursion backwards. We need to start with $J_{n-1}(x)$ and $J_{n-2}(x)$.

rhs = zeros(n,1); rhs(n-1:n) = J((n-2:n-1)',x); format long e v = B\rhs

v = 7.651976865579656e-01 4.400505857449330e-01 1.149034849319004e-01 1.956335398266838e-02 2.476638964109952e-03 2.497577302112342e-04 2.093833800238925e-05 1.502325817436807e-06 9.422344172604491e-08 5.249250179911870e-09 2.630615123687451e-10 1.198006746303136e-11 4.999718179448401e-13 1.925616764480172e-14 6.885408200044221e-16 2.297531532210343e-17 7.186396586807488e-19 2.115375568053260e-20 5.880344573595754e-22 1.548478441211652e-23 3.873503008524655e-25 9.227621982096665e-27 2.098223955943776e-28 4.563424055950103e-30 9.511097932712488e-32 1.902951751891381e-33 3.660826744416801e-35 6.781552053554108e-37 1.211364502417112e-38 2.089159981718163e-40

These values are accurate to almost all the significant digits shown. The backwards recursion suppresses the unwanted $Y_n(x)$ and the process is stable.

#### Center

Can we avoid having to know those two starting values? Yes, by using a matrix formulation of the Miller algorithm. Insert the normalization condition in the first row. Now $B$ is tridiagonal, except for first row, which is

[1 0 2 0 2 0 2 ... ]

Other rows have $1$ s on the super- and subdiagonal and $-2n/x$ on the diagonal of the $(n+1)$ -st row.

```
n = 30;
B = blt(n,x);
B(1,1:2) = [1 0];
B(1,3:2:n) = 2;
spy(B)
title('B')
```

The beauty of this approach is that it does not require any values of the Bessel function *a priori*. The right hand side is simply a unit vector.

rhs = eye(n,1); format long e v = B\rhs semilogy(0:n-1,v,'o-') xlabel('n') title('J_n(1)')

v = 7.651976865579665e-01 4.400505857449335e-01 1.149034849319005e-01 1.956335398266840e-02 2.476638964109954e-03 2.497577302112343e-04 2.093833800238926e-05 1.502325817436807e-06 9.422344172604494e-08 5.249250179911872e-09 2.630615123687451e-10 1.198006746303136e-11 4.999718179448401e-13 1.925616764480171e-14 6.885408200044220e-16 2.297531532210343e-17 7.186396586807487e-19 2.115375568053260e-20 5.880344573595753e-22 1.548478441211652e-23 3.873503008524655e-25 9.227621982096663e-27 2.098223955943775e-28 4.563424055950102e-30 9.511097932712487e-32 1.902951751891381e-33 3.660826744416762e-35 6.781552053355331e-37 1.211364395117367e-38 2.088559301926495e-40

#### Relative error

w = J((0:n-1)',x); relerr = (v - w)./w; semilogy(0:n-1,abs(relerr) + eps*(v==w),'o-') xlabel('n') title('relative error')

The relative error deteriorates for the last few values of `n`. This is not unexpected because we have effectively truncated a semi-infinite matrix.

#### Triangular Factors

The LU decomposition of $B$ shows no pivoting for this value of $x$, but there is fill-in. The normalization condition is mixing in with the recurrence.

[L,U] = lu(B); spy(L) title('L') snapnow spy(U) title('U')

#### References

J. C. P. Miller, "On the choice of standard solutions for a homogeneous linear equation of the second order", *Quart. J. Mech. Appl. Math.* 3, 1950, 225-235.

Walter Gautschi, "Computational Aspects of Three-Term Recurrence Relations", *SIAM Review* 9, 1967, 24-82.

Get the MATLAB code

Published with MATLAB® R2017a

## Recent Comments