Symbolic Analysis with SymPy

This introduction to SymPy was written using the following sources:

Another important source of documentation about SymPy is the official documentation.

Instructions: Create a new directory called SymPy with a notebook called SymPyTour. Give it a heading 1 cell title Symbolic Analysis with SymPy. Read this page, typing in the code in the code cells and executing them as you go.

Do not copy/paste.

Type the commands yourself to get the practice doing it. This will also slow you down so you can think about the commands and what they are doing as you type them.</font>

Save your notebook when you are done, then try the accompanying exercises.

Introduction

SymPy is a symbolic mathematics package for Python. In symbolic mathematics, the mathematical entities we want to manipulate are not necessarily concrete, specific numbers such as $2.0$. Instead, they can be symbols, such as $x$ or $y$ that only sometimes take on concrete values.

Let's start out by importing necessary packages from IPython:



In [ ]:

    
from IPython.html.widgets import interact
from IPython.display import display

It is important to not use %pylab inline when using SymPy. The reason for this is that %pylab imports names that conflict with those in SymPy. Instead, use %matplotlib inline and manually import numpy and matplotlib.



In [ ]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

The sympy namespace imports many things when imported using Python's * notation. You should only do this when you are confident the names imported by sympy will not conflict with other names in your code.



In [ ]:

    
from sympy import *

When using SymPy in the IPython Notebook, it is recommended to initialize $\LaTeX$ printing of expressions:



In [ ]:

    
init_printing(use_latex='mathjax')

Symbolic variables

In SymPy you need to create symbols for the symbolic variables you want to work with. You can create a new symbol using the Symbol class:



In [ ]:

    
x = Symbol('x')

Once you have created symbols, you can use them in mathematical expressions:



In [ ]:

    
(pi + x)**2

The function IPython.display.display works just like print, but also handles rich representations such as $\LaTeX$. This is convenient as you can use display inside loops and other structures.



In [ ]:

    
display((pi + x)**2)

You can create a set of symbolic using the symbols function:



In [ ]:

    
a, b, c = symbols("a, b, c")

You can add mathematical assumptions to symbols when you create them by passing in keyword arguments. This creates the symbol x with the assumption that it is a real number:



In [ ]:

    
x = Symbol('x', real=True)

The set of all assumptions about a symbol can be probed using a set of attributes with the naming convention is_[assumption]:



In [ ]:

    
x.is_imaginary

When an assumption not determined, it will return None:



In [ ]:

    
x.is_positive



In [ ]:

    
x = Symbol('x', positive=True)



In [ ]:

    
x > 0

Assumptions propagate through mathematical expressions and are computed dynamically. Negating a positive number gives a negative one:



In [ ]:

    
(-x).is_positive, (-x).is_negative

Here is a list of all assumptions that symbols know about. What happens if you remove the k == k.lower() test from this code? What do you think those additional is_* attributes correspond to?



In [ ]:

    
for k in dir(x):
    if k.startswith('is_') and k == k.lower():
        print k

Complex numbers

The imaginary unit is denoted I in SymPy:



In [ ]:

    
1+1*I



In [ ]:

    
I**2



In [ ]:

    
I.is_imaginary



In [ ]:

    
(x*I + 1)**2

The abs function will take the complex absolute value:



In [ ]:

    
abs(x*I + 1)

Real, Integer and Rational numbers

There are three basic numerical types in SymPy:

Real
Rational
Integer

The Real and Integer classes for SymPy are different from Python's builtin float and int types. However, you will almost never need to manually create Real and Integer instances. Any time a Python float or int appears in a symbolic expression, they are converted to a SymPy Real or Integer instance:



In [ ]:

    
srepr(2.0*x)

When working with numerical types in SymPy it is important to carefully distinguish SymPy types from the native Python numerical types (float and int). Python's builtin type function can be used to show the type of an object. You can use sympify to create a SymPy type from a Python type. To go in the other direction, using Python's builtin float or int:



In [ ]:

    
print type(1.0)
print type(sympify(1.0))
print type(float(sympify(1.0)))

Rational numbers must be created manually by passing the numerator and denominator to the Rational class:



In [ ]:

    
r1 = Rational(4,5)
r2 = Rational(5,4)

Basic arithmetic operations work with rational numbers:



In [ ]:

    
r1



In [ ]:

    
r1+r2



In [ ]:

    
sqrt(r1/r2)



In [ ]:

    
r1.is_rational, r1.is_Rational

Numerical evaluation

SymPy integers and real numbers support arbitrary precision. To evaluate a symbolic expression numerically to a given precision use the evalf method:



In [ ]:

    
pi.evalf(n=50)

You can also use the N function to accomplish the same thing:



In [ ]:

    
y = (x + pi)**2



In [ ]:

    
N(y, 5)

When you numerically evaluate algebraic expressions you often want to substitute a symbol with a numerical value. You can do that using the subs function:



