11. sympy-examples

Symbolic Computation with SymPy¶

This Notebook closely follows the SymPy tutorial.

In [2]:

import sympy as sym

import math

Basic Symbolic Manipulation¶

In [3]:

print(math.sqrt(2)**2)

2.0000000000000004

In [4]:

print(sym.sqrt(2)**2)

2

In [5]:

print(sym.sqrt(8)**2)

8

We can do symbolic math not just on numbers, but we can tell SymPy what to treat as a variable.

In [6]:

x, y, z = sym.symbols("x y z")

In [7]:

expr = x + 2*y
expr

Out[7]:

$\displaystyle x + 2 y$

In [8]:

expr - 1

Out[8]:

$\displaystyle x + 2 y - 1$

In [9]:

expr - y

Out[9]:

$\displaystyle x + y$

In [10]:

f = x*expr
f

Out[10]:

$\displaystyle x \left(x + 2 y\right)$

In [11]:

g = sym.expand(f)
g

Out[11]:

$\displaystyle x^{2} + 2 x y$

In [12]:

sym.simplify(g**2)

Out[12]:

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

Substitution¶

In [13]:

expr = sym.sin(x*2*sym.pi)
expr

Out[13]:

$\displaystyle \sin{\left(2 \pi x \right)}$

In [14]:

expr.subs(x,1)

Out[14]:

$\displaystyle 0$

In [15]:

a = expr.subs(x,0.125)
a

Out[15]:

$\displaystyle \frac{\sqrt{2}}{2}$

In [16]:

type(a)

Out[16]:

sympy.core.mul.Mul

Note that this is not a floating point number -- it is still a SymPy object. To make it floating point, we can use evalf().

In [17]:

b = a.evalf()
print(b, type(b))

0.707106781186548 <class 'sympy.core.numbers.Float'>

This is still a SymPy object, because SymPy can do arbitrary precision.

In [18]:

sym.pi.evalf(100)

Out[18]:

$\displaystyle 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$

Want regular python types?

In [19]:

c = float(b)
print(c, type(c))

0.7071067811865476 <class 'float'>

Using SymPy Interactively¶

For convenience, we will import everything from sympy and use it interatively in a calculator-like environment.

In [20]:

from sympy import *

The following statements in SymPy tutorials are meant to make output look pretty in Jupyter Notebook. However, they are not needed for the latest versions of iPython and Jupyter Notebook, so we omit them.

from sympy.interactive import printing
printing.init_printing()

Python and SymPy¶

SymPy symbols are just objects and when you do operations on two sympy objects the result is a sympy object. When you combine a sympy and python object, the result is also a sympy object. But now we need to be careful when doing fractions. For instance doing x + 1/3 will first compute 1/3 in python (giving 0 or 0.333... depending on the division operator) and then add it to the sympy x symbol. The Rational() function makes this all happen in sympy.

In [21]:

f = expr + Rational(1,3)
f

Out[21]:

$\displaystyle \sin{\left(2 \pi x \right)} + \frac{1}{3}$

In [22]:

expr + 1/3

Out[22]:

$\displaystyle \sin{\left(2 \pi x \right)} + 0.333333333333333$

Equality¶

= is still the assignment operator of python (it does not mean symbolic equality), and == is still the logical test (exact structural equality). There is a separate object, Eq() to specify symbolic equality

In [23]:

x + 1 == 4

Out[23]:

False

In [24]:

Eq(x + 1, 4)

Out[24]:

$\displaystyle x + 1 = 4$

In [25]:

solve(Eq(x+1, 4))

Out[25]:

[3]

In [26]:

a = (x + 1)**2
b = x**2 + 2*x + 1    # these are algebraically equal

In [27]:

a == b

Out[27]:

False

In [28]:

simplify(a - b)   # this will test equality algebraically

Out[28]:

$\displaystyle 0$

In [29]:

a = cos(x) + sym.I*sin(x)
a

Out[29]:

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

In [30]:

simplify(a)

Out[30]:

$\displaystyle e^{i x}$

More Substitution¶

Note that substitution returns a new expression: SymPy objects are immutable.

In [31]:

expr = cos(x)
expr.subs(x, 0)

Out[31]:

$\displaystyle 1$

In [32]:

expr

Out[32]:

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

In [33]:

x

Out[33]:

$\displaystyle x$

Multiple substitutions by passing a list of tuples.

In [34]:

expr = x**3 + 4*x*y - z
expr

Out[34]:

$\displaystyle x^{3} + 4 x y - z$

In [35]:

expr.subs([(x, 2), (y, 4)])

