#### Modeling Physical Systems

Modeling physical systems is simpler than you may think, so do not make it too hard. In general, you just break up a hard problem into a series of easy ones. For example, you solve most introductory problems by doing the following:

- Drawing a free body diagram of all the forces on your object of interest.
- Writing an expression for the acceleration of the object, in terms of the forces on it, using Newton’s 2nd law.
- Solving for the final kinematic conditions, such as position and velocity, as a function of initial conditions, the forces involved, and the time interval.

The problem you run into, however, is that the forces on the object are often a function of the other kinematic equations—usually the velocity or position. The trick, therefore, is to program a computer to follow the same easy method, but do it in small time steps where you can assume the acceleration is constant.

#### Python Programming

I like the Python programming language, as it is free, easy to use, and fast becoming a new standard for computational physics. I use the Anaconda interactive development environment by Continuum Analytics. The book shown is both free and an excellent introduction to using Python for computations.

#### Example **Problem:** A Skydiver

A skydiver jumps out of a plane holding his body in a *cannonball* shape (terminal velocity = 200 ft/s) for 1 minute, then he changes to a *spread eagle* shape (terminal velocity = 150 ft/s) for 2 minutes. At this time he pulls a cord that starts a series of streamers and parachutes that unfurl in 12 stages, over the course of one minute, lowering his terminal velocity by 10 ft/s every 5 seconds. After one more minute, he pulls his main parachute, which smoothly over the course of one second, lowers his terminal velocity to 15 ft/s. He then gently floats until he hits the ground 10 minutes after he jumped out of the plane.

a) What was the elevation from which the skydiver was dropped?

b) Plot the skydiver’s height above the ground, downward velocity, and downward acceleration, as functions of time.

**Solution using Python 3 code:**

First I import some packages from the numpy and matplotlib packages. Alternatively, you can use the pylab package which contains the most used versions of both.

from numpy import linspace,zeros

import matplotlib.pyplot as plt

Next, I define a function that applies Newton’s 2nd law.

def acceleration(v,vt):

'''Returns the net force per mass on the skydiver.

The acceleration is a function of the skydiver's

velocity (v) and terminal velocity (vt), and uses the

global variable g, which represents the gravitational

field of Earth. Down is positive.'''

global g

a = g * (1 - (v/vt)**2)

return a

And another function that simply defines the terminal velocity based on the problem.

def terminal_velocity(t):

'''Returns the terminal velocity in ft/s of the skydiver, as

a function of time, in seconds, based on the problem's

specifications.'''

if t < 60.0:

vt = 200.0 # cannonball

elif t < 3*60.0:

vt = 150.0 # spread eagle

elif t < 4*60.0:

dt = t - 3*60.0

vt = 150.0 - 10.0*int(1 + dt/5)

elif t < 5*60:

vt = 30 # 150 ft/s - 12 * 10 ft/s = 30 ft/s

elif t < 5*60 + 1:

dt = t - 5*60

vt = 30 - dt*15 #

else:

vt = 15.0 ## at landing

return vt

Now is the main program. First I define the initial conditions and define a set of numpy arrays that will hold the position, velocity, and acceleration. Remember that I have defined down to be positive.

g = 32 # Acceleration of gravity in ft/s**2

N = 6000 # Number of time steps

tmax = 10.0*60 # The total time flight in seconds

t = linspace(0.0,tmax,N) # This initiates an array for the time

x,v,a = zeros(N), zeros(N), zeros(N) # This initiates arrays

dt = t[1] – t[0] # We are using a constant time step

x[0] = 0.0; v[0] = 0.0; a[0] = g # sets initial conditions

Next we calculate the position, velocity, and acceleration during the first time step. Now [0] is initial and [1] is final. The hard part, however, is guessing the average acceleration between time [0] and time [1]. I first guess it is the initial acceleration. Then after calculating the final acceleration, I iterate until convergence.

a_avg = a[0] # guess the average acceleration

a_old = 10000 # give some dummy value to get into loop

while abs(a_avg – a_old) > 1E-15:

v[1] = v[0] + a_avg*dt

x[1] = x[0] + v[0]*dt + a_avg*(dt**2)/2

vt = terminal_velocity(t[1])

a[1] = acceleration(v[1],vt)

a_old = a_avg

a_avg = (v[0] + v[1])/2

Now that we know two prior accelerations, we can use linear

extrapolation to estimate the average acceleration based on the last two accelerations. Notice that we will call the rate of change

of the acceleration the “jerk.”

for i in range(2,N): ## all but the first one!!!

#Estimate the acceleration half way between points i and i+1

jerk = (a[i-1]-a[i-2])/dt

a_avg = a[i-1] + jerk*(dt/2)

# Kinematics

v[i] = v[i-1] + a_avg * dt

x[i] = x[i-1] + v[i-1] * dt + 0.5 * a_avg * dt**2

#Newton's 2nd Law

vt = terminal_velocity(t[i])

a[i] = acceleration(v[i],vt)

Now I set the height of the plane to be the final position, and print it out. Notice in Python 3 print is a function, but old Python 2 code will make it a statement, so watch out for parentheses.

h = x[-1]

print("The altitude of the plane was "+ str(int(round(h/1000))) + ",000 ft.")

Now I make a three panel plot. Pay attention to the little things I did, as I made it more complicated than I had to show you how to plot things. Also notice how we can just do regular algebra using numpy arrays.

plt.close()

plt.subplot(3,1,1)

plt.ylabel(r"$\bf \frac{a}{g}$", fontsize=10)

plt.plot(t/60,a/g)

TITLE = "A Skydiver"

plt.title(TITLE, weight='bold', fontsize = 16)

plt.xlim(0,10)

if abs(max(a)) > abs(min(a)):

plt.ylim(-max(a/g),max(a/g))

else:

plt.ylim(min(a/g),-min(a/g))

plt.subplot(3,1,2)

plt.xlim(0,10)

plt.ylabel(r”$\bf v \rm{(ft/s)}$”, fontsize=10)

plt.plot(t/60,v)

plt.subplot(3,1,3)

plt.ylabel(r”$\bf height \ \rm{(‘000 ft)}$”, fontsize=10)

plt.xlabel(r’$\bf time \ \rm{(min)}$’, fontsize=12)

plt.xlim(0,10)

plt.plot(t/60,(h-x)/1000)

plt.tight_layout()

plt.savefig(“skydiver_figure.jpg”) # saves a jpg

plt.savefig(“skydiver_figure.pdf”) # saves a pdf

plt.show()

And, that is that. Notice how most of this is simply doing a good high school physics problem, and just making the final values become the initial values after every iteration.