# Symbolic Computation with SymPy
Symbolic computation deals with the computation of mathematical objects symbolically. This means that the mathematical objects are represented exactly, not approximately, and mathematical expressions with unevaluated variables are left in symbolic form.


```python
from sympy import *
```

Here is a brief example, if you wanted to compute a square root in python you would do:


```python
import math
math.sqrt(8)
```




    2.8284271247461903



This is an approximate solution because $\sqrt{8}$ is an irrational number. But if you were to use the sympy version of of the square root function:


```python
sqrt(8)
```




$\displaystyle 2 \sqrt{2}$



Notice that this is an exact solution, where the irrational part is left unevaluated by default.

Let's move to a more interesting examples

---
## The power of symbolic computation


```python
x, y = symbols('x y')
type(x)
```




    sympy.core.symbol.Symbol




```python
y
```




$\displaystyle y$



Now the code understands that `x` and `y` are symbols, not variables. These symbols can be used to build expressions:


```python
expr = x + 2*y
expr
```




$\displaystyle x + 2 y$




```python
expr + 1
```




$\displaystyle x + 2 y + 1$




```python
(x**2)*expr
```




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



Notice that in this last expression, the $x^2$ term was not distributed, but that can be requested.


```python
expand((x**2)*expr)
```




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



The reverse operation of `expand` is `factor`:


```python
factor(x**2 + x)
```




$\displaystyle x \left(x + 1\right)$



Here are a few more interesting examples:


```python
sp.diff(x*sin(x), x)
```




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




```python
sp.integrate(x*cos(x) + sin(x), x)
```




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




```python
sp.integrate(sin(x**2), (x, -oo, oo))
```




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




```python
sp.limit(sin(x)/x, x, 0)
```




$\displaystyle 1$




```python
sp.solve(x**2-4, x)
```




    [-2, 2]




```python
t = symbols('t')
y = Function('y')
ode = dsolve(Eq(y(t).diff(t, t) - y(t), exp(t)), y(t))
ode
```




$\displaystyle y{\left(t \right)} = C_{2} e^{- t} + \left(C_{1} + \frac{t}{2}\right) e^{t}$




```python
latex(ode)
```




    'y{\\left(t \\right)} = C_{2} e^{- t} + \\left(C_{1} + \\frac{t}{2}\\right) e^{t}'



---
## Symbols VS Variables

Let's understand a bit of details of SymPy's implementation.

`symbols` takes a string of variable names separated by spaces or commas, and creates Symbols out of them. We can then assign these to variable names.


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

The important thing to notice here is that the name of a Symbol and the name of the variable it is assigned to need not have anything to do with one another. In other words, **the variable and the symbol are not the same thing**.


```python
a, b, c = symbols('x, y, z')
```


```python
a
```




$\displaystyle x$




```python
crazy = symbols('unralated')
crazy + 1
```




$\displaystyle unralated + 1$



Take the following example. Changing `x` to 2 had no effect on `expr`. This is because `x = 2` changes the Python variable `x` to 2, but has no effect on the SymPy Symbol x, which was what we used in creating `expr`.


```python
x = symbols('x')
expr = x + 1
x = 2
display(expr)
```


$\displaystyle x + 1$


---
## Basic Operations
Here we discuss some of the most basic operations needed for expression manipulation in SymPy.


```python
x, y, z = symbols("x y z")
```

One of the most common things you might want to do with a mathematical expression is **substitution**. Substitution replaces all instances of something in an expression with something else. It is done using the `subs` method. For example:


```python
expr = cos(x) + 1
expr.subs(x, y)
```




$\displaystyle \cos{\left(y \right)} + 1$




```python
expr.subs(x, 0)
```




$\displaystyle 2$




```python
expr.subs(x, 2)
```




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



