Solving linear recurrences

Lesson 5 solved first-order recurrences by unrolling them. That method breaks down for recurrences such as the Fibonacci recurrence, where each term depends on two earlier ones: substituting the recurrence into itself produces a growing tangle of terms. This lesson develops a general method for the most important class of recurrences, the linear recurrences with constant coefficients. The method reduces the problem to solving a polynomial equation, and it explains why so many counting sequences grow like powers of a fixed number. It also produces one of the most surprising formulas in elementary mathematics, an exact expression for the Fibonacci numbers involving the irrational number \(\sqrt{5}\).

Linear homogeneous recurrences

A linear homogeneous recurrence of order \(k\) with constant coefficients has the form

\[ a_n = c_1 a_{n-1} + c_2 a_{n-2} + \cdots + c_k a_{n-k}, \]

where \(c_1, \dots, c_k\) are fixed numbers with \(c_k \ne 0\). Linear means the earlier terms appear only to the first power, with no products such as \(a_{n-1} a_{n-2}\). Homogeneous means there is no extra term that depends only on \(n\); the Tower of Hanoi recurrence \(H_n = 2H_{n-1} + 1\) is not homogeneous because of the \(+1\). Constant coefficients means the \(c_i\) do not depend on \(n\). The Fibonacci recurrence \(F_n = F_{n-1} + F_{n-2}\) is the standard example, of order 2 with \(c_1 = c_2 = 1\).

Two properties make these recurrences tractable. First, if two sequences satisfy the recurrence, so does any combination \(\alpha a_n + \beta b_n\), because the recurrence is linear: substituting the combination into each side gives the same combination of the two sides. In the language of linear algebra, the solutions form a vector space. Second, a solution is determined by its first \(k\) terms. So the space of solutions has dimension \(k\), and if we can find \(k\) suitably different solutions, every solution is a combination of them.

The characteristic equation

Look for solutions of the simplest possible form, \(a_n = r^n\) for a non-zero number \(r\). Substituting into a second-order recurrence \(a_n = c_1 a_{n-1} + c_2 a_{n-2}\) gives \(r^n = c_1 r^{n-1} + c_2 r^{n-2}\), and dividing by \(r^{n-2}\),

\[ r^2 = c_1 r + c_2, \qquad\text{that is}\qquad r^2 - c_1 r - c_2 = 0. \]

This is the characteristic equation of the recurrence. The sequence \(r^n\) is a solution exactly when \(r\) is a root of it. For an order-\(k\) recurrence the characteristic equation is \(r^k - c_1 r^{k-1} - \cdots - c_k = 0\), a polynomial of degree \(k\).

Theorem (distinct roots). Suppose the characteristic equation of \(a_n = c_1 a_{n-1} + c_2 a_{n-2}\) has two distinct roots \(r_1\) and \(r_2\). Then every solution has the form

\[ a_n = \alpha r_1^n + \beta r_2^n \]

for constants \(\alpha\) and \(\beta\), which are determined by the initial conditions.

Proof. Each of \(r_1^n\) and \(r_2^n\) is a solution, so by linearity every combination \(\alpha r_1^n + \beta r_2^n\) is a solution. Now take any solution \(a_n\). Choose \(\alpha\) and \(\beta\) so that the combination agrees with \(a_n\) at \(n = 0\) and \(n = 1\): we need \(\alpha + \beta = a_0\) and \(\alpha r_1 + \beta r_2 = a_1\). This pair of linear equations has a unique solution because its determinant, \(r_2 - r_1\), is not zero. Then \(a_n\) and the combination are two solutions of the recurrence with the same first two terms. Since a solution is determined by its first two terms, they agree for every \(n\).

Binet's formula for the Fibonacci numbers

For \(F_n = F_{n-1} + F_{n-2}\), the characteristic equation is \(r^2 - r - 1 = 0\), with roots

\[ \varphi = \frac{1 + \sqrt{5}}{2} \approx 1.618, \qquad \psi = \frac{1 - \sqrt{5}}{2} \approx -0.618. \]

The number \(\varphi\) is the golden ratio. The general solution is \(F_n = \alpha\varphi^n + \beta\psi^n\). The initial conditions \(F_0 = 0\) and \(F_1 = 1\) give \(\alpha + \beta = 0\) and \(\alpha\varphi + \beta\psi = 1\), so \(\alpha = -\beta\) and \(\alpha(\varphi - \psi) = 1\). Since \(\varphi - \psi = \sqrt{5}\), we get \(\alpha = 1/\sqrt{5}\) and

\[ F_n = \frac{1}{\sqrt{5}}\Bigl(\varphi^n - \psi^n\Bigr) = \frac{1}{\sqrt{5}}\biggl(\Bigl(\frac{1 + \sqrt{5}}{2}\Bigr)^n - \Bigl(\frac{1 - \sqrt{5}}{2}\Bigr)^n\biggr). \]

This is Binet's formula, named after Jacques Binet, who published it in 1843, although Abraham de Moivre and others knew it earlier. It is remarkable that a formula full of square roots of 5 always produces a whole number. Because \(|\psi| < 1\), the term \(\psi^n/\sqrt{5}\) is less than \(1/2\) in absolute value for every \(n \ge 0\) and shrinks rapidly, so \(F_n\) is simply the nearest integer to \(\varphi^n / \sqrt{5}\). It follows that the Fibonacci numbers grow exponentially, by a factor approaching \(\varphi\) at each step: \(F_{n+1}/F_n \to \varphi\).

The same phenomenon holds for any linear recurrence: the root of largest absolute value, the dominant root, controls the growth. In the example above, \(a_n = 2 \cdot 3^n - 2^n\) grows like \(2 \cdot 3^n\), and the ratio \(a_{n+1}/a_n\) tends to 3. This is how the domino tilings and bit strings of lesson 5 grow: by a factor of about 1.618 for each unit of length.

Repeated roots

If the characteristic equation has a repeated root \(r\), then \(r^n\) gives only one solution, and we need a second. The answer is \(n r^n\).

Theorem (repeated root). If \(r^2 - c_1 r - c_2 = (r - s)^2\), then every solution of \(a_n = c_1 a_{n-1} + c_2 a_{n-2}\) has the form \(a_n = (\alpha + \beta n)s^n\).

Proof sketch. Expanding \((r - s)^2\) shows \(c_1 = 2s\) and \(c_2 = -s^2\). Check that \(b_n = n s^n\) satisfies the recurrence: \(c_1 b_{n-1} + c_2 b_{n-2} = 2s(n - 1)s^{n-1} - s^2(n - 2)s^{n-2} = s^n(2n - 2 - n + 2) = n s^n = b_n\). The solutions \(s^n\) and \(n s^n\) can match any two initial values, because the equations \(\alpha = a_0\) and \((\alpha + \beta)s = a_1\) always have a solution when \(s \ne 0\). The rest of the argument is as for distinct roots.

For recurrences of higher order the same pattern holds: each distinct root \(r\) of multiplicity \(m\) contributes solutions \(r^n, n r^n, \dots, n^{m-1} r^n\), and the general solution is a combination of all of them. The roots may also be complex numbers. For example \(a_n = -a_{n-2}\) has characteristic equation \(r^2 + 1 = 0\) with roots \(\pm i\), and its real solutions are periodic with period 4, such as 1, 0, \(-1\), 0, 1, 0, ... Complex roots always produce oscillation of this kind.

Non-homogeneous recurrences

A non-homogeneous linear recurrence has an extra term: \(a_n = c_1 a_{n-1} + c_2 a_{n-2} + f(n)\). The Tower of Hanoi recurrence, with \(f(n) = 1\), is an example. The key fact, again from linearity, is that the difference of two solutions solves the homogeneous recurrence. So the general solution is one particular solution plus the general solution of the homogeneous recurrence.

