# Modeling and Simulation in Python


Copyright 2017 Allen Downey

License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)



```python
import numpy as np

from sympy import *

init_printing()
```

This notebook uses SymPy to derive the equations of motion for a springy pendulum and a rigid pendulum.

Here are the symbols we need:


```python
x, y, vx, vy, ax, ay = symbols('x, y, vx, vy, ax, ay')
length0, k, m, g, R, t = symbols('length0, k, m, g, R, t')
```

We'll use vectors P, V, and A to represent position, velocity and acceleration.  The Vector class in modsim.py doesn't play nicely with SymPy, so I'll use NumPy arrays instead:


```python
P = np.array([x, y])
P
```




    array([x, y], dtype=object)




```python
V = np.array([vx, vy])
V
```




    array([vx, vy], dtype=object)




```python
A = np.array([ax, ay])
A
```




    array([ax, ay], dtype=object)



The only vector operations we need are `mag` and `hat`:


```python
def mag(P):
    return sqrt(P[0]**2 + P[1]**2)
```




```python
mag(P)
```


```python
def hat(P):
    return P / mag(P)
```


```python
hat(P)
```




    array([x/sqrt(x**2 + y**2), y/sqrt(x**2 + y**2)], dtype=object)



For convenience, I'll define intermediate variables like length:


```python
length = mag(P)
length
```

`f_spring` is the force on the particle due to the spring


```python
f_spring = -k * (length - length0) * hat(P)
f_spring
```




    array([-k*x*(-length0 + sqrt(x**2 + y**2))/sqrt(x**2 + y**2),
           -k*y*(-length0 + sqrt(x**2 + y**2))/sqrt(x**2 + y**2)], dtype=object)



xhat and yhat are unit vectors along the x and y axes: 


```python
xhat = np.array([1, 0])
yhat = np.array([0, 1])
```

Now I can write force due to gravity as a Vector


```python
f_grav = -m * g * yhat
f_grav
```




    array([0, -g*m], dtype=object)



To write $\sum F = ma$, I'll define the left-hand side and right-hand sides of the equations separately.


```python
lhs = f_spring + f_grav
```


```python
rhs = m * A
```

Now I can make two equations, one for each component of F and A:


```python
eq1 = Eq(lhs[0], rhs[0])
eq1
```


```python
eq2 = Eq(lhs[1], rhs[1])
eq2
```

Now I want equations that are explicit in ax and ay.  In this case I can get them easily by dividing through by m.  But for the rigid pendulum we will need to use solve, so I want to demonstrate that here.


```python
soln = solve([eq1, eq2], [ax, ay])
```

Now we can extract the expressions for ax and ay


```python
soln[ax]
```


```python
soln[ay]
```

And we can get SymPy to format the result for LaTeX:


```python
print(latex(soln[ax]))
```

    \frac{k length_{0} x}{m \sqrt{x^{2} + y^{2}}} - \frac{k x}{m}



```python
print(latex(soln[ay]))
```

    - g + \frac{k length_{0} y}{m \sqrt{x^{2} + y^{2}}} - \frac{k y}{m}


Or generate Python code we can paste into a slope function:


```python
print(python(soln[ax]))
```

    k = Symbol('k')
    length0 = Symbol('length0')
    x = Symbol('x')
    m = Symbol('m')
    y = Symbol('y')
    e = k*length0*x/(m*sqrt(x**2 + y**2)) - k*x/m



```python
print(python(soln[ay]))
```

    g = Symbol('g')
    k = Symbol('k')
    length0 = Symbol('length0')
    y = Symbol('y')
    m = Symbol('m')
    x = Symbol('x')
    e = -g + k*length0*y/(m*sqrt(x**2 + y**2)) - k*y/m


To see these equations run, see pendulum.ipynb

### Rigid pendulum

Solving the rigid pendulum is almost the same, except we need a third equation to represent the geometric constraint.  The simplest form of the constraint is:

$ x^2 + y ^2 = length$

But this equation doesn't involve vx, vy, ax, and ay, so it's not much help.  However, if we take the time derivative of both sides, we get

$ 2 x vx + 2 y vy = 0 $

And if we take the time derivative one more time (and divide through by 2), we get

$ x ax + y ay + vx^2 + vy^2 = 0 $

And that's just what we need.


```python
eq3 = Eq(x*ax + y*ay + vx**2 + vy**2, 0)
eq3
```

Now we can represent the force due to tension as a vector with unknown magnitude, R, and direction opposite P. 


```python
f_tension = -R * hat(P)
```

Again, we can write $\sum F = ma$


```python
lhs = f_grav + f_tension
```


```python
rhs = m * A
```

And make one equation for each dimension:


```python
eq4 = Eq(lhs[0], rhs[0])
eq4
```




```python
eq5 = Eq(lhs[1], rhs[1])
eq5
```

Now we have three equations in three unknowns:


```python
soln = solve([eq3, eq4, eq5], [ax, ay, R])
```

And we can get explicit expressions for ax and ay


```python
soln[ax]
```




```python
soln[ay]
```

Again, we can get the results in LaTeX or Python.


```python
print(latex(soln[ax]))
```

    \frac{x}{x^{2} + y^{2}} \left(g y - vx^{2} - vy^{2}\right)



```python
print(latex(soln[ay]))
```

    - \frac{1}{x^{2} + y^{2}} \left(g x^{2} + y \left(vx^{2} + vy^{2}\right)\right)


To see these equations run, see pendulum2.ipynb


```python
print(python(soln[ax]))
```

    x = Symbol('x')
    g = Symbol('g')
    y = Symbol('y')
    vx = Symbol('vx')
    vy = Symbol('vy')
    e = x*(g*y - vx**2 - vy**2)/(x**2 + y**2)



```python
print(python(soln[ay]))
```

    g = Symbol('g')
    x = Symbol('x')
    y = Symbol('y')
    vx = Symbol('vx')
    vy = Symbol('vy')
    e = -(g*x**2 + y*(vx**2 + vy**2))/(x**2 + y**2)



```python

```

2c029682db528178b3b9ac4e00de2780ce7a022b
ipynb
Jupyter Notebook
code/pendulum_sympy.ipynb
leilamerz/ModSimPy
b877b6053a461c179643f55d0f1c1d929af0aef8
MIT
2019-04-27T22:43:12.000Z
2019-11-11T15:12:23.000Z
code/pendulum_sympy.ipynb
leilamerz/ModSimPy
b877b6053a461c179643f55d0f1c1d929af0aef8
MIT
2019-10-09T18:50:22.000Z
2022-03-21T01:39:48.000Z
code/pendulum_sympy.ipynb
leilamerz/ModSimPy
b877b6053a461c179643f55d0f1c1d929af0aef8
MIT
2021-04-13T01:10:41.000Z
2021-04-13T01:10:41.000Z
Qwen/Qwen-72B
1. YES
2. YES