Introduction to SymPy

This is to present a short introduction to symbolic mathematics in Python using SymPy

Following the flow of : http://www.cfm.brown.edu/people/dobrush/am33/SymPy/index.html



In [10]:

    
import sympy as sym
init_printing()  # to get nice \LaTeX output

Symbols

There are two ways in which a mathematical symbol can be defined:



In [11]:

    
x = sym.Symbol('x')



In [12]:

    
x+1









    Out[12]:





$$x + 1$$

For multiple symbols one can do



In [13]:

    
y, z = sym.symbols('y z')
y+z+1









    Out[13]:





$$y + z + 1$$



In [14]:

    
x='abc'



In [16]:

    
expr = x + 'def' ; expr









    Out[16]:





'abcdef'

Substitute



In [17]:

    
x = sym.symbols('x')
expr = x+1
expr









    Out[17]:





$$x + 1$$



In [20]:

    
expr.subs(x,1) # should be 2









    Out[20]:





$$2$$



In [21]:

    
expr.subs(x,z)









    Out[21]:





$$z + 1$$

The importance of substition is the following

We can simply evaluate an expression at a point. e.g. for the expression cos(x)+1 we want to know its value when x=0:



In [22]:

    
expr = cos(x)+1



In [23]:

    
expr.subs(x,0)









    Out[23]:





$$2$$

But also is very important that we can substitute a variable with a new expression! This way we can build sequences...e.g. $x^{x^{x^{x}}}$



In [24]:

    
expr.subs(x, x+y)









    Out[24]:





$$\cos{\left (x + y \right )} + 1$$



In [27]:

    
expr.subs(x, x**y)









    Out[27]:





$$\cos{\left (x^{y} \right )} + 1$$

N.B. subs() returns a new expression, i.e. the expr is not changed!

subs is immutable! To perform multiple substitution we use a list of tuples:



In [28]:

    
expr = cos(x)+2*sin(y)+z



In [31]:

    
expr.subs([(x,1),(y,3), (z,0)])









    Out[31]:





$$2 \sin{\left (3 \right )} + \cos{\left (1 \right )}$$

This is very usefull because we can use list comprehension for the substitutions!



In [33]:

    
expr = x**4 - 4*x**3 + 4*x**2 - 2*x + 3
expr









    Out[33]:





$$x^{4} - 4 x^{3} + 4 x^{2} - 2 x + 3$$



In [36]:

    
replacements = [(x**i, y**i) for i in range(5) if i%2==0]  #replace all instances of x that 
                                                           # have even powers with y



In [37]:

    
expr.subs(replacements)









    Out[37]:





$$- 4 x^{3} - 2 x + y^{4} + 4 y^{2} + 3$$

Strings -> Sympy Expressions

We can use the sympify() ( NOT SIMPLIFY ) to convert string to sympy expressions



In [40]:

    
str_expr = "x**2+3*x+2"



In [42]:

    
expr = sym.sympify(str_expr)



In [44]:

    
expr









    Out[44]:





$$x^{2} + 3 x + 2$$



In [45]:

    
expr.subs(x,2)









    Out[45]:





$$12$$

evalf()

Use evalf() to evaluate a numerical expression into a floating point number:



In [46]:

    
expr= sqrt(8)



In [47]:

    
expr









    Out[47]:





$$2 \sqrt{2}$$



In [51]:

    
expr.evalf()









    Out[51]:





$$2.82842712474619$$



In [53]:

    
expr.evalf(50) #first 50 digits









    Out[53]:





$$2.8284271247461900976033774484193961571393437507539$$

To evaluate an expression with a symbol to a point we can use the subs flag in evalf() using a dictionary : Symbol: point



In [55]:

    
expr = cos(2*x)
expr









    Out[55]:





$$\cos{\left (2 x \right )}$$



In [56]:

    
expr.evalf(subs={x:0.2})









    Out[56]:





$$0.921060994002885$$



In [58]:

    
cos(2*0.2) #it works









    Out[58]:





$$0.921060994002885$$

There are roundoff errors smaller than the desired precision. Such numbers can be cut using the chop flag in eval():



In [59]:

    
one = cos(1)**2 + sin(1)**2



In [60]:

    
(one-1).evalf()









    Out[60]:





$$-4.0 \cdot 10^{-124}$$



In [61]:

    
(one-1).evalf(chop=True)









    Out[61]:





$$0$$

Lambdify

lambdify() is used to convert a SymPy expression to another library (usually NumPy). lambdify() acts as a lambda function but also converting the SymPy expressions to a another numerical library.

For example if we want to use the subs flag of the evalf() function to evaluate the function at 1000 points, then we can use numpy and lambdify...



In [62]:

    
import numpy as np



In [68]:

    
a = np.arange(100)



In [69]:

    
expr = sin(x)



In [74]:

    
f = lambdify(x, expr, "numpy")  # takes an argument x=x0 and returns the expr(x0)



In [75]:

    
f(a)









    Out[75]:





