# Tutorial - SymPy in 10 minutes
    
    

## Introduction

In this tutorial we will learn the basics of [Jupyter](http://jupyter.org/) and [SymPy](https://es.wikipedia.org/wiki/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](http://live.sympy.org/) or [SymPy Gamma](http://www.sympygamma.com/). 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](http://docs.sympy.org/latest/tutorial/index.html).

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

```python
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

```python
import sympy as sym
```

and then call the functions from SymPy in the following manner

```python
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](https://athena.brynmawr.edu/jupyter/hub/dblank/public/Jupyter%20Notebook%20Users%20Manual.ipynb).


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


```python
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.


```python
1 + 3
```

Let us make some calculations


```python
factorial(5)
```


```python
Out[4]* 10
```


```python
1 // 3
```


```python
1 / 3
```


```python
S(1) / 3
```

We can evaluate an expression to its floating point version


```python
sqrt(2*pi)
```


```python
float(Out[9])
```

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.


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


```python
volume
```


```python
float(volume)
```

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



Definite and indefinite integrals


```python
integrate(sin(x), x)
```


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

We can define a function and integrate it


```python
f = lambda x: x**2 + 5
```


```python
f(5)
```


```python
integrate(f(z), z)
```


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

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



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

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 ``??``


```python
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
    



```python
integrate??
```

## Examples

### Solving algebraic equations

Right now, there are two main options for solution of algebraic systems, namely: 
[``solveset`` and ``solve``](http://docs.sympy.org/latest/tutorial/solvers.html).
The preferred method is ``solveset`` (see this
[explanation](http://docs.sympy.org/latest/modules/solvers/solveset.html), although
there are systems that can be solved using ``solve`` and not ``solveset``.


To solve equations using ``solveset``:


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

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


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

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:


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


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

### 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.


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


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



```python
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]$$



```python
A * B
```




$$\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``.


```python
simplify(A * B)
```




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



### Plotting

We can make 2D and 3D plots


```python
%matplotlib inline
from sympy.plotting import plot3d
```


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


```python
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.


```python
f = lambda x: x**2
```


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


```python
f(x).diff(x)
```


```python
g = lambda x: sin(x)
```


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

Yes, SymPy knows about the chain rule.

To finish, let us solve a second order ODE

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


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


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

## 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


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


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




    array([1, 4, 9])



We can try a more difficult one


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


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

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


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


```python
plt.plot(pts, fun_pts);
```

## References

- Žiga Lenarčič. ["10 minute (wx)Maxima tutorial:"](https://andrejv.github.io/wxmaxima/tutorials/10minute.zip), (2008). Accessed on: July 25, 2017.
- Borland, David, and Russell M. Taylor II. ["Rainbow color map (still) considered harmful."](https://data3.mprog.nl/course/15%20Readings/40%20Reading%204/Borland_Rainbow_Color_Map.pdf) IEEE computer graphics and applications 27.2 (2007): 14-17.

The following cell change the style of the notebook.


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





<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>
<style>
/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}	
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}	
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}
</style>








```python

```

d49f5d102021494cc8f5e7314c43f115de66821a
ipynb
Jupyter Notebook
notebooks/sympy/sympy_in_10_minutes.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2017-06-29T17:45:20.000Z
2022-02-06T20:14:29.000Z
notebooks/sympy/sympy_in_10_minutes.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
notebooks/sympy/sympy_in_10_minutes.ipynb
nicoguaro/AdvancedMath
2749068de442f67b89d3f57827367193ce61a09c
MIT
2019-04-22T08:08:56.000Z
2022-01-27T08:15:53.000Z
Qwen/Qwen-72B
1. YES
2. YES