# Python - Symbolic Mathematics (`sympy`)


```python
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
```


```python
sp.init_printing()
```

### `sympy` treats stuff fundementally different than `numpy`


```python
np.sqrt(8)
```


```python
type(np.sqrt(8))
```


```python
sp.sqrt(8)
```


```python
type(sp.sqrt(8))
```


```python
np.sqrt(8) == sp.sqrt(8)
```


```python
np.pi
```


```python
np.pi + 4
```


```python
sp.pi
```


```python
sp.pi + 4
```

#### You can use `np.float64()` to turn a `sympy` float to a `numpy` float


```python
my_sympy_to_numpy_pi = np.float64(sp.pi)
```


```python
type(my_sympy_to_numpy_pi)
```


```python
my_sympy_to_numpy_pi + 4
```

### sympy has its own way to handle rational numbers


```python
sp.Rational(3,5)
```


```python
sp.Rational(3,5) + sp.Rational(1,3)
```

#### Adding `.n()` to the end of a sympy expression will `evaluate` expression


```python
sp.Rational(3,5).n()
```


```python
sp.pi.n()
```

#### You can add a value to `.n(value)` to set the number of significant figures


```python
sp.Rational(222,7).n()
```


```python
sp.Rational(222,7).n(1)
```


```python
sp.Rational(222,7).n(3)
```

#### `nsimplify()` will sort-of do the reverse of `.n()`


```python
sp.nsimplify(0.6)
```


```python
sp.nsimplify(4.242640687119286)
```


```python
sp.nsimplify(sp.pi, tolerance=1e-2)
```


```python
sp.nsimplify(sp.pi, tolerance=1e-5)
```


```python
sp.nsimplify(sp.pi, tolerance=1e-6)
```

### ... to $\infty$ and beyond


```python
sp.oo
```


```python
sp.oo + 3
```


```python
1e199 < sp.oo
```

---
## Matrix

* There are a couple of different ways to create a Matrix
* `Matrix(rows of lists)`
* `Matrix(shape, list)`

#### `rows of lists` is pretty straightforward


```python
sp.Matrix([[1,2],[3,4]])
```

#### `shape, list` gives you more flexibility


```python
sp.Matrix(2,2,[1,2,3,4])
```


```python
sp.Matrix(1,4,[1,2,3,4])
```


```python
sp.Matrix(4,1,[1,2,3,4])
```


```python
my_matrix = sp.Matrix(2,2,[1,2,3,4])

my_matrix
```


```python
my_matrix.det()
```


```python
my_matrix ** -1
```


```python
(my_matrix ** -1) * my_matrix
```

---
# Symbolic

### You have to explicitly tell `SymPy` what symbols you want to use.

* Once you declare symbols to use in `sympy` they are unavaliable for other packages


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

### Expressions are then able use these symbols


```python
my_equation = (a * x**2 * y)  + (b * x * y) + (c * x * y**2) + 6

my_equation
```


```python
my_equation + 3
```


```python
my_equation / x
```


```python
sp.simplify(my_equation / x)
```


```python
sp.collect(my_equation,x)
```


```python
sp.collect(my_equation,y)
```

### You can evaluate equations for specific values


```python
my_equation_x = sp.collect(my_equation,x)

my_equation_x
```


```python
my_equation_x.subs({y:sp.Rational(1,2), a:4, b:2, c:8})
```

### You can evaluate equations sybolically


```python
my_equation_x
```


```python
my_y = (2*x + 3)

my_y
```

#### Replace $y$ with $2x + 3$ using `.subs()`


```python
my_equation_x.subs(y, my_y)
```

#### Multiply everything through using `.expand()`


```python
sp.expand(my_equation_x.subs(y, my_y))
```

#### Collect the terms of $x$ using `.collect()`


```python
sp.collect(sp.expand(my_equation_x.subs(y, my_y)),x)
```

#### Evaluate the expression for some values of $a, b, c$ using `.subs()`


```python
sp.collect(sp.expand(my_equation_x.subs(y, my_y)),x).subs({a:4, b:2, c:8})
```

---
## Calculus


