# Python - Symbolic Mathematics (`sympy`)


```python
%matplotlib inline

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
sp.sqrt(8)
```


```python
np.pi
```


```python
sp.pi
```

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


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


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

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


```python
sp.nsimplify(0.125)
```


```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
1e9 < sp.oo
```

### Primes


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


```python
sp.nextprime(2018)
```


```python
sp.factorint(11192018)
```

# Symbolic

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


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

### Expressions are then able use these symbols


```python
my_equation = 2*x + y

my_equation
```


```python
my_equation + 3
```


```python
my_equation - x
```


```python
my_equation / x
```


```python
my_greek_equation = mu**2 / rho * (a + b)

my_greek_equation
```

### `SymPy` has all sorts of ways to manipulates symbolic equations


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


```python
another_equation = (x + 2) * (x - 3)

another_equation
```


```python
sp.expand(another_equation)
```


```python
long_equation = 2*y*x**3 + 12*x**2 - x + 3 - 8*x**2 + 4*x + x**3 + 5 + 2*y*x**2 + x*y

long_equation
```


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


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

### You can evaluate equations for specific values


```python
trig_equation = a*sp.sin(2*x + y) + b*sp.cos(x + 2*y)

trig_equation
```


```python
trig_equation.subs({a:2, b:3, x:4, y:5})
```


```python
trig_equation.subs({a:2, b:3, x:4, y:5}).n()
```


```python
sp.expand(trig_equation, trig=True)
```


```python
sp.collect(sp.expand(trig_equation, trig=True),sp.cos(x))
```

#### You can evaluate/simplify equations sybolically


```python
my_equation_xyz = sp.sqrt((x * (y - 4*x)) / (z * (y - 3*x)))

my_equation_xyz
```


```python
my_equation_x = (3 * a * y) / (9 * a - y)

my_equation_x
```


```python
my_new_xyz = my_equation_xyz.subs(x, my_equation_x)

my_new_xyz
```


```python
sp.simplify(my_new_xyz)
```

# System of equations

$$
\begin{array}{c}
9x - 2y = 5 \\
-2x + 6y = 10 \\
\end{array}
\hspace{3cm}
\left[
\begin{array}{cc}
9 & -2 \\
-2 & 6 \\
\end{array}
\right]
\left[
\begin{array}{c}
x\\
y
\end{array}
\right]
=
\left[
\begin{array}{c}
5\\
10
\end{array}
\right]
$$


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

b_matrix = sp.Matrix([[5],
                      [10]])
```


```python
a_matrix, b_matrix
```


```python
a_matrix.inv()
```


```python
a_matrix.inv() * a_matrix
```


```python
a_matrix.inv() * b_matrix
```

# Solving equations - `solve`


```python
yet_another_equation = x**3 + x + 10

yet_another_equation
```


```python
sp.solve(yet_another_equation,x)
```

#### ... complex numbers


```python
sp.I
```


```python
a_complex_number = 2 + 3 * sp.I

a_complex_number
```


```python
sp.re(a_complex_number), sp.im(a_complex_number)
```

### ... solving more symbolically


```python
symbolic_equation = a*x**2 + b*x +c

symbolic_equation
```


```python
sp.solve(symbolic_equation, x)
```

## Calculus


```python
symbolic_equation
```


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


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


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


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


```python
sp.integrate(symbolic_equation,(x,0,5)).subs({a:2, b:7, c:3}).n()
```


```python
trig_equation
```


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


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

### Ordinary differential equation - `dsolve`


```python
f = sp.Function('f')
```


```python
f(x)
```


```python
sp.Derivative(f(x),x,x)
```


```python
equation_ode = sp.Derivative(f(x), x, x) + 9*f(x)

equation_ode
```


```python
sp.dsolve(equation_ode, f(x))
```

### Limits


```python
limit_equation = (1 + (1 / x)) ** x

limit_equation
```

$$\lim _{x\to 5 }\left(1+{\frac {1}{x}}\right)^{x}$$


```python
sp.limit(limit_equation, x, 5)
```


```python
sp.limit(limit_equation, x, 5).n()
```

$$\lim _{x\to \infty }\left(1+{\frac {1}{x}}\right)^{x}$$


```python
sp.limit(limit_equation, x, sp.oo)
```


```python
sp.limit(limit_equation, x, sp.oo).n()
```

### Summation

$$ \sum{\frac {x^{a}}{a!}} $$


```python
sum_equation = x**a / sp.factorial(a)

