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.

In [1]:
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 ($x_i$ and $v_i$). 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:

  1. 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).
  2. 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}$.
  3. 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$).
  4. 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.

In [2]:
g = 9.8  # N/kg         # define constant gravity field
jhat = array((0,1))     # define constant vertical unit vector
def gravity_force(m):   # function requires one input parameter (mass m)
    Fg = -m*g*jhat      #          compute force vector
    return Fg           #          return the force vector to the loop

2. Define time domain.

In [3]:
t0,tf = 0, 2.0          # initial and final time (seconds)
dt = 0.01               # delta t (seconds)

3. Define intial conditions.

In [4]:
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.

In [5]:
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:

In [6]:
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:

In [7]:
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.

In [8]:
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()
No description has been provided for this image

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.

In [9]:
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=f'{g:.1f}')
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
No description has been provided for this image

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:

2D spring

Or observe elliptical gravitational orbits:

gravitational orbit