```python
my_equation
```


```python
sp.diff(my_equation,x)
```


```python
sp.diff(my_equation,x,2)
```


```python
sp.integrate(my_equation,x)
```


```python
sp.integrate(my_equation,(x,0,5))   # limits x = 0 to 5
```


```python
sp.integrate(my_equation,(x,0,5)).n()
```


```python
sp.integrate(my_equation,(x,0,5)).subs({y:sp.Rational(1,2), a:4, b:2, c:8})
```


```python
sp.integrate(my_equation,(x,0,5)).subs({y:sp.Rational(1,2), a:4, b:2, c:8}).n()
```

---
# Solving equations - `solve`

* Need to rearrange equation so that it is set to equal zero 

### One equation

$$ \large
3x - 3 =  30 \hspace{1cm} \rightarrow \hspace{1cm} 3x - 33 = 0
$$


```python
equation_in_x = 3 * x - 33

equation_in_x
```


```python
sp.solve([equation_in_x], [x])
```

### System of equations

$$ \large
\begin{array}{c}
9x - 2y = 5\\
-2x + 6y = 10\\
\end{array}
$$


```python
equation_a = 9*x - 2*y - 5
equation_b = -2*x + 6*y - 10
```


```python
sp.solve([equation_a, equation_b], [x,y])
```

---

### We can also do it the Matrix way

$$ \large
\begin{bmatrix}
9 & -2 \\
-2 & 6 \\
\end{bmatrix}
\begin{bmatrix}
x \\
y
\end{bmatrix}
=
\begin{bmatrix}
5 \\
10
\end{bmatrix}
$$

$$ \large
A \vec x = b
$$


```python
my_A = sp.Matrix(2,2,[9, -2,-2, 6])

my_A
```


```python
my_b = sp.Matrix(2,1,[5,10])

my_b
```

$$ \large
\vec x = A^{-1}b
$$


```python
my_A ** -1 * my_b
```

---
## Let's do some graphing stuff ...

### In the following examples - notice the difference between `numpy` variables and `sympy` variables!


&nbsp;

$$
\large f(x) = 2\,\cos(5x) \ e^{-x}
$$

&nbsp;

### Need to create a `numpy` arrays to do the graphing


```python
# 200 points between -pi and pi

my_np_x = np.linspace(-np.pi, np.pi, 200)
```


```python
my_np_fx = 2 * np.cos(5 * my_np_x) * np.exp(-my_np_x)
```


```python
plt.style.use('ggplot')
```


```python
fig,ax = plt.subplots(1,1)
fig.set_size_inches(10,6)

fig.tight_layout()

ax.set_ylim(-11,11)
ax.set_xlim(-np.pi,np.pi)

ax.set_xlabel("This is X", fontsize = 14)
ax.set_ylabel("This is Y", fontsize = 14)

ax.plot(my_np_x, my_np_fx, 
        color = (1.0, 0.0, 0.0, 0.5), 
        marker='None', 
        linestyle='-', 
        linewidth = 6,
        label = "f(x)")

ax.legend(loc = 0, fontsize = 24, shadow = True);
```

---
## Taylor Series

* Taylor series are polynomial approximations at a specific point


```python
my_sp_fx = 2 * sp.cos(5 * x) * sp.exp(-x)
my_sp_fx
```

### Taylor series of f(x) at f(x) = 0


```python
my_taylor = sp.series(my_sp_fx, x, x0 = 0)

my_taylor
```

### If you want differnt number of terms

* n = magnitude of the highest term
* n = 4 means all terms up to x$^{4}$ or $\mathcal{O}(4)$


```python
my_taylor = sp.series(my_sp_fx, x, x0 = 0, n=4)

my_taylor
```


```python
my_taylor.removeO()
```


```python
my_taylor.removeO().n(3)
```


```python
# Make NumPy versions of the term to plot

my_np_1term = -2.0 * my_np_x + 2.0
my_np_2term = -24.0 * my_np_x**2 - 2.0 * my_np_x + 2.0
my_np_3term = 24.7 * my_np_x**3 - 24.0 * my_np_x**2 - 2.0 * my_np_x + 2.0
```


