Recurrence Relations

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

Many functions which arise as solutions to second-order differential equations satisfy recurrence relations of the form

(Subscript[f, n])[
        x]  ==  α[n, x] (Subscript[f, n - 1])[x]  +  β[n, x] 
          (Subscript[f, n - 2])[x]

where α[n, x] and β[n, x] are simple coefficients that depend upon the form of the underlying differential equation.  For example, the associated Legendre functions (Subscript[P, ℓ,m])[x] satisfy the three-term recurrence relation

(ℓ - m) 
        (Subscript[P, ℓ, m])[x]  ==  x(2  ℓ - 1) 
          (Subscript[P, ℓ - 1, m])[x]  -  (ℓ + m - 1) 
          (Subscript[P, ℓ - 2, m])[x]

with respect to ℓ.  Thus. given the closed-form initial values:

(Subscript[P, m, m])[
        x]  =  (((-1))^m) (2  m - 
              1)!! (((1 - (x^2)))^(m/2))
(Subscript[P, m + 1, m])[
        x]  =  x(2  m + 1) (Subscript[P, m, m])[x]

we can compute (Subscript[P, ℓ, m]) for ℓ > m + 1 by applying the recursion relation ℓ - m - 1 times beginning with the initial values (Subscript[P, m, m]) and (Subscript[P, m + 1, m]).   This method can be used to obtain either symbolic or numerical values of the function.

Traditional procedural programming languages, such as Fortran or C, rely on recurrence relations for numerical evaluation of many functions.  The implementation of this type algorithm always involves some form of loop statement.  Mathematica provides similar tools, but also provides more transparent recursive methods which appear to be closer to the more elegant formal statement of the algorithm.  In this notebook we compare several methods for implementing recursive algorithms for Legendre functions.  We also investigate the speed, stability, and accuracy of these methods.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1];

In order to compare our results for associated Legendre functions with the built-in function SphericalHarmonicY, it is useful to define a normalization function as follows.

norm[l_, 
        m_]  :=   Sqrt[ ((2  l + 1)/(4  π)) (((l - m)!)/((l + m)!))]

Loop method

Recurrence relations can be implemented in several ways.  We begin by presenting a Legendre function pleg based upon a Do loop.  This implementation should look quite familiar to anyone with experience programming in C or Fortran.  The procedure begins by evaluating two starting values.  If ℓ = m or ℓ = m + 1 it simply returns one of the starting values; otherwise, the recursion relation is developed in a Do loop.

Note: it is good practice to document your functions and to check for legal arguments.  In this case, no values are returned for illegal arguments.

ClearAll[pleg]

\(\(\(pleg::usage\  = 
	\*""pleg[ℓ,m,x] computes the Legendre polynomial (Subscript[P, ℓ, m])[x] using the two-term recursion relation."";\)\(
\)
    \)\)
pleg[ℓ_Integer  /; ℓ ≥ 0,  
          m_Integer  /; m ≥ 0,  x_]  :=  
        Block[{pm2, pm1, pm, 
            k}, 
pm2  =  (((-1))^m) (2  m - 
                    1)!! (((1 - (x^2)))^(m/2)); 
pm1  =  x (2  m + 
                  1) pm2; 
Which[
NumberQ[x]  &&  Abs[x] > 1,  
            Return[], 
m > ℓ,  
            Return[],  
ℓ == m,  
            Return[pm2],  
ℓ == m + 1, 
            Return[pm1], 
ℓ > 
              m + 1,  
Do[ 
pm  =  (x (2  k - 1) pm1  -  (k + m - 1) pm2) / (k - 
                      m); 
pm2  =  
                pm1; 
pm1  =  pm, {k, 
                m + 2,  ℓ}]; Return[pm]];];
pleg[ℓ_Integer  /; ℓ ≥ 
            0,   m_Integer  /; m < 0,  
        x_]  :=  (((-1))^m) (((ℓ + m)!)/((ℓ - m)!)) pleg[ℓ, -m, x]

