misc

Sympy

Following the Sympy tutorial

SymPy Tutorial

This tutorial assumes a decent mathematical background.

What is Symbolic Computation?

the mathematical objects are represented exactly, not approximately, and mathematical expressions with unevaluated variables are left in symbolic form.

1
2
3
4
5
6
7
8
9
10
>>> import sympy as s
>>> import math as m
>>> m.sqrt(8)
2.8284271247461903
>>> s.sqrt(8)
2*sqrt(2)
>>> m.sqrt(8)*m.sqrt(8) == 8.0
False
>>> s.sqrt(8)*s.sqrt(8) == 8.0
True

Unlike many symbolic manipulation systems, variables in SymPy must be defined before they are used

Create symbols

1
2
3
4
5
6
7
>>> x,y=s.symbols("x y")
>>> x,y
(x, y)
>>> x
x
>>> type(x)
<class 'sympy.core.symbol.Symbol'>

Use symbols

symbols takes a string of variable names separated by spaces or commas, and creates Symbols out of them.

1
2
3
4
5
6
7
8
9
>>> expr = x + 2*y
>>> expr
x + 2*y
>>> expr * x
x*(x + 2*y)
>>> expr - x
2*y
>>> expr / x
(x + 2*y)/x

most simplifications are not performed automatically.

Both forms are useful in different circumstances.

In SymPy, there are functions to go from one form to the other

Simplification

1
2
3
4
5
6
7
8
9
10
>>> simplify = s.expand
>>> factor = s.factor
>>> expr
x + 2*y
>>> expr*x
x*(x + 2*y)
>>> simplify(expr*x)
x**2 + 2*x*y
>>> factor(_)
x*(x + 2*y)

The Power of Symbolic Computation

The real power of a symbolic computation system such as SymPy is the ability to do all sorts of computations symbolically. SymPy can simplify expressions, compute derivatives, integrals, and limits, solve equations, work with matrices, and much, much more, and do it all symbolically. It includes modules for plotting, printing (like 2D pretty printed output of math formulas, or LATEX), code generation, physics, statistics, combinatorics, number theory, geometry, logic, and more.

Prettify the output

1
2
>>> help(s.init_printing)
>>> s.init_printing() # This will automatically enable the best printer available in your environment.
1
2
3
4
init_printing(pretty_print=True) # default argument
>>> s.init_printing()
>>> s.sqrt(8)
2⋅√2
1
2
3
4
5
6
>>> s.init_printing(use_unicode=True)
>>> s.sqrt(8)
2⋅√2
>>> s.diff(s.sin(x)*s.exp(x) , x)
 x           x
ℯ ⋅sin(x) + ℯ ⋅cos(x)
1
s.init_printing(use_latex='mathjax')

Get the latex expression

1
2
>>> s.latex(s.diff(s.sin(x)*s.exp(x) , x))
e^{x} \sin{\left(x \right)} + e^{x} \cos{\left(x \right)}

Derivative

1
s.diff(expression , x)

Integral

1
2
3
4
5
s.integrate(expression  , x)
s.integrate(exprssion , (x,-oo,+oo))

>>> s.integrate(s.sin(x) , (x,s.pi/4,s.pi*3/4))
√2

Limit

1
2
3
4
s.limit(expression,x,to_what)

>>> s.limit(s.sin(x)/x , x , 0)
1

Normal equation

1
2
3
4
s.solve(expression , x) # expression = 0

>>> s.solve(x**2 - 2, x)
[-√2, √2]

Differential equation

1
2
3
4
5
6
7
8
y = s.Function("y")
s.dsovle(s.Eq(expression , target ) , y(x))
# expression = target

>>> s.dsolve(s.Eq(y(x).diff(x,x)-y(x) , s.exp(x)) , y(x))
           -x   ⎛     x⎞  x
y(x) = C₂⋅ℯ   + ⎜C₁ + ─⎟⋅ℯ
                ⎝     2⎠

Gotcha

To begin, we should make something about SymPy clear.

SymPy is nothing more than a Python library, like NumPy, Django, or even modules in the Python standard library sys or re.

What this means is that SymPy does not add anything to the Python language.

Gotcha1

Symbols can have names longer than one character if we want.

1
crazy = sympy.symbols("unrelated")

==So it’s a good idea to let symbol names and Python variable names coincide !!==

Update value in a expression
  1. Evaluating an expression at a point. For example, if our expression is cos(x) + 1 and we want to evaluate it at the point x = 0, so that we get cos(0) + 1, which is 2.
  2. Replacing a subexpression with another subexpression.
  3. To perform multiple substitutions at once, pass a list of (old, new) pairs to subs.
