# Symbolic Mathematics in Python 

There are times when you need to solve a difficult problem symbollically or analytically. If you have ever used Wolfram Alpha, then you have already done this. Sympy is a python library that allows you to do symbolic mathematics in python. 


```
import sympy as sym
```

## 1. Introduction

### Example 1.1

If you try to write the follwing in python by itself, you will get an error telling you x is undefined:
$$x-x$$


```
x-x
```

(The error above is on purpose). Variables in python need to be defined before you can say something specific about them


```
x=102
x-x
```

If you are trying to show that $x-x=0$ is true for any $x$, the above answer would not be valid. Instead you can use a symbolic expression to show that it is true

**First we define the variable as a symmbolic expression**


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


```

```

**Now we can use the variable in a symbolic expression**


```
x-x
```

### Example 1.2

Sympy can be used to perform algebreic operations (among other things). Consider the following expression: $$(3a-4b)^3$$


We can use symppy to expand the expression algebraically.

**First we need to define the variables as symbolic expressions**


```
a,b = sym.symbols('a,b')
```

**Side note** Notice that the left hand side of the epression has two variables being defined. Python can define more than one variable at a time:


```
x1,y1 =10,20
print(x1)
print(y1)
```

    10
    20


**Back to the expression** We can define an expression using the variables $a$ and $b$.


```
expr = (3*a-4*b)**3
print(expr)
```

    (3*a - 4*b)**3


We can also make it look nicer in our notebook. This doesn't affect the math, but it makes our notebook more readable.


```
sym.init_printing()
```


```
expr
```

**Now we expand the function algebreically**


```
expr.expand()
```

Sympy can also factor the equation


```
sym.factor(26*a**3-108*a**2*b+144*a*b**2-64*b**3)
```

If you want to copy and paste a result, you print the result.


```
print(sym.factor(26*a**3-108*a**2*b+144*a*b**2-64*b**3))
```

    2*(a - 2*b)*(13*a**2 - 28*a*b + 16*b**2)


You can also chain together functions


```
expr.expand().factor()
```

### Exercise 1.1

Show that the following two expressions are true.
$$(2w-3z)(2w+3z)=4w^2-9z^2$$
$$(2w-3z)^2\ne4w^2-9z^2$$


```
w,z = sym.symbols('w,z')
expr = (2*𝑤-3*𝑧)*(2*𝑤+3*𝑧)
sym.init_printing()
expr
```


```
expr.expand()
```


```
expr = (2*𝑤-3*𝑧)**2
expr
```


```
expr.expand()
```

## 2. Solving Equations

Sympy can be used to symbolilically solve equations. As before, you need to define which variables are symbols

### Example 2.1

Use sympy to solve the following equation
$$ax^3+bx^2+cx+d=0$$


```
# Define the variables
a,b,c,d,x = sym.symbols('a,b,c,d,x')
```


```
# Define the expression
expr=a*x**3+b*x**2+c*x+d
expr
```

We can use the `solvset` function to solve this equation


```
solutions=sym.solveset(expr,x)
```


