This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.

Summary

Chapter topics

This document has introduced computing with arrays and plotting curve data stored in arrays. The Numerical Python package contains lots of functions for array computing, including the ones listed in the table below. Plotting has been done with tools that closely resemble the syntax of MATLAB.

Construction	Meaning
`array(ld)`	copy list data `ld` to a `numpy` array
`asarray(d)`	make array of data `d` (no data copy if already array)
`zeros(n)`	make a `float` vector/array of length `n`, with zeros
`zeros(n, int)`	make an `int` vector/array of length `n` with zeros
`zeros((m,n))`	make a two-dimensional `float` array with shape (`m`,`n`)
`zeros_like(x)`	make array of same shape and element type as `x`
`linspace(a,b,m)`	uniform sequence of `m` numbers in $ [a,b] $
`a.shape`	tuple containing `a`'s shape
`a.size`	total no of elements in `a`
`len(a)`	length of a one-dim. array `a` (same as `a.shape[0]`)
`a.dtype`	the type of elements in `a`
`a.reshape(3,2)`	return `a` reshaped as $ 3\times 2 $ array
`a[i]`	vector indexing
`a[i,j]`	two-dim. array indexing
`a[1:k]`	slice: reference data with indices `1,2,...,k-1`
`a[1:10:3]`	slice: reference data with indices `1,4,7`
`b = a.copy()`	copy an array
`sin(a), exp(a), ...`	`numpy` functions applicable to arrays
`c = concatenate((a, b))`	`c` contains `a` with `b` appended
`c = where(cond, a1, a2)`	`c[i] = a1[i]` if `cond[i]`, else `c[i] = a2[i]`
`isinstance(a, ndarray)`	is `True` if `a` is an array

Array computing

When we apply a Python function f(x) to a Numerical Python array x, the result is the same as if we apply f to each element in x separately. However, when f contains if statements, these are in general invalid if an array x enters the boolean expression. We then have to rewrite the function, often by applying the where function from Numerical Python.

Plotting curves

The sections MATLAB-style plotting with Matplotlib and Matplotlib; pyplot prefix provide a quick overview of how to plot curves with the aid of Matplotlib. The same examples coded with the Easyviz plotting interface appear in the section SciTools and Easyviz.

Making movies

Each frame in a movie must be a hardcopy of a plot in PNG format. These plot files should have names containing a counter padded with leading zeros. One example may be tmp_0000.png, tmp_0001.png, tmp_0002.png. Having the plot files with names on this form, we can make an animated GIF movie in the file movie.gif, with two frames per second, by

os.system('convert -delay 50 tmp_*.png movie.gif')

Alternatively, we may combine the plot files to a Flash video:

os.system('ffmpeg -r 5 -i tmp_%04d.png -vcodec flv movie.flv')

Other formats can be made using other codecs, see the section Making videos.

Terminology

The important topics in this document are

array computing
vectorization
plotting
animations

Example: Animating a function

Problem

In this document's summarizing example we shall visualize how the temperature varies downward in the earth as the surface temperature oscillates between high day and low night values. One question may be: What is the temperature change 10 m down in the ground if the surface temperature varies between 2 C in the night and 15 C in the day?

Let the $ z $ axis point downwards, towards the center of the earth, and let $ z=0 $ correspond to the earth's surface. The temperature at some depth $ z $ in the ground at time $ t $ is denoted by $ T(z,t) $. If the surface temperature has a periodic variation around some mean value $ T_0 $, according to $$ \begin{equation*} T(0,t)=T_0+ A\cos (\omega t),\end{equation*} $$ one can find, from a mathematical model for heat conduction, that the temperature at an arbitrary depth is $$ \begin{equation} T(z,t) = T_0 + Ae^{-az}\cos (\omega t - az),\quad a =\sqrt{\omega\over 2k}\tp \tag{18} \end{equation} $$ The parameter $ k $ reflects the ground's ability to conduct heat ($ k $ is called the thermal diffusivity or the heat conduction coefficient).

The task is to make an animation of how the temperature profile in the ground, i.e., $ T $ as a function of $ z $, varies in time. Let $ \omega $ correspond to a time period of 24 hours. The mean temperature $ T_0 $ is taken as 10 C, and the maximum variation $ A $ is assumed to be 10 C. The heat conduction coefficient $ k $ may be set as $ 1\hbox{ mm}^2/\hbox{s} $ (which is $ 10^{-6} \hbox{ m}^2/\hbox{s} $ in proper SI units).

Solution

To animate $ T(z,t) $ in time, we need to make a loop over points in time, and in each pass in the loop we must save a plot of $ T $, as a function of $ z $, to file. The plot files can then be combined to a movie. The algorithm becomes

for $ t_i=i\Delta t $, $ i=0,1,2\ldots,n $:

plot the curve $ y(z)=T(z,t_i) $
store the plot in a file

combine all the plot files into a movie