Out[35]:

$\displaystyle 40 - z$

Simplifying¶

There is not unique definition of what the simplest form of an expression is.

simplify() tries lots of methods for simplification.

In [36]:

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

Out[36]:

$\displaystyle 1$

In [37]:

simplify( (x**3 + x**2 - x - 1)/(x**2 + 2*x + 1) )

Out[37]:

$\displaystyle x - 1$

In [38]:

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

Out[38]:

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

However, sometimes it doesn't have your idea of what the simplest form is.

In [39]:

simplify(x**2 + 2*x + 1)

Out[39]:

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

Instead, factor may be what you want.

In [40]:

factor(x**2 + 2*x + 1)

Out[40]:

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

Polynomial Simplification¶

In [41]:

expand((x + 1)**2)

Out[41]:

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

In [42]:

expand((x + 2)*(x - 3))

Out[42]:

$\displaystyle x^{2} - x - 6$

In [43]:

expand( (x + 1)*(x - 2) - (x - 1)*x)

Out[43]:

$\displaystyle -2$

In [44]:

factor(x**2*z + 4*x*y*z + 4*y**2*z)

Out[44]:

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

In [45]:

factor_list(x**2*z + 4*x*y*z + 4*y**2*z)

Out[45]:

(1, [(z, 1), (x + 2*y, 2)])

collect collects common powers.

In [46]:

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

Out[46]:

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

In [47]:

collected_expr = collect(expr, x)
collected_expr

Out[47]:

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

cancel cancels.

In [48]:

a = (x**2 + 2*x + 1)/(x**2 + x)
a

Out[48]:

$\displaystyle \frac{x^{2} + 2 x + 1}{x^{2} + x}$

In [49]:

cancel(a)

Out[49]:

$\displaystyle \frac{x + 1}{x}$

trigsimp simplifies trigonometric identities.

In [50]:

trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)

Out[50]:

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

In [51]:

trigsimp(sin(x)*tan(x)/sec(x))

Out[51]:

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

The tutorial discusses some of the nuances of simplification of powers and special functions.

Calculus¶

Derivatives.

In [52]:

x = symbols('x')
diff(cos(x), x)

Out[52]:

$\displaystyle - \sin{\left(x \right)}$

In [53]:

diff(exp(x**2), x, 2).subs(x, 5).evalf()

Out[53]:

$\displaystyle 7344499732413.36$

Third derivative:

In [54]:

diff(x**4, x, 3)

Out[54]:

$\displaystyle 24 x$

In [55]:

diff?

Signature: diff(f, *symbols, **kwargs)
Docstring:
Differentiate f with respect to symbols.

Explanation
===========

This is just a wrapper to unify .diff() and the Derivative class; its
interface is similar to that of integrate().  You can use the same
shortcuts for multiple variables as with Derivative.  For example,
diff(f(x), x, x, x) and diff(f(x), x, 3) both return the third derivative
of f(x).