Note that although this pleg can be used for either symbolic or numerical values of x, we anticipate that it is most useful for numerical arguments and decline to require simplification, which would inflate the evaluation time considerably.  However, to compare pleg and SphericalHarmonicY symbolically, it is helpful to simplify first; that step is superfluous for numerical evaluation.

compare[f_, l_, x_]  :=  
      Table[{ f[l, m, Cos[x]] // Simplify, 
            SphericalHarmonicY[ l, m, x, 0]/norm[l, m]}, {m, -l, l}]  // 
        TableForm

First, we verify that pleg reproduces the symbolic expressions for spherical harmonics.

compare[pleg, 4, x]

Next, we compare numerical values with SphericalHarmonicY.

compare[pleg, 20, 0.1]  // N[#, 16] &

Recursive implementation

Direct recursion

Mathematica permits recursive definitions to be made directly.  Thus, expressing the recurrence relation in terms of SetDelayed functions, the recurrence relation can be implemented without using loop counters.  Hence, the recursive implementation is much more compact.   However, it remains good practice to verify that the arguments are legal.

ClearAll[ pleg2 ]

\(\(pleg2::usage\  = \*""pleg2[ℓ,m,x] computes the Legendre polynomial (Subscript[P, ℓ, m])[x] using the two-term recursion relation in a recursive procedure."";\)
\)
pleg2[ℓ_Integer  /; ℓ ≥ 0, 
          m_Integer  /; m ≥ 0, x_]  :=  
        Which[
NumberQ[x]  &&  Abs[x] > 1,  
          Return[], 
m > ℓ,  
          Return[], 
ℓ == 
            m, (((-1))^m) (2  m - 
                  1)!! (((1 - (x^2)))^(m/2)), 
ℓ == m + 1,  
          x (2  m + 1) pleg2[m, m, x], 
ℓ > 
            m + 1, (x (2  ℓ - 1) pleg2[ℓ - 1, m, 
                    x]  -  (ℓ + m - 1) 
                  pleg2[ℓ - 2, m, x]) / (ℓ - 
                m)];
pleg2[ℓ_Integer  /; ℓ ≥ 
            0,   m_Integer  /; m < 0,  
        x_]  :=  (((-1))^m) (((ℓ + m)!)/((ℓ - m)!)) pleg2[ℓ, -m, x]

Once again we check this version of the Legendre function for both symbolic and numerical accuracy.

LegendreP[20, -10, x]  // Simplify

pleg2[20, -10, x] // Simplify

compare[pleg2, 20, 0.1]  // N[#, 16] &

Although this algorithm is quite simple to state, it appears to be prohibitively slow — much slower than the direct do loop.  To verify this observation, we use the Timing function.

Timing[pleg[20, 2, 0.99]]

Timing[pleg2[20, 2, 0.99]]

The problem is that recursive evaluation can require many function calls.  Fortunately, there is a simple solution to this problem.

Dynamic programming

An improved algorithm for recursive functions uses the construction
	f[x_] := f[x] = rhs
known as dynamic programming.  With this procedure each evaluation of the function is stored as a rule associated with the function.  Therefore, the recursive algorithm needs to evaluate each occurrence of the same function call only once rather than every time that it is encountered, which often produces dramatic improvements in speed.

ClearAll[pleg3]

\(\(pleg3::usage\  = \*""pleg3[ℓ,m,x] computes the Legendre polynomial (Subscript[P, ℓ, m])[x] using the two-term recursion relation in a recursive procedure."";\)
\)
\(\(pleg3[ℓ_Integer\  /; ℓ ≥ 0, 
          m_Integer\  /; m ≥ 0, 
          x_]\  := \ \(pleg3[ℓ, m, x]\  = 
          Which[
	NumberQ[x]\  && \ Abs[x] > 1, \ Return[], 
	
            m > ℓ, \ 
            Return[], 
	ℓ == 
              m, \(\((\(-1\))\)\^m\) \(\((2  m - 
                    1)\)!!\)\ \((1 - x\^2)\)\^\(m/2\), 
	ℓ == m + 1, \ 
            x\ \((2  m + 1)\)\ pleg3[m, m, x], 
	ℓ > 
              m + 1, \((x\ \((2  ℓ - 1)\)\ pleg3[ℓ - 1, m, 
                      x]\  - \ \((ℓ + m - 1)\) 
                    pleg3[ℓ - 2, m, x])\)\ /\ \((ℓ - 
                  m)\)]\);\)\)