In [ ]:

    
y.subs(x, 1.5)



In [ ]:

    
N(y.subs(x, 1.5))

The subs function can also be used to substitute Symbols and other expressions:



In [ ]:

    
y.subs(x, a+pi)

You can also pass subs a dict of keys and values to be substituted:



In [ ]:

    
e = x*y



In [ ]:

    
e.subs({x:4, y:pi})

Often times, it is useful to convert from a SymPy expression to a Python function that can be used with NumPy arrays. This can be done using SymPy's lambdify function in the following manner:



In [ ]:

    
z = Symbol('z')
f = lambdify([z], cos(pi*z)**2, modules='numpy')
x = np.linspace(0.0, 2*np.pi, 100)
y = f(x)
plt.plot(x, y)

Algebra

SymPy has a number of powerful functions and methods for performing algebraic manipulations.



In [ ]:

    
x, y, z = symbols('x, y, z')

expand and factor

The expand and factor functions are roughly inverses of each other. The expand function will expand a product of expressions into a sum of expression:



In [ ]:

    
e = (x+1)*(x+2)*(x+3)
e



In [ ]:

    
expand(e)

factor does the inverse operation:



In [ ]:

    
factor(_)

expand also has a number of keyword arguments that enable it to expand other types of mathematical expressions, such as trigonometric function or complex numbers:



In [ ]:

    
e = sin(a+b)



In [ ]:

    
expand(e)



In [ ]:

    
expand(e, trig=True)



In [ ]:

    
expand((x+y)**3)



In [ ]:

    
expand(x+y, complex=True)

simplify

A common task in symbolic mathematics is to "simplify" an expression. SymPy has a simplify function that applies a set of transformations using heuristics in an attempt to produce the most simplified version of an expression:



In [ ]:

    
simplify((x+x*y)/x)



In [ ]:

    
simplify(sin(a)**2 + cos(a)**2)



In [ ]:

    
simplify(cos(x)/sin(x))

collect and coeff

When working with polynomials, it is often useful to pick out terms that are powers of a given expression. This can be done with collect:



In [ ]:

    
e = x*y + x - 3 + 2*x**2 - z*x**2 + x**3



In [ ]:

    
collect(e, x)

The coeff method will return the coefficient of a given term:



In [ ]:

    
collect(e, x).coeff(x, 2)

cancel

The cancel function will try to cancel common factors in fractional expressions:



In [ ]:

    
cancel((x**2 + 2*x + 1)/(x**2 + x))



In [ ]:

    
cancel((exp(x)+2.0*exp(x))/(exp(x)))

Other simplification functions

SymPy has a number of other functions for performing algebraic manipulations:

trigsimp
expand_trig
powsimp
expand_power_exp
expand_power_base
powdenest
expand_log
logcombine
rewrite
expand_func
hyperexpand
combsimp

Calculus

SymPy is able to perform many operations from calculus using symbolic expressions. This includes differentiation, integration, Taylor series, limits, sums and products.



In [ ]:

    
x, y, z = symbols('x, y, z')
a, b, c = symbols('a, b, c')

Differentiation

You can take a deriviate using the diff function:



In [ ]:

    
diff(a*x**2, x)



In [ ]:

    
diff(cos(a*x+b), x)

Derivatives of higher orders are specified by a third argument to diff. Here is a second derivative:



In [ ]:

    
diff(cos(a*x+b), x, 2)

diff can also be used to take partial derivatives:



In [ ]:

    
f = sin(x*y) + cos(y*z)



In [ ]:

    
diff(f, x)



In [ ]:

    
diff(f, y)



In [ ]:

    
diff(f, z)

Note that diff takes a derivative immediately. To construct an unevaluated derivative, use the Derivative objects, which accepts the same arguments as diff. Here is a mixed second partial derivative:



In [ ]:

    
d = Derivative(f, x, 1, y, 1)

An unevaluated derivative can be evaluated using the doit() method:



In [ ]:

    
d.doit()

Integration

Symbolic integrals can be done using the integrate function:



In [ ]:

    
integrate(6*x**5, x)



In [ ]:

    
integrate(sin(x), x)



In [ ]:

    
integrate(log(x), x)



In [ ]:

    
integrate(2*x + sinh(x), x)

SymPy also knows many integrals involving special functions:



In [ ]:

    
integrate(exp(-x**2)*erf(x), x)

Definite integrals can be done by specifying limits in the form (variable, min, max):



In [ ]:

    
integrate(x**3, (x, -1, 1))



In [ ]:

    
integrate(sin(x), (x, 0, pi/2))



In [ ]:

    
integrate(cos(x), (x, -pi/2, pi/2))

In SymPy the symbol for $\infty$ is oo (two lowercase "o" characters). This allows integrals over infinite ranges:



