# 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. 


```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$$


```python
x-x
```

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


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




    0



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**


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

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


```python
x-x
```




$\displaystyle 0$



### 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**


```python
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:


```python
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$.


```python
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.


```python
sym.init_printing()
```


```python
expr
```

**Now we expand the function algebreically**


```python
expr.expand()
```

Sympy can also factor the equation


```python
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.


```python
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


```python
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$$


```python
# Your code here
w,z=sym.symbols('w,z')
expr=(2*w-3*z)*(2*w+3*z)
print(expr,'=',expr.expand())
print("this shows that the first equation is true")
expr2=(2*w-3*z)**2
print(expr2,'=',expr2.expand())
print("this shows that the second equation is true, because (2*w-3*z)**2 is not equal to (4*w**2 - 9*z**2)")
```

    (2*w - 3*z)*(2*w + 3*z) = 4*w**2 - 9*z**2
    this shows that the first equation is true
    (2*w - 3*z)**2 = 4*w**2 - 12*w*z + 9*z**2
    this shows that the second equation is true, because (2*w-3*z)**2 is not equal to (4*w**2 - 9*z**2)


## 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$$


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


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

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


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


```python
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)}



```python
solutions
```

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


```python
# 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.


```python
# Your code here
#kc=(products/reactants)
#8.05=(x^2/(0.250-x)). Use ICE table to solve for this
#0=(x^2/(0.250-x)-8.05)
sym.solveset?
x = sym.symbols('x')
sym.solveset((((x**2)/(0.250-x))-8.05),x)
```

Why did you pick your answer?

The concentration should be positive, and so the answer is that each species is ~0.24M. 

## 3. Calculus

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

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


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

Similarly we can integrate:


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

We can also carry out definite integrals:


```python
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$


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


```python
from sympy import *
sym.integrate(exp(2*x), (x, 0, 5))
```

### Exercise 3.2
Solve the following definite 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`


```python
# Your code here
from sympy import *
σ,μ,x=sym.symbols('σ,μ,x')
sym.integrate((1/(σ*sqrt(2*sym.pi)))*exp(-(1/2)*((x-μ)/σ)**2), (x, -oo, oo))
```




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



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$:


```python
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}}}}$


```python
# Your code here
expr = x + cos(x)
p = sym.plot(expr)
```


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

# Lecture

## L1. Hydrogen Atom

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


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

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

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


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

We can use a loop to print out many energies


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


```python
r = sympy.symbols('r')
sympy.physics.hydrogen.R_nl(1, 0, r, 1)
#principle number, angular momentum, variable, Z (charge)
```


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

And the probablity distribution function


```python
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
```


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

Plot a 2s orbital


```python
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
```


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

We can change the range to see the node better.


```python
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. 


```python
r = sympy.symbols('r')
sympy.physics.hydrogen.R_nl(2, 1, r, 1)
#principle number, angular momentum, variable, Z (charge)
```


```python
sym.plot(sympy.physics.hydrogen.R_nl(2, 1, r, 1),(r,0,15))
```


```python
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
```


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


```python
r = sympy.symbols('r')
sympy.physics.hydrogen.R_nl(3, 0, r, 1)
#principle number, angular momentum, variable, Z (charge)
```


```python
sym.plot(sympy.physics.hydrogen.R_nl(3, 0, r, 1),(r,0,20))
```


```python
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
```


```python
sym.plot(prob_3s,(r,1,20))
```


```python
r = sympy.symbols('r')
sympy.physics.hydrogen.R_nl(4, 0, r, 1)
#principle number, angular momentum, variable, Z (charge)
```


```python
sym.plot(sympy.physics.hydrogen.R_nl(4, 0, r, 1),(r,0,25))
```


```python
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
```


```python
sym.plot(prob_4s,(r,1,25))
```


```python
r = sympy.symbols('r')
sympy.physics.hydrogen.R_nl(3, 2, r, 1)
#principle number, angular momentum, variable, Z (charge)
```


```python
sym.plot(sympy.physics.hydrogen.R_nl(3, 2, r, 1),(r,0,25))
```


```python
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
```




$\displaystyle \frac{8 r^{4} e^{- \frac{2 r}{3}}}{98415}$




```python
sym.plot(prob_3d,(r,0,25))
```


```python

```

1248163b29cd1049392c4b29b98859390a300fd2
ipynb
Jupyter Notebook
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-daliahassan98
1aff940ba2d3075e006242b57512d953655b563e
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-daliahassan98
1aff940ba2d3075e006242b57512d953655b563e
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-daliahassan98
1aff940ba2d3075e006242b57512d953655b563e
MIT
Qwen/Qwen-72B
1. YES
2. YES