# SymPy - Symbolic algebra in Python

J.R. Johansson (jrjohansson at gmail.com)

The latest version of this [IPython notebook](http://ipython.org/notebook.html) lecture is available at [http://github.com/jrjohansson/scientific-python-lectures](http://github.com/jrjohansson/scientific-python-lectures).

The other notebooks in this lecture series are indexed at [http://jrjohansson.github.io](http://jrjohansson.github.io).


```python
%matplotlib inline
import matplotlib.pyplot as plt
```

## Introduction

There are two notable Computer Algebra Systems (CAS) for Python:

* [SymPy](http://sympy.org/en/index.html) - A python module that can be used in any Python program, or in an IPython session, that provides powerful CAS features. 
* [Sage](http://www.sagemath.org/) - Sage is a full-featured and very powerful CAS enviroment that aims to provide an open source system that competes with Mathematica and Maple. Sage is not a regular Python module, but rather a CAS environment that uses Python as its programming language.

Sage is in some aspects more powerful than SymPy, but both offer very comprehensive CAS functionality. The advantage of SymPy is that it is a regular Python module and integrates well with the IPython notebook. 

In this lecture we will therefore look at how to use SymPy with IPython notebooks. If you are interested in an open source CAS environment I also recommend to read more about Sage.

To get started using SymPy in a Python program or notebook, import the module `sympy`:


```python
from sympy import *
```

To get nice-looking $\LaTeX$ formatted output run:


```python
init_printing()

# or with older versions of sympy/ipython, load the IPython extension
#%load_ext sympy.interactive.ipythonprinting
# or
#%load_ext sympyprinting
```

## Symbolic variables

In SymPy we need to create symbols for the variables we want to work with. We can create a new symbol using the `Symbol` class:


```python
x = Symbol('x')
```


```python
(pi + x)**2
```


```python
# alternative way of defining symbols
a, b, c = symbols("a, b, c")
```


```python
type(a)
```




    sympy.core.symbol.Symbol



We can add assumptions to symbols when we create them:


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


```python
x.is_imaginary
```




    False




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


```python
x > 0
```

### Complex numbers

The imaginary unit is denoted `I` in SymPy. 


```python
1+1*I
```


```python
I**2
```


```python
(x * I + 1)**2
```

### Rational numbers

There are three different numerical types in SymPy: `Real`, `Rational`, `Integer`: 


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


```python
r1
```


```python
r1+r2
```


```python
r1/r2
```

## Numerical evaluation

SymPy uses a library for arbitrary precision as numerical backend, and has predefined SymPy expressions for a number of mathematical constants, such as: `pi`, `e`, `oo` for infinity.

To evaluate an expression numerically we can use the `evalf` function (or `N`). It takes an argument `n` which specifies the number of significant digits.


```python
pi.evalf(n=50)
```


```python
y = (x + pi)**2
```


```python
N(y, 5) # same as evalf
```

When we numerically evaluate algebraic expressions we often want to substitute a symbol with a numerical value. In SymPy we do that using the `subs` function:


```python
y.subs(x, 1.5)
```


```python
N(y.subs(x, 1.5))
```

The `subs` function can of course also be used to substitute Symbols and expressions:


```python
y.subs(x, a+pi)
```

We can also combine numerical evaluation of expressions with NumPy arrays:


```python
import numpy
```


```python
x_vec = numpy.arange(0, 10, 0.1)
```


```python
y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
```


```python
fig, ax = plt.subplots()
ax.plot(x_vec, y_vec);
```

However, this kind of numerical evaluation can be very slow, and there is a much more efficient way to do it: Use the function `lambdify` to "compile" a SymPy expression into a function that is much more efficient to evaluate numerically:


```python
f = lambdify([x], (x + pi)**2, 'numpy')  # the first argument is a list of variables that
                                         # f will be a function of: in this case only x -> f(x)
```


```python
y_vec = f(x_vec)  # now we can directly pass a numpy array and f(x) is efficiently evaluated
```

The speedup when using "lambdified" functions instead of direct numerical evaluation can be significant, often several orders of magnitude. Even in this simple example we get a significant speed up:


```python
%%timeit

y_vec = numpy.array([N(((x + pi)**2).subs(x, xx)) for xx in x_vec])
```

    10 loops, best of 3: 28.2 ms per loop



```python
%%timeit

y_vec = f(x_vec)
```

    The slowest run took 8.86 times longer than the fastest. This could mean that an intermediate result is being cached 
    100000 loops, best of 3: 2.93 µs per loop


## Algebraic manipulations

One of the main uses of an CAS is to perform algebraic manipulations of expressions. For example, we might want to expand a product, factor an expression, or simplify an expression. The functions for doing these basic operations in SymPy are demonstrated in this section.

### Expand and factor

The first steps in an algebraic manipulation 


```python
(x+1)*(x+2)*(x+3)
```


```python
expand((x+1)*(x+2)*(x+3))
```

The `expand` function takes a number of keywords arguments which we can tell the functions what kind of expansions we want to have performed. For example, to expand trigonometric expressions, use the `trig=True` keyword argument:


```python
sin(a+b)
```


```python
expand(sin(a+b), trig=True)
```

See `help(expand)` for a detailed explanation of the various types of expansions the `expand` functions can perform.

The opposite of product expansion is of course factoring. To factor an expression in SymPy use the `factor` function: 


```python
factor(x**3 + 6 * x**2 + 11*x + 6)
```

### Simplify

The `simplify` tries to simplify an expression into a nice looking expression, using various techniques. More specific alternatives to the `simplify` functions also exists: `trigsimp`, `powsimp`, `logcombine`, etc. 

The basic usages of these functions are as follows:


```python
# simplify (sometimes) expands a product
simplify((x+1)*(x+2)*(x+3))
```


```python
# simplify uses trigonometric identities
simplify(sin(a)**2 + cos(a)**2)
```


```python
simplify(cos(x)/sin(x))
```

### apart and together

To manipulate symbolic expressions of fractions, we can use the `apart` and `together` functions:


```python
f1 = 1/((a+1)*(a+2))
```


```python
f1
```


```python
apart(f1)
```


```python
f2 = 1/(a+2) + 1/(a+3)
```


```python
f2
```


```python
together(f2)
```

Simplify usually combines fractions but does not factor: 


```python
simplify(f2)
```

## Calculus

In addition to algebraic manipulations, the other main use of CAS is to do calculus, like derivatives and integrals of algebraic expressions.

### Differentiation

Differentiation is usually simple. Use the `diff` function. The first argument is the expression to take the derivative of, and the second argument is the symbol by which to take the derivative:


```python
y
```


```python
diff(y**2, x)
```

For higher order derivatives we can do:


```python
diff(y**2, x, x)
```


```python
diff(y**2, x, 2) # same as above
```

To calculate the derivative of a multivariate expression, we can do:


```python
x, y, z = symbols("x,y,z")
```


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

$\frac{d^3f}{dxdy^2}$


```python
diff(f, x, 1, y, 2)
```

## Integration

Integration is done in a similar fashion:


```python
f
```


```python
integrate(f, x)
```

By providing limits for the integration variable we can evaluate definite integrals:


```python
integrate(f, (x, -1, 1))
```

and also improper integrals


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

Remember, `oo` is the SymPy notation for inifinity.

### Sums and products

We can evaluate sums using the function `Sum`:


```python
n = Symbol("n")
```


```python
Sum(1/n**2, (n, 1, 10))
```


```python
Sum(1/n**2, (n,1, 10)).evalf()
```


```python
Sum(1/n**2, (n, 1, oo)).evalf()
```

Products work much the same way:


```python
Product(n, (n, 1, 10)) # 10!
```

## Limits

Limits can be evaluated using the `limit` function. For example, 


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

We can use 'limit' to check the result of derivation using the `diff` function:


```python
f
```


```python
diff(f, x)
```

$\displaystyle \frac{\mathrm{d}f(x,y)}{\mathrm{d}x} = \frac{f(x+h,y)-f(x,y)}{h}$


```python
h = Symbol("h")
```


```python
limit((f.subs(x, x+h) - f)/h, h, 0)
```

OK!

We can change the direction from which we approach the limiting point using the `dir` keywork argument:


```python
limit(1/x, x, 0, dir="+")
```


```python
limit(1/x, x, 0, dir="-")
```

## Series

Series expansion is also one of the most useful features of a CAS. In SymPy we can perform a series expansion of an expression using the `series` function:


```python
series(exp(x), x)
```

By default it expands the expression around $x=0$, but we can expand around any value of $x$ by explicitly include a value in the function call:


```python
series(exp(x), x, 1)
```

And we can explicitly define to which order the series expansion should be carried out:


```python
series(exp(x), x, 1, 10)
```

The series expansion includes the order of the approximation, which is very useful for keeping track of the order of validity when we do calculations with series expansions of different order:


```python
s1 = cos(x).series(x, 0, 5)
s1
```


```python
s2 = sin(x).series(x, 0, 2)
s2
```


```python
expand(s1 * s2)
```

If we want to get rid of the order information we can use the `removeO` method:


```python
expand(s1.removeO() * s2.removeO())
```

But note that this is not the correct expansion of $\cos(x)\sin(x)$ to $5$th order:


```python
(cos(x)*sin(x)).series(x, 0, 6)
```

## Linear algebra

### Matrices

Matrices are defined using the `Matrix` class:


```python
m11, m12, m21, m22 = symbols("m11, m12, m21, m22")
b1, b2 = symbols("b1, b2")
```


```python
A = Matrix([[m11, m12],[m21, m22]])
A
```


```python
b = Matrix([[b1], [b2]])
b
```

With `Matrix` class instances we can do the usual matrix algebra operations:


```python
A**2
```


```python
A * b
```

And calculate determinants and inverses, and the like:


```python
A.det()
```


```python
A.inv()
```

## Solving equations

For solving equations and systems of equations we can use the `solve` function:


```python
solve(x**2 - 1, x)
```


```python
solve(x**4 - x**2 - 1, x)
```

System of equations:


```python
solve([x + y - 1, x - y - 1], [x,y])
```

In terms of other symbolic expressions:


```python
solve([x + y - a, x - y - c], [x,y])
```

## Further reading

* http://sympy.org/en/index.html - The SymPy projects web page.
* https://github.com/sympy/sympy - The source code of SymPy.
* http://live.sympy.org - Online version of SymPy for testing and demonstrations.

## Versions


```python
%reload_ext version_information

%version_information numpy, matplotlib, sympy
```




<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>2.7.10 64bit [GCC 4.2.1 (Apple Inc. build 5577)]</td></tr><tr><td>IPython</td><td>3.2.1</td></tr><tr><td>OS</td><td>Darwin 14.1.0 x86_64 i386 64bit</td></tr><tr><td>numpy</td><td>1.9.2</td></tr><tr><td>matplotlib</td><td>1.4.3</td></tr><tr><td>sympy</td><td>0.7.6</td></tr><tr><td colspan='2'>Sat Aug 15 11:37:37 2015 JST</td></tr></table>



b8ceade190497b6f196ba0f059e9671e256eed37
ipynb
Jupyter Notebook
Lecture-5-Sympy.ipynb
hbayraktaroglu/scientific-python-lectures
e9221b58153aad15c323e3af5f2fde999bf8b7a1
CC-BY-3.0
2015-01-03T05:28:29.000Z
2022-03-31T21:06:24.000Z
Lecture-5-Sympy.ipynb
hbayraktaroglu/scientific-python-lectures
e9221b58153aad15c323e3af5f2fde999bf8b7a1
CC-BY-3.0
2015-04-12T06:17:06.000Z
2021-02-02T09:57:19.000Z
Lecture-5-Sympy.ipynb
hbayraktaroglu/scientific-python-lectures
e9221b58153aad15c323e3af5f2fde999bf8b7a1
CC-BY-3.0
2015-01-05T00:48:26.000Z
2022-03-21T05:51:01.000Z
Qwen/Qwen-72B
1. YES
2. YES