# 4. Introduction to Sympy for Mathematicians

## Table of Contents
- __Introduction to Sympy__
- __Plots__
- __Polynomials__
- Series
    - Taylor Series
    - Summing Series
- Linear Algebra
- Differential Equations
- __Limits, Derivatives, Integrals__
    - __Limits__
    - __Derivatives__
    - __Integrals__

Next week? Choose among:

- Stochastics/statistics: `random`, `scipy.statistics` , `statsmodels`
- Machine Learning 2/Deep Learning: `scikit-learn`, `keras`, `tensorflow`
- SAGE
- Other: ___________?



## 4.1 Introduction to SymPy

`sympy` is a package for symbolic calculations in Python, similar to *Wolfram Alpha/Mathematica*. It works with expressions containing symbols. Its core symbolic capabilities include:

- Elementary algebra
- Pattern matching and substitution
- Symbolic transformations and simplifications
- All the common mathematical functions and many special ones
- Calculus: limits, derivatives, integrals, Taylor series, and so on
- Solving many kinds of equations, including differential equations
- Linear algebra
- Polynomials
- Infinite sums and products
- Planar geometry
- Physical constants and units
- Numeric computation in arbitrary precision

Today we will import everything contained in the package:


```python
# import sympy as sym
from sympy import *
init_printing()
```

SymPy is also available at http://live.sympy.org, an online Python shell, with `sympy` already imported, which
allows you to perform simple calculations. The full documentation of this package can be found here: https://docs.sympy.org/latest/index.html. We can also watch a YouTube video:


```python
from IPython.display import YouTubeVideo
YouTubeVideo('kx2GzBeGPco', width=533, height=300)
```










Symbols are basic bricks used to construct expressions.
Each symbol has a name used for printing expressions. Objects of the class `Symbol` should be created and assigned to Python variables in order to be used in expressions. The symbol name and the name of the variable to which this symbol is assigned are two independent things, and one may write `abc=Symbol('xyz')`. But then one has to write `abc` in input expressions, while `SymPy` will write `xyz` in output ones, producing unnecessary confusion. The Python variable name should better be the same as the symbol name.

In languages specifically designed for symbolic calculations, such as *Mathematica*, if a variable to which nothing has been assigned is used, it automatically means a symbol with the same name. Python has not been designed for symbolic calculations. If you use a variable to which nothing has been assigned, you will get an error message. Symbol objects have to be created explicitly.


```python
x = Symbol('x')
x
```


```python
a = x**2-1
a
```

However, we can also use the command `sympify`, then we don't need to define $x$ manually:


```python
sympify('x^2 - 1')
```




$\displaystyle x^{2} - 1$



Internally expressions are trees. The function `srepr` returns a string representing this tree.


```python
srepr(a)
```




    "Add(Pow(Symbol('x'), Integer(2)), Integer(-1))"



One may use the functions `Add`, `Mul`, `Pow`, etc. instead of the binary operations `+`, `*`, `**`, etc.


```python
Add(Pow(Symbol('x'), Integer(2)), Integer(-1))
```

Some packages are available that provide Greek letters:


```python
from sympy.abc import _clash2

f = sympify('beta^alpha * x^(alpha - 1) * e^(-beta * x)/ Gamma(alpha)', locals=_clash2)
f
```


```python
srepr(f)
```




    "Mul(Pow(Symbol('beta'), Symbol('alpha')), Pow(Symbol('e'), Mul(Integer(-1), Symbol('beta'), Symbol('x'))), Pow(Symbol('x'), Add(Symbol('alpha'), Integer(-1))), Pow(Function('Gamma')(Symbol('alpha')), Integer(-1)))"



Several symbols can be defined at once. The string is split at spaces.


```python
y,z = symbols('y z')
```

Let's substitute $y+1$ for $x$.


```python
b = a.subs(x,y+1)
b
```

`sympy` automatically applies simplifications of elementary functions which are correct everywhere.


```python
cos(-x)
```


```python
cos(pi/4)
```

The method `rewrite` tries to rewrite an expression in terms of a given function.


```python
cos(x).rewrite(exp)
```

Be careful that some functions cannot always be simplified by default:


```python
sqrt(x)**4, sqrt(x**4)
```

Why not $x^2$? Try $x=i$.


```python
x = symbols('x',real=True)
sqrt(x**4)
```

## 4.2 Plots

Let's plot `a = x**2-1` and `b = a.subs(x,y+1)`. To do this, we first change the variable $y$ to $x$:


```python
a = x**2-1
b = b.subs(y,x)
b
```

`sympy` uses `matplotlib`. However, it distributes $x$ points adaptively, not uniformly.


