The code above for simulating vertical vibrations in a vehicle is naturally implemented as a Python function. This function can take the most important physical parameters of the problem as input, along with information about the file with road shapes. We allow for defining road shapes either through a file on a web site or a local file.
def bumpy_road(url=None, m=60, b=80, k=60, v=5):
"""
Simulate vertical vehicle vibrations.
========= ==============================================
variable description
========= ==============================================
url either URL of file with excitation force data,
or name of a local file
m mass of system
b friction parameter
k spring parameter
v (constant) velocity of vehicle
Return data (list) holding input and output data
[x, t, [h,F,u], [h,F,u], ...]
========= ==============================================
"""
# Download file (if url is not the name of a local file)
if url.startswith('http://') or url.startswith('file://'):
import urllib
filename = os.path.basename(url) # strip off path
urllib.urlretrieve(url, filename)
else:
# Check if url is the name of a local file
filename = url
if not os.path.isfile(filename):
print url, 'must be a URL or a filename'
sys.exit(1) # abort program
# else: ok
h_data = np.loadtxt(filename) # read numpy array from file
x = h_data[0,:] # 1st column: x coordinates
h_data = h_data[1:,:] # other columns: h shapes
t = x/v # time corresponding to x
dt = t[1] - t[0]
def f(u):
return k*u
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])
return data
Note that function arguments can be given default values (known as keyword arguments in Python). Python has a lot of operating system functionality, such as checking if a file, directory, or link exist, creating or removing files and directories, running stand-alone applications, etc.
After calling
road_url = 'http://hplbit.bitbucket.org/data/bumpy/bumpy.dat.gz'
data = solve(url=road_url, m=60, b=200, k=60, v=6)
the data array contains single arrays and triplets of arrays,
[x, t, [h,F,u], [h,F,u], ..., [h,F,u]]
This list, or any Python object, can be stored on file for later retrieval of the results, using the pickling functionality in Python:
import cPickle as pickle
# Or using a more advanced module
import dill as pickle
outfile = open('bumpy.res', 'w')
pickle.dump(data, outfile)
outfile.close()
The advantage of the dill module over cPickle is that it can
store a wider range of Python objects (e.g., lambda functions).
The code above and the bumpy_road function are found in the
file
bumpy.py.
Since the roads have a quite noise shape, the force \( F=-ma \) looks
very noisy, while the response \( u(t) \) to this excitation is
significantly less noisy,
see the bottom plot in Figure 5 for an example.
It may be useful to compute the root mean square value of \( u \) to get
a number for the typical amplitude of the vibrations:
$$ u_{\mbox{rms}} =
\sqrt{T^{-1}\int_0^T u^2dt} \approx \sqrt{\frac{1}{N+1}\sum_{i=0}^N (u^n)^2}
\thinspace .
$$
The last expression can be computed as follows using a vectorized
sum (np.sum(u**2)):
u_rms = []
for h, F, u in data[2:]:
u_rms.append(np.sqrt((1./len(u))*np.sum(u**2))
Very often in numerical computing we have some list/array and want to compute a new list/array where each element is some expression involving an element of the first list/array, typically
v = []
for element in u:
v.append(expression(element))
This code can more compactly be written as a list comprehension:
v = [expression(element) for element in u]
We may use the list comprehension construction for computing the u_rms
values:
u_rms = [np.sqrt((1./len(u))*np.sum(u**2))
for h, F, u in data[2:]]