It can be wise to make a general animate function where we just feed in some $ f(x,t) $ function and make all the plot files. If animate has arguments for setting the labels on the axis and the extent of the $ y $ axis, we can easily use animate also for a function $ T(z,t) $ (we just use $ z $ as the name for the $ x $ axis and $ T $ as the name for the $ y $ axis in the plot). Recall that it is important to fix the extent of the $ y $ axis in a plot when we make animations, otherwise most plotting programs will automatically fit the extent of the axis to the current data, and the tick marks on the $ y $ axis will jump up and down during the movie. The result is a wrong visual impression of the function.

The names of the plot files must have a common stem appended with some frame number, and the frame number should have a fixed number of digits, such as 0001, 0002, etc. (if not, the sequence of the plot files will not be correct when we specify the collection of files with an asterisk for the frame numbers, e.g., as in tmp*.png). We therefore include an argument to animate for setting the name stem of the plot files. By default, the stem is tmp_, resulting in the filenames tmp_0000.png, tmp_0001.png, tmp_0002.png, and so forth. Other convenient arguments for the animate function are the initial time in the plot, the time lag $ \Delta t $ between the plot frames, and the coordinates along the $ x $ axis. The animate function then takes the form

def animate(tmax, dt, x, function, ymin, ymax, t0=0,
            xlabel='x', ylabel='y', filename='tmp_'):
    t = t0
    counter = 0
    while t <= tmax:
        y = function(x, t)
        plot(x, y, '-',
             axis=[x[0], x[-1], ymin, ymax],
             title='time=%2d h' % (t/3600.0),
             xlabel=xlabel, ylabel=ylabel,
             savefig=filename + '%04d.png' % counter)
        savefig('tmp_%04d.pdf' % counter)
        t += dt
        counter += 1

The $ T(z,t) $ function is easy to implement, but we need to decide whether the parameters $ A $, $ \omega $, $ T_0 $, and $ k $ shall be arguments to the Python implementation of $ T(z,t) $ or if they shall be global variables. Since the animate function expects that the function to be plotted has only two arguments, we must implement $ T(z,t) $ as T(z,t) in Python and let the other parameters be global variables (or one can let T be a class with the parameters as attributes and a method __call__(self, x, t)). The T(z,t) implementation then reads

def T(z, t):
    # T0, A, k, and omega are global variables
    a = sqrt(omega/(2*k))
    return T0 + A*exp(-a*z)*cos(omega*t - a*z)

Suppose we plot $ T(z,t) $ at $ n $ points for $ z\in [0,D] $. We make such plots for $ t\in [0,t_{\max}] $ with a time lag $ \Delta t $ between the them. The frames in the movie are now made by

# set T0, A, k, omega, D, n, tmax, dt
z = linspace(0, D, n)
animate(tmax, dt, z, T, T0-A, T0+A, 0, 'z', 'T')

We have here set the extent of the $ y $ axis in the plot as $ [T_0-A,T_0+A] $, which is in accordance with the $ T(z,t) $ function.

The call to animate above creates a set of files with names of the form tmp_*.png. Out of these files we can create an animated GIF movie or a video in, e.g., Flash format by running operating systems commands with convert and avconv (or ffmpeg):

os.system('convert -delay 50 tmp_*.png movie.gif')
os.system('avconv -i tmp_%04d.png -r 5 -vcodec flv movie.flv')

See the section Making videos for how to create videos in other formats.

It now remains to assign proper values to all the global variables in the program: n, D, T0, A, omega, dt, tmax, and k. The oscillation period is 24 hours, and $ \omega $ is related to the period $ P $ of the cosine function by $ \omega = 2\pi /P $ (realize that $ \cos (t2\pi/P) $ has period $ P $). We then express $ P=24 $ h as $ 24\cdot 60\cdot 60 $ s and compute $ \omega $ as $ 2\pi/P \approx 7\cdot 10^{-5}\hbox{ s}^{-1} $. The total simulation time can be 3 periods, i.e., $ t_{\max} = 3P $. The $ T(z,t) $ function decreases exponentially with the depth $ z $ so there is no point in having the maximum depth $ D $ larger than the depth where $ T $ is visually zero, say 0.001. We have that $ e^{-aD}=0.001 $ when $ D=-a^{-1}\ln 0.001 $, so we can use this estimate in the program. The proper initialization of all parameters can then be expressed as follows:

k = 1E-6        # thermal diffusivity (in m*m/s)
P = 24*60*60.   # oscillation period of 24 h (in seconds)
omega = 2*pi/P
dt = P/24       # time lag: 1 h
tmax = 3*P      # 3 day/night simulation
T0 = 10         # mean surface temperature in Celsius
A = 10          # amplitude of the temperature variations in Celsius
a = sqrt(omega/(2*k))
D = -(1/a)*log(0.001) # max depth
n = 501         # no of points in the z direction

Note that it is very important to use consistent units. Here we express all units in terms of meter, second, and Kelvin or Celsius.