array([ 0.        ,  0.84147098,  0.90929743,  0.14112001, -0.7568025 ,
       -0.95892427, -0.2794155 ,  0.6569866 ,  0.98935825,  0.41211849,
       -0.54402111, -0.99999021, -0.53657292,  0.42016704,  0.99060736,
        0.65028784, -0.28790332, -0.96139749, -0.75098725,  0.14987721,
        0.91294525,  0.83665564, -0.00885131, -0.8462204 , -0.90557836,
       -0.13235175,  0.76255845,  0.95637593,  0.27090579, -0.66363388,
       -0.98803162, -0.40403765,  0.55142668,  0.99991186,  0.52908269,
       -0.42818267, -0.99177885, -0.64353813,  0.29636858,  0.96379539,
        0.74511316, -0.15862267, -0.91652155, -0.83177474,  0.01770193,
        0.85090352,  0.90178835,  0.12357312, -0.76825466, -0.95375265,
       -0.26237485,  0.67022918,  0.98662759,  0.39592515, -0.55878905,
       -0.99975517, -0.521551  ,  0.43616476,  0.99287265,  0.63673801,
       -0.30481062, -0.96611777, -0.7391807 ,  0.1673557 ,  0.92002604,
        0.82682868, -0.02655115, -0.85551998, -0.89792768, -0.11478481,
        0.77389068,  0.95105465,  0.25382336, -0.67677196, -0.98514626,
       -0.38778164,  0.56610764,  0.99952016,  0.51397846, -0.44411267,
       -0.99388865, -0.62988799,  0.31322878,  0.96836446,  0.73319032,
       -0.17607562, -0.92345845, -0.82181784,  0.0353983 ,  0.86006941,
        0.89399666,  0.10598751, -0.77946607, -0.94828214, -0.24525199,
        0.68326171,  0.98358775,  0.37960774, -0.57338187, -0.99920683])

Simplification

Functions that are used to simplify expressions

simplify()



In [77]:

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









    Out[77]:





$$1$$



In [79]:

    
expr = (x**3 + x**2 - x - 1)/(x**2 + 2*x + 1);
expr









    Out[79]:





$$\frac{x^{3} + x^{2} - x - 1}{x^{2} + 2 x + 1}$$



In [80]:

    
simplify(expr)









    Out[80]:





$$x - 1$$



In [83]:

    
expr = gamma(x)/gamma(x-2); expr   ## gamma(x) := Gamma functions









    Out[83]:





$$\frac{\Gamma{\left(x \right)}}{\Gamma{\left(x - 2 \right)}}$$



In [84]:

    
simplify (expr)









    Out[84]:





$$\left(x - 2\right) \left(x - 1\right)$$

factor() & factor_list()

factor is used on polynomials with rational coefficients and factors it inot irreducible factors over the rational numbers



In [88]:

    
expr = x**3 - x**2 + x- 1;
expr









    Out[88]:





$$x^{3} - x^{2} + x - 1$$



In [89]:

    
factor(expr)









    Out[89]:





$$\left(x - 1\right) \left(x^{2} + 1\right)$$



In [90]:

    
expr = x**2*z + 4*x*y*z + 4*y**2*z; expr









    Out[90]:





$$x^{2} z + 4 x y z + 4 y^{2} z$$



In [91]:

    
factor (expr)









    Out[91]:





$$z \left(x + 2 y\right)^{2}$$



In [93]:

    
factor_list(expr)









    Out[93]:





$$\left ( 1, \quad \left [ \left ( z, \quad 1\right ), \quad \left ( x + 2 y, \quad 2\right )\right ]\right )$$

expand()

expand() is used to make an expression bigger. The opposite of factor().



In [96]:

    
expr = (x+1)**2;
expr









    Out[96]:





$$\left(x + 1\right)^{2}$$



In [97]:

    
expand(expr)









    Out[97]:





$$x^{2} + 2 x + 1$$



In [ ]:

collect() and coeff()

collect() is used to collect common powers of a term



In [98]:

    
expr = x*y +x-3+2*x**2 - z*x**2 + x**3; expr









    Out[98]:





$$x^{3} - x^{2} z + 2 x^{2} + x y + x - 3$$



In [99]:

    
cl = collect(expr, x)



In [100]:

    
cl









    Out[100]:





$$x^{3} + x^{2} \left(- z + 2\right) + x \left(y + 1\right) - 3$$

Very usefull when combined with the coeff() function that fives the coefficients of a symbol in the expression:



In [101]:

    
expr









    Out[101]:





$$x^{3} - x^{2} z + 2 x^{2} + x y + x - 3$$



In [102]:

    
expr.coeff(x,2)









    Out[102]:





$$- z + 2$$



In [103]:

    
expr.coeff(y, 1)









    Out[103]:





$$x$$



In [105]:

    
cl.coeff(x,3)









    Out[105]:





$$1$$



In [106]:

    
cl.coeff(x,2)









    Out[106]:





$$- z + 2$$

cancel()

cancel() is used to take any rational function and put it into standard canonical form.



In [107]:

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









    Out[107]:





$$\frac{1}{x} \left(x + 1\right)$$



