Preface

This section gives an introduction to recurrences; most of them we will meet later in this chapter.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340

Glossary

Recurrences

Our object of interest is a sequence (either finite or infinite) of numbers rather than continuous functions. Generally speaking, a sequence is a particular example of a piecewise function that has constant value during some fixed interval. Traditionally, a sequence is denoted by \( {\bf a} = \left\{ a_n \right\}_{n \ge 0} = \left\{ a_0 , a_1 , a_2 , \ldots \right\} \) or simply { an }. We will consider sequences where the index runs from zero to infinity (or some finite number. However, there are known cases when context requires sequences with subscripts from minus infinity to infinity. More often than not, the elements of a sequence are defined implicitly via some equation.

A recurrence or recurrence relation is an equation that relates different members of a sequence of numbers \( {\bf a} = \left\{ a_n \right\}_{n \ge 0} = \left\{ a_0 , a_1 , a_2 , \ldots \right\} , \) where an are the values to be determined. A solution of a recurrence is any sequence that satisfies the recurrence throughout its range.
The order of a recurrence relation is the difference between the largest and smallest subscripts of the members of the sequence that appear in the equation. The general form of a recurrence relation of order p is \( a_n = f(n, a_{n-1} , a_{n-2} , \ldots , a_{n-p}) \) for some function f. A recurrence of a finite order is usually referred to as a difference equation. A full history recurrence is one that depends on all the previous functions.

Our primary interest in recurrences comes from two sources: discretization of differential equations and representation of their solutions in the form of infinite series. Since computers love recurrences, computer science and numerical analysis courses cover this topic in more detail and from different prospective than we do (and need). Now it is right time to classify recurrences further and give some examples.

If in the difference equation \( x_n = f\left( n, x_{n-1} , x_{n-2} , \ldots , x_{n-p} \right) \) of order p, the function f is linear in all its arguments, then the equation is called linear. The first and second order linear difference equations have the following form:
  • \( a_n x_n + b_n x_{n-1} = f_n \)                    first order linear equation,
  • \( a_n x_n + b_n x_{n-1} + c_n x_{n-2} = f_n \)     second order linear equation,
where {fn}, {an}, {bn}, and {cn} are given (known) sequences. When all members of these sequences except {fn} do not depend on n, the difference equation is said to have constant coefficients; otherwise, these difference equations are said to have variable coefficients.

The sequence {fn} is called a nonhomogeneous sequence, or forcing/input sequence of the difference equation. If all elements of {fn} are zero, then the linear difference equation is called a homogeneous equation; otherwise, we call it nonhomogeneous or inhomogeneous.

A difference equation usually has infinitely many solutions. In order to pin down a solution, we have to know one or more of its elements. If we have a difference equation of order p, we need to specify p sequential values of the sequence, called the initial conditions. Do, for first order recurrences, we have to specify only one element of the sequence {xn}, say the first one x0=a; for the second order difference equations we have to know two elements: x0=a and x1=b; and so forth. It should be noted that the unique solution of a recurrence may be specified by imposing restrictions other than initial conditions.

Recurrences, although a very tedious computation method by hand, is very simple to do in Mathematica. The best way to learn how to do recurrences in Mathematica are by examples, and a perfect example for this topic is the Fibonacci integer sequence.

Example: We can define the Fibonacci recurrence (which is a second order difference equation \( F_{n+2} = F_{n+1} + F_n , \quad F_0 = 0, \ F_1 =1 \) ) as

fib = (If[#1 == 1 || #1 == 2, 1, #0[#1 - 1] + #0[#1 - 2]]) &
f[n_]:=Union @@ NestList[{{0,1},{1,1}}.# &, {1, 1}, n]
fibonacciList[n_] := Module[{x = 0}, NestList[x + (x = #) &, 1, n - 1]]
fibonacciList[10]
{1, 1, 2, 3, 5, 8, 13, 21, 34, 55}
Next example from the Mathematica docs shows one way down values are more expressive (dynamic programming):
fib[n_] :=
Module[{f},
f[1] = f[2] = 1; f[i_] := f[i] = f[i - 1] + f[i - 2];
f[n] ]
fibonacciSequence[n_] :=
Module[{fPrev = 0, fNext = 1, i = 0},
While[i++ < n, {fPrev, fNext} = {fNext, fPrev + fNext}];
fNext]
Another efficient way to define the Fibonacci recurrence is to use the explicit formula (which is not available for most other nonlinear recurrences): `
\[ f_n = \frac{1}{\sqrt{5}} \left[ \left( \frac{1 + \sqrt{5}}{2} \right)^n - \left( \frac{1 - \sqrt{5}}{2} \right)^n \right] , \qquad n=0,1,2,\ldots . \]
fibos[n_] := RootReduce@(((1 + Sqrt[5])/2)^n - ((1 - Sqrt[5])/2)^n)/Sqrt[5]
Mathematica also has a special command to evaluate n-th Fibonacci number: Fibonacci[n] .

Tail recursive Fibonacci sequence generator

fiboSequence[n_, a_, b_] := fiboSequence[n - 1, b, Sow[a] + b]
fiboSequence[0, a_, b_] := Sow[a]
fiboSequence[n_] := Reap[fiboSequence[n, 0, 1]][[2, 1]]

fiboSequence[15]
{0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610}
Pattern matching Fibonacci sequence generator
fiboSequence2[n_] :=
Quiet@ReplaceRepeated[{0, 1}, {x___, a_, b_} :> {x, a, b, a + b}, MaxIterations -> n - 1]

fiboSequence2[15]
{0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610}
Matrix method:
fn[n_]:=First[MatrixPower[{{1,1},{1,0},n-1].{1,0}]
Another useful way uses LinearRecurrence:
LinearRecurrence[{1, 1}, {1, 1}, 10]
For instance, imagine that we started with a single pair of male and female rabbits. If these rabbits and their descendants reproduced at top speed ("like bunnies") for 7777 years, without any deaths, we would have enough rabbits to cover the entire state of Rhode Island.    ■

Example: The number of Dyck words with n parenthesis pairs satisfies the full history recurrence:

\[ C_{n+1} = \sum_{i=0}^n C_i C_{n-i} , \qquad C_0 =1 . \]
Dyck words and language are named after the German mathematician Walther von Dyck.    ■

There are known two basic methods to find a solution of homogeneous constant coefficient difference equation \( x_{n+p} = a_p x_{n+p-1} + a_{p-1} x_{n+p-2} + \cdots + a_1 x_n : \) test substitution \( x_{n} = r^n \) and method of generating functions. The former reduces the problem to a polynomial equation, called the characteristic equation, and the solution depends on its roots: whether they are real, multiple, or complex. The latter is more powerful and could be applied to nonhomogeneous equations as well. Till some extend, the generating function method resembles application of the Laplace transform for solving differential equations. We illustrate their applications in the following examples.

Example: Consider a homogeneous second order difference equation subject to some initial conditions:

\[ y_{n+2} = 6\,y_{n+1} - y_n , \qquad n=0,1,2,\ldots , \qquad y_0 = 1, \quad y_1 =3. \]
Mathematica has two commands to solve the above recurrence:
LinearRecurrence[{6, -1}, {1, 3}, 10]
{1, 3, 17, 99, 577, 3363, 19601, 114243, 665857, 3880899}
Here the first pair represents the coefficients of the difference equation of order 2; next pair tells Mathematica the first two initial values of the recurrence. The final number (10) is the number of terms the user wants to see as output. Another option:
Transpose[ NestList[Flatten[{Rest[#], ListCorrelate[{-1, 6}, #]}] &, {1, 3}, 10]][[1]]
{1, 3, 17, 99, 577, 3363, 19601, 114243, 665857, 3880899, 22619537}
So we get exactly the same answer for the first 10 terms in the required sequence. On the other hand, Mathematica has a special command to determine the general solution:
RSolve[{x[n+2]-6*x[n+1]+x[n]==0, x[0]==1, x[1]==3},x[n],n]
{{x[n] -> 1/2 ((3 - 2 Sqrt[2])^n + (3 + 2 Sqrt[2])^n)}}
So now we have the general formula:
\[ y_n = \frac{1}{2} \left( 3- 2\sqrt{2} \right)^n + \frac{1}{2} \left( 3+ 2\sqrt{2} \right)^n , \qquad n=0,1,2,\ldots . \]

Now we show how to apply two theoretical methods. Substituting \( y_n = r^n into the given recurrence, we get

\[ r^{n+2} = 6\,r^{n+1} - r^n \qquad \Longrightarrow \qquad r^2 -6r +1 =0. \]
The algebraic equation \( r^2 -6r +1 =0 \) is called the characteristic equation for the given recurrence. Correspondingly, the polynomial \( r^2 -6r +1 \) is referred to as the characteristic polynomial. Since the characteristic equation has two distinct real roots
Solve[r^2 - 6*r + 1 == 0, r]
{{r -> 3 - 2 Sqrt[2]}, {r -> 3 + 2 Sqrt[2]}}
the general solution of the given difference equation of the second order is
\[ y_n = A\,\left( 3+2\sqrt{2} \right)^n + B \,\left( 3-2\sqrt{2} \right)^n , \qquad n=0,1,2,\ldots ; \]
where A and B are some constants to be determined. Substitution of the general form into the initial conditions yields
\[ y_0 = A+ B =1, \qquad y_1 = A\,\left( 3+2\sqrt{2} \right) + B \,\left( 3-2\sqrt{2} \right) =3 . \]
This system of algebraic equations is easy to solve: \( A = B = \frac{1}{2} . \)

Now we apply the generating function method. Let Y(z) be the generating function for the (unknown) sequence yn:
\[ Y(z) = \sum_{n\ge 0} y_n z^n . \]
Multiplying both sides of the recurrence \( y_{n+2} = 6\,y_{n+1} - y_n \) by zn+2 and summing the results, we obtain
\[ \sum_{n\ge 0} y_{n+2} z^{n+2} = 6 \sum_{n\ge 0} y_{n+1} z^{n+2} - \sum_{n\ge 0} y_{n} z^{n+2} . \]
In the first sum at left hand side, we shift the index of summation by setting k = n+2, which gives
\[ \sum_{n\ge 0} y_{n+2} z^{n+2} = \sum_{k\ge 2} y_k z^{k} = \sum_{k\ge 0} y_k z^{k} - y_0 - y_1 z = Y(z) - 1-3z. \]
Similar for two other sums:
\begin{eqnarray*} \sum_{n\ge 0} y_{n+1} z^{n+2} &=& \sum_{k\ge 1} y_k z^{k+1} = z\,\sum_{k\ge 1} y_k z^{k} = z\,\sum_{k\ge 0} y_k z^{k} - z\,y_0 = z\,Y(z) - z , \\ sum_{n\ge 0} y_{n} z^{n+2} &=& z^2 sum_{n\ge 0} y_{n} z^{n} = z^2 Y(z) . \end{eqnarray*}
These formulas allow us to rewrite the original equation for sums as
\[ Y(z) - 1-3z = 6 \left( z\,Y(z) - z \right) - z^2 Y(z) \qquad\Longrightarrow \qquad Y(z) \left[ 1 -6z + z^2 \right] = 1+3z -6z = 1 -3z . \]
Therefore, the generating function for sequence of numbers that satisfy the given recurrence and initial conditions becomes
\[ Y(z) = \frac{1-3z}{1 -6z + z^2} = \frac{1}{1 -6z + z^2} - \frac{3z}{1 -6z + z^2} . \]
To find the value of yn, one needs to extract the n-th coefficient:
\[ y_n = \left[ z^n \right] Y(z) = \left[ z^n \right] \frac{1}{1 -6z + z^2} - \left[ z^n \right] \frac{3z}{1 -6z + z^2} . \]
Recalling the generating function for the Chebyshev polynomial of the second kind
\[ U(z) = \frac{1}{1 -6z + z^2} = \sum_{n\ge 0} U_n (3)\,z^n , \]
we see that our solution can be represented as a linear combination of Chebyshev polynomials evaluated at x = 3:
\[ y_n = U_n (3) -3\, U_{n-1} (3) , \qquad n=1,2,\ldots . \]
Now we check with Mathematica
CoefficientList[Series[(1-3*z)/(1-6 z + z^2), {z, 0, 10}],z]
{1, 3, 17, 99, 577, 3363, 19601, 114243, 665857, 3880899, 22619537}
Compare with Chebyshev polynomials
Table[ChebyshevU[n,3] -3*ChebyshevU[n-1,3], {n,1,10}]
{1, 3, 17, 99, 577, 3363, 19601, 114243, 665857, 3880899, 22619537}
So we see that our formula is correct. ■

Example: Consider nonhomogeneous second order difference equation

\[ y_{n+2} = 6\,y_{n+1} - y_n +n , \qquad n=0,1,2,\ldots , \qquad y_0 = 1, \quad y_1 =3. \]
■

Example: Consider a homogeneous third order difference equation

\[ y_{n+3} = 6\,y_{n+2} - 12\,y_{n+1} + 8\,y_n , \qquad n=0,1,2,\ldots , \qquad y_0 = 1, \quad y_1 =6 , \quad y_2 =8. \]
The characteristic polynomial of the given recurrence relation is \( r^3 -6\,r^2 +12\,r -8 = (r-2)^3 . \) So it has only one triple root (of multiplicity 3) r = 2. Now we seek the general solution in the form
\[ y_{n} = a\,2^n + bn\,2^n + cn^2 2^n , \qquad n=0,1,2,\ldots , \]
with some unknown yet coefficients a, b, and c. To determine them, we use the initial conditions:
\[ y_{0} = a =1, \quad y_{1} = 2a + 2b + 2c =6, \quad y_{2} = 4a +8b +16c = 8. \]
Solving with Mathematica, we obtain
Solve[{b + c == 3, 8*b + 16 c == 8}, {b, c}]
{{b -> 5, c -> -2}}
\[ y_{n} = 2^n + 5n\,2^n -2n^2 2^n = 2^n \left( 1 + 5n -2n^2 \right) , \qquad n=0,1,2,\ldots . \]
■

Example: Consider a homogeneous second order difference equation

\[ y_{n+2} = 2\,y_{n+1} - 10\,y_n , \qquad n=0,1,2,\ldots , \qquad y_0 = 1, \quad y_1 =3. \]
The characteristic polynomial of the given recurrence relation is \( r^2 -2r +10 =0, \) which has two complex conjugate roots: \( r = 1 \pm 3{\bf j} . \) We represent these wo roots in polar form:
\[ r = \sqrt{10} \,e^{\pm {\bf j} \theta} , \qquad \mbox{where } \theta = \arctan 3 \approx 1.24905 . \]
Extracting real and imaginary parts from \( r^n , \) we obtain two linearly independent real real-valued solutions:
\[ 10^{n/2} \,\cos n\theta \qquad \mbox{and} \qquad 10^{n/2} \,\sin n\theta . \]
Therefore, the general solution is their linear combination
\[ y_n = a\,10^{n/2} \,\cos n\theta + b\,10^{n/2} \,\sin n\theta , \]
where coefficients a and b are determined from the initial conditions:
\[ y_0 = a = 1, \qquad y_1 = 10^{1/2} \left( \cos \theta + b\,\sin \theta \right) =3 . \]
So
\[ a=1, \qquad b= 3\cdot 10^{-1/2} \left( \sin\theta \right)^{-1} - \cot \theta = \frac{2}{3} . \]
Hence, the required solution is
\[ y_n = 10^{n/2} \left( \cos n\theta + \frac{2}[3} \,\sin n\theta \right) . s \]
■

First order Recurrences

Consider some population when its size is measured by possible fixed interval in time. Thus, we can describe the population size by a sequence { xn }, with x0 denoting the initial population size, x1 the population size at the next generation measured at time t1, x2 the population size at the next generation measured at time t2, and so on. Usually the time interval between generations is taken to be a constant.

For example, suppose the population changes only through births and deaths, so that xn+1 - xn is the number of births minus the number of deaths over the time interval from tn to tn+1. Suppose further that the birth and death rates are constants b and d, respectively. Then

\[ x_{n+1} - x_n = (b-d)\,x_n \qquad\mbox{or}\qquad x_{n+1} = \left( 1+b-d \right) x_n = r\,x_n , \]
where r=1+b-d. This is the linear homogeneous difference equation of order 1. If initial population x0 is know, then population at any time tn is: \( x_n = r^n x_0 . \)

Consider another recurrence that generates the Pell numbers:

y_{n+2} = 2 y_{n+1} + y_n ,         y_0 =0, y_1 =1,     n=0,1,2,.... .

Then these numbers can be generated as follows:

a = 1; b = 0; c = 0; lst = {b}; Do[c = a + b + c; AppendTo[lst, c];
�� a = b; b = c, {n, 20}];
lst
CoefficientList[Series[x/(1 - 2*x - x^2), {x, 0, 20}], x] (* 20 numbers *)
Expand[Table[((1 + Sqrt[2])^n - (1 - Sqrt[2])^n)/(2 Sqrt[2]), {n, 0, 20}]]

  2. Rivera-Figueroa, A. and Rivera-Rebolledo, J.M., A new method to solve the second-order linear difference equations with constant coefficients, International Journal of Mathematical Education in Science and Technology, Volume 47, 2016 - Issue 4, pp. 636--649. https://doi.org/10.1080/0020739X.2015.1091517
  3. Tisdell, C.C., Critical perspectives of the ‘new’ difference equation solution method of Rivera-Figueroa and Rivera-Rebolledo, International Journal of Mathematical Education in Science and Technology, 2018, https://doi.org/10.1080/0020739X.2018.1469796

 

Return to Mathematica page

Return to the main page (APMA0330)
Return to the Part 1 (Plotting)
Return to the Part 2 (First Order ODEs)
Return to the Part 3 (Numerical Methods)
Return to the Part 4 (Second and Higher Order ODEs)
Return to the Part 5 (Series and Recurrences)
Return to the Part 6 (Laplace Transform)
Return to the Part 7 (Boundary Value Problems)