In [ ]:

    
integrate(exp(-x), (x, 0, oo))



In [ ]:

    
integrate(exp(-x**2), (x, -oo, oo))

Double integrals can be performed by giving a sequence of variables and/or ranges to integrate over. In this case, the Integral class is used to construct an unevaluated integral:



In [ ]:

    
i = Integral(exp(-x**2 - y**2), (x,-oo,oo), (y,-oo,oo))
i

Again, the doit() method performs the integrals:



In [ ]:

    
i.doit()

Sums and products

Finite and infinite sums and products can be construted, in both evaluated and unevaluated forms.



In [ ]:

    
n = Symbol('n')

Use the summation function to directly perform a sum:



In [ ]:

    
s = summation(1/n**2, (n, 1, 10))
s

Here is an unevaluated infinite sum:



In [ ]:

    
s = Sum(1/n**2, (n, 1, oo))

The doit() method will attempt to produce an exact, symbolic answer:



In [ ]:

    
s.doit()

If that is not possible, you can always use the evalf method to obtain a numerical answer to a given precision:



In [ ]:

    
s.evalf()

Here is an unevaluated product:



In [ ]:

    
p = Product(n, (n, 1, 10))
p



In [ ]:

    
p.doit()

Right now evalf doesn't work with Product, which looks like a bug:



In [ ]:

    
p.evalf()

Limits

SymPy is capable of taking limits symbolically using the limit function:



In [ ]:

    
limit(sin(x)/x, x, 0)

Let's use the definition of the derivative as a limit to check the derivative of this function:



In [ ]:

    
f



In [ ]:

    
diff(f, x)



In [ ]:

    
h = Symbol('h')



In [ ]:

    
df_dx = (f.subs(x, x+h) - f)/h
df_dx

Let's construct an unevaluated limit using the Limit class and then take the limit $h\rightarrow0$ to recover the exact answer above:



In [ ]:

    
l = Limit(df_dx, h, 0)
l



In [ ]:

    
l.doit()

Some other examples of limits:



In [ ]:

    
limit(x, x, oo)



In [ ]:

    
limit(1/x, x, oo)

limit and Limit also accept a dir keyword argument that can be used to specify the direction of the limit:



In [ ]:

    
limit(1/x, x, 0, dir='+')



In [ ]:

    
limit(1/x, x, 0, dir='-')

Taylor series expansion

SymPy can perform Taylor series about arbitrary points using the series function.

Expand the function $e^x$ about the default point ($x=0$) to the default order (5):



In [ ]:

    
series(exp(x), x)

Here is the Taylor series about zero, up to $7th$ order:



In [ ]:

    
series(exp(x), x, 0, 7)

Let's plot the series approximation alongside the exact function to get a sense of how well the approximate works. Note how we use the removeO method to remove the $\mathcal{O}(x^7)$ term:



In [ ]:

    
plot(exp(x), _.removeO(), (x,0-5,0+5))

Solving equations

SymPy can solve algebraic equations symbolically using the solve function.



In [ ]:

    
x, y, a, b = symbols('x, y, a, b')

If solve is passed an expression, it is assumed that expression is equal to $0$ when converting it to an equation.



In [ ]:

    
solve(x**4 - 1, x)

Note how in this form, solve returns a list of solutions. In general, the number of solutions will correspond to the order of the equation.

If solve is passed dict=True, it will return a dictionary of solutions. This can be useful if you want to substitute the solutions later using subs.



In [ ]:

    
solve(x**4 - 1, x, dict=True)

You can also pass solve an equation object constructed using the Eq class. This let's you specify non-zero right and left hand sides of the equation:



In [ ]:

    
eq = Eq(x**4, 1)
eq



In [ ]:

    
soln = solve(eq, dict=True)
soln

We can then substitute the solutions to verify that they solve the orginal equation:



In [ ]:

    
[eq.subs(s) for s in soln]

You can also pass a list of equations and variables to solve for:



In [ ]:

    
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])

Linear Algebra

In most cases, you will use NumPy for linear algebra. This works well if each element of your arrays is a concrete number. However, in some cases you will want to work with matrices of symbolic expressions. That can be done using the SymPy Matrix class. The API presented by the Matrix class is very similar to that of NumPy.



In [ ]:

    
x = Symbol('x')
y = Symbol('y')



In [ ]:

    
A = Matrix([[1,x], [y,1]])
A



In [ ]:

    
A[0,0]



In [ ]:

    
A[:,1]

You can create symbolic expressions using these matrices:



In [ ]:

    
A**2

Here is a symbolic matrix inverse:



In [ ]:

    
A.inv()

For small matrices, you can even find eigenvalues symbolically:



In [ ]:

    
A.eigenvals()

The subs method works as expected:



In [ ]:

    
A.subs({x:0,y:I})

All content is under a modified MIT License, and can be freely used and adapted. See the full license text here.