# Introduction to Computational Dynamics¶

In this notebook we'll do some actual computer programming, using functions and loops to compute the motion of an object subjected to a force.

If you're new to programming (or just new to Python), you should first review the notebook titled "Functional Programming".

## A. Set up¶

You will need to import the following packages.

```
from numpy import *
import matplotlib.pyplot as plt
```

## B. Theory¶

A complex physical system may be converted to a computational Physics problem when a direct solution is difficult or impossible. With a computational model, a computer can work out how the system evolves by adding together many small time steps.

To simulate the continuous functions $x(t)$ and $v(t)$, we will convert integral equations, such as

$$x(t) = x_0 + \int_{0}^t v(t) \; dt$$

$$v(t) = v_0 + \int_{0}^t a(t) \; dt$$

into approximate finite difference equations, such as

$$x_{i+1} = x_i + v_i \Delta t$$

$$v_{i+1} = v_i + a_i \Delta t$$

where $i$ labels the many small time steps taken by the computer. The result is a sequence of $x$ positions $\left \{ x_0,x_1,x_2,... \right \}$ at times $\left \{t_0, t_1,t_2,... \right \}$.

At each time step $i$, the current value of the force is found, which may depend on the current position and velocity. From this, the current value of the acceleration is obtained using Newton's Second Law: $a_i = F/m$. With acceleration $a_i$ we find the next velocity $v_{i+1}$. With the velocity $v_i$ we find the next position $x_{i+1}$. This process is then repeated in a loop for as many time steps as we need.

For motion in 2 or 3 dimensions, we repeat this independently for

- $x$, $v_x$, $a_x$ and $F_x$
- $y$, $v_y$, $a_y$ and $F_y$
- $z$, $v_z$, $a_z$ and $F_z$

In the sections below, we find the position of a body as follows:

**Write a force function**: Program a function that returns the force on the body. It may depend on the body's location, its mass, and/or on other quantities (such as velocity or some material properties).**Define the time domain**: initial time $t_i$, final time $t_f$, and the time interval $\Delta t$. The time interval should be small enough that the velocity of the body does not change significantly over any step from $t_i$ to $t_{i+1}$.**Define the initial conditions**: provide values for the initial position and velocity of the body ($\vec{x}_0$ and $\vec{v}_0$) and the body's mass ($m$).**Perform a computational loop**over many small time steps. At each new time, find the new force, then the new acceleration, then the new velocity and position. Repeat the loop every $\Delta t$ until you arrive at the final time of interest.

If you save the position and velocity every time step, you may then plot the trajectory or otherwise study the motion.

## C. Projectile Motion Example¶

Our first computational model is the motion of an object of mass $m$ subject to a uniform gravity force, $\vec{F} = -mg \hat{\jmath}$.

Let's follow the procedure listed above for a 1 kg body launched from the origin with initial speed and direction $v_0=10$ m/s and $\theta_0=60^\circ$.

**1. Write a force function.**

This force is very simple since it is the same at every point in space, and does not change in time.

```
g = 9.8 # N/kg # define constant gravity field
jhat = array((0,1)) # define constant vertical unit vector
def gravity_force(m):
Fg = -m*g*jhat
return Fg
```

**2. Define time domain.**

```
t0,tf = 0, 2.0 # initial and final time (seconds)
dt = 0.01 # delta t (seconds)
```

**3. Define intial conditions.**

```
v0 = 10 # initial speed (m/s)
Q0 = 60 # initial direction (degrees)
rpd = pi/180 # radians per degree conversion factor
m = 1.0 # mass of body (kg)
r = array(( 0, 0 )) # initial position vector (m)
v = v0*array(( cos(Q0*rpd), sin(Q0*rpd) )) # initial velocity vector (m/s)
t = t0 # initial time
```

We want to keep track of the position $\vec{r}$ and the velocity $\vec{v}$ at
all times. One way to do this is with a
Python list, which uses the square
brackets. We'll make a list of vectors, adding a new one at each time step in a
loop. We'll use `r_list`

, `v_list`

and `t_list`

for these list
quantities.

```
r_list = [r,] # start a list of position vectors
v_list = [v,] # start a list of velocity vectors
t_list = [t,] # start a list of times
```

**4. Perform the loop.**

There are several ways to proceed. Here I've used a `while`

loop:

```
while t <= tf: # keep repeating until time tf
F = gravity_force(m) # get current force
a = F/m # current acceleration
t = t + dt # increment the time
r = r + v*dt # increment the position
v = v + a*dt # increment the velocity
r_list.append(r) # append the new r to list r_list
v_list.append(v) # .. v v_list
t_list.append(t) # .. t t_list
```

