APM_2F005 - Algorithms For Discrete Mathematics Bachelor 2
Lecturer: Lucas Gerin

Symbolic computation 2: Linear recurrences

Table of contents

• Warm-up
• Solving linear recurrences with Sympy
  • The function rsolve
  • Application: Fibonacci
• Two-dimensional recurrence: matrices with SymPy
• Rational fractions

## loading python libraries

# necessary to display plots inline:
%matplotlib inline   

# load the libraries
import matplotlib.pyplot as plt # 2D plotting library
import numpy as np              # package for scientific computing  
from pylab import *

from math import *              # package for mathematics (pi, arctan, sqrt, factorial ...)
import sympy as sympy             # package for symbolic computation
from sympy import *

Warm-up

The aim of this Lab session is to use SymPy to solve automatically some simple recurrences. We first solve a particular case "by hand" and then use the SymPy function rsolve.

Exercise 1. Solving a recurrence (almost) by hand

We want to find an explicit formula for the sequence $(u_n)$ defined by \begin{equation} \begin{cases} u_0&=1,\\ u_{n}&=2u_{n-1} +3n^2. \qquad (\forall n\geq 1). \end{cases} \tag{$\star$} \end{equation}

1. . Let $\alpha,a,b,c$ be four real parameters. Let $(w_n)_{n\geq 0}$ be the sequence defined by $$ w_n=\alpha 2^n +an^2+bn+c. $$ Use `SymPy` to find $\alpha,a,b,c$ such that for every $0\leq n\leq 3$, $$ w_n=u_n. $$ To solve a system of equations with `SymPy` with unknowns $x,y$, write something like ``` solve([x-y-2,3*y+x],[x,y]) ```
2. . Prove with `SymPy` that $(w_n)$ defined as above is the unique solution of ($\star$). (It might be helpful to define a function `w(n,alpha,a,b,c)` which computes $w_n$.)

# ------------- Question 1 --------------
var('n',integer=True)
var('alpha a b c')
def w(n,alpha,a,b,c):
    return alpha*2**n +a*n**2 +b*n +c

def u(n):
    if n==0:
        return 1
    else:
        return 2*u(n-1)+3*n**2

FourEquations=[w(n,alpha,a,b,c)-u(n) for n in range(4)]

Solutions=solve(
# Quantities which are to be zero:
FourEquations,
# Unknowns:
[alpha,a,b,c]
)
print('--------------')
print('Question 1:', Solutions)

--------------
Question 1: {alpha: 19, a: -3, b: -12, c: -18}

1. We have four unknowns $\alpha,a,b,c$ and four equations \begin{align*} u_0&=1,\\ u_1&=2u_0+3\times 1^2,\\ u_2&=2u_1+3\times 2^2,\\ u_3&=2u_2+3\times 3^2,\\ \end{align*} With the above script we obtain that if the formula is true then necessarily $$ u_n=19\times 2^n -3n^2-12n-18. $$

#-------------- Question 2 ----------------
print('--------------')
c_s=Solutions[c]
alpha_s=Solutions[alpha]
b_s=Solutions[b]
a_s=Solutions[a]
Expr=w(n,alpha_s,a_s,b_s,c_s)-2*w(n-1,alpha_s,a_s,b_s,c_s)-3*n**2
print('Question 2:')
print('With these coefficients the simplificaion of w(n)-2w(n-1)-3n**2 gives for every n : '+str(simplify(Expr)))

#--------------- We check the solutions
print('--------------')
print('Safety check:')
print('First values of w_n: '+str([w(i,alpha_s,a_s,b_s,c_s) for i in range(10)]))

# In order to check our results:
def Recurrence(n):
    if n==0:
        return 1
    else:
        return 2*Recurrence(n-1)+3*n**2

print('First values (with recurrence formula) : '+str([Recurrence(n) for n in range(10)]))

--------------
Question 2:
With these coefficients the simplificaion of w(n)-2w(n-1)-3n**2 gives for every n : 0
--------------
Safety check:
First values of w_n: [1, 5, 22, 71, 190, 455, 1018, 2183, 4558, 9359]
First values (with recurrence formula) : [1, 5, 22, 71, 190, 455, 1018, 2183, 4558, 9359]