pleg3[ℓ_Integer  /; ℓ ≥ 
            0,   m_Integer  /; m < 0,  
        x_]  :=  (((-1))^m) (((ℓ + m)!)/((ℓ - m)!)) pleg3[ℓ, -m, x]

?pleg3

The evaluation of this function requires a time intermediate between pleg and pleg2.

Timing[pleg3[20, 2, 0.99]]

The second evaluation of the same function is practically instantaneous because Mathematica remembers its previous result and does not repeat the work that was needed to generate it.

Timing[pleg3[20, 2, 0.99]]

Actually, not only does it remember the result we requested, it also remembers all of its intermediate results as well.  Therefore, any of those can be recalled without computational effort.  Of course, one now trades memory against speed and may find that the kernel eventually exits due to lack of memory.

?? pleg3

Timing[ 
        pleg[20, 2, 
          x] ][[1]]

Timing[ 
        pleg2[20, 2, 
          x] ][[1]]

Timing[ 
        pleg3[20, 2, 
          x] ][[1]]

Direct iteration is clearly the fastest method, and nondynamic recursion the slowest, with dynamic programming providing intermediate speed and better clarity.  Obviously, when speed is at a premium, the traditional algorithms are often best.

Counting recursive function calls

The function Trace can be used to investigate how an expression is evaluated by Mathematica.  When used with a single argument expr, it produces a nested list of the expressions it encountered during the evaluation of expr.  Each of these intermediate expression is wrapped in HoldForm to prevent further evaluation of the result of Trace.  For example, a trace for the simple recursive evaluation of pleg2 gives the following result.

stuff  =  Trace[pleg2[3, 2, x]]

FullForm shows that HoldForm prevents the various expressions in this trace from being evaluated themselves.

stuff // FullForm

Although the complete trace is instructive, it is also rather indigestible.  We can determine the total number of times that pleg2 is evaluating by giving a second argument to Trace[expr,pattern], which provides a list of all expressions involving pattern encountered during the evaluation of expr.

stuff  =  
      Trace[ pleg2[3, 2, x],  pleg2[__] ] // Flatten

Therefore, the following function counts the total number of calls to pleg2 needed for a particular recursive evaluation.

totalCalls2[ℓ_, m_]  :=  
      Length[ Trace[pleg2[ℓ, m, x],  pleg2[__]] // 
          Flatten ]

We now find that the total number of function calls can become quite large for recursive evaluations.

totalCalls2[20, 2]

By contrast, the number of calls used with dynamic programming is much smaller.

totalCalls3[ℓ_, m_]  :=  
      Length[ Trace[pleg3[ℓ, m, x],  pleg3[__]] // 
          Flatten ]

totalCalls3[30, 2]

It is also instructive to tabulate the calls to functions of various degree.  For this we use a replacement rule in the argument list of Trace to produce a simple list of degree values.

partialCalls2[ℓ_, m_]  :=  
      Block[{k, kList}, 
        kList  =  
          Flatten[Trace[pleg2[ℓ, m, x],  
                pleg2[k_, __] -> k ] /. HoldForm -> List]; 
        Table[{k, Count[kList, k]}, {k, 0, ℓ}]]

partialCalls2[20, 2]  // TableForm

Similarly, we can produce a list of partial calls for pleg3.  However, to obtain a nontrivial result, it is necessary to clear all of the rules associated with pleg3 and to execute again its definition cell.

partialCalls3[ℓ_, m_]  :=  
      Block[{k, kList}, 
        kList  =  
          Flatten[Trace[pleg3[ℓ, m, x],  
                pleg3[k_, __] -> k ] /. HoldForm -> List]; 
        Table[{k, Count[kList, k]}, {k, 0, ℓ}]]

partialCalls3[20, 2]  // TableForm