You can pass evaluate=False to get an unevaluated Derivative class.  Note
that if there are 0 symbols (such as diff(f(x), x, 0), then the result will
be the function (the zeroth derivative), even if evaluate=False.

Examples
========

>>> from sympy import sin, cos, Function, diff
>>> from sympy.abc import x, y
>>> f = Function('f')

>>> diff(sin(x), x)
cos(x)
>>> diff(f(x), x, x, x)
Derivative(f(x), (x, 3))
>>> diff(f(x), x, 3)
Derivative(f(x), (x, 3))
>>> diff(sin(x)*cos(y), x, 2, y, 2)
sin(x)*cos(y)

>>> type(diff(sin(x), x))
cos
>>> type(diff(sin(x), x, evaluate=False))
<class 'sympy.core.function.Derivative'>
>>> type(diff(sin(x), x, 0))
sin
>>> type(diff(sin(x), x, 0, evaluate=False))
sin

>>> diff(sin(x))
cos(x)
>>> diff(sin(x*y))
Traceback (most recent call last):
...
ValueError: specify differentiation variables to differentiate sin(x*y)

Note that ``diff(sin(x))`` syntax is meant only for convenience
in interactive sessions and should be avoided in library code.

References
==========

.. [1] http://reference.wolfram.com/legacy/v5_2/Built-inFunctions/AlgebraicComputation/Calculus/D.html

See Also
========

Derivative
idiff: computes the derivative implicitly
File:      c:\users\win10hx64\appdata\local\programs\python\python310\lib\site-packages\sympy\core\function.py
Type:      function

Differentiate different variables:

In [56]:

expr = exp(x*y*z)
diff(expr, y, x, z)

Out[56]:

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

Unevaluated derivatied:

In [57]:

deriv = Derivative(expr, x, y, z)
deriv

Out[57]:

$\displaystyle \frac{\partial^{3}}{\partial z\partial y\partial x} e^{x y z}$

In [58]:

deriv.doit()

Out[58]:

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

Integrals:

In [59]:

integrate(cos(x), x)

Out[59]:

$\displaystyle \sin{\left(x \right)}$

Definite integral:

In [60]:

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

Out[60]:

$\displaystyle 1$

Double integral:

In [61]:

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

Out[61]:

$\displaystyle \pi$

If it is unable to do the integral, it returns an Integral object:

In [62]:

expr = integrate(x**x, x)
print(expr)
expr

Integral(x**x, x)

Out[62]:

$\displaystyle \int x^{x}\, dx$

In [94]:

a = x / sqrt(x**4 + 10*x**2 - 96*x - 71)   # example from Wikipedia Risch algorithm page)
a
f = cos(x*2*pi)
f

Out[94]:

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

In [93]:

integrate(f,(x,0,1))

Out[93]:

$\displaystyle 0$

In [84]:

sym.printing.python(a)

Out[84]:

"x = Symbol('x')\ne = 1 + x + x**2/2 - x**4/8 - x**5/15 - x**6/240 + x**7/90 + 31*x**8/5760 + x**9/5670 + O(x**10)"

In [95]:

integrate(a, (x, -2, -1))  
   # this has a known solution, but SymPy fails to find it

Out[95]:

$\displaystyle \int\limits_{-2}^{-1} \frac{x}{\sqrt{x^{4} + 10 x^{2} - 96 x - 71}}\, dx$

Limits:

In [66]:

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

Out[66]:

$\displaystyle 1$

Series expansions:

In [67]:

expr = exp(sin(x))
a = expr.series(x, 0, 10)    # expansion about x=0 up to 10th order

Note: For SymPy 0.7.5 and earlier, the order term O() is only supported for expansions about 0 or oo. So when you compute the series about a point other than 0 or oo, the result will be shifted to 0. You need to shift back by first removing the O() term.

In [68]:

print(a)

1 + x + x**2/2 - x**4/8 - x**5/15 - x**6/240 + x**7/90 + 31*x**8/5760 + x**9/5670 + O(x**10)

For SymPy 0.7.6, this is not an issue:

In [69]:

c = log(x).series(x, x0=1, n=6)
c

Out[69]:

$\displaystyle -1 - \frac{\left(x - 1\right)^{2}}{2} + \frac{\left(x - 1\right)^{3}}{3} - \frac{\left(x - 1\right)^{4}}{4} + \frac{\left(x - 1\right)^{5}}{5} + x + O\left(\left(x - 1\right)^{6}; x\rightarrow 1\right)$

In [70]:

sym.printing.ccode(simplify(c.removeO()))

Out[70]:

'(1.0/5.0)*pow(x, 5) - 5.0/4.0*pow(x, 4) + (10.0/3.0)*pow(x, 3) - 5*pow(x, 2) + 5*x - 137.0/60.0'

Solvers¶

If no Eq() is done, then it is assumed to be equal to 0.

In [71]:

solve(Eq(x**2, 1), x)

Out[71]:

[-1, 1]

In [72]:

solve(x**2 - 1, x)

Out[72]:

[-1, 1]

In [73]:

solve([x - y + 2, x + y - 3], [x, y])

Out[73]:

{x: 1/2, y: 5/2}

Roots will report if a solution is multiple by listing it multiple times.

In [74]:

roots(x**3 - 6*x**2 + 9*x, x)

Out[74]:

{3: 2, 0: 1}

0 is 1 root, and 3 is 2 more roots:

Differential Equations:¶

You need an undefined function (f and g already are by our init_session() above, but we've probably reset these:

In [75]:

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

In [76]:

f(x)

Out[76]:

$\displaystyle f{\left(x \right)}$

In [77]:

f(x).diff(x)

Out[77]:

$\displaystyle \frac{d}{d x} f{\left(x \right)}$

In [78]:

diffeq = Eq(f(x).diff(x, 2) - 2*f(x).diff(x) + f(x), sin(x))

In [79]:

diffeq

Out[79]:

$\displaystyle f{\left(x \right)} - 2 \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = \sin{\left(x \right)}$

In [80]:

dsolve(diffeq, f(x))

Out[80]:

$\displaystyle f{\left(x \right)} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{\cos{\left(x \right)}}{2}$