# SymPy demo

If you need to install SymPy:

    pip install sympy


```python
import sympy
sympy.__version__
```




    '1.6.2'




```python
sympy.init_printing(use_latex='mathjax')
```

## Rearrange and simplify equations

We have an expression for Vp in terms of Young's modulus, $E$, and shear modulus, $\mu$:

$$ V_\mathrm{P} = \sqrt{\frac{\mu\,(E-4\mu)}{\rho\,(E-3\mu)}} $$

We need single symbols for quantities, so I'll use $\alpha$ for $V_\mathrm{P}$ and $\beta$ for $V_\mathrm{S}$, and $\gamma$ for their ratio.


```python
alpha, beta, gamma = sympy.symbols("alpha, beta, gamma")
lamda, mu, E, K, M, rho = sympy.symbols("lamda, mu, E, K, M, rho")
```

Now we can use these symbols to make an expression.


```python
from sympy import sqrt

alpha_expr = sqrt((mu * (E - 4*mu)) / (rho * (E - 3*mu)))
alpha_expr
```




$\displaystyle \sqrt{\frac{\mu \left(E - 4 \mu\right)}{\rho \left(E - 3 \mu\right)}}$




```python
print(sympy.latex(alpha_expr))
```

    \sqrt{\frac{\mu \left(E - 4 \mu\right)}{\rho \left(E - 3 \mu\right)}}


We also know that

$$ \mu = \frac{3KE}{9K-E} $$


```python
mu_expr = (3 * K * E) / (9 * K - E)
```

Now we can substitute this into the first expression.


```python
subs = alpha_expr.subs(mu, mu_expr)
subs
```




$\displaystyle \sqrt{3} \sqrt{\frac{E K \left(- \frac{12 E K}{- E + 9 K} + E\right)}{\rho \left(- E + 9 K\right) \left(- \frac{9 E K}{- E + 9 K} + E\right)}}$



This is a bit ugly! Let's simplify it:


```python
from sympy import simplify

simplify(subs)
```




$\displaystyle \sqrt{3} \sqrt{\frac{K \left(E + 3 K\right)}{\rho \left(- E + 9 K\right)}}$



We can get this as LaTeX plain-text if we want.


```python
print(sympy.latex(simplify(subs)))
```

    \sqrt{3} \sqrt{\frac{K \left(E + 3 K\right)}{\rho \left(- E + 9 K\right)}}


## Solve an equation

We'll solve

 $$ x + 2y = 3 $$
 $$ 3x + 4y = 17 $$ 
 
We'll re-write this using $x_1$ and $x_2$ instead of $x$ and $y$. They are just the names for the unknowns, but it's easier if we just think of there being one, multi-valued unknown. 

 $$ x_1 + 2x_2 = 3 $$
 $$ 3x_1 + 4x_2 = 17 $$ 
 
Multi-valued quantities are _vectors_, usually written in bold face: $\mathbf{x} = [x_1, x_2]$.
 
Now we can rewrite this in the form $\mathbf{A}\mathbf{x} = \mathbf{b}$ (which is analogous to the $\mathbf{G}\mathbf{m} = \mathbf{d}$ form of many geophysical problems). $\mathbf{A}$ is the matrix containing all the parameters, or coefficients, of the variables (the unknowns in this case) on the left-hand side of the equation. And $\mathbf{b}$ is a vector containing the known outputs &mdash; the right-hand side of the equation.
 
 $$  \mathbf{A}\mathbf{x} = \mathbf{b} $$
 
 $$ \begin{bmatrix} 1, 2 \\ 3, 4 \end{bmatrix}  \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} =  \begin{bmatrix} 3 \\ 17 \end{bmatrix}$$

To multiply a matrix by a vector, we multiply the first row by the vector and add the result. So the first row yields $(1 x_1 + 2 x_2)$ &mdash; the first of the two equations we started with. And the second row gives us $(3 x_1 + 4 x_2)$ &mdash; the second equation.

So now 'all' we have to do is find the vector $\mathbf{x}$ that satisfies this new algebraic equation:

 $$ \mathbf{A}\mathbf{x} = \mathbf{b} $$
 $$ \Rightarrow \mathbf{x} = \mathbf{A}^{-1} \mathbf{b}$$
 
There's one catch... That symbol $\mathrm{A}^{-1}$ doesn't mean the reciprocal. It means the inverse. And that's where the fun starts.

We won't go into it now, but the inverse can be hard to compute. Sometimes it's impossible. So mathematicians have come up with lots of other ways to solve this kind of equation. Welcome to the world of **linear algebra**!

### Solve with SymPy


```python
from sympy.solvers import solve
from sympy import symbols

# Define the symbols we are going to use.
x_1, x_2 = symbols('x_1, x_2')

# Define the equations, making them equal zero.
equations = [x_1 + 2*x_2 - 3,
             3*x_1 + 4*x_2 - 17]

# Solve for x_1 and x_2.
solve(equations, (x_1, x_2))
```




$\displaystyle \left\{ x_{1} : 11, \  x_{2} : -4\right\}$



We can also solve it as a linear system.

This requires us to formulate the equations as a single matrix, called an _augmented matrix_. It might look a bit funny, but it's a standard way to solve this kind of problem:


```python
from sympy import Matrix, solve_linear_system

x_1, x_2 = symbols('x_1, x_2')

system = Matrix([[1, 2,  3],
                 [3, 4, 17]])

solve_linear_system(system, x_1, x_2)
```




$\displaystyle \left\{ x_{1} : 11, \  x_{2} : -4\right\}$



### Solve with `np.linalg`

We can solve linear systems without SymPy.


```python
import numpy as np

A = np.array([[1, 2],
              [3, 4]])

b = np.array([3, 17])
```


```python
x = np.linalg.inv(A) @ b

x
```




    array([11., -4.])



We can check that this is actually a solution:


```python
A @ x
```




    array([ 3., 17.])



Indeed, this is the `b` we started with. However, not all matrices are invertible, so this method won't always work. Sometimes we need `solve`:


```python
np.linalg.solve(A, b)
```




    array([11., -4.])



Least squares is another option:


```python
x, *_ = np.linalg.lstsq(A, b, rcond=-1)

x
```




    array([11., -4.])



----

&copy; 2020 Agile Scientific &mdash; licensed CC-BY

78b473dcd0fac3e17fcc8e508f65fa9994282993
ipynb
Jupyter Notebook
SymPy_demo.ipynb
EvanBianco/geocomputing_demos
96dd8e0d7293077d39099d13f102709fb6812961
Apache-2.0
2021-07-29T15:17:08.000Z
2021-07-29T15:17:15.000Z
SymPy_demo.ipynb
EvanBianco/geocomputing_demos
96dd8e0d7293077d39099d13f102709fb6812961
Apache-2.0
SymPy_demo.ipynb
EvanBianco/geocomputing_demos
96dd8e0d7293077d39099d13f102709fb6812961
Apache-2.0
2021-07-30T18:03:32.000Z
2021-07-30T18:03:32.000Z
Qwen/Qwen-72B
1. YES
2. YES