1.32 Lesson 32

```# 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.
#

```