2. We first observe that eq.$(\star)$ has a unique solution since the sequence $(u_n)$ is a recurrence of order one: it is uniquely determined by $u_0$. Besides, set
   $$ w_n=19\times 2^n -3n^2-12n-18. $$
   Parameters are chosen so that $w_0=u_0=1$. The above script shows that for every $n$ we have that $w_n-2w_{n-1}-3n^2=0$. Therefore $(w_n)$ satisfies the recurrence $(\star)$ and starts with the same initial value $w_0=1$. Therefore $w_n=u_n$ for every $n$.

Solving recurrences with SymPy: rsolve

The function rsolve

We will use Sympy to obtain explicit formulas for some sequences defined by linear recurrences. More precisely, we will see how to obtain an explicit formula for $u_n$ in two cases:

1. Linear recurrence of order one: this is a sequence $(u_n)_{n\geq 0}$ is defined by $$ \begin{cases} u_0&=a, \\ u_{n}&=\alpha u_{n-1}+f(n),\qquad (n\geq 1), \end{cases} $$ where $a,\alpha$ are some given constants and $f$ is an arbitrary function.
2. Linear recurrence of order two: this is a sequence $(u_n)_{n\geq 0}$ is defined by $$ \begin{cases} u_0&=a, \\ u_1&=b,\\ u_{n}&=\alpha u_{n-1}+\beta u_{n-2}+f(n),\qquad (n\geq 2), \end{cases} $$ where $a,,b,\alpha,\beta$ are some given constants and $f$ is an arbitrary function.

Some known examples fit in this settings:

1. Geometric sequences: $u_0=a,u_n=ru_{n-1}$.
2. Arithmetic sequences: $u_0=a,u_n=u_{n-1}+r$.
3. The Fibonacci sequence: $F_1=1,F_2=1,F_{n}=F_{n-1}+F_{n-2}$.

The following script shows how to solve the two recurrences \begin{align*} u_0=5,\qquad u_{n}&=3u_{n-1},\\ v_0=1,\qquad v_{n}&=2v_{n-1}+n,\\ \end{align*} using function rsolve.

# First example: a geometric sequence
u = Function('u')                 # declares the name of the sequence
n = symbols('n',integer=True)     # declares the variable
f = u(n)-3*u(n-1)                 # the expression which is to be zero
ExplicitFormula1=rsolve(f,u(n),
                        {u(0):5}) # initial condition

print('The formula for u(n) is '+str(ExplicitFormula1)+'')

# Second example with a non-linear term
print('-------')
v = Function('v')                 # declares the name of the sequence
n = symbols('n',integer=True)     # declares the variable
f = v(n)-2*v(n-1)-n           # the expression which is to be zero
ExplicitFormula2=rsolve(f,v(n),
                        {v(0):1}) # initial condition

print('The formula for v(n) is '+str(ExplicitFormula2)+'')

# Third example 
print('-------')
w = Function('w')                 # declares the name of the sequence
n = symbols('n',integer=True)     # declares the variable
g = w(n)-w(n-2)-8*n-2                 # the expression which is to be zero
ExplicitFormula3=simplify(rsolve(g,w(n),{w(0):2,w(1):10}))

print('The formula for w(n) is '+str(ExplicitFormula3)+'')

The formula for u(n) is 5*3**n
-------
The formula for v(n) is 2**n/2 + n/2 + 1/2
-------
The formula for w(n) is -(-1)**n/2 + 2*n**2 + 5*n + 5/2

Exercise 2. A first example with rsolve

1. Use `SymPy` to solve the recurrence of the Fibonacci sequence and find an explicit formula for $F_n$.
   (To set up two initial conditions you must write `{u(1):1,u(2):1}`.)
2. (Theory) Use the formula obtained at Question 1 to prove that $$ \lim_{n\to +\infty} \frac{F_n}{F_{n-1}}= \varphi, $$ where $\varphi=\frac{1+\sqrt{5}}{2}$ is the golden ratio.

