# SymPy


```python
import sympy as sp
sp.init_printing(use_unicode=True)
from scipy import constants as cons # physical constants
import numpy as np
import matplotlib.pyplot as plt
```

## Basics

Initialization


```python
x, y, z = sp.symbols('x y z')
```

Definition of Sympy symbols (variables).


```python
expr1 = (2*x + 3) * (x - 6)
expr1b = (x - 6) * (2*x + 3)
expr2 = 2*x**2 -9*x - 18
```


```python
expr1 == expr1b, expr1 == expr2
```




    (True, False)



The first 2 expressions are identical, despite different term order.
The 3rd expression is syntactically different, but mathematically identical.


```python
expr1 - expr2
```


```python
sp.simplify(expr1 - expr2)
```

Simplify reveals that expr1 and expr2 are identical.


```python
expr1.subs(x, sp.pi)
```


```python
sp.simplify(expr1.subs(x, expr1))
```

Substitution of variables


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


```python
sp.integrate(expr1, (x, 2, 10))
```

Integration without and with given interval


```python
1/2, sp.S(1) / 2
```

Defition of numbers as symbols

## Example 1: Kinetic Energy in Special Relativity


```python
beta, c, m = sp.symbols(['beta', 'c', 'm'])
```


```python
gamma = 1/sp.sqrt(1-beta**2)
gamma
```


```python
Ekin = (gamma - 1)*m*c**2
Ekin
```


```python
Ekin.subs(c, cons.c).subs(beta, 0.1).subs(m, 1)
```

Kinetic energy (in J) of 1kg moving at 10% speed of light.


```python
v = beta * c
Ekin_nr = sp.S(1)/2*m*v**2
Ekin_nr
```

Non-relativistic formula for kinetic energy.


```python
Ekin.series(beta)
```

The non-relativistic formula is actually the first term of the Taylor expansion of the relativistic kinetic energy equation.


```python
err = sp.simplify(Ekin / Ekin_nr)
err
```

Ratio of actual kinetic energy to non-relativistic approximation

## Plotting

### Plotting with Sympy


```python
sp.plotting.plot(err, (beta, 0, 0.9))
```

### Plotting using Numpy


```python
f = sp.lambdify(beta, err) # converts Sympy expression into function
f
```




    <function _lambdifygenerated(beta)>




```python
x = np.linspace(0.0001, 1, 100)
```


```python
plt.plot(x, f(x))
```

## Example 2: Tsiolkovsky's Rocket Equation

https://en.wikipedia.org/wiki/Tsiolkovsky_rocket_equation

### Derivation of the Equation from Newton's Law


```python
mdot, ve, m0, mdry, t, Dv = sp.symbols('mdot v_e m_0 m_dry t \Delta{V}', 
                                positive=True, real=True) # assumptions given here
                                            # for formula simplification by Sympy
```


```python
F_exh = mdot * ve # thrust force is mass flow times exhaust velocity
F_exh
```


```python
m = m0 - mdot * t # current mass
m
```


```python
a = F_exh / m # current acceleration
a
```


```python
t_max = (m0 - mdry) / mdot # time the thrust is applied until reaction mass is consumed
t_max
```


```python
Dv = sp.integrate(a, (t, 0, t_max)).simplify()
Dv
```


```python
sp.expand_log(Dv)
```

One can trivially reformulate this equation into the textbook one using

$$
    \log{a} - log{b} = \log{\frac{a}{b}}
$$


```python
Dv_textbook = ve * sp.log(m0/mdry)
Dv_textbook
```

### Calculating Mass Ratio for Earth Orbit


```python
v_leo = 9700 # low earth orbit, in m/s
ve_HO = 4400 # exhaust velocity for liquid hydrogen+ oxygen rockets
```


```python
m0_result = sp.solve(Dv.subs(mdry, 1).subs(ve, ve_HO) - v_leo, m0)
m0_result[0].evalf() # the solver returns a list of results, even if it is only 1 element like here
```

According to the rocket equation, a 1 stage hydrogen/oxygen rocket must must carry more than 8 times the dry mass (i.e. payload, tanks, engines, etc.) in fuel to reach low-earth orbit.

In practice this ratio is even higher because air friction and gravitational pull has been neglected here.

a25e248ca7f938f833e43a802d741adaa4b83008
ipynb
Jupyter Notebook
notebooks/symbolic_computation.ipynb
lungben/python_tutorial
b5cab0bee30cdebe6db2d671cce0c9230896b402
MIT
notebooks/symbolic_computation.ipynb
lungben/python_tutorial
b5cab0bee30cdebe6db2d671cce0c9230896b402
MIT
2019-07-30T16:59:31.000Z
2019-07-30T16:59:31.000Z
notebooks/symbolic_computation.ipynb
lungben/python_tutorial
b5cab0bee30cdebe6db2d671cce0c9230896b402
MIT
2019-12-25T14:41:27.000Z
2019-12-25T14:41:27.000Z
Qwen/Qwen-72B
1. YES
2. YES