Finding a particular solution is usually done by guessing its form from \(f(n)\). If \(f(n)\) is a constant, try a constant; if it is a polynomial of degree \(d\), try a polynomial of degree \(d\); if it is \(t^n\), try a multiple of \(t^n\). If the guess happens to solve the homogeneous recurrence, multiply it by \(n\) and try again. For the Tower of Hanoi, try the constant \(p\): \(p = 2p + 1\) gives \(p = -1\). The homogeneous recurrence \(H_n = 2H_{n-1}\) has solutions \(\alpha 2^n\). So \(H_n = \alpha 2^n - 1\), and \(H_1 = 1\) gives \(\alpha = 1\): \(H_n = 2^n - 1\), as before.

A computational check

The method is easy to check by computer. The code below solves any second-order recurrence with distinct real roots from its coefficients and initial values, compares the closed form with direct iteration, and shows the ratio of successive Fibonacci numbers approaching the golden ratio. It also computes \(F_n\) with exact integer arithmetic through the matrix identity \(\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}\), which with repeated squaring needs only about \(\log_2 n\) matrix multiplications.

from math import sqrt

def iterate(c1, c2, a0, a1, n):
    seq = [a0, a1]
    while len(seq) <= n:
        seq.append(c1 * seq[-1] + c2 * seq[-2])
    return seq

def closed_form(c1, c2, a0, a1):
    disc = c1 * c1 + 4 * c2
    r1, r2 = (c1 + sqrt(disc)) / 2, (c1 - sqrt(disc)) / 2
    beta = (a1 - a0 * r1) / (r2 - r1)
    alpha = a0 - beta
    return lambda n: alpha * r1 ** n + beta * r2 ** n, (r1, r2)

for (c1, c2, a0, a1) in [(5, -6, 1, 4), (1, 1, 0, 1), (2, 3, 1, 1)]:
    f, roots = closed_form(c1, c2, a0, a1)
    direct = iterate(c1, c2, a0, a1, 10)
    sign = "+" if c2 >= 0 else "-"
    print(f"a_n = {c1}a_(n-1) {sign} {abs(c2)}a_(n-2), roots {roots[0]:.4f}, {roots[1]:.4f}")
    print("   iteration  :", direct)
    print("   closed form:", [round(f(n)) for n in range(11)])

fib = iterate(1, 1, 0, 1, 30)
print("F(30)/F(29) =", fib[30] / fib[29], " golden ratio =", (1 + sqrt(5)) / 2)

def mat_mult(A, B):
    return [[A[0][0]*B[0][0] + A[0][1]*B[1][0], A[0][0]*B[0][1] + A[0][1]*B[1][1]],
            [A[1][0]*B[0][0] + A[1][1]*B[1][0], A[1][0]*B[0][1] + A[1][1]*B[1][1]]]

def fib_matrix(n):
    result, base = [[1, 0], [0, 1]], [[1, 1], [1, 0]]
    while n:
        if n & 1:
            result = mat_mult(result, base)
        base = mat_mult(base, base)
        n >>= 1
    return result[0][1]

print("F(200) =", fib_matrix(200))

a_n = 5a_(n-1) - 6a_(n-2), roots 3.0000, 2.0000
   iteration  : [1, 4, 14, 46, 146, 454, 1394, 4246, 12866, 38854, 117074]
   closed form: [1, 4, 14, 46, 146, 454, 1394, 4246, 12866, 38854, 117074]
a_n = 1a_(n-1) + 1a_(n-2), roots 1.6180, -0.6180
   iteration  : [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
   closed form: [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
a_n = 2a_(n-1) + 3a_(n-2), roots 3.0000, -1.0000
   iteration  : [1, 1, 5, 13, 41, 121, 365, 1093, 3281, 9841, 29525]
   closed form: [1, 1, 5, 13, 41, 121, 365, 1093, 3281, 9841, 29525]
F(30)/F(29) = 1.6180339887482036  golden ratio = 1.618033988749895
F(200) = 280571172992510140037611932413038677189525

Floating-point arithmetic loses accuracy for large \(n\): Binet's formula in double precision gives correct Fibonacci numbers only up to about \(F_{70}\). For exact values the integer matrix method, or simple iteration with Python's unlimited integers, is the right tool. Closed forms are for understanding growth; exact large values come from integer computation.