Scientific Programming with Python

---

Acceleration and Force

Recall the law of gravity near the Earth’s surface, describing the position of an object travelling vertically:

y = y₀ + v_y0t − ½ gt²

If we take the derivative of this expression we get one for the velocity:

v = dy/dt = v_y0 − gt

which shows how it becomes more and more negative with time in response to the downward acceleration of gravity.

If we take another derivative we get an expression for the acceleration:

a = dv/dt = −g

We can write this in vector notation using bold face to represent vector quantities, for example g = (0, −g) since it’s a constant value pointing downward:

a = dv/dt = g

Experimental examination of gravity and other forces in nature such as friction tell us that the important kinematic quantity here is the vector acceleration, which is zero except in the presence of such forces:

m a = Σ_i F_i

where m is the mass of the moving object. This is Newton’s second law of motion.

The force of gravity is therefore

F_G = m g

---

Air Resistance, Analytically

Another force experienced by an object moving through the air is air resistance, which can be approximated by

F_D = −½ ρAC_D|v|v

The minus sign indicates that the force is opposite to the direction of the velocity; ρ is the density of air, 1.225 kg/m³ at sea level, A is the projected area of the object perpendicular to the velocity (i.e. πr² for a sphere), and C_D is a dimensionless quantity known as the the drag coefficient.

An object traveling through the air therefore has an acceleration given by

m a = m g − ½ ρAC_D|v|v

From this expression we can see that a larger cross-sectional area increases the force of air resistance, so that a feather will drop more slowly than a hammer — but not on the Moon, where there is no air!

At low speeds C_D ~ b/|v|, so the magnitude of the force varies directly with speed, while at high speeds C_D is approximately constant (in the range of 0.1 – 1) and the magnitude of the force varies as the square of the speed. The former is the result of the displacement of air mass by the object with laminar (smooth) flow of air around it; the latter reflects turbulent flow produced by friction around the ball.

A ball thrown through the air, which has relatively low speed, therefore has an acceleration given by

dv/dt = g − ½ (ρAb/m)v

The factor −½ (ρAb/m) has units of inverse time, so we can define a characteristic time:

τ = 2m/(ρAb)

and

dv/dt = g − v/τ

This differential equation can be integrated relatively easily:

∫_v0^v dv/(g − v/τ) = ∫₀^t dt

−τ ln(v − gτ)|_v0^v = t

v = gτ + (v₀ − gτ)e^-t/τ

The transient exponential term will go away in several characteristic times, and the ball will eventually reach a terminal velocity vf = gτ. The greater the density of air and/or the projected area of the ball, the smaller the terminal velocity. Conversely, the greater the ball’s mass the greater its terminal velocity.

A typical value of b = 17 m/s, and for a baseball, A = π (3.7 cm)² = 0.0043 m², m = 0.15 kg, then τ = 3.4 s, so the terminal velocity is –33 m/s.

For a human in free fall, the terminal speed will be 50 m/s to 90 m/s, depending on how they hold themselves (flat or bullet-like).

We can now define a drag velocity function for any given time and initial velocity:

import numpy as np

g = 9.81
tau = 3.4
vf = g * tau

def vd(t, v0):
    ext = np.exp(-t/tau)
    (vx0, vy0) = v0
    return np.array( (t, vx0 * ext, -vf + (vy0 + vf)*ext) )

As usual, there’re a few things going on in this function:

• The argument t can be an array of values, which when divided by tau will result in an array of each element divided by tau:

np.array( (0, 1, 2) )/tau

⇒ array([ 0. , 0.29411765, 0.58823529])

Such element-by-element operations on an array are known as broadcasting.

• Similarly, we’re using the exponential function exp() from numpy, which operates element-by-element on arrays:

np.exp( -np.array( (0, 1, 2) )/tau )

⇒ array([ -1. , -1.34194177, -1.80080771])


Note: there’s also an exp() function in the math module, but it will only accept scalar (non-sequence) values.

• We are using the np.array() function to return an array of arrays, as the three calculations in the return statement are each arrays of length n, producing a 3 x n array: [ t, vx, vy ].

An array of time values can be generated with

t = np.arange(0,20,0.5)

This function is like the built-in range() but it will allow non-integer steps! More importantly, it returns an array of values, rather than a range, so we can hand it to vdt() to process :

v0 = (10, 15)
vdt = vd(t,v0)

The result is a triplet of arrays [ t, vx, vy ]:

print(vdt)

⇒
[[ 0.00000000e+00 5.00000000e-01 ...  1.85000000e+01  1.90000000e+01  1.95000000e+01]
 [ 1.00000000e+01 8.63243197e+00 ...  4.33438100e-02  3.74162491e-02  3.22993225e-02]
 [ 1.50000000e+01 8.38726154e+00 ... -3.31444153e+01 -3.31730775e+01 -3.31978199e+01]]

And we’ll plot this to see what it looks like:

def vdplot(vdt):
    from matplotlib.pyplot import plot, \
        title, xlabel, ylabel, legend, \
        show, close
    (t, vdx, vdy) = vdt
    plot(t,vdx, label='vx', linestyle='dashed')
    plot(t,vdy, label='vy')
    title('Velocity of a Ball through the Air.')
    xlabel('t')
    ylabel('v')
    legend()
    show()
    close()

vdplot(vdt)

There are new python and matplotlib features here:

• The plot() functions here uses a combination of positional and keyword arguments like label= and linestyle=.
  
  Keyword arguments don’t have to be in a particular position, unlike regular arguments like the x and y arrays, which must be first and second, respectively.
• The label= options for plot() are used by the legend() function to produce a box with the data names and styles.
• The linestyle= option lets you choose a variety of common ways to distinguish lines, such as 'solid' (the default), 'dotted', etc.
• In addition to the dashed line, matplotlib automatically assigns different colors to the different curves, but you can set your own with an option like color='red'.
• If you want to plot just the data points, use the function scatter().

---

Air Resistance, Numerically

Numerical integration of a differential equation is also possible with the SciPy package:

from scipy.integrate import quad

Here we’re importing the function quad(integrand, initial, final) which performs integration by quadrature, approximating the area under a curve by rectangles or trapezoids, as you probably learned about in calculus:

In this case we are using the integrand in the integral equation above to calculate the time t corresponding to each speed v in the y direction (Remember that in one dimension v becomes v, a signed value that is negative when an object is falling):

∫_v0^v dv/(-g − v/τ) = ∫₀^t dt = t

Coded in a function:

g = 9.81
tau = 3.4
vf = g * tau
integrand = lambda v: -1/(g*(1 + v/vf))

The keyword lambda is an alternative way to define one-line functions. It’s known as an anonymous function because no name is defined for it, though it can be assigned to a name as above.

We want to integrate this function between v0[1] and various values of v before it reaches the terminal velocity -vf = -g * tau, the result being times corresponding to those velocities:

vs = vf * np.arange(v0[1]/vf, -1, -0.05)
# starts at v = v0[1] and excludes -vf (why?)

It’s helpful to plot this function to get an idea of its behavior:

import matplotlib.pyplot as plt
plt.plot(vs, integrand(vs))
plt.xlabel('v')
plt.ylabel('integrand')
plt.show()
plt.close()

Although the integrand is negative, integrating from a positive v0[1] to a negative v = -vf reverses the sign to produce a positive time.

The function quad() can’t handle numpy arrays as input, so we’ll use a list comprehension:

events = [ quad(integrand, v0[1], v) + (v,) for v in vs ]

The returned values have the form (t, error, v) (the last one appended as a tuple). To print them more nicely, we can use a formatting string:

for event in events: print('%5.2f %5.2e %6.2f' % tuple(event))

⇒
-0.00 0.00e+00  15.00
 0.12 1.32e-15  13.33
 0.24 2.70e-15  11.66
 0.37 4.12e-15  10.00
 0.50 5.60e-15   8.33
 0.64 7.14e-15   6.66
 0.79 8.75e-15   4.99
 0.94 1.04e-14   3.33
 1.10 1.22e-14   1.66
 1.26 1.40e-14  -0.00
 1.44 1.60e-14  -1.67
 1.62 1.80e-14  -3.34
 1.82 2.02e-14  -5.00
 2.02 2.24e-14  -6.67
 2.24 2.49e-14  -8.34
 2.48 2.75e-14 -10.01
 2.73 3.03e-14 -11.67
 3.00 3.33e-14 -13.34
 3.30 3.66e-14 -15.01
 3.62 1.03e-13 -16.68
 3.98 1.58e-12 -18.34
 4.38 2.45e-11 -20.01
 4.83 3.97e-10 -21.68
 5.36 6.92e-09 -23.35
 5.98 5.74e-12 -25.02
 6.73 3.65e-10 -26.68
 7.71 3.62e-08 -28.35
 9.09 1.14e-09 -30.02
11.45 1.92e-09 -31.69

We can see the error is very small so we can be confident in the result. Comparing to the result from the exact solution which is calculated for specific times:

# convert previous result to array and grab first column of times
ts = np.array(events)[:,0]
vdts = vd(ts,v0)
# vdts is a 3 x n array; we need n x 3 for the next step,
# so use the array method transpose():
for vdt in vdts.transpose(): print('%5.2f %5.2f %6.2f' % tuple(vdt))

 1.26 6.90   0.00
 1.44 6.55  -1.67
 1.62 6.21  -3.34
 1.82 5.86  -5.00
 2.02 5.52  -6.67
 2.24 5.17  -8.34
 2.48 4.83 -10.01
 2.73 4.48 -11.67
 3.00 4.14 -13.34
 3.30 3.79 -15.01
 3.62 3.45 -16.68
 3.98 3.10 -18.34
 4.38 2.76 -20.01
 4.83 2.41 -21.68
 5.36 2.07 -23.35
 5.98 1.72 -25.02
 6.73 1.38 -26.68
 7.71 1.03 -28.35
 9.09 0.69 -30.02
11.45 0.34 -31.69

The velocities in the third columns match! Remember, in the first case they were input to a calculation, in the second they were output from a calculation.

What’s the point other than confirming that integration by quadrature works? Well, now we can use the more general formula for drag which is important for large velocities:

m a = m g − ½ ρAC_Dvv

dv/dt = g − (1 + |v|/c)v/τ

∫_v0^v dv/(g − (1 + |v|/c)v/τ) = ∫₀^t dt = t

The constant c = b / C_D is very large; for b = 17 m/s and C_D = 0.5, c = 34 m/s. This is comparable to the terminal velocity, ~gτ ~ 33 m/s, so the extra term becomes significant on that order of speed.