```python
%matplotlib inline
```


```python
p = plot(a, b, show=False)
p[1].line_color = "r"
p.legend = True
p.show()
```

We can also create plots in which functions are defined implicitly:


```python
equation = Eq(x**2,y)
equation
```


```python
plot_implicit(equation,(x,-1,1),(y,0,1))
```


```python
plot_implicit(x**2+y**2>1,(x,-1,1),(y,-1,1))
```

Some additional plotting functions can be imported from `sympy.plotting`.


```python
from sympy.plotting import plot_parametric, plot3d, plot3d_parametric_line, plot3d_parametric_surface
```


```python
t = Symbol('t')
plot_parametric(sin(2*t),cos(3*t),(t,0,2*pi),
                title='Lissajous',xlabel='x',ylabel='y')
```

We can also create some 3d plots using `matplotlib` via `sympy`:


```python
# show plots in new window
%matplotlib qt
```

First we plot the surface $f(x,y) = x^2 + y^2$:


```python
plot3d(x**2+y**2,(x,-2,2),(y,-2,2))
```




    <sympy.plotting.plot.Plot at 0x20cf7d1c5c8>



Now a parametric space curve - a spiral.


```python
plot3d_parametric_line(cos(t),sin(t),0.1*t,(t,0,4*pi))
```




    <sympy.plotting.plot.Plot at 0x20cf859c708>



And finally a torus:


```python
u,v = symbols('u v')
a = .5
plot3d_parametric_surface((1+a*cos(u))*cos(v),(1+a*cos(u))*sin(v),a*sin(u),(u,0,2*pi),(v,0,2*pi))
```




    <sympy.plotting.plot.Plot at 0x20cf7ea4f08>



## 4.3 Polynomials

`SymPy` does not expand brackets automatically. The function `expand` is used for this.


```python
a = (x+y-z)**6
a
```


```python
a = expand(a)
a
```

Degree of the polynomial $a$ in $x$:


```python
degree(a,x)
```

Roots of the polynomial $a$:


```python
roots(a,x)
```

Let's collect terms with same power of $x$ together.


```python
collect(a,y)
```

Any polynomial with integer coefficients can be factorized into polynomials with integer coefficients (which cannot be factorized further). There exist efficient algorithms to do this.


```python
a = factor(a)
a
```

`sympy` does not automatically cancel ratios of polynomials by their greatest common divisor. The function `cancel` is used for this.


```python
a = (x**3-y**3)/(x**2-y**2)
a
```


```python
cancel(a)
```

`sympy` does not automatically bring sums of rational expressions to common denominator. The function `together` is used for this.


```python
a = y/(x-y)+x/(x+y)
a
```


```python
together(a)
```

The function `simplify` tries to rewrite an expression *in a simplest form*. This concept is not well defined (different forms may be considered simplest in different contexts), and there exists no algorithm for such simplification. The function `simplify` works heuristically, and it is not possible to guess in advance what simplifications it will try. It is very convenient in interactive sessions in order to check if it will succeed in rewriting an expression in some reasonably good form. But it is not desirable to use it in programs. There one should better use more specialized functions which perform well defined expression transformations.


```python
simplify(a)
```

We can solve the equation $a = 0$ for $x$:


```python
solve(Eq(a,0),x)[1] # solve(Eq(a,0),x)
```

Partial fraction decomposition with respect to $x$:


```python
apart(a,x)
```

Let's substitute some values for the symbols $x$ and $y$.


```python
evaluation = a.subs({x:1,y:2})
evaluation
```

And how much is it numerically?


```python
evaluation.n()
```

`subs` can also be used to interchange two variables. However, we need to be careful. Compare the difference:


```python
a.subs([(x,y),(y,x)])
# a.subs([(y,x),(x,y)])
```


```python
a.subs([(x,y),(y,x)],simultaneous=True)
```

<div class="alert alert-info"><h4>Homework</h4><p>Read Section 4.4-4.6.</p></div>

## 4.4 Series

### 4.4.1 Taylor Series

The Taylor series can be computed using the function `series`:


```python
exp(x).series(x,0,n=5)
```

Let's prepare 3 series.


```python
sinx = series(sin(x),x,0,n=8)
sinx
```


```python
cosx = series(cos(x),x,n=8)
cosx
```


```python
tanx = series(tan(x),x,n=8)
tanx
```

Products and ratios of series are not calculated automatically, the function `series` should be applied to them.


```python
series(sinx / cosx,n=8) # sinx / cosx
```

The following series should be equal to 1. But since `sinx` and `cosx` are known only with a limited accuracy, we obtain 1 with the same accuracy.