1
2
3
4
5
>>> expr = x + y(x)
>>> expr
x + y(x)
>>> expr.subs(x,s.sqrt(2))
y(√2) + √2
1
2
3
4
>>> expr
3*x + y
>>> expr.subs([(x,5),(y,1)])
16

There are two important things to note about subs. First, it returns a new expression. SymPy objects are immutable. That means that subs does not modify it in-place.

==In fact, since SymPy expressions are immutable, no function will change them in-place. All functions will return new expressions !==

Gotcha2

SymPy does not extend Python syntax is that = does not represent equality in SymPy.

There is a separate object, called Eq, which can be used to create symbolic equalities

1
2
3
>>> s.Eq(y(x) , s.exp(x))
        x
y(x) = ℯ

Gotcha3

Do not use == to test if two expressions are equal, use a.equals(b) instead

1
2
3
4
>>> a = cos(x)**2 - sin(x)**2
>>> b = cos(2*x)
>>> a.equals(b)
True

Gotcha4

Combination of different types

Whenever you combine a SymPy object and a SymPy object, or a SymPy object and a Python object, you get a SymPy object, but whenever you combine two Python objects, SymPy never comes into play, and so you get a Python object.

SymPy Features

Converting Strings to SymPy Expressions

The sympify function (that’s sympify, not to be confused with simplify) can be used to convert strings into SymPy expressions.

Warning : sympify uses eval. Don’t use it on unsanitized input.

1
2
3
>>> new = sympify('z**4 + z**2 + z')
>>> new
z**4 + z**2 + z

Evaluate an expression

To evaluate a numerical expression into a floating point number, use evalf.

SymPy can evaluate floating point expressions to arbitrary precision. By default, 15 digits of precision are used, but you can pass any number as the argument to evalf.

1
2
3
4
5
6
>>> expr.subs(y , x**2)
x**2 + 3*x
>>> _.subs(x,sqrt(2))
2 + 3*sqrt(2)
>>> _.evalf(10)
6.242640687

it is more efficient and numerically stable to pass the substitution to evalf using the subs flag, which takes a dictionary of Symbol: point pairs.

1
2
3
4
>>> new.evalf( 3, subs={z:5})
655.
>>> new
z**4 + z**2 + z

Sometimes there are roundoff errors smaller than the desired precision that remain after an expression is evaluated. Such numbers can be removed at the user’s discretion by setting the chop flag to True.

1
2
3
4
5
>>> one = cos(1)**2 + sin(1)**2
>>> (one - 1).evalf()
-0.e-124
>>> (one - 1).evalf(chop=True)
0

if you wanted to evaluate an expression at a thousand points, using SymPy would be far slower than it needs to be, especially if you only care about machine precision. Instead, you should use libraries like NumPy and SciPy.

Lambdify

lambdify acts like a lambda function, except it converts the SymPy names to the names of the given numerical library, usually NumPy. 

1
2
3
4
5
6
7
>>> import numpy
>>> domain = numpy.arange(10)
>>> domain
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> function = lambdify(z,new,"numpy")
>>> function(domain)
array([   0,    3,   22,   93,  276,  655, 1338, 2457, 4168, 6651])

Warning : lambdify uses eval. Don’t use it on unsanitized input.

Printing

There are several printers available in SymPy. The most common ones are

  • str
  • srepr
  • ASCII pretty printer
  • Unicode pretty printer
  • LaTeX
  • MathML
  • Dot

In order to see latex output, you should in an GUI Python console, for example, JupyterNoteBook

latex display

Simplification

SymPy has dozens of functions to perform various kinds of simplification. There is also one general function called simplify() that attempts to apply all of these functions in an intelligent way to arrive at the simplest form of an expression.

the gamma function
$$
\Gamma(n) = (n-1)!
$$

1
2
>>> simplify(gamma(x)/gamma(x-1))
x - 1

simplify() is best when used interactively, when you just want to whittle down an expression to a simpler form. You may then choose to apply specific functions once you see what simplify() returns, to get a more precise result. It is also useful when you have no idea what form an expression will take, and you need a catchall function to simplify it

Polynomial/Rational Function Simplification

Given a polynomial, expand() will put it into a canonical form of a sum of monomials.

1
2
3
>>> expand( (x-2)**2 + (x+sqrt(2))**2)
   2
2⋅x  - 4⋅x + 2⋅√2⋅x + 6

factor() takes a polynomial and factors it into irreducible factors over the rational numbers. For polynomials, factor() is the opposite of expand().

1
2
3
>>> factor( x**2 - 2*x + 1)
       2
(x - 1)

