%md <h1>An Introduction to Recursion using SageMath</h1>

%md <p>Here are a few tips on how to get started using Sage to work with recursion.



# An Introduction to Recursion using SageMath

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

Consider the recurrence relation $a_{n+2}=-a_{n+1}+2 a_n$ with the initial conditions $a_0=1$, $a_1=5$.

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")')

"/projects/sage/sage-7.5/local/share/maxima/5.35.1/share/solve_rec/solve_rec.mac"

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.

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


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 advice 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)





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)) 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')




For a recursive version of the binary search algorithm, see our introduction the SageMath implementation of Binary Tree Sorting.

Derangements

m=4/3

m.n(digits=8)

1.3333333