Stability

In general one must carefully consider the stability  of recursion relations.  Consider a three-term linear recurrence relation with respect to a degree parameter n, such as the one used above to generate associated Legendre polynomials (with n represented by ℓ above).  Such a relation will have two linearly independent solutions, corresponding to the two solutions to the second-order differential equation, but only one of those solutions corresponds to the desired function while the other may be described as spurious.  Suppose that we propose to employ the recursion relation in the direction of increasing n.  Given that machine numbers have finite precision, the starting values are not truly exact but rather can be represented as two linear combinations of the true solution with a small admixture of the spurious solution that originates due to truncation.  In the best of all possible worlds, the spurious solution will be exponentially damped as the recursion relation is applied so that only the true solution emerges and we obtain an accurate result.  However, Murphy's Law dictates that it is likely that the spurious admixture grows exponentially as the recurrence relation is applied so that soon roundoff error is amplified and overwhelms the true solution we seek.  Under those unfortunate circumstances the recurrence relation cannot be applied in that direction.  Fortunately, if we apply the relation in the opposite direction instead, the spurious solution will be exponentially damped and the method will succeed.
	Thus, we can classify recurrence relations as either stable, unstable, or neutral with respect to increasing n according to the behavior of roundoff errors as the relation is iterated.  The relation is stable if the spurious solution is exponentially damped or unstable if it grows exponentially.  The relation is considered neutral if the spurious solution has a power-law behavior, which is sometimes described as exponentially neutral.  The behavior of many common recurrence relations can be found in standard references, such as Abramowitz and Stegun, which report that upward recursion is stable for Legendre functions (Subscript[P, n])[x] but unstable for Bessel functions (Subscript[J, n])[x].  Conversely, upward recursion is unstable for (Subscript[Q, n])[x] and stable for (Subscript[Y, n])[x], which are the complementary solutions to the Legendre or Bessel equations.  Hence, we can confidently employ upward recursion for Legendre functions of type (Subscript[P, n]), as done herein, but must use downward recursion for Bessel functions of type (Subscript[J, n]).  However, downward recursion is a little more complicated because starting values are not immediately available.

We can easily show that starting errors remain small as the Legendre recurrence relation is iterated by comparing results with and without deliberately introducing a small error term.

ClearAll[testLegendre];
testLegendre[0, x_, eps_]  :=  1  +  eps;
testLegendre[1, x_, eps_]  :=  x;
testLegendre[n_Integer  /; n > 1, x_, 
        eps_]  := testLegendre[n, x, 
          eps]  = ( (2  n - 1) 
                x testLegendre[n - 1, x, 
                  eps]  -  (n - 1) testLegendre[n - 2, x, 
                  eps ]) / n

Table[{k, (testLegendre[k, 0.99, 0. ]  -  
              testLegendre[k, 0.99, 10*$MachineEpsilon])},  {k, 0, 20}]  // 
      TableForm

Upward recursion is clearly stable for Legendre functions.  By contrast, upward recurrence for Bessel functions is unstable.  Although it is not immediately obvious how we should generate starting values, we can use the built-in function for this purpose.  We also check below that our recurrence method reproduces the built-in function for relatively small n.  However, the errors grow as n increases.  Fortunately, Mathematica provides superior algorithms for evaluating Bessel functions.  Nevertheless, the distinction between stable and unstable recursion should now be clear.

ClearAll[testBessel];
testBessel[0, x_, eps_]  :=  BesselJ[0, x]  +  eps;
testBessel[1, x_, eps_]  :=  BesselJ[1, x];
testBessel[n_Integer  /; n > 1, x_, 
        eps_]  :=  testBessel[n, x, 
          eps]  =  ((2  n - 2)/x) testBessel[n - 1, x, eps]  -  
          testBessel[n - 2, x, eps ]

Table[{k, (testBessel[k, 0.99, 0. ]  -  
              BesselJ[k, 0.99]), (testBessel[k, 0.99, 0. ]  -  
              testBessel[k, 0.99, 10*$MachineEpsilon])}, {k, 0, 20}]  // 
      TableForm