Notice that this last example the value of the expression was not evaluated. That is because $\cos\left(2\right)$ is irrational (proof [here](https://math.stackexchange.com/questions/94478/sin-1-circ-is-irrational-but-how-do-i-prove-it-in-a-slick-way-and-tan1)). If we wanted to evaluate the expression, we would use the `evalf` method.


```python
expr.subs(x, 2).evalf()
```




$\displaystyle 0.583853163452858$



Notice that the output has the traditional 15 decimal places. But the `evalf` method lets you choose how many digits you want.


```python
expr.subs(x, 2).evalf(30)
```




$\displaystyle 0.583853163452857613002431770499$



Take a look at the first 100 digits of $\pi$


```python
pi.evalf(100)
```




$\displaystyle 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$



To perform multiple substitutions at once, pass a list of `(old, new)` pairs to `subs`.


```python
expr = x**3 + 4*x*y - z
expr.subs([(x, 2), (y, 4), (z, 0)])
```




$\displaystyle 40$



You can use python's list comprehensions to generate substitutions. In the example below we replace all instances of $x$ that have an even power with $y$.


```python
expr = x**4 - 4*x**3 + 4*x**2 - 2*x + 3
replacements = [(x**i, y**i) for i in range(5) if i % 2 == 0]
expr.subs(replacements)
```




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



---
## Printing
As you may have noticed, the output of symbolyc expressions in the jupyter notebook environment automatically prints the rendered $\LaTeX$ expression. But you can print outputs in other formats as well.


```python
expr = Integral(sqrt(1/x), x)
expr
```




$\displaystyle \int \sqrt{\frac{1}{x}}\, dx$



You can ask for the latex source code.


```python
latex(expr)
```




    '\\int \\sqrt{\\frac{1}{x}}\\, dx'



Use the `pprint` to print the expression using Unicode characters.


```python
pprint(expr)
```

    ⌠           
    ⎮     ___   
    ⎮    ╱ 1    
    ⎮   ╱  ─  dx
    ⎮ ╲╱   x    
    ⌡           


Or non-unicode charcters


```python
pprint(expr, use_unicode=False)
```

      /          
     |           
     |     ___   
     |    / 1    
     |   /  -  dx
     | \/   x    
     |           
    /            


---
## Simplification
One of the most useful features of a symbolic manipulation system is the ability to simplify mathematical expressions. SymPy has dozens of functions to perform various kinds of simplification. There is also one general function called `simplify()` that attempts to apply all of these functions in an intelligent way to arrive at the simplest form of an expression.


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

SymPy knows trigonometrical identities


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




$\displaystyle 1$



Facotoring out commom factors in polynomials


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




$\displaystyle x - 1$



And even understands identities from big formulas of statitistical distributions, like the gamma function.


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




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



There is a caveat: "simplify" does not have a well-defined target, so the `simplify()` method might not be what you need.


```python
expand((x + 1)**2)
```




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




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




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



Another powerfull tool is the `rewrite` method. The example below rewrites the tangent functions in terms of the sine function


```python
tan(x).rewrite(sin)
```




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



And the factorial faction using the gamma function.


```python
factorial(x).rewrite(gamma)
```




$\displaystyle \Gamma\left(x + 1\right)$



There are other types of simplifications as well, for powers, more elaborate trigonometric identities, logarithims, etc. Check out the documentations [here](https://docs.sympy.org/latest/tutorial/simplification.html).

---
## Calculus


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


```python
expr = exp(x**2 * y**2)
expr
```




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



The first steps we usually take in calculus usually envolves limits.


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




$\displaystyle 1$



and we learn that limits might be different from the right and from the left.


```python
limit(1/x, x, 0, '+')
```




$\displaystyle \infty$




```python
limit(1/x, x, 0, '-')
```




$\displaystyle -\infty$



You can take analytical derivatives with the `diff` method


```python
diff(expr, x)
```




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



You can also take multiple derivatives by chaining input arguments


```python
diff(expr, x, x, y)
```




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



The same idea can be used for integrals.


```python
integrate(cos(x), x)
```




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




```python
integrate(cos(x), y)
```




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



To evaluate an integral, pass variable and their limits as tuples


```python
integrate(cos(x), (x, 0, pi/2))
```




$\displaystyle 1$




```python
integrate(cos(x), (y, 0, pi/2))
```




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



You can even integrate to infinity


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




$\displaystyle \pi$




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




$\displaystyle \pi$



SymPy can also do taylor expansions aroundo a given point. The example below make a 4th degree approximation to the sine function aroun $x=0$.


```python
sin(x).series(x, 0, 4)
```




$\displaystyle x - \frac{x^{3}}{6} + O\left(x^{4}\right)$



And if you are an economist and do not care the $O(x^n)$ notation, you can remove it.


```python
sin(x).series(x, 0, 4).removeO()
```




$\displaystyle - \frac{x^{3}}{6} + x$



---
## Solvers

Writing equations in SymPy assumes $=0$ at the end, unless you are explicit about the right side.


```python
Eq(x**2 -1)
```




$\displaystyle x^{2} - 1 = 0$




```python
Eq(x**2, 1)
```




$\displaystyle x^{2} = 1$




```python
Eq(x**2, x-1)
```




$\displaystyle x^{2} = x - 1$



We can solve equations


```python
expr = Eq(x**2 - 5*x + 6)
solve(expr, x)
```




    [2, 3]



Sympy can even give elaborate answers to solutions. In this case, use `solveset` instead of `solve` 


```python
solveset(x - x, x, domain=S.Reals)
```




$\displaystyle \mathbb{R}$




```python
solveset(x - x, x)
```




$\displaystyle \mathbb{C}$




```python
solveset(sin(x) - 1, x)
```




$\displaystyle \left\{2 n \pi + \frac{\pi}{2}\; |\; n \in \mathbb{Z}\right\}$




```python
solveset(x**2 + 1, x)
```




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



---
## Solving Differential Equations
The first step is to declare functions


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

`f` and `g` are now undefined functions. We can call `f(x)`, and it will represent an unknown function.


```python
f(x)
```




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




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




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




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




$\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)}$



To solve the equation, use the `dsolve` method.


```python
dsolve(diffeq, f(x))
```




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



Here, $C_{i}$ denotes arbitrary constants.

---
## Matrices
Matrix operations in SymPy are the same as in numpy. The only diffenrence is that now we can use symbols.


```python
x, y = symbols('x, y')
M = Matrix([[1, 2, 3], [x, 3, 1], [4, y, 6]])
M
```




$\displaystyle \left[\begin{matrix}1 & 2 & 3\\x & 3 & 1\\4 & y & 6\end{matrix}\right]$




```python
M.det()
```




$\displaystyle 3 x y - 12 x - y - 10$



Now take an example of a mtrix with incomplete rank


```python
M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])
M
```




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




```python
M.nullspace()  # or "kernel of the linear transformation"
```




    [Matrix([
     [-1],
     [-1],
     [ 1]])]




```python
M.columnspace()  # or "span of the linear transformation"
```




    [Matrix([
     [1],
     [2],
     [3]]), Matrix([
     [1],
     [1],
     [1]])]



Notice that each output is a sympy object


```python
display(M.diagonalize()[0])
display(M.diagonalize()[1])
```


$\displaystyle \left[\begin{matrix}-1 & - \frac{2 \sqrt{11}}{7} - \frac{3}{7} & - \frac{3}{7} + \frac{2 \sqrt{11}}{7}\\-1 & \frac{2}{7} - \frac{\sqrt{11}}{7} & \frac{2}{7} + \frac{\sqrt{11}}{7}\\1 & 1 & 1\end{matrix}\right]$



$\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & 3 - \sqrt{11} & 0\\0 & 0 & 3 + \sqrt{11}\end{matrix}\right]$


