Considering the application where the present mathematical model describes
the vibrations of a vehicle driving along a bumpy road, we need
to establish the force array F from the shape of the road \( h(x) \).
Various shapes are
available as a file with web address http://hplbit.bitbucket.org/data/bumpy/bumpy.dat.gz. The Python functionality for downloading this gzip
compressed file as a local file bumpy.dat.gz
and reading it into a numpy array goes as follows:
filename = 'bumpy.dat.gz'
url = 'http://hplbit.bitbucket.org/data/bumpy/bumpy.dat.gz'
import urllib
urllib.urlretrieve(url, filename)
h_data = np.loadtxt(filename) # read numpy array from file
The h_data object is a rectangular numpy array where the first
column contains the \( x \) coordinates along the road and the next columns
contain various road shapes \( h(x) \). We can extract the \( x \) data and
redefine h_data to contain solely the \( h(x) \) shapes:
x = h_data[0,:] # 1st column: x coordinates
h_data = h_data[1:,:] # other columns: h shapes
In general, the syntax a[s:t:i,2] gives a view (not a copy) to the part of
the array a where the first index goes from s to t,
but not including the t value, in increments of i, and the
second index is fixed at 2. Just writing : for an index means all
legal values of this index.
Given \( h(x) \), the corresponding acceleration \( a(t) \) needed in the force \( F(t)=-ma(t) \), follows from \( a(t)=h''(vt)v^2 \), where \( v \) is the velocity of the vehicle. The computation may utilize a finite difference approximation for the second-order derivative \( h'' \) and be encapsulated in a Python function:
def acceleration(h, x, v):
"""Compute 2nd-order derivative of h."""
# Method: standard finite difference aproximation
d2h = np.zeros(h.size)
dx = x[1] - x[0]
for i in range(1, h.size-1, 1):
d2h[i] = (h[i-1] - 2*h[i] + h[i+1])/dx**2
# Extraplolate end values from first interior value
d2h[0] = d2h[1]
d2h[-1] = d2h[-2]
a = d2h*v**2
return a
Note that here, h is a one-dimensional array containing the \( h(x) \)
values corresponding to a given coordinate array x. Also note that
we for mathematical simplicity set \( h''(x) \) at the end points equal to
\( h''(x) \) at the closest interior point.
The computations of d2h above was done array element by array element.
This loop can be a slow process in Python for long arrays. To speed
up computations dramatically, we can invoke a vectorization of
the above algorithm. This means that we get rid of the loops and
perform arithmetics on complete (or almost complete) arrays.
The vectorized form of the
acceleration function goes like
def acceleration_vectorized(h, x, v):
"""Compute 2nd-order derivative of h. Vectorized version."""
d2h = np.zeros(h.size)
dx = x[1] - x[0]
d2h[1:-1] = (h[:-2] - 2*h[1:-1] + h[2:])/dx**2
# Extraplolate end values from first interior value
d2h[0] = d2h[1]
d2h[-1] = d2h[-2]
a = d2h*v**2
return a
For each shape \( h(x) \) we want to compute the corresponding
vertical displacement \( u(t) \) using the mathematical model (1).
This can be accomplished by looping over the columns of h_data and
calling solver for each column, i.e., each realization of
the force \( F \).
The major arrays from the computations
are collected in a list data. The two first elements in
data are x and t. The next elements are
3-lists [h, a, u] for each road shape.
Note that some elements in data are arrays while others are list of
arrays. This composition is convenient when analyzing and visualizing
key quantities in the problem.
The computations of u for each road shape can be done as follows:
data = [x, t] # key input and output data (arrays)
for i in range(h_data.shape[0]):
h = h_data[i,:] # extract a column
a = acceleration(h, x, v)
F = -m*a
u = solver(t=t, I=0, m=m, b=b, f=f, F=F)
data.append([h, F, u])
A parameter choice \( m=60 \) kg, \( v=5 \) m/s, \( k=60 \) N/m, and \( b=80 \) Ns/m corresponds to a velocity of 18 km/h and a mass of 60 kg, i.e., bicycle conditions.