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

  1. What happens if you ask SageMath or sympy to compute a Taylor series expansion of sqrt(x) at x = 0? Explain the result.

  2. 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))\).