```python
fig,ax = plt.subplots(1,1)
fig.set_size_inches(10,8)

fig.tight_layout()

ax.set_ylim(-4,4)
ax.set_xlim(-1,1)

ax.set_xlabel("This is X", fontsize = 14)
ax.set_ylabel("This is Y", fontsize = 14)

ax.plot(my_np_x, my_np_fx, 
        color = (1.0, 0.0, 0.0, 0.5), 
        marker='None', 
        linestyle='-', 
        linewidth = 10,
        label = "f(x)")

ax.plot(my_np_x, my_np_1term, color='b', marker='None', linestyle='--', label="1-term")
ax.plot(my_np_x, my_np_2term, color='g', marker='None', linestyle='--', label="2-term")
ax.plot(my_np_x, my_np_3term, color='k', marker='None', linestyle='--', label="3-term")

ax.legend(loc = 0, fontsize = 18);
```

---
## General Equation Solving - `nsolve`

$$
\Large f(x) = 2\,\cos(5x) \ e^{-x} \\[10pt]
\Large g(x) = \frac{3}{2} \left [\frac{x^3}{\pi} - \pi x \right]\\
$$

### Make a `numpy` version of g(x)


```python
my_np_gx = 3/2 * (my_np_x ** 3 / np.pi - np.pi * my_np_x)
```

### Where do they cross? - The graph


```python
fig,ax = plt.subplots(1,1)
fig.set_size_inches(10,6)

fig.tight_layout()

ax.set_ylim(-7,7)
ax.set_xlim(-np.pi,np.pi)

ax.set_xlabel("This is X", fontsize = 14)
ax.set_ylabel("This is Y", fontsize = 14)

ax.plot(my_np_x, my_np_fx, 
        color = (1.0, 0.0, 0.0, 0.5), 
        marker='None', 
        linestyle='-', 
        linewidth = 6,
        label = "f(x)")

ax.plot(my_np_x, my_np_gx, 
        color = (0.25, 0.0, 0.75, 0.5), 
        marker='None', 
        linestyle='-', 
        linewidth = 6,
        label = "g(x)")


ax.legend(loc = 0, fontsize = 24);
```

### Where do they cross? - The `sympy` solution

### Make a `sympy` version of g(x)


```python
my_sp_gx = sp.Rational(3,2) * (x ** 3 / sp.pi - sp.pi * x)

my_sp_gx
```


```python
my_sp_fx, my_sp_gx
```

### Need to provide an initial guess


```python
my_guess = 3.0

sp.nsolve(my_sp_fx - my_sp_gx, x, my_guess)
```


```python
all_guesses = (3.0, 0, -1.0)

for val in all_guesses:
    result = sp.nsolve(my_sp_fx - my_sp_gx, x, val)
    print(result)
```

### Your guess has to be (somewhat) close or the solution will not converge:


```python
my_guess = -40

sp.nsolve(my_sp_fx - my_sp_gx, x, my_guess)
```

---
## Primes


```python
# List of primes in the range 0 -> 100

list(sp.primerange(0,100))
```


```python
# The 100th prime number

sp.prime(100)
```


```python
# The next prime after 2020

sp.nextprime(2020)
```


```python
# The prime factors of this year

sp.factorint(2022)
```

# `SymPy` can do *so* much more. It really is magic. 

## Complete documentation can be found [here](http://docs.sympy.org/latest/index.html)

e37e741db22ba6ce9ce42f7639effaa4399fefd9
ipynb
Jupyter Notebook
Python_SymPy.ipynb
UWashington-Astro300/Astro300-W22
371eb704030a104cb8e826bb14c353e5a863b0f3
MIT
Python_SymPy.ipynb
UWashington-Astro300/Astro300-W22
371eb704030a104cb8e826bb14c353e5a863b0f3
MIT
Python_SymPy.ipynb
UWashington-Astro300/Astro300-W22
371eb704030a104cb8e826bb14c353e5a863b0f3
MIT
Qwen/Qwen-72B
1. YES
2. YES