Linear recurrence with constant coefficients

For order 1, the recurrence

$${\displaystyle a_{n}=ra_{n-1}}$$

has the solution ${\displaystyle a_{n}=r^{n}}$ with ${\displaystyle a_{0}=1}$ and the most general solution is ${\displaystyle a_{n}=kr^{n}}$ with ${\displaystyle a_{0}=k}$. The characteristic polynomial equated to zero (the characteristic equation) is simply ${\displaystyle t-r=0}$.

Solutions to such recurrence relations of higher order are found by systematic means, often using the fact that ${\displaystyle a_{n}=r^{n}}$ is a solution for the recurrence exactly when ${\displaystyle t=r}$ is a root of the characteristic polynomial. This can be approached directly or using generating functions (formal power series) or matrices.

Consider, for example, a recurrence relation of the form

$${\displaystyle a_{n}=Aa_{n-1}+Ba_{n-2}.}$$

When does it have a solution of the same general form as ${\displaystyle a_{n}=r^{n}}$? Substituting this guess (ansatz) in the recurrence relation, we find that

$${\displaystyle r^{n}=Ar^{n-1}+Br^{n-2}}$$

must be true for all ${\displaystyle n>1}$.

Dividing through by ${\displaystyle r^{n-2}}$, we get that all these equations reduce to the same thing:

$${\displaystyle {\begin{aligned}r^{2}&=Ar+B,\\r^{2}-Ar-B&=0,\end{aligned}}}$$

which is the characteristic equation of the recurrence relation. Solve for ${\displaystyle r}$ to obtain the two roots ${\displaystyle \lambda _{1}}$, ${\displaystyle \lambda _{2}}$: these roots are known as the characteristic roots or eigenvalues of the characteristic equation. Different solutions are obtained depending on the nature of the roots: If these roots are distinct, we have the general solution

$${\displaystyle a_{n}=C\lambda _{1}^{n}+D\lambda _{2}^{n}}$$

while if they are identical (when ${\displaystyle A^{2}+4B=0}$), we have

$${\displaystyle a_{n}=C\lambda ^{n}+Dn\lambda ^{n}}$$

This is the most general solution; the two constants ${\displaystyle C}$ and ${\displaystyle D}$ can be chosen based on two given initial conditions ${\displaystyle a_{0}}$ and ${\displaystyle a_{1}}$ to produce a specific solution.

In the case of complex eigenvalues (which also gives rise to complex values for the solution parameters ${\displaystyle C}$ and ${\displaystyle D}$), the use of complex numbers can be eliminated by rewriting the solution in trigonometric form. In this case we can write the eigenvalues as ${\displaystyle \lambda _{1},\lambda _{2}=\alpha \pm \beta i.}$ Then it can be shown that

$${\displaystyle a_{n}=C\lambda _{1}^{n}+D\lambda _{2}^{n}}$$

can be rewritten as[4]^{: 576–585}

$${\displaystyle a_{n}=2M^{n}\left(E\cos(\theta n)+F\sin(\theta n)\right)=2GM^{n}\cos(\theta n-\delta ),}$$

where

$${\displaystyle {\begin{array}{lcl}M={\sqrt {\alpha ^{2}+\beta ^{2}}}&\cos(\theta )={\tfrac {\alpha }{M}}&\sin(\theta )={\tfrac {\beta }{M}}\\C,D=E\mp Fi&&\\G={\sqrt {E^{2}+F^{2}}}&\cos(\delta )={\tfrac {E}{G}}&\sin(\delta )={\tfrac {F}{G}}\end{array}}}$$

Here ${\displaystyle E}$ and ${\displaystyle F}$ (or equivalently, ${\displaystyle G}$ and ${\displaystyle \delta }$) are real constants which depend on the initial conditions. Using

$${\displaystyle \lambda _{1}+\lambda _{2}=2\alpha =A,}$$