```python
series(sinx**2+cosx**2,n=8)
```

Series can be differentiated and integrated using the functions `diff` and `integrate`, respectively.


```python
diff(cosx,x)
```


```python
integrate(cosx,x)
```

It is not possible to substitute a numerical value for the expansion variable in a series (and hence to plot it). To this end one has to remove the $\mathcal{O}$ term first, transforming a series into a polynomial.


```python
sinx.removeO()
```

### 4.4.2 Summing Series

In `scipy`, `oo` means $\infty$.


```python
n = Symbol('n')
summation(1/n**2,(n,1,oo))
```


```python
summation(1/n**4,(n,1,oo))
```

An unevaluated sum is denoted as `Sum`.


```python
a = Sum(x**n/factorial(n),(n,0,oo))
Eq(a,a.doit())
```

## 4.5 Linear Algebra


```python
a,b,c,d,e,f = symbols('a b c d e f')
```

A matrix can be constructed from a list of lists.


```python
M=Matrix([[a,b,c],[d,e,f]])
M
```




$\displaystyle \left[\begin{matrix}a & b & c\\d & e & f\end{matrix}\right]$




```python
M.shape
```

A matrix can be constructed from a function.


```python
def g(i,j):
    return Rational(1,i+j+1)
Matrix(3,3,g)
```




$\displaystyle \left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5}\end{matrix}\right]$



Or from an undefined function.


```python
g = Function('g')
M = Matrix(3,3,g)
M
```




$\displaystyle \left[\begin{matrix}g{\left(0,0 \right)} & g{\left(0,1 \right)} & g{\left(0,2 \right)}\\g{\left(1,0 \right)} & g{\left(1,1 \right)} & g{\left(1,2 \right)}\\g{\left(2,0 \right)} & g{\left(2,1 \right)} & g{\left(2,2 \right)}\end{matrix}\right]$




```python
M[1,2] = 0
M
```




$\displaystyle \left[\begin{matrix}g{\left(0,0 \right)} & g{\left(0,1 \right)} & g{\left(0,2 \right)}\\g{\left(1,0 \right)} & g{\left(1,1 \right)} & 0\\g{\left(2,0 \right)} & g{\left(2,1 \right)} & g{\left(2,2 \right)}\end{matrix}\right]$




```python
M[2,:]
```




$\displaystyle \left[\begin{matrix}g{\left(2,0 \right)} & g{\left(2,1 \right)} & g{\left(2,2 \right)}\end{matrix}\right]$



Diagonal matrices can be constructed using the function `diag`:


```python
M = Matrix([[a,1],[0,a]])
diag(1,M,2)
```




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



With `sympy`, we can just use `*` for matrix multiplication:


```python
A = Matrix([[a,b],[c,d]])
B = Matrix([[1,2],[3,4]])
A*B
```




$\displaystyle \left[\begin{matrix}a + 3 b & 2 a + 4 b\\c + 3 d & 2 c + 4 d\end{matrix}\right]$




```python
B*A
```




$\displaystyle \left[\begin{matrix}a + 2 c & b + 2 d\\3 a + 4 c & 3 b + 4 d\end{matrix}\right]$




```python
simplify(A**(-1))
```




$\displaystyle \left[\begin{matrix}\frac{d}{a d - b c} & - \frac{b}{a d - b c}\\- \frac{c}{a d - b c} & \frac{a}{a d - b c}\end{matrix}\right]$



Consider the following matrix:


```python
x = Symbol('x',real=True)
M = Matrix([[(1-x)**3*(3+x),4*x*(1-x**2),-2*(1-x**2)*(3-x)],
          [4*x*(1-x**2),-(1+x)**3*(3-x),2*(1-x**2)*(3+x)],
          [-2*(1-x**2)*(3-x),2*(1-x**2)*(3+x),16*x]])
M
```




$\displaystyle \left[\begin{matrix}\left(1 - x\right)^{3} \left(x + 3\right) & 4 x \left(1 - x^{2}\right) & \left(3 - x\right) \left(2 x^{2} - 2\right)\\4 x \left(1 - x^{2}\right) & - \left(3 - x\right) \left(x + 1\right)^{3} & \left(2 - 2 x^{2}\right) \left(x + 3\right)\\\left(3 - x\right) \left(2 x^{2} - 2\right) & \left(2 - 2 x^{2}\right) \left(x + 3\right) & 16 x\end{matrix}\right]$



The determinant.


```python
det(M)
```

The eigenvalues and their multiplicities.


```python
M.eigenvals()
```

The eigenvector corresponding to the second eigenvalue.


