# 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)
sym.init_printing()
expr
```


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


```python
expr = (2*w-3*z)**2
sym.init_printing()
expr
```


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


```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
x = sym.symbols('x')

expr=((x*x)/(0.250-x))-8.05
expr
```


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

    {-8.29268379803378, 0.242683798033777}


Why did you pick your answer?


```python
#concentrations cannot be negative
#[CO]=[Cl2]=0.243 M
#[COCl2]=0.25 M-0.243 M
```



## 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
print("1.")
x,y,xy = sym.symbols('x,y,xy')

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


```python
print("2.")

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`


```python
# Your code here
o,u,oo = sym.symbols('o,u,oo')

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

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

### 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
print("1.")
expr = x+sym.cos(x)
p = sym.plot(expr)
```


```python
print("2.")
expr = (1/sym.sqrt(2*sym.pi))*sym.E**(-(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)
```


```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
# Your code here
print("2p orbtial")
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
print("3s orbtial")
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,8))
```


```python
print("4s orbtial")
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,8))
```


```python
print("3d orbtial")
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
```


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


```python

```

dadb124184ac40e4fc77f5c46f5daee002a8e7af
ipynb
Jupyter Notebook
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-m-jacobo
a90a9d42b9be1deb936d6ecfec312b5bca2e018e
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-m-jacobo
a90a9d42b9be1deb936d6ecfec312b5bca2e018e
MIT
symbolic_math.ipynb
sju-chem264-2019/9-26-2019-symbolic-math-m-jacobo
a90a9d42b9be1deb936d6ecfec312b5bca2e018e
MIT
Qwen/Qwen-72B
1. YES
2. YES