u = Function('u')                 # declares the name of the sequence
n = symbols('n',integer=True)     # declares the variable
f = u(n)-u(n-1)-u(n-2)                 # the expression which is to be zero
ExplicitFormula=simplify(rsolve(f,u(n),{u(1):1,u(2):1}))

Formula2=latex(simplify(ExplicitFormula))
print('The formula for u(n) is '+str(Formula2)+'')

The formula for u(n) is \frac{2^{- n} \sqrt{5} \left(- \left(1 - \sqrt{5}\right)^{n} + \left(1 + \sqrt{5}\right)^{n}\right)}{5}

1. We put the above answer in $\texttt{LateX}$: $$ F_n=\frac{2^{- n}}{5} \sqrt{5} \left(\left(1 + \sqrt{5}\right)^{n} - \left(- \sqrt{5} + 1\right)^{n}\right). $$
2. We obtain (by setting $\psi=(1-\sqrt{5})/2$): \begin{align*} \frac{F_{n}}{F_{n-1}}&= \frac{\tfrac{1}{\sqrt{5}}\varphi^{n} - \tfrac{1}{\sqrt{5}}\psi^n }{\tfrac{1}{\sqrt{5}}\varphi^{n-1} - \tfrac{1}{\sqrt{5}}\psi^{n-1} },\\ &= \frac{\tfrac{1}{\sqrt{5}}\varphi - \psi \tfrac{1}{\sqrt{5}}(\psi/\varphi)^{n-1} }{\tfrac{1}{\sqrt{5}} - \tfrac{1}{\sqrt{5}}(\psi/\varphi)^{n-1} },\qquad \text{(we divide everything by $\varphi^{n-1}$.)}\\ &\to \varphi, \end{align*} since $|\psi/\varphi|<1$.

The output of `rsolve` is an expression which depends on the symbolic variable `n`. If we want to evaluate this expression (for instance for $n=10$) we must write:

Value=ExplicitFormula.subs(n,10)
print('10th Fibonacci number = '+str(Value))
print('After simplification : '+str(simplify(Value)))

10th Fibonacci number = sqrt(5)*(-(1 - sqrt(5))**10 + (1 + sqrt(5))**10)/5120
After simplification : 55

Two-dimensional recurrence: matrices with SymPy

Exercise 3. A double linear recurrence

Theory

1. Prove by induction that there exist integers $a_n,b_n$ such that for every $n\geq 1$ $$ (1+\sqrt{2})^n=a_n+b_n\sqrt{2}. $$
2. Find a $2\times 2$ matrix $A$ such that $$ \binom{a_{n+1}}{b_{n+1}}=A\times \binom{a_{n}}{b_{n}}. $$
3. Use `SymPy` to find an explicit formula for $a_n$.
   (In `SymPy` matrices are defined by `A=Matrix([[a,b],[c,d]])`. To raise $A$ to some power just write `A**n`.)
4. Find (with pen/paper) real numbers $c,R$ such that $$ a_n \stackrel{n\to +\infty}{\sim } cR^n. $$

1. For $n=1$ this is true with $a_1=b_1=1$. Assume this holds for some $n\geq 1$. \begin{align*} (1+\sqrt{2})^{n+1}&=(1+\sqrt{2})\times(1+\sqrt{2})^{n}\\ &=(1+\sqrt{2})\times(a_n+b_n\sqrt{2})\\ &=a_n+b_n\sqrt{2}+a_n\sqrt{2}+2b_n\\ &=(a_n+2b_n)+(a_n+b_n)\sqrt{2}. \end{align*} Finally, one can write $$ (1+\sqrt{2})^{n+1}=a_{n+1}+b_{n+1}\sqrt{2}, $$ where $$ a_{n+1}=a_n+2b_n,\qquad b_{n+1}=a_n+b_n, $$ which are easily proved to be integers by induction.
2. By the above computation we have $$ \begin{cases} a_{n+1}&=a_n+2b_n,\\ b_{n+1}&=a_n+b_n, \end{cases} $$ which can be written as: $$ \binom{a_{n+1}}{b_{n+1}}= \begin{pmatrix} 1 & 2\\ 1 & 1 \end{pmatrix} \times \binom{a_{n}}{b_{n}}. $$ By induction this implies that $$ \binom{a_{n}}{b_{n}}= \begin{pmatrix} 1 & 2\\ 1 & 1 \end{pmatrix}^{n-1} \times \binom{a_1}{b_1} =\begin{pmatrix} 1 & 2\\ 1 & 1 \end{pmatrix}^{n-1} \times \binom{1}{1}. $$

