Lecture 36: Symbolic Computation with sympy¶
SymPy is pure Python package for symbolic computation. The pure means that, unlike in SageMath, no modifications to the Python scripting language have been made. SymPy is integrated in SageMath and at times we need to be aware of how to use its functionality properly. SageMath is licensed under the GNU GPL, whereas the license of SymPy is BSD.
sympy outside SageMath¶
Outside SageMath, we can compute with SymPy as in the Sage cell server, point your browser to the URL <http://live.sympy.org>.
A comparison of SymPy versus SageMath is given at <https://github.com/sympy/sympy/wiki/SymPy-vs.-Sage>.
SymPy is one of the core packages of the SciPy stack <http://www.scipy.org>.
sympy inside SageMath¶
We explicitly applied sympy to define a generator for the terms in a Taylor series.
If we use sympy in SageMath, then we may have to be explicit
about the use of its functions.
For example, the var
in SageMath is different
from the var
in sympy.
import sympy
x = sympy.var('x')
type(x)
We see that the type of x
is sympy.core.symbol.Symbol.
Observe the difference with the var
that we are normally using in SageMath:
y = var('y')
type(y)
and the type is
sage.symbolic.expression.Expression.
If we use the sqrt()
on a variable that is a sympy symbol,
then we get a symbolic expression of SageMath on return.
print(type(sqrt(x)))
print(type(sqrt(y)))
In both cases we see the same type.
We can make a generator object on the sqrt(1-x).
That works.
tsqrtx = sympy.series(sqrt(1-x),x,0,n=None)
print(tsqrtx)
print([next(tsqrtx) for k in range(10)])
But it will not work on sqrt(1-y).
tsqrty = sympy.series(sqrt(1-y),y,0,n=None)
print(tsqrty)
next(tsqrty)
next(tsqrty)
Although tsqrty
is also some type of generator,
the second application of the generator triggers
the StopIteration
exception.
So for the generating function for the Catalan numbers,
we must use a sympy symbol x
and not a SageMath variable y.
g = (1 - sqrt(1-4*x))/(2*x)
cg = sympy.series(g,x,0,n=None)
print([next(cg) for k in range(11)])
print(catalan_number(10))
Pattern Matching¶
Another feature better in SymPy than in SageMath is pattern matching, which can be useful when combined with substitution. To run this code we can also open a Terminal session, select Misc, and select python.
reset()
from sympy import *
x = Symbol('x')
y = Wild('y')
d = (10*x**3).match(y*x**3)
d
On return is a dictionary that has y_
as key
with 10
as corresponding value, which expresses
the match between y
and 10
in the expressions 10*x**3
and y*x**3
.
Solving Recurrence Relations¶
SymPy can solve recurrence relations.
Consider the Fibonacci numbers,
defined as \(f(n) = f(n-1) + f(n-2)\).
We import rsolve
and declare f
as a function of n
.
reset()
from sympy import Function, rsolve
from sympy.abc import n
f = Function('f')
The equation \(f(n) = f(n-1) + f(n-2)\)
is entered as f(n+2) - f(n+1) - f(n)
and assigned to the variable Fib
.
This variable is the first argument of rsolve
and we solve for f(n)
.
Fib = f(n+2) - f(n+1) - f(n)
print('solving', Fib)
s = rsolve(Fib, f(n))
s
The solution s
has the following form:
C0*(1/2 + sqrt(5)/2)**n + C1*(-sqrt(5)/2 + 1/2)**n
.
Because we have a two terms relation,
for a unique solution, we need to initial conditions,
as an additional dictionary argument.
s2 = rsolve(Fib, f(n), {f(0): 0, f(1): 1})
s2
and the solution s2
equals
sqrt(5)*(1/2 + sqrt(5)/2)**n/5 - sqrt(5)*(-sqrt(5)/2 + 1/2)**n/5
.
As a verification, let us apply the formula to compute the
first 10 Fibonacci numbers. But first, we need to convert the SymPy s2
into a SageMath expression:
fs2 = SR(s2)
Now we can evaluate and expand the expressions involving sqrt(5)
,
which are then automatically simplified.
[fs2(n=k).expand() for k in range(10)]
The Nearest Algebraic Number¶
SymPy depends only on mpmath
and mpmath
exports the pslq
function which allows to compute the nearest algebraic number,
given a floating-point approximation. Let us illustrate this
with \(\sqrt{2}\) which is a root of x^2 - 2
.
from mpmath import pslq
s2 = RR(sqrt(2))
pwrs = [s2**k for k in range(3)]
In pwrs
are the first three powers of s2
starting at the zeroth power.
cff = pslq(pwrs, verbose=True)
cff
The output in cff
is the list [2, 0, -1]
which are the
coefficients of the polynomial \(2 - x^2\).
This is the polynomial for which s2
is an approximate roots
with 53 bits of precision.
Assignments¶
What happens if you ask SageMath or sympy to compute a Taylor series expansion of
sqrt(x)
atx = 0
? Explain the result.The cost \(T(n)\) of a merge sort on a list of \(n\) numbers is governed by the following recurrence relation:
\[T(n) = 2 T\left(\frac{n}{2}\right) + n-1, \quad T(1) = 0.\]Use rsolve to find an explicit solution for \(T(n)\). Simplify so you can see \(T(n)\) is \(O(n \log(n))\).