collect() collects common powers of a term in an expression. 

1
2
3
>>> collect( 2*x**2 + 4*x**4 + 16*x**16 + 7*x**7 , x)
    16      7      4      2
16⋅x   + 7⋅x  + 4⋅x  + 2⋅x

expr.coeff(x, n) gives the coefficient of x**n in expr

1
2
>>> collect( 2*x**2 + 4*x**4 + 16*x**16 + 7*x**7 , x).coeff(x,16)
16

cancel() will take any rational function and put it into the standard canonical form

1
2
3
>>> cancel( (x-2)**3 / (x-2) )
 2
x  - 4⋅x + 4

apart() performs a partial fraction decomposition on a rational function

Trigonometric Simplification

Power

power

This is important to remember, because by default, SymPy will not perform simplifications if they are not true in general.

1
2
3
# so we can specify type of number when we defining they
>>> x, y = symbols('x y', positive=True)
>>> a, b = symbols('a b', real=True)

Exponentials and logarithms

Logarithms have similar issues as powers.

Special Functions

SymPy implements dozens of special functions, ranging from functions in combinatorics to mathematical physics.

  • factorial(n)
  • binomial(n, k)
  • ….

Rewrite

To rewrite an expression in terms of a function,useexpr.rewrite(function)

1
2
>>> gamma(x).rewrite(factorial)
(x - 1)!

Piecewise Functions

1
formula = Piecewise( (x**2, x>2), (1+x, x<=2) )

Combine Functions

1
2
3
4
5
6
7
8
9
10
11
12
13
>>> f
⎧0.5⋅log(x)  for x ≥ 1
⎨
⎩ 2⋅x - 1    otherwise
>>> ff = s.simplify(f.subs(x,f))
>>> ff
⎧⎧0.5⋅log(log(x)) - 0.346573590279973  for x ≥ 1
⎪⎨                                                for x ≥ 7.38905609893065
⎪⎩         0.5⋅log(2⋅x - 1)            otherwise
⎨
⎪          ⎧1.0⋅log(x) - 1  for x ≥ 1
⎪          ⎨                                             otherwise
⎩          ⎩   4⋅x - 3      otherwise

$$
f = \begin{cases} \ln{\sqrt{x}} , & x \ge 1\
2x-1 , & x<1
\end{cases}
$$

$$
ff = f(f(x))
$$

Calculus

Tips

Do differential n times

1
2
3
4
5
6
7
8
9
10
11
12
>>> def d(f,n):
...     return f if n == 0 else d(s.diff(f,x) , n-1)
>>> d(x**4 , 3)
24⋅x
>>> d(x**4 , 1)
   3
4⋅x
>>> d(x**4 , 2)
    2
12⋅x
>>> d(x**4 , 3)
24⋅x

Taylor expansion

Taylor expansion with Peano Reminder at x0

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
>>> x = s.symbols('x')
>>> f = s.sin(x)
>>> f.series(x,x0=0,n=3)
     ⎛ 3⎞
x + O⎝x ⎠
>>> f.series(x,x0=0,n=5)
     3
    x     ⎛ 5⎞
x - ── + O⎝x ⎠
    6
>>> f.series(x,x0=0,n=7)
     3     5
    x     x     ⎛ 7⎞
x - ── + ─── + O⎝x ⎠
    6    120
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# For abstract f=f(x)
>>> f = s.Functions('f')
>>> s.series(f(x) , x , x0=s.symbols('x0') , n = 3)
                                                ⎛  2        ⎞│
                                              2 ⎜ d         ⎟│
                                      (x - x₀) ⋅⎜────(f(ξ₁))⎟│
                                                ⎜   2       ⎟│
                 ⎛ d        ⎞│                  ⎝dξ₁        ⎠│ξ₁=x₀    ⎛        3        ⎞
f(x₀) + (x - x₀)⋅⎜───(f(ξ₁))⎟│      + ───────────────────────────── + O⎝(x - x₀) ; x → x₀⎠
                 ⎝dξ₁       ⎠│ξ₁=x₀                 2
# get rid of O(x)
>>> s.series(f(x) , x , x0=s.symbols('x0') , n = 3).removeO()
          ⎛  2        ⎞│
        2 ⎜ d         ⎟│
(x - x₀) ⋅⎜────(f(ξ₁))⎟│
          ⎜   2       ⎟│
          ⎝dξ₁        ⎠│ξ₁=x₀            ⎛ d        ⎞│
───────────────────────────── + (x - x₀)⋅⎜───(f(ξ₁))⎟│      + f(x₀)
              2                          ⎝dξ₁       ⎠│ξ₁=x₀