#Question 3
A=Matrix([[1,2],[1,1]])
var('n', integer=True)
Power=A**(n-1)
#print(Power)

an=latex(simplify(Power[0,0]+Power[0,1]))
print(an)
bn=latex(simplify(Power[1,0]+Power[1,1]))
print(bn)

\frac{\left(1 - \sqrt{2}\right)^{n}}{2} + \frac{\left(1 + \sqrt{2}\right)^{n}}{2}
\frac{\sqrt{2} \left(- \left(1 - \sqrt{2}\right)^{n} + \left(1 + \sqrt{2}\right)^{n}\right)}{4}

3. We export the result in LateX: \begin{align*} a_n&=\frac{\left(1 - \sqrt{2}\right)^{n}}{2} + \frac{\left(1 + \sqrt{2}\right)^{n}}{2}\\ b_n&= \frac{\sqrt{2} \left(- \left(1 - \sqrt{2}\right)^{n} + \left(1 + \sqrt{2}\right)^{n}\right)}{4} \end{align*} (You can indeed check that $a_n+\sqrt{2}b_n=(1+\sqrt{2})^n$.)
4. As $\left|-\sqrt{2}+1\right|<1$, we have that $(-\sqrt{2}+1)^n$ tends to zero. It follows that \begin{align*} a_n&=\frac{1}{2} \left(1 + \sqrt{2}\right)^{n} +\mathrm{o}(1)\\ &\sim \frac{1}{2} \left(1 + \sqrt{2}\right)^{n}. \end{align*}

Exercise 4. Rational fractions. (Taken from BX2023's test)

We set $$ u_1=\frac{1}{1+1}, \quad u_2=\frac{1}{1+\frac{1}{2+1}},\quad u_3=\frac{1}{1+\frac{1}{2+\frac{1}{3+1}}}, \quad u_4=\frac{1}{1+\frac{1}{2+\frac{1}{3+\frac{1}{4+1}}}}, \dots $$ Using SymPy, write $u_{20}$ as a rational fraction $a/b$. (To check your result: the approximate value is $0.69777...$.)

# Tests
u1=Rational(1,2)
u2=Rational(1,1+Rational(1,3))
u3=Rational(1,1+Rational(1,2+Rational(1,4)))

def u_n(n):
    u=n+1
    for k in range(1,n):
        u=n-k+Rational(1,u)
    return Rational(1,u)

for n in range(1,21):
    print(n,'---',u_n(n))
    
    

1 --- 1/2
2 --- 3/4
3 --- 9/13
4 --- 37/53
5 --- 187/268
6 --- 1129/1618
7 --- 7933/11369
8 --- 63621/91177
9 --- 573561/821986
10 --- 5742571/8229836
11 --- 63224941/90609397
12 --- 759216193/1088053549
13 --- 9875036179/14152185188
14 --- 138308505777/198213712978
15 --- 2075328803577/2974210627873
16 --- 33214434676489/47600517297953
17 --- 564774524186833/809393860526194
18 --- 10167887629480051/14571878633638372
19 --- 193221133200680401/276910505412260141
20 --- 3864956170297235421/5538974690732597941

We set \begin{align*} v_0&=21\\ v_1&=20+\frac{1}{v_0}=20+\frac{1}{20+1}\\ v_2&=19+\frac{1}{v_1}=19+\frac{1}{20+\frac{1}{20+1}}\\ & \dots\\ v_{20}&=1+\frac{1}{v_{19}}=1+\frac{1}{2+\frac{1}{\dots+\frac{1}{{\dots+\frac{1}{20+1}}}}} \end{align*} And we finally put $u_{20}=\frac{1}{v_{20}}$. Therefore the above script computes: $$ u_{20}=\frac{3864956170297235421}{5538974690732597941}. $$

Exercise 5. Multivariate Polynomials (Taken from BX2024's Test)

Use SymPy to find the coefficient of $x^3y^2z^7$ in the polynomial $(x+2y+z)^{12}$.

var('x y z')

Expression=(x+2*y+z)**12
display(expand(Expression))

$\displaystyle x^{12} + 24 x^{11} y + 12 x^{11} z + 264 x^{10} y^{2} + 264 x^{10} y z + 66 x^{10} z^{2} + 1760 x^{9} y^{3} + 2640 x^{9} y^{2} z + 1320 x^{9} y z^{2} + 220 x^{9} z^{3} + 7920 x^{8} y^{4} + 15840 x^{8} y^{3} z + 11880 x^{8} y^{2} z^{2} + 3960 x^{8} y z^{3} + 495 x^{8} z^{4} + 25344 x^{7} y^{5} + 63360 x^{7} y^{4} z + 63360 x^{7} y^{3} z^{2} + 31680 x^{7} y^{2} z^{3} + 7920 x^{7} y z^{4} + 792 x^{7} z^{5} + 59136 x^{6} y^{6} + 177408 x^{6} y^{5} z + 221760 x^{6} y^{4} z^{2} + 147840 x^{6} y^{3} z^{3} + 55440 x^{6} y^{2} z^{4} + 11088 x^{6} y z^{5} + 924 x^{6} z^{6} + 101376 x^{5} y^{7} + 354816 x^{5} y^{6} z + 532224 x^{5} y^{5} z^{2} + 443520 x^{5} y^{4} z^{3} + 221760 x^{5} y^{3} z^{4} + 66528 x^{5} y^{2} z^{5} + 11088 x^{5} y z^{6} + 792 x^{5} z^{7} + 126720 x^{4} y^{8} + 506880 x^{4} y^{7} z + 887040 x^{4} y^{6} z^{2} + 887040 x^{4} y^{5} z^{3} + 554400 x^{4} y^{4} z^{4} + 221760 x^{4} y^{3} z^{5} + 55440 x^{4} y^{2} z^{6} + 7920 x^{4} y z^{7} + 495 x^{4} z^{8} + 112640 x^{3} y^{9} + 506880 x^{3} y^{8} z + 1013760 x^{3} y^{7} z^{2} + 1182720 x^{3} y^{6} z^{3} + 887040 x^{3} y^{5} z^{4} + 443520 x^{3} y^{4} z^{5} + 147840 x^{3} y^{3} z^{6} + 31680 x^{3} y^{2} z^{7} + 3960 x^{3} y z^{8} + 220 x^{3} z^{9} + 67584 x^{2} y^{10} + 337920 x^{2} y^{9} z + 760320 x^{2} y^{8} z^{2} + 1013760 x^{2} y^{7} z^{3} + 887040 x^{2} y^{6} z^{4} + 532224 x^{2} y^{5} z^{5} + 221760 x^{2} y^{4} z^{6} + 63360 x^{2} y^{3} z^{7} + 11880 x^{2} y^{2} z^{8} + 1320 x^{2} y z^{9} + 66 x^{2} z^{10} + 24576 x y^{11} + 135168 x y^{10} z + 337920 x y^{9} z^{2} + 506880 x y^{8} z^{3} + 506880 x y^{7} z^{4} + 354816 x y^{6} z^{5} + 177408 x y^{5} z^{6} + 63360 x y^{4} z^{7} + 15840 x y^{3} z^{8} + 2640 x y^{2} z^{9} + 264 x y z^{10} + 12 x z^{11} + 4096 y^{12} + 24576 y^{11} z + 67584 y^{10} z^{2} + 112640 y^{9} z^{3} + 126720 y^{8} z^{4} + 101376 y^{7} z^{5} + 59136 y^{6} z^{6} + 25344 y^{5} z^{7} + 7920 y^{4} z^{8} + 1760 y^{3} z^{9} + 264 y^{2} z^{10} + 24 y z^{11} + z^{12}$

Sympy finds $$ (x+2y+z)^{12}= x^{12}+ \dots + 31680x^3y^2z^7+\dots +z^{12}. $$