$${\displaystyle \lambda _{1}\cdot \lambda _{2}=\alpha ^{2}+\beta ^{2}=-B,}$$

one may simplify the solution given above as

$${\displaystyle a_{n}=(-B)^{\frac {n}{2}}\left(E\cos(\theta n)+F\sin(\theta n)\right),}$$

where ${\displaystyle a_{1}}$ and ${\displaystyle a_{2}}$ are the initial conditions and

$${\displaystyle {\begin{aligned}E&={\frac {-Aa_{1}+a_{2}}{B}}\\F&=-i{\frac {A^{2}a_{1}-Aa_{2}+2a_{1}B}{B{\sqrt {A^{2}+4B}}}}\\\theta &=\arccos \left({\frac {A}{2{\sqrt {-B}}}}\right)\end{aligned}}}$$

In this way there is no need to solve for ${\displaystyle \lambda _{1}}$ and ${\displaystyle \lambda _{2}}$.

In all cases—real distinct eigenvalues, real duplicated eigenvalues, and complex conjugate eigenvalues—the equation is stable (that is, the variable ${\displaystyle a}$ converges to a fixed value [specifically, zero]) if and only if both eigenvalues are smaller than one in absolute value. In this second-order case, this condition on the eigenvalues can be shown[5] to be equivalent to ${\displaystyle |A|<1-B<2}$, which is equivalent to ${\displaystyle |B|<1}$ and ${\displaystyle |A|<1-B}$.

Solving the homogeneous equation

$${\displaystyle x_{t}=a_{1}x_{t-1}+\cdots +a_{n}x_{t-n}}$$

involves first solving its characteristic polynomial

$${\displaystyle \lambda ^{n}=a_{1}\lambda ^{n-1}+\cdots +a_{n-2}\lambda ^{2}+a_{n-1}\lambda +a_{n}}$$

for its characteristic roots λ₁, ..., λ_n. These roots can be solved for algebraically if n ≤ 4, but not necessarily otherwise. If the solution is to be used numerically, all the roots of this characteristic equation can be found by numerical methods. However, for use in a theoretical context it may be that the only information required about the roots is whether any of them are greater than or equal to 1 in absolute value.

It may be that all the roots are real or instead there may be some that are complex numbers. In the latter case, all the complex roots come in complex conjugate pairs.

If all the characteristic roots are distinct, the solution of the homogeneous linear recurrence

$${\displaystyle x_{t}=a_{1}x_{t-1}+\cdots +a_{n}x_{t-n}}$$

can be written in terms of the characteristic roots as

$${\displaystyle x_{t}=c_{1}\lambda _{1}^{t}+\cdots +c_{n}\lambda _{n}^{t}}$$

where the coefficients c_i can be found by invoking the initial conditions. Specifically, for each time period for which an iterate value is known, this value and its corresponding value of t can be substituted into the solution equation to obtain a linear equation in the n as-yet-unknown parameters; n such equations, one for each initial condition, can be solved simultaneously for the n parameter values. If all characteristic roots are real, then all the coefficient values c_i will also be real; but with non-real complex roots, in general some of these coefficients will also be non-real.

If there are complex roots, they come in conjugate pairs and so do the complex terms in the solution equation. If two of these complex terms are c_jλ^t
_j and c_j+1λ^t
_j+1, the roots λ_j can be written as

$${\displaystyle \lambda _{j},\lambda _{j+1}=\alpha \pm \beta i=M\left({\frac {\alpha }{M}}\pm {\frac {\beta }{M}}i\right)}$$

where i is the imaginary unit and M is the modulus of the roots:

$${\displaystyle M={\sqrt {\alpha ^{2}+\beta ^{2}}}.}$$

Then the two complex terms in the solution equation can be written as

