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.
#(Note: The html font may not show the appropriate Python spacing)
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(1,N-1): ## all but the first one!!!
#Estimate the acceleration half way between points i and i+1
jerk = (a[i]-a[i-1])/dt
a_avg = a[i] + jerk*(dt/2)
# Kinematics
v[i+1] = v[i] + a_avg * dt
x[i+1] = x[i] + v[i] * dt + 0.5 * a_avg * dt**2
#Newton's 2nd Law
vt = terminal_velocity(t[i+1])
a[i] = acceleration(v[i+1],vt)
# Now I set the height of the plane to be the final position, and print it out. This is because I defined down to the positive in the program, and zero position as the plane itself.
# 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, just to show you how to plot things. Also notice how we can 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.
USAYPT Problem Master 