# SymPy tutorial

**SymPy** is a Python package for performing **symbolic mathematics**
which can perform algebra, integrate and differentiate equations, 
find solutions to differential equations, and *numerically solve
messy equations* -- along other uses.

CHANGE LOG
    
    2017-06-12  First revision since 2015-12-26.

Let's import sympy and initialize its pretty print functionality 
which will print equations using LaTeX.
Jupyter notebooks uses Mathjax to render equations
so we specify that option.


```python
import sympy as sym
sym.init_printing(use_latex='mathjax')

#  If you were not in a notebook environment,
#  but working within a terminal, use:
#
#  sym.init_printing(use_unicode=True)
```

## Usage

These sections are illustrated with examples drawn from
[rlabbe](https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/Appendix-A-Installation.ipynb) from his appendix for Kalman Filters.

It is important to distinguish a Python variable
from a **declared symbol** in sympy.


```python
phi, x = sym.symbols('\phi, x')

#  x here is a sympy symbol, and we form a list:
[ phi, x ]
```




$$\left [ \phi, \quad x\right ]$$



Notice how we used a LaTeX expression for the symbol `phi`.
This is not necessary, but if you do the output will render nicely as LaTeX.

Also notice how $x$ did not have a numerical value for the list to evaluate.

So what is the **derivative** of $\sqrt{\phi}$ ?


```python
sym.diff('sqrt(phi)')
```




$$\frac{1}{2 \sqrt{\phi}}$$



We can **factor** equations:


```python
sym.factor( phi**3 - phi**2 + phi - 1 )
```




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



and we can **expand** them:


```python
((phi+1)*(phi-4)).expand()
```




$$\phi^{2} - 3 \phi - 4$$



You can also use strings for equations that use symbols that you have not defined:


```python
x = sym.expand('(t+1)*2')
x
```




$$2 t + 2$$



## Symbolic solution

Now let's use sympy to compute the **Jacobian** of a matrix. 
Suppose we have a function,

$$h=\sqrt{(x^2 + z^2)}$$

for which we want to find the Jacobian with respect to x, y, and z.


```python
x, y, z = sym.symbols('x y z')

H = sym.Matrix([sym.sqrt(x**2 + z**2)])

state = sym.Matrix([x, y, z])

H.jacobian(state)
```




$$\left[\begin{matrix}\frac{x}{\sqrt{x^{2} + z^{2}}} & 0 & \frac{z}{\sqrt{x^{2} + z^{2}}}\end{matrix}\right]$$



Now let's compute the discrete process noise matrix $\mathbf{Q}_k$ given the continuous process noise matrix 
$$\mathbf{Q} = \Phi_s \begin{bmatrix}0&0&0\\0&0&0\\0&0&1\end{bmatrix}$$

and the equation

$$\mathbf{Q} = \int_0^{\Delta t} \Phi(t)\mathbf{Q}\Phi^T(t) dt$$

where 
$$\Phi(t) = \begin{bmatrix}1 & \Delta t & {\Delta t}^2/2 \\ 0 & 1 & \Delta t\\ 0& 0& 1\end{bmatrix}$$


```python
dt = sym.symbols('\Delta{t}')

F_k = sym.Matrix([[1, dt, dt**2/2],
                  [0,  1,      dt],
                  [0,  0,      1]])

Q = sym.Matrix([[0,0,0],
                [0,0,0],
                [0,0,1]])

sym.integrate(F_k*Q*F_k.T,(dt, 0, dt))
```




$$\left[\begin{matrix}\frac{\Delta{t}^{5}}{20} & \frac{\Delta{t}^{4}}{8} & \frac{\Delta{t}^{3}}{6}\\\frac{\Delta{t}^{4}}{8} & \frac{\Delta{t}^{3}}{3} & \frac{\Delta{t}^{2}}{2}\\\frac{\Delta{t}^{3}}{6} & \frac{\Delta{t}^{2}}{2} & \Delta{t}\end{matrix}\right]$$



## Numerical solution

You can find the *numerical value* of an equation by substituting in a value for a variable:


```python
x = sym.symbols('x')

w = (x**2) - (3*x) + 4
w.subs(x, 4)
```




$$8$$



Typically we want a numerical solution where the analytic solution is messy,
that is, we want a **solver**.
This is done by specifying a sympy equation, for example:


```python
LHS = (x**2) - (8*x) + 15
RHS = 0
#  where both RHS and LHS can be complicated expressions.

solved = sym.solveset( sym.Eq(LHS, RHS), x, domain=sym.S.Reals )
#  Notice how the domain solution can be specified.

solved
#  A set of solution(s) is returned.
```




$$\left\{3, 5\right\}$$




```python
#  Testing whether any solution(s) were found:
if solved != sym.S.EmptySet:
    print("Solution set was not empty.")
```

    Solution set was not empty.



```python
#  sympy sets are not like the usual Python sets...
type(solved)
```




    sympy.sets.sets.FiniteSet




```python
#  ... but can easily to converted to a Python list:
l = list(solved)
print( l, type(l) )
```

    ([3, 5], <type 'list'>)



```python
LHS = (x**2)
RHS = -4
#  where both RHS and LHS can be complicated expressions.

solved = sym.solveset( sym.Eq(LHS, RHS), x )
#  Leaving out the domain will include the complex domain.

solved
```




$$\left\{- 2 i, 2 i\right\}$$



## Application to financial economics

We used sympy to deduce parameters of Gaussian mixtures
in module `lib/ys_gauss_mix.py` and the explanatory notebook
is rendered at https://git.io/gmix 

5e9266c19ae2a3149b23ef250a800397c9b532df
ipynb
Jupyter Notebook
docs/fecon235-08-sympy.ipynb
maidenlane/five
bf14dd37b0f14d6998893c2b0478275a0fc55a82
BSD-3-Clause
2020-04-24T05:29:26.000Z
2020-04-24T05:29:26.000Z
fecon235-master/docs/fecon235-08-sympy.ipynb
maidenlane/five
bf14dd37b0f14d6998893c2b0478275a0fc55a82
BSD-3-Clause
fecon235-master/docs/fecon235-08-sympy.ipynb
maidenlane/five
bf14dd37b0f14d6998893c2b0478275a0fc55a82
BSD-3-Clause
2020-04-24T05:34:06.000Z
2020-04-24T05:34:06.000Z
Qwen/Qwen-72B
1. YES
2. YES