$${\displaystyle {\begin{aligned}c_{j}\lambda _{j}^{t}+c_{j+1}\lambda _{j+1}^{t}&=M^{t}\left(c_{j}\left({\frac {\alpha }{M}}+{\frac {\beta }{M}}i\right)^{t}+c_{j+1}\left({\frac {\alpha }{M}}-{\frac {\beta }{M}}i\right)^{t}\right)\\[6pt]&=M^{t}\left(c_{j}\left(\cos \theta +i\sin \theta \right)^{t}+c_{j+1}\left(\cos \theta -i\sin \theta \right)^{t}\right)\\[6pt]&=M^{t}{\bigl (}c_{j}\left(\cos \theta t+i\sin \theta t\right)+c_{j+1}\left(\cos \theta t-i\sin \theta t\right){\bigr )}\end{aligned}}}$$

where θ is the angle whose cosine is α/M and whose sine is β/M; the last equality here made use of de Moivre's formula.

Now the process of finding the coefficients c_j and c_j+1 guarantees that they are also complex conjugates, which can be written as γ ± δi. Using this in the last equation gives this expression for the two complex terms in the solution equation:

$${\displaystyle 2M^{t}\left(\gamma \cos \theta t-\delta \sin \theta t\right)}$$

which can also be written as

$${\displaystyle 2{\sqrt {\gamma ^{2}+\delta ^{2}}}M^{t}\cos(\theta t+\psi )}$$

where ψ is the angle whose cosine is γ/√γ² + δ² and whose sine is δ/√γ² + δ².

Depending on the initial conditions, even with all roots real the iterates can experience a transitory tendency to go above and below the steady state value. But true cyclicity involves a permanent tendency to fluctuate, and this occurs if there is at least one pair of complex conjugate characteristic roots. This can be seen in the trigonometric form of their contribution to the solution equation, involving cos θt and sin θt.

In the second-order case, if the two roots are identical (λ₁ = λ₂), they can both be denoted as λ and a solution may be of the form

$${\displaystyle x_{t}=c_{1}\lambda ^{t}+c_{2}t\lambda ^{t}.}$$

An alternative solution method involves converting the nth order difference equation to a first-order matrix difference equation. This is accomplished by writing w_1,t = y_t, w_2,t = y_t−1 = w_1,t−1, w_3,t = y_t−2 = w_2,t−1, and so on. Then the original single nth-order equation

$${\displaystyle y_{t}=a_{1}y_{t-1}+a_{2}y_{t-2}+\cdots +a_{n}y_{t-n}+b}$$

can be replaced by the following n first-order equations:

$${\displaystyle {\begin{aligned}w_{1,t}&=a_{1}w_{1,t-1}+a_{2}w_{2,t-1}+\cdots +a_{n}w_{n,t-1}+b\\w_{2,t}&=w_{1,t-1}\\&\,\,\,\vdots \\w_{n,t}&=w_{n-1,t-1}.\end{aligned}}}$$

Defining the vector w_i as

$${\displaystyle \mathbf {w} _{i}={\begin{bmatrix}w_{1,i}\\w_{2,i}\\\vdots \\w_{n,i}\end{bmatrix}}}$$

this can be put in matrix form as

$${\displaystyle \mathbf {w} _{t}=\mathbf {A} \mathbf {w} _{t-1}+\mathbf {b} }$$

Here A is an n × n matrix in which the first row contains a₁, ..., a_n and all other rows have a single 1 with all other elements being 0, and b is a column vector with first element b and with the rest of its elements being 0.

This matrix equation can be solved using the methods in the article Matrix difference equation. In the homogeneous case y_i is a para-permanent of a lower triangular matrix [6]

The recurrence

$${\displaystyle y_{t}=a_{1}y_{t-1}+\cdots +a_{n}y_{t-n}+b,}$$

can be solved using the theory of generating functions. First, we write ${\textstyle Y(x)=\sum _{t\geq 0}y_{t}x^{t}}$. The recurrence is then equivalent to the following generating function equation:

$${\displaystyle Y(x)=a_{1}xY(x)+a_{2}x^{2}Y(x)+\cdots +a_{n}x^{n}Y(x)+{\frac {b}{1-x}}+p(x)}$$

