Author: Ken Levasseur
Description: Worksheets related to Applied Discrete Structures

# An Introduction to Sets and Combinatorics using Sage

Here are a few tips on how to get started using Sage to work with recursion.

## Solution of Linear Recurrence Relations

Sage cannot currently solve this system, but maxima does. So we first load maxima’s solve_rec function. If this changes, I’d appreciate it if anyone would let me know!

maxima('load("solve_rec")')


Sage’s maxima function sends the argument provided as a string to maxima and returns the maxima output. Here is a second order linear recurrence relation with constant coefficients.

a[n]=(-2)^(n+1)/3+5/3

This tells us that the solution is  $a(n) = (1/3) (-2)^{n+1} + 5/3$, but to compute with this result we need to extract the right side and use the expression to define a function called $a$.

a(n)=sol.rhs()
a

n |--> 1/3*(-2)^(n + 1) + 5/3

Note:  At one point, I got advise to use the somewhat more complicated expression

a(n)=  maxima(sol.rhs()).sage()

which also works.  It may be possible that certain solutions where the maxima output has a form that is incompatible with Sage and the more complicated expression is needed.

Here are the first few terms of the sequence and a plot.

map(a,range(10))

[1, 3, -1, 7, -9, 23, -41, 87, -169, 343]
list_plot(map(lambda k:[k,a(k)],range(10)))


## A recursive definition

Here is how to define a sequence recursively. For example suppose $a_{0}=4$ and for positive $n$, $a_{n}=3 a_{n-1}$.

var('a')

a
def a(n):
if n==0:
return 4
else:
return 3*a(n-1)


What follows are some calculations and a plot of the sequence.

a(5)

972
data=[a(n) for n in range(20)]

data


[4, 12, 36, 108, 324, 972, 2916, 8748, 26244, 78732, 236196, 708588, 2125764, 6377292, 19131876, 57395628, 172186884, 516560652, 1549681956, 4649045868]
var('pts')

pts
pts=[(n,a(n)/n) for n in range(10)]

pts

[(0, 4), (1, 12), (2, 36), (3, 108), (4, 324), (5, 972), (6, 2916), (7, 8748), (8, 26244), (9, 78732)]
list_plot(pts,plotjoined=True, color='purple')


## Non-constant coefficients

maxima('solve_rec(a[n+2]+a[n+1]-2*(n+4)*a[n]=0,a[n],a[0]=1,a[1]=3)')

false