A scientific application

Physical problem and mathematical model

The task is to make a simulation program that can predict how a (simple) mechanical system oscillates in response to environmental forces. Introducing \(u(t)\) as some displacement of the system at time \(t\), application of Newton’s second law of motion to such a mechanical system often results in the following type of equation for \(u\):

\[\tag{1} mu'' + f(u') + s(u) = F(t),\]

The prime, as in \(u'\), denotes differentiation with respect to time (\(u'(t)\) or \(du/dt\)). Furthermore, \(m\) is the mass of the system, \(f(u')\) is a friction force that gives rise to a damping of the motion, \(s(u)\) represents a restoring force, such as a spring, and \(F(t)\) models the external environmental forces on the system. Equation (1) must be accompanied by two initial conditions: \(u(0)=I\) and \(u'(0)=V\). The values of these have no effect on the steady state behavior of \(u(t)\) for large values of \(t\), since this behavior is determined by the force \(F(t)\) and the system parameters \(m\), \(f(u')\), and \(s(u)\).

There are two types of the friction force \(f\): linear damping \(f(u')=bu'\) and quadratic damping \(f(u')=bu'|u'|\). The input data consists of \(m\), \(b\), \(s(u)\), \(F(t)\), \(I\), \(V\), and specification of linear or quadratic damping. The unknown quantity to be computed is \(u(t)\) for \(t\in (0,T]\).

One example where the model above has relevance, is the vertical vibration of a vehicle in response to a bumpy road. Let \(h(x)\) be the height of the road at some coordinate \(x\) along the road. When driving along this road with constant velocity \(v\), the vehicle is moved up and down in time according to \(h(vt)\), resulting in an external vertical force \(F(t)=-mh''(vt)v^2\). We assume that the vehicle has springs and dampers that here are modeled as \(bu'\) (damper) and \(s(s)=ku\) (spring), with given damping parameter \(b\) and spring constant \(k\). The unknown \(u(t)\) is the vertical displacement of the vehicle relative to the road. Figure Vehicle on a bumpy road illustrates the situation. You may view an animation of such a motion on the web (by the way, this animation was created by coupling the Pysketcher tool for figure drawings with the differential equation solver presented later, see the script bumpy_road_fig.py for all details).

Vehicle on a bumpy road

Another example regards the vertical shaking of a building due to earthquake-induced movement of the ground. If the vertical displacement of the ground is recorded as a function \(d(t)\), this results in a vertical force \(F(t)=-md''(t)\). The soil foundation acts as a spring and damper on the building, modeled through the damping parameter \(b\) and normally a linear spring term \(s(u)=ku\).

In both cases we drop the effect of gravity, which is just a constant compression of the spring.

Our task is to compute and analyze the vertical \(u(t)\) vibrations of a vehicle, given the shape \(h(x)\) of the road and some velocity \(v\).

Numerical model

The differential equation problem (1) can be solved by introducing finite difference approximations for \(u''\) and \(u'\). In case of quadratic damping one can use a geometric mean to approximate \(u'|u'|\) and thereby linearize the equations. The result of using such numerical methods is an algorithm for computing \(u(t)\) at discrete points in time. Let \(u^n\) be the approximation to \(u\) at time \(t_n=n\Delta t\), \(n=1,2,\ldots\), where \(\Delta t\) is a (small) time interval. For example, if \(\Delta t = 0.1\), we find approximations \(u^1\) to \(u\) at \(t=0.1\), \(u^2\) at \(t=0.2\), \(u^3\) to \(t=0.3\), and so forth. Any value \(u^{n+1}\) can be computed if \(u^n\) and \(u^{n-1}\) are known (i.e., previously computed). The formula for \(u^{n+1}\) is, in case of linear damping \(f(u')=bu'\),

\[\tag{2} u^{n+1} = \left(2mu^n + (\frac{b}{2}\Delta t - m)u^{n-1} + \Delta t^2(F^n - s(u^n)) \right)(m + \frac{b}{2}\Delta t)^{-1},\]

where \(F^n\) means \(F(t)\) evaluated for \(t=t_n\). A special formula must be applied for \(n=0\):

\[\tag{3} u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) \thinspace .\]

For quadratic damping we have a slightly different formula,

\[u^{n+1} = \left( m + b|u^n-u^{n-1}|\right)^{-1}\times\nonumber\]

\[\tag{4} \quad \left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \Delta t^2 (F^n - s(u^n)) \right),\]

and again a special formula for \(u^1\):

\[\tag{5} u^1 = u^0 + \Delta t V + \frac{\Delta t^2}{2m}\left(-bV|V| - s(u^0) + F^0\right) \thinspace .\]

The implementation of the computational algorithm can make use of an array u to represent \(u^n\) as u[n]. The force \(F(t_n)\) is assumed to be available as an array element F[n]. The following Python function computes u given an array t with time points \(t_0,t_1,\ldots\), the initial displacement I, mass m, damping parameter b, restoring force s(u), environmental forces F as an array (corresponding to t).

Simple implementation

Let us first implement the computational formulas for the linear damping case in a short and compact Python function:

from numpy import *

def solver_linear_damping(I, V, m, b, s, F, t):
    N = t.size - 1              # No of time intervals
    dt = t[1] - t[0]            # Time step
    u = zeros(N+1)              # Result array
    u[0] = I
    u[1] = u[0] + dt*V + dt**2/(2*m)*(-b*V - s(u[0]) + F[0])

    for n in range(1,N):
        u[n+1] = 1./(m + b*dt/2)*(2*m*u[n] + \
                 (b*dt/2 - m)*u[n-1] + dt**2*(F[n] - s(u[n])))
    return u

Dissection of the code

Functions in Python start with def, followed by the function name and the list of input objects separated by comma. The function body is indented, and the first non-indented line signifies the end of the function body block. Output objects are returned to the calling code by a return statement.

The arguments to this function and the variables created inside the function are not declared with type. We therefore need to know what the variables are supposed to be: I, V, m, and b are real numbers, while F and t are one-dimensional arrays of the same length, where F holds \(F(t_n)\) and t holds \(t_n\), \(n=0,1,\ldots,N+1\).

The number of elements in an array t is given by t.size (or len(t), but t.size works for multi-dimensional arrays too). Arrays are indexed by square brackets, and indices always start at 0. Numerical codes frequently needs to loop over array indices, i.e., a set of integers. Such a set is produced by range(start, stop, increment), which returns a list of integers start, start+increment, start+2*increment, and so on, up to but not including stop. Writing just range(stop) means range(0, stop, 1). The particular call range(1, N) used in the code above results in a list of integers: 1, 2, ..., N-1.

Array functionality is enabled by the numpy package, which offers functions such as zeros and linspace, as known from Matlab. Here we import all objects in the numpy package by the statement

from numpy import *

Comments start with the character # and the rest of the line is then ignored by Python.

Every variable in Python is an object. In particular, the s function above is a function object, transferred to the function as any other object, and called as any other function. Transferring a function as argument to another function is therefore simpler and cleaner in Python than in, e.g., C, C++, Java, C#, and Matlab.

How can we use the solver_linear_damping function? We need to call it with relevant values for the arguments. Suppose we want to solve a vibration problem with \(I=1\), \(V=0\), \(F=0\), \(m=2\), \(b=0.2\), \(s(u)=2u\), \(\Delta t=0.2\), for \(t\in [0,10\pi]\). This will be a damped sinusoidal solution (setting \(b=0\) will result in \(u(t)=\cos t\)). The test code becomes

from solver import solver_linear_damping
from numpy import *

def s(u):
    return 2*u

T = 10*pi      # simulate for t in [0,T]
dt = 0.2
N = int(round(T/dt))
t = linspace(0, T, N+1)
F = zeros(t.size)
I = 1; V = 0
m = 2; b = 0.2
u = solver_linear_damping(I, V, m, b, s, F, t)

from matplotlib.pyplot import *
plot(t, u)
savefig('tmp.pdf')   # save plot to PDF file
savefig('tmp.png')   # save plot to PNG file
show()

The solver_linear_damping function resides in the file solver.py, so if we make the call to this function in separate file (assumed above), we have to import the function as shown. We also need functions and variables from the numpy package, here pi, linspace, and zeros. Normally, there is one statement per line in Python programs, and there is no need to end the statement with a semicolon. However, if we want multiple statements on a line, they must be separated by semi colon as demonstrated in the initialization of I and V. Plotting the computed curve makes use of the matplotlib package, where we call its plot, savefig, and show functions. Figure Plot of computed curve shows the result.

Plot of computed curve

More advanced implementation

Let us extend the function above to also include the quadratic damping formulas. The \(F(t)\) function in (1) is required to be available as an array in the solver_linear_damping function, but now we will allow for more flexibility: the F argument may either be a Python function F(t) or a Python array. In addition, we add some checks that variables have correct values or are of correct type.

import numpy as np

def solver(I, V, m, b, s, F, t, damping='linear'):
    """
    Solve m*u'' + f(u') + s(u) = F for time points in t.
    u(0)=I and u'(0)=V,
    by a central finite difference method with time step dt.
    If damping is 'linear', f(u')=b*u, while if damping is
    'quadratic', we have f(u')=b*u'*abs(u').
    s(u) is a Python function, while F may be a function
    or an array (then F[i] corresponds to F at t[i]).
    """
    N = t.size - 1              # No of time intervals
    dt = t[1] - t[0]            # Time step
    u = np.zeros(N+1)           # Result array
    b = float(b); m = float(m)  # Avoid integer division

    # Convert F to array
    if callable(F):
        F = F(t)
    elif isinstance(F, (list,tuple,np.ndarray)):
        F = np.asarray(F)
    else:
        raise TypeError(
            'F must be function or array, not %s' % type(F))

    u[0] = I
    if damping == 'linear':
        u[1] = u[0] + dt*V + dt**2/(2*m)*(-b*V - s(u[0]) + F[0])
    elif damping == 'quadratic':
        u[1] = u[0] + dt*V + \
               dt**2/(2*m)*(-b*V*abs(V) - s(u[0]) + F[0])
    else:
        raise ValueError('Wrong value: damping="%s"' % damping)

    for n in range(1,N):
        if damping == 'linear':
            u[n+1] = (2*m*u[n] + (b*dt/2 - m)*u[n-1] +
                      dt**2*(F[n] - s(u[n])))/(m + b*dt/2)
        elif damping == 'quadratic':
            u[n+1] = (2*m*u[n] - m*u[n-1] + b*u[n]*abs(u[n] - u[n-1])
                      - dt**2*(s(u[n]) - F[n]))/\
                      (m + b*abs(u[n] - u[n-1]))
    return u, t

Dissection of the code

Two types of import. This time we replace from numpy import *, which imports over 500 variables and functions, by import numpy as np, which just imports one variable, the package numpy, here under the nickname np. All variables and functions in numpy must now be reached via the prefix np., as in np.zeros, and np.linspace, and np.pi (for \(\pi\)). The advantage of the prefix is that we clearly see where functionality comes from. The disadvantage is that mathematical formulas like \(\sin (\pi x)\) must be written np.sin(np.pi*x). We can always perform an explicit import of some names, like sin and pi, to write the formula as sin(pi*x):

import numpy as np
from numpy import sin, pi

def myfunction(x):
    return sin(pi*x)

Another disadvantage with the np. prefix is that names are no longer (almost) the same as in Matlab. Many will therefore prefer to do from numpy import * and skip the prefix.

Doc strings. The string, enclosed in triple double-quotes, right after the function definition, is a doc string used for documenting the function. Various tools can extract function definitions and doc strings to automatically produce nicely typeset manuals.

Avoiding integer division. At one line we explicitly convert b and m to float variables. This is not strictly necessary, but if we supply b=2 and m=4, a computation like b/m will give zero as result because both b and m are then integers and b/m implies integer division, not the mathematical division of real numbers (this is not a special feature of Python - all languages with a strong heritage from C invoke integer division if both operands in a division are integers). We need to make sure that at least one of the operands in a division is a real number (float) to ensure the intended mathematical operation. Looking at the formulas, there is never a problem with integer division in our implementation, because dt is computed from t, which has float elements (linspace makes float elements by default), and all divisions in our code involve dt as one of the operands. However, future edits may alter the way formulas are written, so to be on the safe side we ensure that real input parameters are float objects.

Flexible variable type. As mentioned, we allow F to be either a function or an array. In the former case, we convert F to an array such that the rest of the code can assume that F is indeed an array. It is easy to check the type of variables in Python. The test if callable(F) is true if the object F can be called as a function.

Checking correct variable type. To test that F is an array, we can use isinstance(F, np.ndarray) (ndarray is the name of the array type in numpy). In the code above we allow that F can also be a list or a tuple. Running F through the asarray function makes an array out of F in the cases where F is a list or tuple. The final else: clause takes care of the situation where F is neither a function, nor an object that can easily be converted to an array, and an error message is issued. More precisely, we raise a TypeError exception indicating that we have encountered a wrong type. The error message will contain the type of F as obtained from type(F). Instead of using the syntax isinstance(F, list) we may test type(F) == list or even type(F) in (list,tuple,np.ndarray).

We could simplify the if-else test involving F to just the two lines

if callable(F):
    F = F(t)

if we are sure that F is either a function or a numpy array. Should the user send in something else for F, Python will encounter a run-time error when trying to index F as in F[0] (in the statement computing u[1]).

To be absolutely safe, we should test that the other arguments are of right type as well. For example,

if not isinstance(I, (float,int)):
    raise TypeError('V must be float or int, not %s' % type(V))

and similar for V, m, b, while s must be tested by callable(s), t must be np.ndarray, and damping must be str.

Calling the function. Below is a simple example on how the solver function can be called to solve the differential equation problem \(mu'' + bu' + ku = A\sin\pi t\), \(u(0)=I\), \(u'(0)=0\):

import numpy as np
from numpy import sin, pi  # for nice math

def F(t):
    # Sinusoidal bumpy road
    return A*sin(pi*t)

def s(u):
    return k*u

A = 0.25
k = 2
t = np.linspace(0, 20, 2001)
u, t = solver(I=0, V=0, m=2, b=0.05, s=s, F=F, t=t)

# Show u(t) as a curve plot
import matplotlib.pyplot as plt
plt.plot(t, u)
plt.show()

Local and global variables. We use in this example a linear spring function \(s(u)\),

def f(u):
    return k*u

Here, u is a local variable, which is accessible just inside in the function, while k is a global variable, which must be initialized outside the function prior to calling f.

Advanced programming of functions with parameters

The best way to implement mathematical functions that has a set of parameters in addition some independent variables is to create a class where the parameters are attributes and a __call__ method evaluates the function formula given the independent variables as arguments. This requires, of course, knowledge of classes and special methods like __call__.

As an example, the f(u) function above can be implemented as

class Spring:
    def __init__(self, k):
        self.k = k
    def __call__(self, u):
        return self.k*u

f = Spring(k)

Because of the __call__ method, we can make calls f(u). Note that the k parameter is bundled with the function in the f object.

The excitation force

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.

A high-level solve function

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.

Storing Python objects in files

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.

Computing the root mean square value

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 First realization of a bumpy road, with corresponding excitation of the wheel and resulting vertical vibrations 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:]]