sum_equation
```

$$ \sum _{a=0}^{3}{\frac {x^{a}}{a!}} $$


```python
sp.summation(sum_equation, [a, 0, 3])
```


```python
sp.summation(sum_equation.subs({x:1}), [a, 0, 3])
```


```python
sp.summation(sum_equation.subs({x:1}), [a, 0, 3]).n()
```

$$ \sum _{a=0}^{10}{\frac {x^{a}}{a!}} $$


```python
sp.summation(sum_equation.subs({x:1}), [a, 0, 10]).n()
```

$$ \sum _{a=0}^{\infty}{\frac {x^{a}}{a!}} $$


```python
sp.summation(sum_equation, [a, 0, sp.oo])
```

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

$$
\large y_1 = \frac{x^3}{4} - 3x
$$


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


```python
my_np_y1 = my_np_x ** 3 / 4 - 3 * my_np_x
```


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

fig.tight_layout()

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

ax.set_xlabel("This is X")
ax.set_ylabel("This is Y")

ax.plot(my_np_x, my_np_y1, color='r', marker='None', linestyle='-', linewidth=4);
```

### Fourier Series


```python
my_sp_y1 = x ** 3 / 4 - 3 * x

my_sp_y1
```


```python
my_fourier = sp.fourier_series(my_sp_y1, (x, -sp.pi, sp.pi))

my_fourier
```


```python
my_fourier.truncate(3).n(2)
```


```python
my_np_1term = -4.1 * np.sin(my_np_x)
my_np_2term = -4.1 * np.sin(my_np_x) + 0.91 * np.sin(2*my_np_x)
my_np_3term = -4.1 * np.sin(my_np_x) + 0.91 * np.sin(2*my_np_x) - 0.47 * np.sin(3*my_np_x)
```


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

fig.tight_layout()

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

ax.set_xlabel("This is X")
ax.set_ylabel("This is Y")

ax.plot(my_np_x, my_np_y1, color='r', marker='None', linestyle='-', linewidth=8)

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);
```

### Another Function

$$
\large y_2 = 2\,\sin(5x) \ e^{-x}
$$


```python
my_np_y2 = 2 * np.sin(5 * my_np_x) * np.exp(-my_np_x)
```


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

fig.tight_layout()

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

ax.set_xlabel("This is X")
ax.set_ylabel("This is Y")

ax.plot(my_np_x, my_np_y2, color='r', marker='None', linestyle='-', linewidth=4);
```

### Taylor Expansions


```python
my_sp_y2 = 2 * sp.sin(5 * x) * sp.exp(-x)
my_sp_y2
```


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

my_taylor
```


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


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

## General Equation Solving - `nsolve`

$$
\large y_1 = \frac{x^3}{4} - 3x\\
\large y_2 = 2\,\sin(5x) \ e^{-x}
$$

### Where do they cross? - The graph


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

fig.tight_layout()

ax.set_ylim(-5,5)
ax.set_xlim(-np.pi,4)

ax.set_xlabel("This is X")
ax.set_ylabel("This is Y")

ax.plot(my_np_x, my_np_y1, color='b', marker='None', linestyle='--', linewidth = 4)
ax.plot(my_np_x, my_np_y2, color='r', marker='None', linestyle='-', linewidth = 4);
```

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


```python
my_sp_y1, my_sp_y2
```


```python
my_guess = 3.3

sp.nsolve(my_sp_y1 - my_sp_y2, x, my_guess)
```


```python
all_guesses = (3.3, 0, -0.75)

for val in all_guesses:
    result = sp.nsolve(my_sp_y1 - my_sp_y2, 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_y1 - my_sp_y2, x, my_guess)
```

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

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

6e6cc23e0bd9acb5a0530935a31ede1834c2b51d
ipynb
Jupyter Notebook
Python_SymPy.ipynb
UWashington-Astro300/Astro300-A18
f03b1967195aae65f46729aa285948d1e8648a0c
MIT
Python_SymPy.ipynb
UWashington-Astro300/Astro300-A18
f03b1967195aae65f46729aa285948d1e8648a0c
MIT
Python_SymPy.ipynb
UWashington-Astro300/Astro300-A18
f03b1967195aae65f46729aa285948d1e8648a0c
MIT
Qwen/Qwen-72B
1. YES
2. YES