```python
M.eigenvects()[1][2][0]
```




$\displaystyle \left[\begin{matrix}\frac{x}{2} + \frac{1}{2}\\\frac{x + 1}{x - 1}\\1\end{matrix}\right]$



## 4.6 Differential Equations


```python
t = Symbol('t')
x = Function('x')
p = Function('p')
```

First we solve the first-order ODE $\frac{\mathrm{d}}{\mathrm{d}t}x(t) = -x(t)$.


```python
dsolve(Eq(diff(x(t),t),-x(t)),x(t))
```

Now the second-order ODE $\frac{\mathrm{d}^2}{\mathrm{d}t^2}x(t) = -x(t)$.


```python
dsolve(Eq(diff(x(t),t,2),x(t)),x(t))
```

Finally a system of first-order equations:
\begin{align*}
\begin{cases}
\frac{\mathrm{d}}{\mathrm{d}t}x(t) = p(t), \\
\frac{\mathrm{d}}{\mathrm{d}t}p(t) = -x(t).
\end{cases}
\end{align*}


```python
dsolve((diff(x(t),t)-p(t),diff(p(t),t)+x(t)))
```

<div class="alert alert-success"><h4>Done reading?</h4><p>Then read the final section:</p></div>

## 4.7 Limits, Derivatives, Integrals

### 4.7.1 Limits

Limits can be computed by using the function `limit`. We will do this for the function

\begin{align*}
f(x) = \frac{\tan(\sin x)-\sin(\tan x)}{x^7}.
\end{align*}


```python
%matplotlib inline
x = Symbol('x')
f = (tan(sin(x))-sin(tan(x)))/x**7
plot(f, (x,-2,2))
```


```python
limit(f,x,0)
```

This limit is easy to compute as it is a two-sided limit. Now let's calculate one-sided limits for the function

\begin{align*}
f(x) = \frac{\tan(\sin x)-\sin(\tan x)}{x^7 + \exp\left(-\frac{1}{x}\right)}.
\end{align*}


```python
a = (tan(sin(x))-sin(tan(x)))/(x**7 + exp(-1/x))
plot(a, xlim=[-.25,.25], ylim=[-0.5,0.5]) # does not always work
```


```python
limit(a,x,0,'+')
```


```python
limit(a,x,0,'-')
```

### 4.7.2 Derivatives

We compute derivatives using the function `diff`.


```python
a = x*sin(x+y)
diff(a,x,1)
```


```python
diff(a,y)
```

Let's take the second derivative in $x$ and the first one in $y$.


```python
diff(a,x,2,y)
```

Expressions with undefined functions can also be differentiated.


```python
f = Function('f')
b = diff(x*f(x**2),x)
```

What's this?


```python
b
```




$\displaystyle 2 x^{2} \left. \frac{d}{d \xi_{1}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=x^{2} }} + f{\left(x^{2} \right)}$



The function `Derivative` represents an unevaluated derivative. It can be evaluated by the method `doit`.


```python
a = Derivative(sin(x),x)
Eq(a,a.doit())
# diff(sin(x),x)
```

### 4.7.3 Integrals

Integrals can be computed using the function `integrate`.


```python
a = 1/(x*(x**2-2)**2)
integrate(a,x)
```


```python
integrate(1/(exp(x)+1),x)
```


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


```python
a = integrate(x**x,x)
a.doit()
```

This is an unevaluated integral (so it is written as `Integral`).


```python
print(a)
```

    Integral(x**x, x)



```python
a = Integral(sin(x),x)
Eq(a,a.doit())
```

### What's next? SAGE.

### 🔴 *Next Week:*


```python
import numpy as np
np.random.choice(['Statistics','Something else'], p=[0.99,0.01])
```




    'Statistics'



886b8ee9d57b7581ff6b615b2ae841bdad3b103d
ipynb
Jupyter Notebook
Notebooks/4 - Sympy for Mathematicians.ipynb
Roshanmahes/Python-for-Teachers
bc82948cb02305fc944caea0ece0941643d344a6
MIT
2020-08-15T12:45:39.000Z
2022-01-29T13:46:15.000Z
Notebooks/4 - Sympy for Mathematicians.ipynb
Roshanmahes/Python-for-Teachers
bc82948cb02305fc944caea0ece0941643d344a6
MIT
Notebooks/4 - Sympy for Mathematicians.ipynb
Roshanmahes/Python-for-Teachers
bc82948cb02305fc944caea0ece0941643d344a6
MIT
2020-10-02T22:20:42.000Z
2021-06-03T10:30:32.000Z
Qwen/Qwen-72B
1. YES
2. YES