Tutorial - SymPy in 10 minutes

Introduction

In this tutorial we will learn the basics of Jupyter and SymPy. SymPy is a Python library for symbolic computation. It provides computer algebra capabilities either as a standalone application, as a library to other applications, or live on the web as SymPy Live or SymPy Gamma. Sympy is similar to other CAS (Computer Algebra Software) like Mathematica, Maple or Maxima.

A more complete tutorial can be found at http://docs.sympy.org/latest/tutorial/index.html.

Before using SymPy we should load it, like any other Python libary. We will use

init_session()

to make some imports, this will help us in its interactive use.

For scripting it would be better to do the imports diffferently, for example

import sympy as sym

and then call the functions from SymPy in the following manner

x = sym.Symbols("x")
expr = sym.cos(x)**2 + 3*x
deriv = expr.diff(x)

where we computed the derivative of $cos^2(x) + 3x$, that should be $-2\sin(x)\cos(x) + 3$.

For further information on the Jupyter Notebook you can refer to the User Manual.



In [1]:

    
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from sympy import *



In [2]:

    
init_session()
plt.style.use("seaborn-notebook")









    



IPython console for SymPy 1.2 (Python 3.6.5-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.2/

Let us start with some simple calculations. Belowe we have a code cell with an addition. Place the cursor on it and press SHIFT + ENTER to evaluate it.



In [3]:

    
1 + 3









    Out[3]:





$$4$$

Let us make some calculations



In [4]:

    
factorial(5)









    Out[4]:





$$120$$



In [5]:

    
Out[4]* 10









    Out[5]:





$$1200$$



In [6]:

    
1 // 3









    Out[6]:





$$0$$



In [7]:

    
1 / 3









    Out[7]:





$$0.3333333333333333$$



In [8]:

    
S(1) / 3









    Out[8]:





$$\frac{1}{3}$$

We can evaluate an expression to its floating point version



In [9]:

    
sqrt(2*pi)









    Out[9]:





$$\sqrt{2} \sqrt{\pi}$$



In [11]:

    
float(Out[9])









    Out[11]:





$$2.5066282746310007$$

In the previous line we used the expression Out[8] that stores the output from the evaluation in cell 8 (In[8]). We can also store expressions as variables, just as any Python variable.



In [12]:

    
radius = 10
height = 100
area = pi * radius**2
volume = area * height



In [13]:

    
volume









    Out[13]:





$$10000 \pi$$



In [14]:

    
float(volume)









    Out[14]:





$$31415.926535897932$$

So far, we have just used SymPy as a calculator. Let us try more advanced calculations

Definite and indefinite integrals



In [15]:

    
integrate(sin(x), x)









    Out[15]:





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



In [16]:

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









    Out[16]:





$$2$$

We can define a function and integrate it



In [17]:

    
f = lambda x: x**2 + 5



In [18]:

    
f(5)









    Out[18]:





$$30$$



In [19]:

    
integrate(f(z), z)









    Out[19]:





$$\frac{z^{3}}{3} + 5 z$$



In [20]:

    
integrate(1/(x**2 + y), x)









    Out[20]:





$$- \frac{\sqrt{- \frac{1}{y}} \log{\left (x - y \sqrt{- \frac{1}{y}} \right )}}{2} + \frac{\sqrt{- \frac{1}{y}} \log{\left (x + y \sqrt{- \frac{1}{y}} \right )}}{2}$$

If we assume that the denominator is positive, the expression can be factorized further



In [21]:

    
a = symbols("a", positive=True)
integrate(1/(x**2 + a), x)









    Out[21]:





$$\frac{\operatorname{atan}{\left (\frac{x}{\sqrt{a}} \right )}}{\sqrt{a}}$$

We just learnt the basics, we can try some examples now.

Note: If you want to know more about a specific function you can use help() or the IPython magic command ??



In [22]:

    
help(integrate)









    



Help on function integrate in module sympy.integrals.integrals:

integrate(*args, **kwargs)
    integrate(f, var, ...)
    
    Compute definite or indefinite integral of one or more variables
    using Risch-Norman algorithm and table lookup. This procedure is
    able to handle elementary algebraic and transcendental functions
    and also a huge class of special functions, including Airy,
    Bessel, Whittaker and Lambert.
    
    var can be:
    
    - a symbol                   -- indefinite integration
    - a tuple (symbol, a)        -- indefinite integration with result
                                    given with `a` replacing `symbol`
    - a tuple (symbol, a, b)     -- definite integration
    
    Several variables can be specified, in which case the result is
    multiple integration. (If var is omitted and the integrand is
    univariate, the indefinite integral in that variable will be performed.)
    
    Indefinite integrals are returned without terms that are independent
    of the integration variables. (see examples)
    
    Definite improper integrals often entail delicate convergence
    conditions. Pass conds='piecewise', 'separate' or 'none' to have
    these returned, respectively, as a Piecewise function, as a separate
    result (i.e. result will be a tuple), or not at all (default is
    'piecewise').
    
    **Strategy**
    
    SymPy uses various approaches to definite integration. One method is to
    find an antiderivative for the integrand, and then use the fundamental
    theorem of calculus. Various functions are implemented to integrate
    polynomial, rational and trigonometric functions, and integrands
    containing DiracDelta terms.
    
    SymPy also implements the part of the Risch algorithm, which is a decision
    procedure for integrating elementary functions, i.e., the algorithm can
    either find an elementary antiderivative, or prove that one does not
    exist.  There is also a (very successful, albeit somewhat slow) general
    implementation of the heuristic Risch algorithm.  This algorithm will
    eventually be phased out as more of the full Risch algorithm is
    implemented. See the docstring of Integral._eval_integral() for more
    details on computing the antiderivative using algebraic methods.
    
    The option risch=True can be used to use only the (full) Risch algorithm.
    This is useful if you want to know if an elementary function has an
    elementary antiderivative.  If the indefinite Integral returned by this
    function is an instance of NonElementaryIntegral, that means that the
    Risch algorithm has proven that integral to be non-elementary.  Note that
    by default, additional methods (such as the Meijer G method outlined
    below) are tried on these integrals, as they may be expressible in terms
    of special functions, so if you only care about elementary answers, use
    risch=True.  Also note that an unevaluated Integral returned by this
    function is not necessarily a NonElementaryIntegral, even with risch=True,
    as it may just be an indication that the particular part of the Risch
    algorithm needed to integrate that function is not yet implemented.
    
    Another family of strategies comes from re-writing the integrand in
    terms of so-called Meijer G-functions. Indefinite integrals of a
    single G-function can always be computed, and the definite integral
    of a product of two G-functions can be computed from zero to
    infinity. Various strategies are implemented to rewrite integrands
    as G-functions, and use this information to compute integrals (see
    the ``meijerint`` module).
    
    The option manual=True can be used to use only an algorithm that tries
    to mimic integration by hand. This algorithm does not handle as many
    integrands as the other algorithms implemented but may return results in
    a more familiar form. The ``manualintegrate`` module has functions that
    return the steps used (see the module docstring for more information).
    
    In general, the algebraic methods work best for computing
    antiderivatives of (possibly complicated) combinations of elementary
    functions. The G-function methods work best for computing definite
    integrals from zero to infinity of moderately complicated
    combinations of special functions, or indefinite integrals of very
    simple combinations of special functions.
    
    The strategy employed by the integration code is as follows:
    
    - If computing a definite integral, and both limits are real,
      and at least one limit is +- oo, try the G-function method of
      definite integration first.
    
    - Try to find an antiderivative, using all available methods, ordered
      by performance (that is try fastest method first, slowest last; in
      particular polynomial integration is tried first, Meijer
      G-functions second to last, and heuristic Risch last).
    
    - If still not successful, try G-functions irrespective of the
      limits.
    
    The option meijerg=True, False, None can be used to, respectively:
    always use G-function methods and no others, never use G-function
    methods, or use all available methods (in order as described above).
    It defaults to None.
    
    Examples
    ========
    
    >>> from sympy import integrate, log, exp, oo
    >>> from sympy.abc import a, x, y
    
    >>> integrate(x*y, x)
    x**2*y/2
    
    >>> integrate(log(x), x)
    x*log(x) - x
    
    >>> integrate(log(x), (x, 1, a))
    a*log(a) - a + 1
    
    >>> integrate(x)
    x**2/2
    
    Terms that are independent of x are dropped by indefinite integration:
    
    >>> from sympy import sqrt
    >>> integrate(sqrt(1 + x), (x, 0, x))
    2*(x + 1)**(3/2)/3 - 2/3
    >>> integrate(sqrt(1 + x), x)
    2*(x + 1)**(3/2)/3
    
    >>> integrate(x*y)
    Traceback (most recent call last):
    ...
    ValueError: specify integration variables to integrate x*y
    
    Note that ``integrate(x)`` syntax is meant only for convenience
    in interactive sessions and should be avoided in library code.
    
    >>> integrate(x**a*exp(-x), (x, 0, oo)) # same as conds='piecewise'
    Piecewise((gamma(a + 1), -re(a) < 1),
        (Integral(x**a*exp(-x), (x, 0, oo)), True))
    
    >>> integrate(x**a*exp(-x), (x, 0, oo), conds='none')
    gamma(a + 1)
    
    >>> integrate(x**a*exp(-x), (x, 0, oo), conds='separate')
    (gamma(a + 1), -re(a) < 1)
    
    See Also
    ========
    
    Integral, Integral.doit



In [23]:

    
integrate??

Examples

Solving algebraic equations

Right now, there are two main options for solution of algebraic systems, namely: solveset and solve. The preferred method is solveset (see this explanation, although there are systems that can be solved using solve and not solveset.

To solve equations using solveset:



In [24]:

    
a, b, c = symbols("a b c")
solveset(a*x**2 + b*x + c, x)









    Out[24]:





$$\left\{- \frac{b}{2 a} - \frac{\sqrt{- 4 a c + b^{2}}}{2 a}, - \frac{b}{2 a} + \frac{\sqrt{- 4 a c + b^{2}}}{2 a}\right\}$$

We should enter the equations as equated to 0, or as an equation



In [25]:

    
solveset(Eq(a*x**2 + b*x, -c), x)









    Out[25]:





$$\left\{- \frac{b}{2 a} - \frac{\sqrt{- 4 a c + b^{2}}}{2 a}, - \frac{b}{2 a} + \frac{\sqrt{- 4 a c + b^{2}}}{2 a}\right\}$$

Currently solveset is not capable of solving the following types of equations:

Non-linear multivariate system
Equations solvable by LambertW (Transcendental equation solver).

solve can be used for such cases:



In [26]:

    
solve([x*y - 1, x - 2], x, y)









    Out[26]:





$$\left [ \left ( 2, \quad \frac{1}{2}\right )\right ]$$



In [27]:

    
solve(x*exp(x) - 1, x )









    Out[27]:





$$\left [ \operatorname{LambertW}{\left (1 \right )}\right ]$$

Linear Algebra

We use Matrix to create matrices. Matrices can contain expressions. And we use the method .inv() to compute the inverse and * to multiply matrices.



In [28]:

    
A = Matrix([
        [1, -1],
        [1, sin(c)]
    ])
display(A)









    





$$\left[\begin{matrix}1 & -1\\1 & \sin{\left (c \right )}\end{matrix}\right]$$



In [29]:

    
B = A.inv()
display(B)









    





$$\left[\begin{matrix}\frac{\sin{\left (c \right )}}{\sin{\left (c \right )} + 1} & \frac{1}{\sin{\left (c \right )} + 1}\\- \frac{1}{\sin{\left (c \right )} + 1} & \frac{1}{\sin{\left (c \right )} + 1}\end{matrix}\right]$$



In [30]:

    
A * B









    Out[30]:





$$\left[\begin{matrix}\frac{\sin{\left (c \right )}}{\sin{\left (c \right )} + 1} + \frac{1}{\sin{\left (c \right )} + 1} & 0\\0 & \frac{\sin{\left (c \right )}}{\sin{\left (c \right )} + 1} + \frac{1}{\sin{\left (c \right )} + 1}\end{matrix}\right]$$

This should be the identity matrix, let us simplify the expression. There are several simplification functions, and simplify is the most general one. Simplifying is a complicated matter... if you are uncertain; use simplify.



In [31]:

    
simplify(A * B)









    Out[31]:





$$\left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$$

Plotting

We can make 2D and 3D plots



In [32]:

    
%matplotlib inline
from sympy.plotting import plot3d



In [33]:

    
plot(sin(x), (x, -pi, pi));



In [34]:

    
monkey_saddle = x**3 - 3*x*y**2
p = plot3d(monkey_saddle, (x, -2, 2), (y, -2, 2))

Derivatives and differential equations

We can use the function diff or the method .diff() to compute derivatives.



In [35]:

    
f = lambda x: x**2



In [36]:

    
diff(f(x), x)









    Out[36]:





$$2 x$$



In [37]:

    
f(x).diff(x)









    Out[37]:





$$2 x$$



In [38]:

    
g = lambda x: sin(x)



In [39]:

    
diff(g(f(x)), x)









    Out[39]:





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

Yes, SymPy knows about the chain rule.

To finish, let us solve a second order ODE

$$ u''(t) + \omega^2 u(t) = 0$$



In [40]:

    
u = symbols("u", cls=Function)
omega = symbols("omega", positive=True)



In [41]:

    
ode = u(t).diff(t, 2) + omega**2 * u(t)
dsolve(ode, u(t))









    Out[41]:





$$u{\left (t \right )} = C_{1} \sin{\left (\omega t \right )} + C_{2} \cos{\left (\omega t \right )}$$

Turning Sympy expression into evaluable functions

lambdify provides convenient functions to transform sympy expressions to lambda functions which can be used to calculate numerical values very fast.

Let's try a first example



In [42]:

    
f = lambdify(x, x**2, "numpy")
f(3)









    Out[42]:





$$9$$



In [43]:

    
f(np.array([1, 2, 3]))









    Out[43]:





array([1, 4, 9])

We can try a more difficult one



In [44]:

    
fun = diff(sin(x)*cos(x**3) - sin(x)/x, x)
fun









    Out[44]:





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



In [45]:

    
fun_numpy = lambdify(x, fun, "numpy")

and then evaluate it for some points in, let's say, $[0, 5]$



In [46]:

    
pts = np.linspace(0, 5, 1000)
fun_pts = fun_numpy(pts + 1e-6) # To avoid division by 0



In [47]:

    
plt.plot(pts, fun_pts);

References

Žiga Lenarčič. "10 minute (wx)Maxima tutorial:", (2008). Accessed on: July 25, 2017.
Borland, David, and Russell M. Taylor II. "Rainbow color map (still) considered harmful." IEEE computer graphics and applications 27.2 (2007): 14-17.

The following cell change the style of the notebook.



In [49]:

    
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()









    Out[49]:



In [ ]: