SymPy: Symbolic Python

Video

Motivation

Mathematicians tend to work with symbols more than data. While we are often manipulating them on paper, computers can manipulate them for us in many situations.

In addition, arbitrary position integers, polynomials, more general group elements like modular integers or elliptic curve points, and other similar arithmetic constructs can often be manipulated by computers for us with a bit of symbolic magic.

There are two different Python symbolic computation packages with two starkly different philosophies and separate sets of higher level capabilites. For the simple stuff, they are fairly interchangeable - but depending on your domain, Sage may be the right choice for you.

Sympy tends to be easier to use and learn, however, so we will be teaching it.

Core Object(s)

Expressions

Sympy's core object is the expression. Like in Numpy, they are typically built rather than passed to an explicit constructor.

The simplest kind of expression is the symbol. Sympy has a quick interface to symbols for upper and lowercase roman and greek letters:

import sympy
from sympy.abc import x
example_poly = x**2-1
example_poly

Symbols can also be constructed explicitly, if you need longer ones or custom renders:

x1,x2 = sympy.symbols("x_1 x_2")
x1

From symbols, together with the arithmetic operators and functions like sympy.sin, it is possible to construct complicated expressions:

expr = 1+sympy.sin(x)**2/sympy.cos(x)**2
expr

One neat thing you can do with Sympy is simplify expressions:

sympy.simplify(expr)

Note that Sympy can automatically format pretty-printed output for us!

We can get that output as $\LaTeX$:

print(sympy.latex(expr))

This can be a real boon for accurately constructing and representing objects in your research.

Evaluation

SymPy expressions, even when they contain no variables, are symbolic:

sympy.sqrt(2)

This is really useful, but sometimes we want numerical values (like for plotting!)

Expressions support the .n() method for numerical approximations:

a = sympy.sqrt(2)
a.n()

If you want to evaluate at variables, you can substitute expressions for specific values:

from sympy.abc import x,y
expr = (sympy.sin(x)+sympy.cos(y))
expr.subs({x:x**2})
expr.subs({x:1,y:2}).n()

Parsing

Instead of building expressions from their base components - a very reliable method - you can also parse input.

sympy.parse_expr("sin(x**2+y)")

This is commonly used in programs for taking in user data, and can be a quick way to parse input.

The standard form is the "math string" we are by now largely comfortable with, but it can also (sometimes) parse other formats, like $\LaTeX$. (there are philosophical differences that make such parsing... difficult.)

Matrices

Different Sympy domains revolve around different constructs; for example, the Linear Algebra domain revolves around the sympy.Matrix object:

sympy.Matrix([[1,1],[1,0]])

With these, we can finally realize the absurdly fast Fibonacci algorithm we glimpsed with Numpy for arbitrary precision integers:

def fibonacci(n):
    if n<1:
        return n
    iterator =sympy.Matrix([[1,1],[1,0]])
    base_case = sympy.Matrix([1,1])
    return (iterator**(n-1)*base_case)[0]
fibonacci(10**5)

Sympy Matrixes are not like ndarrays; they respond to all our functions and operators as a mathematician would expect a Matrix to; Because they contain Python objects, they can't take advantage of the same parallel computations as Numpy, so their speed relies on the work of linear algebraists, number theorists, and computer scientists - together with the inherent power of the matrix.

Because their notation is more familiar, there is a bit less to go over - but there are some special constructors worth talking about. Note that Sympy - like mathematicians - is mostly concerned with the creation and manipulation of square matrices.

Block construction ends up being a bit finnicky in SymPy:

A = sympy.zeros(3)
B = sympy.ones(3,2)
C = 2*sympy.ones(2,3)
D = sympy.eye(2)
M_blocks = sympy.BlockMatrix([[A,B],[C,D]])
M = sympy.Matrix(M_blocks)
N_blocks = sympy.BlockDiagMatrix(A,B,C,D)
N = sympy.Matrix(N_blocks)

Modules

SymPy functionality is largely split into various Modules - that is, python submodules - which you can read about in the documentation. As of writing, it is really good at:

• polynomials
• trigonometry
• integral and differential calculus
• Quadratic extensions of $\mathbb R$ and $\mathbb Q$, such as the complex numbers (and their extensions)
• Real and complex symbolic linear algebra
• Parsing and simplifying (reasonably simple) expressions

And experimental functionality in a variety of other areas.

If you need symbolic computations of other sorts, I highly encourage you to explore Sage, a truly incredible feat of software stitching. Sage attempts to perform interfaces and interoperability to almost everything mathematicians use - allowing you to take objects out of Macaulay and into Matlab or Mathematica seamlessly. It could be the subject of its own course - and tends to be relatively complicated to install, operate, and configure.

Polynomials

Polynomial coefficients have a way of getting away from you. When I need polynomial expansions, I tend to do some quick sympy expr.expand() to get the computer to do it for me.

