Sharedsage_worksheets / ADS_Recursion.sagewsOpen in CoCalc
Author: Ken Levasseur
Description: Worksheets related to Applied Discrete Structures

An Introduction to Sets and Combinatorics using Sage

Applied Discrete Structures by Alan Doerr & Kenneth Levasseur is licensed under a Creative Commons Attribution - Noncommercial - No Derivative Works 3.0 United States License.

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/3a(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 aa.


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 a0=4a_{0}=4 and for positive nn, an=3an1a_{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