All of this was to get to the following example, the Jacobian Matrix:


```python
from sympy.abc import rho, phi
X = Matrix([rho*cos(phi), rho*sin(phi), rho**2])
Y = Matrix([rho, phi])
X.jacobian(Y)
```




$\displaystyle \left[\begin{matrix}\cos{\left(\phi \right)} & - \rho \sin{\left(\phi \right)}\\\sin{\left(\phi \right)} & \rho \cos{\left(\phi \right)}\\2 \rho & 0\end{matrix}\right]$



Just to grab a better view of what is happening here, let us build a simpler example


```python
x, y, z = symbols('x, y, z')
eq1 = 2*x**2 + 3*y**3 - 2*z
eq2 = 4*x + y**2 - 2*z**3
eq3 = 3*x + 2*y - 3*z**4
M1 = Matrix([eq1, eq2, eq3])
M2 = Matrix([x, y])
M1.jacobian(M2)
```




$\displaystyle \left[\begin{matrix}4 x & 9 y^{2}\\4 & 2 y\\3 & 2\end{matrix}\right]$



This kind of structure will be perfect for DSGE models

e66a2ec8d8378bd87231748ac5fbd18a9bb87c87
ipynb
Jupyter Notebook
Python Lectures/Section 11 - Symbolic Computation with Sympy.ipynb
Finance-Hub/FinanceHubMaterials
e06cae52ac34873413e946810ad5e6bf79b1c0dc
MIT
2019-11-12T04:52:31.000Z
2021-11-12T09:27:08.000Z
Python Lectures/Section 11 - Symbolic Computation with Sympy.ipynb
antoniosalomao/FinanceHubMaterials
0e2aead9c2a7c92a6826b6b47970afbfa30fb1b2
MIT
2020-04-04T03:03:17.000Z
2020-04-04T03:03:17.000Z
Python Lectures/Section 11 - Symbolic Computation with Sympy.ipynb
antoniosalomao/FinanceHubMaterials
0e2aead9c2a7c92a6826b6b47970afbfa30fb1b2
MIT
2019-06-28T15:35:10.000Z
2021-04-04T02:34:10.000Z
Qwen/Qwen-72B
1. YES
2. YES