We encourage you to run the program heatwave.py to see the movie. The hardcopy of the movie is in the file movie.gif. Figure 21 displays two snapshots in time of the $ T(z,t) $ function.

Figure 21: Plot of the temperature $ T(z,t) $ in the ground for two different $ t $ values.

Scaling

In this example, as in many other scientific problems, it was easier to write the code than to assign proper physical values to the input parameters in the program. To learn about the physical process, here how heat propagates from the surface and down in the ground, it is often advantageous to scale the variables in the problem so that we work with dimensionless variables. Through the scaling procedure we normally end up with much fewer physical parameters that must be assigned values. Let us show how we can take advantage of scaling the present problem.

Consider a variable $ x $ in a problem with some dimension. The idea of scaling is to introduce a new variable $ \bar x = x/x_c $, where $ x_c $ is a characteristic size of $ x $. Since $ x $ and $ x_c $ have the same dimension, the dimension cancels in $ \bar x $ such that $ \bar x $ is dimensionless. Choosing $ x_c $ to be the expected maximum value of $ x $, ensures that $ \bar x\leq 1 $, which is usually considered a good idea. That is, we try to have all dimensionless variables varying between zero and one. For example, we can introduce a dimensionless $ z $ coordinate: $ \bar z= z/D $, and now $ \bar z \in [0,1] $. Doing a proper scaling of a problem is challenging so for now it is sufficient to just follow the steps below - and not worry why we choose a certain scaling.

In the present problem we introduce these dimensionless variables: $$ \begin{align*} \bar z &= z/D\\ \bar T &= {T-T_0\over A}\\ \bar t &= {\omega t} \end{align*} $$ We now insert $ z=\bar z D $ and $ t=\bar t/\omega $ in the expression for $ T(z,t) $ and get $$ \begin{equation*} T = T_0 + Ae^{-b\bar z}\cos(\bar t - b\bar z),\quad b=aD\end{equation*} $$ or $$ \begin{equation*} \bar T(\bar z,\bar t) = {T-T_0\over A} = e^{-b\bar z}\cos(\bar t - b\bar z)\tp\end{equation*} $$ We see that $ \bar T $ depends on only one dimensionless parameter $ b $ in addition to the independent dimensionless variables $ \bar z $ and $ \bar t $. It is common practice at this stage of the scaling to just drop the bars and write $$ \begin{equation} T(z,t) = e^{-bz}\cos(t - bz)\tp \tag{19} \end{equation} $$ This function is much simpler to plot than the one with lots of physical parameters, because now we know that $ T $ varies between $ -1 $ and 1, $ t $ varies between 0 and $ 2\pi $ for one period, and $ z $ varies between 0 and 1. The scaled temperature has only one parameter $ b $ in addition to the independent variable. That is, the shape of the graph is completely determined by $ b $.

In our previous movie example, we used specific values for $ D $, $ \omega $, and $ k $, which then implies a certain $ b=D\sqrt{\omega /(2k)} $ ($ \approx 6.9 $). However, we can now run different $ b $ values and see the effect on the heat propagation. Different $ b $ values will in our problems imply different periods of the surface temperature variation and/or different heat conduction values in the ground's composition of rocks. Note that doubling $ \omega $ and $ k $ leaves the same $ b $ - it is only the fraction $ \omega/k $ that influences the value of $ b $.

We can reuse the animate function also in the scaled case, but we need to make a new $ T(z,t) $ function and, e.g., a main program where $ b $ can be read from the command line:

def T(z, t):
    return exp(-b*z)*cos(t - b*z)  # b is global

b = float(sys.argv[1])
n = 401
z = linspace(0, 1, n)
animate(3*2*pi, 0.05*2*pi, z, T, -1.2, 1.2, 0, 'z', 'T')
movie('tmp_*.png', encoder='convert', fps=2,
      output_file='tmp_heatwave.gif')
os.system('convert -delay 50 tmp_*.png movie.gif')

Running the program, found as the file heatwave_scaled.py, for different $ b $ values shows that $ b $ governs how deep the temperature variations on the surface $ z=0 $ penetrate. A large $ b $ makes the temperature changes confined to a thin layer close to the surface, while a small $ b $ leads to temperature variations also deep down in the ground. You are encouraged to run the program with $ b=2 $ and $ b=20 $ to experience the major difference, or just view the ready-made animations.

We can understand the results from a physical perspective. Think of increasing $ \omega $, which means reducing the oscillation period so we get a more rapid temperature variation. To preserve the value of $ b $ we must increase $ k $ by the same factor. Since a large $ k $ means that heat quickly spreads down in the ground, and a small $ k $ implies the opposite, we see that more rapid variations at the surface requires a larger $ k $ to more quickly conduct the variations down in the ground. Similarly, slow temperature variations on the surface can penetrate deep in the ground even if the ground's ability to conduct ($ k $) is low.