Sharedsage_worksheets / ADS_recursion.sagewsOpen in CoCalc
Author: Ken Levasseur
Views : 27
Description: Worksheets related to Applied Discrete Structures
%md <h1>An Introduction to Recursion using SageMath</h1> %md <b><a href="http://faculty.uml.edu/klevasseur/ads2/">Applied Discrete Structures</a></b> by Alan Doerr & Kenneth Levasseur is licensed under a Creative Commons Attribution - Noncommercial - No Derivative Works 3.0 United States License. %md <p>Here are a few tips on how to get started using Sage to work with recursion.

An Introduction to Recursion using SageMath

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.

Consider the recurrence relation an+2=an+1+2an a_{n+2}=-a_{n+1}+2 a_n with the initial conditions a0=1a_0=1, a1=5a_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/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 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 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)
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