where ${\displaystyle p(x)}$ is a polynomial of degree at most ${\displaystyle n-1}$ correcting the initial terms. From this equation we can solve to get

$${\displaystyle Y(x)=\left({\frac {b}{1-x}}+p(x)\right)\cdot {\frac {1}{1-a_{1}x-a_{2}x^{2}-\cdots -a_{n}x^{n}}}.}$$

In other words, not worrying about the exact coefficients, ${\displaystyle Y(x)}$ can be expressed as a rational function ${\displaystyle Y(x)={\frac {f(x)}{g(x)}}.}$

The closed form can then be derived via partial fraction decomposition. Specifically, if the generating function is written as

$${\displaystyle {\frac {f(x)}{g(x)}}=\sum _{i}{\frac {f_{i}(x)}{(x-r_{i})^{m_{i}}}}}$$

then the polynomial ${\displaystyle p(x)}$ determines the initial set of corrections ${\displaystyle z(n)}$, the denominator ${\displaystyle (x-r_{i})^{m}}$ determines the exponential term ${\displaystyle r_{i}^{n}}$, and the degree ${\displaystyle m}$ together with the numerator ${\displaystyle f_{i}(x)}$ determine the polynomial coefficient ${\displaystyle k_{i}(n)}$.

The method for solving linear differential equations is similar to the method above—the "intelligent guess" (ansatz) for linear differential equations with constant coefficients is ${\displaystyle e^{\lambda x}}$ where ${\displaystyle \lambda }$ is a complex number that is determined by substituting the guess into the differential equation.

This is not a coincidence. Considering the Taylor series of the solution to a linear differential equation:

$${\displaystyle \sum _{n=0}^{\infty }{\frac {f^{(n)}(a)}{n!}}(x-a)^{n}}$$

it can be seen that the coefficients of the series are given by the ${\displaystyle n}$-th derivative of ${\displaystyle f(x)}$ evaluated at the point ${\displaystyle a}$. The differential equation provides a linear difference equation relating these coefficients.

This equivalence can be used to quickly solve for the recurrence relationship for the coefficients in the power series solution of a linear differential equation.

The rule of thumb (for equations in which the polynomial multiplying the first term is non-zero at zero) is that:

$${\displaystyle y^{[k]}\to f[n+k]}$$

and more generally

$${\displaystyle x^{m}*y^{[k]}\to n(n-1)...(n-m+1)f[n+k-m]}$$

Example: The recurrence relationship for the Taylor series coefficients of the equation:

$${\displaystyle (x^{2}+3x-4)y^{[3]}-(3x+1)y^{[2]}+2y=0}$$

is given by

$${\displaystyle n(n-1)f[n+1]+3nf[n+2]-4f[n+3]-3nf[n+1]-f[n+2]+2f[n]=0}$$

or

$${\displaystyle -4f[n+3]+2nf[n+2]+n(n-4)f[n+1]+2f[n]=0.}$$

This example shows how problems generally solved using the power series solution method taught in normal differential equation classes can be solved in a much easier way.

Example: The differential equation

$${\displaystyle ay''+by'+cy=0}$$

has solution

$${\displaystyle y=e^{ax}.}$$

The conversion of the differential equation to a difference equation of the Taylor coefficients is

$${\displaystyle af[n+2]+bf[n+1]+cf[n]=0.}$$

It is easy to see that the ${\displaystyle n}$-th derivative of ${\displaystyle e^{ax}}$ evaluated at ${\displaystyle 0}$ is ${\displaystyle a^{n}}$.

Certain difference equations - in particular, linear constant coefficient difference equations - can be solved using z-transforms. The z-transforms are a class of integral transforms that lead to more convenient algebraic manipulations and more straightforward solutions. There are cases in which obtaining a direct solution would be all but impossible, yet solving the problem via a thoughtfully chosen integral transform is straightforward.