However, we all know that polynomials have a lot more going on.

Consider the following polynomial, which we know has a algebraic solution:

from sympy.abc import x
polynomial = x**3+3*x**2+6*x+4
roots = sympy.solve(polynomial)
sympy.factor(polynomial)
roots

We can set some reasonable assumptions on our values, to get restricted answers:

a = sympy.symbols("a",real=True)
polynomial = a**3+3*a**2+6*a+4
sympy.solve(polynomial)
b = sympy.symbols("b",positive=True)
polynomial = b**3+3*b**2+6*b+4
sympy.solve(polynomial)

Sympy can solve in terms of various variables:

expr = y**2-x**3-x-1
sympy.solve(expr,y)
sympy.solve(expr,x)

And even do polynomial long division and gcd:

expr1 = x**3-2*x-x+2
expr2 = x**2+4*x+4
sympy.gcd(expr1,expr2)
sympy.factor(expr1)

sympy.poly(expr1).div(sympy.poly(expr2))

Note that several of the fancier polynomial tricks are hiding various places, like in the sympy.poly object and sympy.polys submodule.

Trigonometry

Sympy has some fairly efficient tools for symbolic trigonometry:

expr = sympy.sin(x+sympy.pi/3)+sympy.sin(x)
sympy.simplify(expr)

str(expr)

It can solve these equations, even as part of many polynomials:

from sympy.abc import theta
expr = sympy.sin(theta)+sympy.cos(theta)**2+1
solutions = sympy.solve(expr)
expr.subs(theta, solutions[0])

However, it might need to be told to think little deeper when substituting:

evaluated = expr.subs(theta, solutions[1])
sympy.simplify(evaluated)

It favors simplifying into as few $\sin$ and $\cos$ terms as it can, and changing other trig functions into those:

expr = sympy.tan(x)**2+sympy.sec(x)**2
sympy.simplify(expr)

Integral and Differential Calculus

SymPy is very good at integral and differential calculus, though it doesn't know the scientific names of beings animalculous.

With symbolic solvers for differentiation, integration, and many ODEs and PDEs in several variables, it is a great resource - just make sure you sanity-check the answers you get.

To take the derivative of a simple expression, do:

from sympy.abc import x
expr = x**3-sympy.sin(x)*x
expr.diff()

Or, for multivariate,

from sympy.abc import x,y
expr = x**3-sympy.sin(y)*x

expr.diff(y)
expr.diff(x)

With a similar syntax for simple integrals:

expr.integrate(x)
expr.integrate(y)

Note the lack of a +C constant; to get constants, add them in yourself or get the general solution to a simple differential equation, which is more complete.

Differential Equations

We won't get into this much here, but there are many solvers for differential equation systems.

As a simple example,

f,g = sympy.symbols('f g',cls=sympy.Function)

expr = f(x).diff(x,x)-2*f(x).diff(x) + 1

We can now solve the differential equation to get a definition for $f$:

eqn = sympy.dsolve(expr,f(x))

And extract the solution as an expression, and its newly defined symbols;

solution = eqn.rhs

C1,C2,x = solution.free_symbols

Quadratic Extensions

The quadratic field we are all most familiar with is the Gaussian Rationals; for those, we can use sympy.I as the imaginary constant;

a = 1+sympy.I
abs(a)

For other quadratic fields, we can construct them either explicitly - with sympy.sqrt - or by minimal polynomial:

poly = x**2-x-1
a = sympy.solve(poly,x)[0]

b = (2*a-3)**2

We may sometimes need to explicitly expand:

b.expand()

but other than that they behave almost entirely as expected.

Sympy also works fairly well for other extensions by $n$th roots, and extensions of extensions, but (understandably) struggles somewhat on explicit computations in more general algebraic extensions.

Linear Algebra

We have already seen constructing matrices, and can imagine multiplying and adding them.

In SymPy, most of the more common matrix operations - standard decompositions, and so on - are stored as properties of the Matrix class:

M = sympy.Matrix([[1,2],[3,4]])
M.QRdecomposition()

To get a full picture of the available methods, you can look through dir(M) on a matrix. The names have been chosen to try to make it easy for mathematicians to understand what the functions do.

Note that you can do symbolic math on matrices - either matrices of symbols:

a_11,a_12,a_21,a_22 =sympy.symbols("a_{1\,1} a_{1\,2} a_{2\,1} a_{2\,2}")

M = sympy.Matrix([[a_11,a_12],[a_21,a_22]])
M.det()

Or matrices that are symbols:

N = sympy.MatrixSymbol('N',2,2)
M*N*M.inv()

Parsing

We will discuss parsing more next week, but play around with sympy.parse to get a feel for how your representations of math translate.

Worksheet

Today's worksheet has you solving a classic related rates problem with Sympy - using its trigonometry, calculus, and solving functionality.