All the data are now stored in lists of vectors, `r_list`

, `v_list`

and `t_list`

.

It is more convenient to use a two-dimensional array for these quantities, where each row is a different time, and the columns hold the $x$ and $y$ components. For example,

$$\texttt{rr} = \begin{bmatrix} x_0 & y_0 \\ x_1 & y_1 \\ x_2 & y_2 \\ ... & \end{bmatrix}$$

To access the $x$ and $y$ position at time step 5, we use `rr[5,:]`

(row 5, the
`:`

means all columns). To get the entire time-history of the $x$ component we
use `rr[:,0]`

(the `:`

means all rows, the `0`

means column `0`

which is the
$x$ data).

To convert our lists of vectors into two-dimensional arrays:

```
rr = array(r_list) # for convenience later, convert the
vv = array(v_list) # lists into arrays
tt = array(t_list) #
```

That's it. Now `rr`

and `vv`

contain the position and velocity of the object at
all times in `tt`

. It is often helpful to graph the results.

**Graph the result.**

Here is a graph of $y$ versus $x$ position, colored accoring to time $t$. The
`matplotlib`

plot method we use here is `scatter`

. You can see its options in
the `scatter`

documentation.

```
x = rr[:,0] # take the first column for x
y = rr[:,1] # take the second column for y
fig,ax = plt.subplots() # start making a plot
ax.scatter(x, y, c=tt, ) # color the dots by time
ax.grid(True) # show a grid on the plot
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('projectile in a vacuum')
ax.set_aspect('equal') # make x- and y-distance the same on screen
plt.show()
```

You should see a perfect parabola. The computer doesn't know where the ground is so it keeps accelerating downward until the final time is reached (this is fixed below).

## D. A Loop in a Loop¶

Now that we have the basic code, we can execute it over and over again in a variety of circumstances.

For example, perhaps we'd like to compare the behavior of projectiles on
different planets with varying gravitational constants. On the moon $g=1.6
\;\rm{N}/\rm{kg}$; on Jupiter $g=23.5 \;\rm{N}/\rm{kg}$. The code above has
been written a bit more compactly and put in a loop where `g`

is given a few
different values. In this case, we prepare the plot first, then perform the
work above for each `g`

.

```
def gravity_force(m):
Fg = -m*g*array((0,1))
return Fg
fig,ax = plt.subplots() # start the plot
ax.grid(True)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('projectiles, changing g')
ax.set_aspect('equal')
m = 1.0 # mass of body (kg)
t0,tf = 0, 2.0 # initial and final time (seconds)
dt = 0.01 # delta t (seconds)
r0 = array(( 0, 0 )) # initial position vector (m)
Q0 = 60 # initial direction (degrees)
v0 = 10*array(( cos(Q0*pi/180), sin(Q0*pi/180) )) # initial velocity (m/s)
for g in [1.6,9.8,23.5]:
r,v,t = r0,v0,t0 # set initial position, velocity and time
r_list,v_list,t_list = [r],[v],[t] # start lists
while t <= tf: # keep repeating until time tf
F = gravity_force(m) # get current force
a = F/m # current acceleration
t = t + dt # increment the time
r = r + v*dt # increment the position
v = v + a*dt # increment the velocity
r_list.append(r) # append r,v,t to lists
v_list.append(v) #
t_list.append(t) #
if r[1] < 0: break # stop the while-loop if y<0 (hit the ground)
rr = array(r_list) # for convenience later, convert the
vv = array(v_list) # lists into arrays
tt = array(t_list) #
# print and plot results:
print(f'projectile with g = {g:.1f} N/kg:')
print(f' range {rr[-1,0]:0.2f} m, flight time {t:0.2f} s')
ax.plot( rr[:,0], rr[:,1], label='{:.1f}'.format(g) )
ax.legend(loc='upper left') # show a legend on the plot
plt.show()
```

projectile with g = 1.6 N/kg: range 10.00 m, flight time 2.00 s projectile with g = 9.8 N/kg: range 8.90 m, flight time 1.78 s projectile with g = 23.5 N/kg: range 3.75 m, flight time 0.75 s

Note the range is incorrect for the moon projectile ($g=1.6 \;\rm{N}/\rm{kg}$) since it never got a chance to land.

## E. Other Examples¶

With the methods above you can do quite a bit.

In Dynamics Examples there are a couple more demonstrations.

For example, you can make "Lissajous figures" with a 2D spring:

Or observe elliptical gravitational orbits: