CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutSign UpSign In

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place. Commercial Alternative to JupyterHub.

| Download

Recursion with SageMath

Views: 1755
Tweet
Share
Share
%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=3an−1a_{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 



[]
Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_parsing.py", line 538, in introspect O = eval(obj if not preparse else preparse_code(obj), File "<string>", line 1, in <module> NameError: name 'rsolve' is not defined