In [108]:

    
expr = 1/x + (3*x/2 - 2)/(x - 4)



In [109]:

    
expr









    Out[109]:





$$\frac{\frac{3 x}{2} - 2}{x - 4} + \frac{1}{x}$$



In [110]:

    
cancel(expr)









    Out[110]:





$$\frac{3 x^{2} - 2 x - 8}{2 x^{2} - 8 x}$$

apart() - Partial Fraction Decomposition

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



In [112]:

    
expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x); expr









    Out[112]:





$$\frac{4 x^{3} + 21 x^{2} + 10 x + 12}{x^{4} + 5 x^{3} + 5 x^{2} + 4 x}$$



In [113]:

    
apart(expr)









    Out[113]:





$$\frac{2 x - 1}{x^{2} + x + 1} - \frac{1}{x + 4} + \frac{3}{x}$$

Trigonometric Simplification

N.B. Remember that the Python's convention for inverse trig equations is acos(), asin() etc..



In [114]:

    
acos(x)









    Out[114]:





$$\operatorname{acos}{\left (x \right )}$$



In [115]:

    
cos(acos(x))









    Out[115]:





$$x$$

trigsimp()

trigsimp() is used to simplify expressions based on trigonometric identities.



In [116]:

    
expr = sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4; expr









    Out[116]:





$$\sin^{4}{\left (x \right )} - 2 \sin^{2}{\left (x \right )} \cos^{2}{\left (x \right )} + \cos^{4}{\left (x \right )}$$



In [117]:

    
trigsimp(expr)









    Out[117]:





$$\frac{1}{2} \cos{\left (4 x \right )} + \frac{1}{2}$$



In [119]:

    
expr = sin(x)*tan(x)/sec(x); expr









    Out[119]:





$$\frac{\sin{\left (x \right )} \tan{\left (x \right )}}{\sec{\left (x \right )}}$$



In [120]:

    
trigsimp(expr)









    Out[120]:





$$\sin^{2}{\left (x \right )}$$



In [121]:

    
expr = cosh(x)**2 + sinh(x)**2 ; expr









    Out[121]:





$$\sinh^{2}{\left (x \right )} + \cosh^{2}{\left (x \right )}$$



In [122]:

    
trigsimp(expr)









    Out[122]:





$$\cosh{\left (2 x \right )}$$

expand_trig()

expand_trig() is used to expand trigonometric functions (apply the sum or double angle identities)



In [123]:

    
expr = sin(x+y); expr









    Out[123]:





$$\sin{\left (x + y \right )}$$



In [124]:

    
expand_trig(expr)









    Out[124]:





$$\sin{\left (x \right )} \cos{\left (y \right )} + \sin{\left (y \right )} \cos{\left (x \right )}$$



In [127]:

    
expr = tan(2*x); expr









    Out[127]:





$$\tan{\left (2 x \right )}$$



In [128]:

    
expand_trig(expr)









    Out[128]:





$$\frac{2 \tan{\left (x \right )}}{- \tan^{2}{\left (x \right )} + 1}$$

Powers

There are identities for powers that hold under constraints. For example:

$x^a x^b = x^{a+b}$ -> holds always
x^a y^a = (xy)^a -> holds when $x,y>0$ and $a\in R$
(x^a)^b = x^{ab} -> holds when b is integer

Symbols and requirements

For this reason when defining symbols we can put specific requirements as their arguments. For example we define for this section:



In [129]:

    
x, y = symbols('x y', positive=True)



In [130]:

    
a,b = symbols('a b', real=True)



In [131]:

    
z, t, c = symbols('z t c ')

powsimp()

powsimp() is used to simplify an expression with powers (using identities 1,2)



In [133]:

    
expr = (x**a)*(x**b); expr









    Out[133]:





$$x^{a} x^{b}$$



In [134]:

    
powsimp(expr)









    Out[134]:





$$x^{a + b}$$

To force a simplification without messing with the assumptions, while



In [136]:

    
powsimp(t**c*z**c)









    Out[136]:





$$t^{c} z^{c}$$

forcing it with an argument



In [137]:

    
powsimp(t**c*z**c, force=True)









    Out[137]:





$$\left(t z\right)^{c}$$

expand_power_exp() and expand_power_base()

expand_power_exp() and expand_power_base() applies the identities 1,2 from right to left, expanding the exponent or the base of an expression



In [138]:

    
expr = x**(a+b); expr









    Out[138]:





$$x^{a + b}$$



In [139]:

    
expand_power_exp(expr)









    Out[139]:





$$x^{a} x^{b}$$



In [140]:

    
expr = (x*y)**a; expr









    Out[140]:





$$\left(x y\right)^{a}$$



In [141]:

    
expand_power_base(expr)









    Out[141]:





$$x^{a} y^{a}$$

powdenest()

This is to apply identity 3 on nested powers



In [144]:

    
expr = (x**a)**b; expr  # automatically is simplified









    Out[144]:





$$x^{a b}$$



In [143]:

    
powdenest(expr)









    Out[143]:





$$x^{a b}$$