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