# Lesson 32. Solving Recurrence Relations
# ----------------------------------------
#
# Consider a linear constant-coefficient homogeneous recurrence relation:
# y(n) =sum(c[k]*y(n-k), k = 1 .. N)
# As we've seen, this can be analyzed using a matrix formulation:
# X(n) = [y(n), y(n+1), ..., y(n+k-1)]
# X(n+1) = M X(n)
# e.g.
> recurrence:= y(n) = c[1]*y(n-1)+c[2]*y(n-2)+c[3]*y(n-3);
recurrence := y(n) = c[1] y(n - 1) + c[2] y(n - 2) + c[3] y(n - 3)
> with(linalg):
Warning: new definition for norm
Warning: new definition for trace
> M:= matrix([[0,1,0],[0,0,1],[c[3],c[2],c[1]]]);
[ 0 1 0 ]
[ ]
M := [ 0 0 1 ]
[ ]
[ c[3] c[2] c[1] ]
> evalm(M &* [y(n),y(n+1), y(n+2)]);
[ y(n + 1), y(n + 2), c[3] y(n) + c[2] y(n + 1) + c[1] y(n + 2) ]
# The eigenvalues of M are the roots of its characteristic polynomial.
> charpoly(M,t) = det(-M + t * &*());
3 2 3 2
t - t c[1] - t c[2] - c[3] = t - t c[1] - t c[2] - c[3]
# In general the characteristic polynomial of M is t^N - sum(c[j]*t^(N-j), j=1..N). Let's make up an example
# where the eigenvalues are -1, 1 and 2.
> poly:= expand((t+1)*(t-1)*(t-2));
3 2
poly := t - 2 t - t + 2
> cexample:= [2,1,-2];
cexample := [2, 1, -2]
> S:= {c[1]=2, c[2]=1, c[3]=-2};
S := {c[1] = 2, c[2] = 1, c[3] = -2}
> M:=map(t -> eval(subs(S,t)), M);
[ 0 1 0 ]
[ ]
M := [ 0 0 1 ]
[ ]
[ -2 1 2 ]
> eigenvects(M);
[2, 1, {[ 1, 2, 4 ]}], [1, 1, {[ 1, 1, 1 ]}], [-1, 1, {[ 1, -1, 1 ]}]
# The eigenvectors correspond to solutions of the form y(n) = t^n where t is the eigenvalue. Every solution will
# be of the form y(n) = sum( a[j] * t[j]^n, j = 1 .. N) where the a[j] solve the equations expressing the initial
# conditions as linear combinations of the eigenvectors. For example, take the initial conditions y(0)=1, y(1)=0,
# y(2)=0.
> eqns:= { 1 = a[1] + a[2] + a[3],
> 0 = -a[1] + a[2] + 2*a[3],
> 0 = a[1] + a[2] + 4*a[3] };
eqns :=
{1 = a[1] + a[2] + a[3], 0 = - a[1] + a[2] + 2 a[3], 0 = a[1] + a[2] + 4 a[3]}
> solve(eqns, {a[1], a[2], a[3]});
{a[2] = 1, a[3] = -1/3, a[1] = 1/3}
> soln:= n -> 1/3*(-1)^n +1 -1/3*2^n;
n n
soln := n -> 1/3 (-1) + 1 - 1/3 2
> seq(soln(nn), nn = 0 .. 10);
1, 0, 0, -2, -4, -10, -20, -42, -84, -170, -340
# The "dominant" eigenvalue here is 2. As n -> infinity,
# y(n)/2^n -> a[3] = -1/3.
#
# We can also try non-homogeneous recurrence relations, of the form
# y(n) = sum(c[k]*y(n-k), k = 1 .. N) + f(n).
# The basic principle here is that, at least if f(n) has a simple enough form, we can guess at the form of a
# solution. For example, if f(n) = r^n for some constant r, we might look for solutions of the form y(n) = b*r^n
# where b is constant. We find b by plugging this in to the recurrence relation.
> recurrence:= y(n) = 2*y(n-1) + y(n-2) - 2*y(n-3) + r^n;
n
recurrence := y(n) = 2 y(n - 1) + y(n - 2) - 2 y(n - 3) + r
> eval(subs(y = (n -> b*r^n), recurrence));
n (n - 1) (n - 2) (n - 3) n
b r = 2 b r + b r - 2 b r + r
> solve(",b);
n
r
---------------------------------------
n (n - 1) (n - 2) (n - 3)
r - 2 r - r + 2 r
> simplify(");
3
r
-----------------
3 2
r - 2 r - r + 2
# So it will work as long as the denominator is not 0, i.e. unless r itself is one of the eigenvalues. Just as in
# ordinary linear algebra, the general solution is this particular solution plus a solution of the homogeneous
# system. This solution of the homogeneous system can be chosen so that the initial conditions are satisfied.
# For example, with r = -2:
> b:= subs(r=-2,");
b := 2/3
> y:= n -> 2/3*(-2)^n + a[1]*(-1)^n + a[2]*1 + a[3]*2^n;
n n n
y := n -> 2/3 (-2) + a[1] (-1) + a[2] + a[3] 2
> eqns:= { y(0) = 1, y(1) = 0, y(2) = 0 };
eqns := {2/3 + a[1] + a[2] + a[3] = 1, - 4/3 - a[1] + a[2] + 2 a[3] = 0,
8/3 + a[1] + a[2] + 4 a[3] = 0}
> solve(eqns, {a[1],a[2],a[3]});
{a[3] = -1, a[2] = 7/3, a[1] = -1}
> y:= subs(", eval(y));
n n n
y := n -> 2/3 (-2) - (-1) + 7/3 - 2
> subs(r=-2, recurrence);
n n n (n - 1) (n - 1)
2/3 (-2) - (-1) + 7/3 - 2 = 4/3 (-2) - 2 (-1) + 7/3
(n - 1) (n - 2) (n - 2) (n - 2) (n - 3)
- 2 2 + 2/3 (-2) - (-1) - 2 - 4/3 (-2)
(n - 3) (n - 3) n
+ 2 (-1) + 2 2 + (-2)
> simplify(lhs(")-rhs("));
0
> seq(y(nn), nn=0 .. 10);
1, 0, 0, -10, -4, -50, -20, -210, -84, -850, -340
# Asymptotically, this y(n) is dominated by the terms 2/3 (-2)^n - 2^n: as n -> infinity,
# y(n)/2^n is approximately 2/3 (-1)^n - 1, alternating between -1/3 and -5/3.
#
# Now that you know (at least in some simple cases) how the solutions are found, I can show you Maple's own
# command for solving recurrence relations: "rsolve".
> y:= 'y': r:= -2: recurrence;
n
y(n) = 2 y(n - 1) + y(n - 2) - 2 y(n - 3) + (-2)
> rsolve(recurrence, y(n));
n
y(0) + 1/2 y(1) - 1/2 y(2) - (1/3 y(0) - 1/3 y(2)) 2
n n n
+ (1/3 y(0) - 1/2 y(1) + 1/6 y(2)) (-1) + 4/3 - 2/3 2 - 4/3 (-1)
n
+ 2/3 (-2)
> rsolve({recurrence, y(0)=1, y(1)=0, y(2)=0}, y(n));
n n n
2/3 (-2) - (-1) + 7/3 - 2
# "rsolve" can also solve some recurrences that have non-constant coefficients or even are nonlinear. Not every
# one, but some.
> rsolve({ y(n) = n * y(n-1), y(0)=1}, y(n));
GAMMA(n + 1)
> rsolve({ y(n) = n * y(n-1)+1, y(0)=1}, y(n));
exp(1) GAMMA(n + 2) - hypergeom([1], [n + 2], 1)
------------------------------------------------
n + 1
> rsolve({ y(n)*(1+y(n-1)) = y(n-1), y(0)=1}, y(n));
1
-----
n + 1
# Curiously, if we write this one as an equation for y(n) in terms of y(n-1), "rsolve" doesn't handle it.
> rsolve({y(n) = y(n-1)/(1+y(n-1)), y(0)=1}, y(n));
y(n - 1)
rsolve({y(0) = 1, y(n) = ------------}, y(n))
1 + y(n - 1)
# Even when "rsolve" isn't able to solve the recurrence relation, it may be able to give you asymptotic
# information.
> asympt(",n);
2 3
_C _C _C 1
1/n + ---- + --- + --- + O(----)
2 3 4 5
n n n n
# Stability of numerical calculations
# -----------------------------------
#
# Recurrence relations are often used in numerical methods. You want to compute some sequence y(n) which
# satisfies a recurrence relation, so you start with known values for y(0) or the first few y(n), and iterate the
# recurrence formula. This may or may not be a good way to calculate y(n). The main thing that can go wrong
# is that the inevitable roundoff errors can grow with each iteration until they overwhelm the true solution.
# For example, suppose you want to approximate the remainders in the series for E:
> w:= n -> n! * (E - Sum(1/k!, k = 0 .. n)):
> seq(evalf(w(nn)), nn = 0 .. 10);
1.718281828, .718281828, .436563656, .30969097, .23876388, .1938194, .162916,
.14041, .1233, .1095, .098
# Catastrophic cancellation is giving us fewer and fewer digits of accuracy. Can we get around that?
# It is easy to see that w(n) satisfies the following recurrence relation:
> recurrence:= y(n) = n*y(n-1) - 1:
> y(0):= evalf(E-1);
y(0) := 1.718281828
> for nn from 1 to 10 do y(nn):= nn*y(nn-1)-1 od:
> seq(y(nn), nn = 1 .. 10);
.718281828, .436563656, .309690968, .238763872, .193819360, .162916160,
.140413120, .123304960, .109744640, .097446400
# Looks good so far, but ...
> for nn from 11 to 20 do y(nn):= nn*y(nn-1)-1 od:
> seq(y(nn), nn = 11 .. 20);
.071910400, -.137075200, -2.781977600, -39.94768640, -600.2152960,
7 8
-9604.444736, -163276.5605, -.2938979089*10 , -.5584060369*10 ,
10
-.1116812075*10
>
# These results are absurd. The w(n) should all be positive, and in fact less than (n+2)/(n+1)^2. Roundoff error
# has done us in. It all started with the single roundoff error in y(0).
#
# What has happened is the following:
# The w(n) we want is only one solution of the recurrence relation. The general solution is this plus a solution of
# the homogeneous system y(n) = n y(n-1). Since y(n) = n! is a solution of the homogeneous system, the general
# solution is y(n) = w(n) + c n!, where c is an arbitrary constant. This constant is
# determined by y(0) = w(0) + c. If we don't have w(0) exactly right, c will not be 0. This means that, even if we
# had no roundoff error in all the remaining calculations, this c n! component will be present, and as n! grows
# rapidly its presence will eventually be overwhelming.
#
# So this was a bad way to calculate w(n).
# A better way is possible in this case: to use the recurrence relation in reverse.
> recurrence:= y(n-1) = (y(n)+1)/n:
# If we knew y(N) for some large N, we could use this to calculate y(n) for smaller n's.
> y(10):= 'y(10)';\
for kk from 1 to 10 do y(10-kk):= (y(11-kk)+1)/(11-kk) od;
y(10) := y(10)
y(9) := 1/10 y(10) + 1/10
11
y(8) := 1/90 y(10) + ----
90
101
y(7) := 1/720 y(10) + ---
720
821
y(6) := 1/5040 y(10) + ----
5040
5861
y(5) := 1/30240 y(10) + -----
30240
36101
y(4) := 1/151200 y(10) + ------
151200
187301
y(3) := 1/604800 y(10) + ------
604800
792101
y(2) := 1/1814400 y(10) + -------
1814400
2606501
y(1) := 1/3628800 y(10) + -------
3628800
6235301
y(0) := 1/3628800 y(10) + -------
3628800
# The dependence of y(n) on y(N) is less and less as N-n increases, so it doesn't really matter much what we take
# for y(N), as long as it is not extremely large.
# A similar problem occurs with the approximate solution of differential equations. To approximate the solution
# of a differential equation y'(x) = f(x, y(x)) with initial value y(0) = y0, we use a recurrence relation. The
# simplest is the Euler method:
# y(x+h) = y(x) + h*f(x,y(x))
# For example, try the equation y' = -y with y(0) = 1. The actual solution is y(x) = exp(-x). Try
# it with h = 0.2.
> y:= 'y': y(0):= 1: h:= 0.2:\
for kk from 1 to 20 do y(kk * h):= y((kk-1)*h) - h * y((kk-1)*h) od:
> with(plots,display):\
display({plot( exp(-x), x = 0 .. 4),\
plot([seq([kk*h, y(kk*h)], kk = 0 .. 20)], style=POINT)});
** Maple V Graphics **
# So far so good: the approximate solution is not quite the same as the true solution, but decreasing h should
# improve things a bit. At least both have the correct behaviour as x -> infinity.
# But what happens for the equation y' = -20 y ? I'll just do a few of the y values.
> for kk from 1 to 5 do y(kk * h):= y((kk-1)*h) - 20*h * y((kk-1)*h) od;
y(.2) := -3.0
y(.4) := 9.00
y(.6) := -27.000
y(.8) := 81.0000
y(1.0) := -243.00000
# The true solution is y(x) = exp(-20 x), which goes to 0 (even more rapidly) as x -> infinity. The numerical
# method, however, has become unstable. In fact, the recurrence can be written as
# y(k h) = (1-20 h) y((k-1) h)
# which has solutions y(n h) = y(0) (1 - 20 h)^n. When 20 h >= 1, these solutions no longer look like the
# solutions of the differential equation, which decrease (or increase) smoothly to 0. They oscillate above
# and below 0, and if 20 h > 2 they increase in absolute value.
# This is a phenomenon known as "stiffness": when a solution of a differential equation strongly attracts other
# solutions, the step size must be chosen very small in order for the numerical method to be stable.
#