```
print(solutions)
```

    {-(-3*c/a + b**2/a**2)/(3*(sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)) - (sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a), -(-3*c/a + b**2/a**2)/(3*(-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a), -(-3*c/a + b**2/a**2)/(3*(-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)) - (-1/2 + sqrt(3)*I/2)*(sqrt(-4*(-3*c/a + b**2/a**2)**3 + (27*d/a - 9*b*c/a**2 + 2*b**3/a**3)**2)/2 + 27*d/(2*a) - 9*b*c/(2*a**2) + b**3/a**3)**(1/3)/3 - b/(3*a)}



```
solutions
```

What if I need help. You can do this with any python function. `function?`


```
# Run this command to see a help box
sym.solveset?
```

### Exercise 2.1

Use the `solveset` function to solve the following chemical problem.

Phosgene gas, $\text{COCl}_2$, dissociates at high temperatures according to the following equilibrium:
    
$$ \text{COCl}_2 \rightleftharpoons  \text{CO} + \text{Cl}_2 $$

At $\text{400 C}$, the equilibrium constant $K_c=8.05$.  

If you start with a  $\text{0.250 M}$ phosgene sample at $\text{400 C}$, determine the concentrations of all species at equilibrium.


```
x = sym.symbols('x')
expr = ((x**2)/(0.250-x))-8.05
expr
```


```
solutions=sym.solveset(expr,x)
print(solutions)
```

    {-8.29268379803378, 0.242683798033777}


Why did you pick your answer?


```
# Concentrations cannot be negative
```



## 3. Calculus

We can use also Sympy to differentiate and integrate. Let us experiment with differentiating the following expression:

$$x ^ 2 - \cos(x)$$


```
sym.diff(x ** 2 - sym.cos(x), x)
```

Similarly we can integrate:


```
sym.integrate(x ** 2 - sym.cos(x), x)
```

We can also carry out definite integrals:


```
sym.integrate(x ** 2 - sym.cos(x), (x, 0, 5))
```

### Exercise 3.1

Use Sympy to calculate the following:

1. $\frac{d(x ^2 + xy - \ln(y))}{dy}$
1. $\int_0^5 e^{2x}\;dx$


```
x,y = sym.symbols('x,y')
sym.diff(x**2+x*y-sym.ln(y), y)
```


```
sym.integrate(sym.E**(2*x), (x,0,5))
```

### Exercise 3.2
Solve the following definate integral
$$\int\limits_{ - \infty }^\infty  {\frac{1}{{\sigma \sqrt {2\pi } }}{e^{ - \frac{1}{2}{{\left( {\frac{{x - \mu }}{\sigma }} \right)}^2}}}}$$

Hint, the sympy symbol for infinity is `oo`


```
from sympy import oo
from sympy import pi
from sympy import sqrt
from sympy import E
from sympy import Symbol

o,u,x = sym.symbols('o,u,x')

sym.integrate((1/(o*sqrt(2*pi))*(E**(-(1/2)*((x-u)/o)**2))), (x,-oo,oo))
```




$\displaystyle \begin{cases} 0.353553390593274 \sqrt{2} \left(2 - \operatorname{erfc}{\left(\frac{0.707106781186548 u}{o} \right)}\right) + 0.353553390593274 \sqrt{2} \operatorname{erfc}{\left(\frac{0.707106781186548 u}{o} \right)} & \text{for}\: \left(\left(2 \left|{\arg{\left(o \right)}}\right| \leq \frac{\pi}{2} \wedge \left|{4 \arg{\left(o \right)} - 2 \arg{\left(u \right)}}\right| < \pi\right) \vee \left(\left|{4 \arg{\left(o \right)} - 2 \arg{\left(u \right)}}\right| \leq \pi \wedge 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \vee \left(\left|{4 \arg{\left(o \right)} - 2 \arg{\left(u \right)}}\right| < \pi \wedge 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \vee 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \wedge \left(\left(2 \left|{\arg{\left(o \right)}}\right| \leq \frac{\pi}{2} \wedge \left|{- 4 \arg{\left(o \right)} + 2 \arg{\left(u \right)} + 2 \pi}\right| < \pi\right) \vee \left(\left|{- 4 \arg{\left(o \right)} + 2 \arg{\left(u \right)} + 2 \pi}\right| \leq \pi \wedge 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \vee \left(\left|{- 4 \arg{\left(o \right)} + 2 \arg{\left(u \right)} + 2 \pi}\right| < \pi \wedge 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \vee 2 \left|{\arg{\left(o \right)}}\right| < \frac{\pi}{2}\right) \\\int\limits_{-\infty}^{\infty} \frac{\sqrt{2} e^{- \frac{0.5 \left(- u + x\right)^{2}}{o^{2}}}}{2 \sqrt{\pi} o}\, dx & \text{otherwise} \end{cases}$




```
from sympy import oo
from sympy import pi
from sympy import sqrt
from sympy import E
from sympy import Symbol

o,u,x = sym.symbols('o,u,x')

x = Symbol('x', real=True); x
u = Symbol('u', real=True); u
o = Symbol('o', real=True); o

sym.integrate((1/(o*sqrt(2*pi))*(E**(-(1/2)*((x-u)/o)**2))), (x,-oo,oo))

```




$\displaystyle \begin{cases} 0.707106781186547 \sqrt{2} & \text{for}\: 2 \left|{\arg{\left(o \right)}}\right| \leq \frac{\pi}{2} \\\int\limits_{-\infty}^{\infty} \frac{\sqrt{2} e^{- \frac{0.5 \left(- u + x\right)^{2}}{o^{2}}}}{2 \sqrt{\pi} o}\, dx & \text{otherwise} \end{cases}$




```
expr = 0.707106781186547*sqrt(2)
a = (0.707106781186547*sqrt(2).evalf())
print (expr,'=', a, '=', a.round(1))

```

    0.707106781186547*sqrt(2) = 0.999999999999999 = 1.


Lookup Gaussian functions: https://en.wikipedia.org/wiki/Gaussian_function
Does your answer maake sense?

## 4. Plotting with Sympy

Finally Sympy can be used to plot functions. Note that this makes use of [matplotlib](http://matplotlib.org/). 

Let us plot $x^2$:


```
expr = x ** 2
p = sym.plot(expr)
```


    <Figure size 640x480 with 1 Axes>


### Exercise 4.1 Plot the following function:

1. $y=x + cos(x)$
1. ${\frac{1}{{ \sqrt {2\pi } }}{e^{ - \frac{x^2}{2}}}}$


```
expr1 = x + sym.cos(x)
p = sym.plot(expr1)
```


```
expr =  (1/sqrt(2*pi))*E**(-x**2/2)
p = sym.plot(expr)
```

# Lecture

## L1. Hydrogen Atom

Sympy has built in modules for the eigenfunctions of the hydrogen atom. 


```
import sympy.physics.hydrogen
import numpy as np
```

You can caluclate the eigenvalues ($E$) in Hartrees

`sym.physics.hydrogen.E_nl(n,Z)`


```
sym.physics.hydrogen.E_nl(1,1)
```

We can use a loop to print out many energies


```
for n in range(1,5):
    print(sym.physics.hydrogen.E_nl(n,1))
```

    -1/2
    -1/8
    -1/18
    -1/32


We can plot the hydrogen radial wavefunction (1s orbital)


```
r = sym.symbols('r')
sympy.symbols('r')
sympy.physics.hydrogen.R_nl(1, 0, r, 1)
```


```
sym.plot(sympy.physics.hydrogen.R_nl(1, 0, r, 1),(r,0,10.50))
```

And the probablity distribution function


```
sympy.symbols('r')
prob_1s=sympy.physics.hydrogen.R_nl(1, 0, r, 1)*sympy.physics.hydrogen.R_nl(1, 0, r, 1)
prob_1s
```


```
sym.plot(prob_1s,(r,0,10))
```

Plot a 2s orbital


```
sympy.symbols('r')
prob_2s=sympy.physics.hydrogen.R_nl(2, 0, r, 1)*sympy.physics.hydrogen.R_nl(2, 0, r, 1)
prob_2s
```


```
sym.plot(prob_2s,(r,0,10))
```

We can change the range to see the node better.


```
sym.plot(prob_2s,(r,1,8))
```

Notice the node!

### Exercise L1.1

Plot the radial distriubution function for a 2p, 3s, 4s, and 3d orbital. 


```
sympy.symbols('r')
prob_2p=sympy.physics.hydrogen.R_nl(2, 1, r, 1)*sympy.physics.hydrogen.R_nl(2, 1, r, 1)
prob_2p
```


```
sym.plot(prob_2p,(r,0,10))
```


```
sympy.symbols('r')
prob_3s=sympy.physics.hydrogen.R_nl(3, 0, r, 1)*sympy.physics.hydrogen.R_nl(3, 0, r, 1)
prob_3s
```


```
sym.plot(prob_3s,(r,1,8))
```


```
sympy.symbols('r')
prob_4s=sympy.physics.hydrogen.R_nl(4, 0, r, 1)*sympy.physics.hydrogen.R_nl(4, 0, r, 1)
prob_4s
```


```
sym.plot(prob_4s,(r,1,8))
```


```
sympy.symbols('r')
prob_3d=sympy.physics.hydrogen.R_nl(3, 2, r, 1)*sympy.physics.hydrogen.R_nl(3, 2, r, 1)
prob_3d
```


```
sym.plot(prob_3d,(r,0,20))
```


```

```

0c08bbd47e83fd24e2dcd03d9777982da86c21b9
ipynb
Jupyter Notebook
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-laurensruiz
79a4fd3c8ab87f8cde87a92b7738c9c5348fcc75
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-laurensruiz
79a4fd3c8ab87f8cde87a92b7738c9c5348fcc75
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-laurensruiz
79a4fd3c8ab87f8cde87a92b7738c9c5348fcc75
MIT
Qwen/Qwen-72B
1. YES
2. YES