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

Exercises

Exercise 1: Fill lists with function values

Define $$ \begin{equation} h(x) = {1\over\sqrt{2\pi}}e^{-\frac{1}{2}x^2}\tp \tag{20} \end{equation} $$ Fill lists xlist and hlist with \( x \) and \( h(x) \) values for 41 uniformly spaced \( x \) coordinates in \( [-4,4] \).

Hint.

Filename: fill_lists.

Exercise 2: Fill arrays; loop version

The aim is to fill two arrays x and y with \( x \) and \( h(x) \) values, respectively, where \( h(x) \) is defined in (20). Let the \( x \) values be as in Exercise 1: Fill lists with function values. Create empty x and y arrays and compute each element in x and y with a for loop. Filename: fill_arrays_loop.

Exercise 3: Fill arrays; vectorized version

Vectorize the code in Exercise 2: Fill arrays; loop version by creating the \( x \) values using the linspace function from the numpy package and by evaluating \( h(x) \) for an array argument. Filename: fill_arrays_vectorized.

Exercise 4: Plot a function

Make a plot of the function in Exercise 1: Fill lists with function values for \( x\in [-4,4] \). Filename: plot_Gaussian.

Exercise 5: Apply a function to a vector

Given a vector \( v=(2,3,-1) \) and a function \( f(x)=x^3 + xe^x + 1 \), apply \( f \) to each element in \( v \). Then calculate by hand \( f(v) \) as the NumPy expression v**3 + v*exp(v) + 1 using vector computing rules. Demonstrate that the two results are equal. Filename: apply_vecfunc.

Exercise 6: Simulate by hand a vectorized expression

Suppose x and t are two arrays of the same length, entering a vectorized expression

y = cos(sin(x)) + exp(1/t)

If x holds two elements, 0 and 2, and t holds the elements 1 and 1.5, calculate by hand (using a calculator) the y array. Thereafter, write a program that mimics the series of computations you did by hand (typically a sequence of operations of the kind we listed in the section Vector arithmetics and vector functions - use explicit loops, but at the end you can use Numerical Python functionality to check the results). Filename: simulate_vector_computing.

Exercise 7: Demonstrate array slicing

Create an array w with values \( 0,0.1,0.2,\ldots,3 \). Write out w[:], w[:-2], w[::5], w[2:-2:6]. Convince yourself in each case that you understand which elements of the array that are printed. Filename: slicing.

Exercise 8: Replace list operations by array computing

The data analysis problem in

the document Loops and lists [4] is solved by list operations. Convert the list to a two-dimensional array and perform the computations using array operations (i.e., no explicit loops, but you need a loop to make the printout). Filename: sun_data_vec.

Exercise 9: Plot a formula

Make a plot of the function \( y(t)=v_0t - \frac{1}{2}gt^2 \) for \( v_0=10 \), \( g=9.81 \), and \( t\in [0,2v_0/g] \). Set the axes labels as time (s) and height (m). Filename: plot_ball1.

Exercise 10: Plot a formula for several parameters

Make a program that reads a set of \( v_0 \) values from the command line and plots the corresponding curves \( y(t)=v_0t - \frac{1}{2}gt^2 \) in the same figure, with \( t\in [0,2v_0/g] \) for each curve. Set \( g=9.81 \).

Hint.

Filename: plot_ball2.

Exercise 11: Specify the extent of the axes in a plot

Extend the program from Exercise 10: Plot a formula for several parameters such that the minimum and maximum \( t \) and \( y \) values are computed, and use the extreme values to specify the extent of the axes. Add some space above the highest curve to make the plot look better. Filename: plot_ball3.

Exercise 12: Plot exact and inexact Fahrenheit-Celsius conversion formulas

A simple rule to quickly compute the Celsius temperature from the Fahrenheit degrees is to subtract 30 and then divide by 2: \( C = (F-30)/2 \). Compare this curve against the exact curve \( C = (F-32)5/9 \) in a plot. Let \( F \) vary between \( -20 \) and \( 120 \). Filename: f2c_shortcut_plot.

Exercise 13: Plot the trajectory of a ball

The formula for the trajectory of a ball is given by $$ \begin{equation} f(x) = x\tan\theta - {1\over 2v_0^2}{gx^2\over\cos^2\theta} + y_0, \tag{21} \end{equation} $$ where \( x \) is a coordinate along the ground, \( g \) is the acceleration of gravity, \( v_0 \) is the size of the initial velocity, which makes an angle \( \theta \) with the \( x \) axis, and \( (0,y_0) \) is the initial position of the ball.

In a program, first read the input data \( y_0 \), \( \theta \), and \( v_0 \) from the command line. Then plot the trajectory \( y=f(x) \) for \( y\geq 0 \). Filename: plot_trajectory.

Exercise 14: Plot data in a two-column file

The file src/plot/xy.dat contains two columns of numbers, corresponding to \( x \) and \( y \) coordinates on a curve. The start of the file looks as this:

     -1.0000       -0.0000
     -0.9933       -0.0087
     -0.9867       -0.0179
     -0.9800       -0.0274
     -0.9733       -0.0374

Make a program that reads the first column into a list x and the second column into a list y. Plot the curve. Print out the mean \( y \) value as well as the maximum and minimum \( y \) values.

Hint.

Filename: read_2columns.

Remarks

The function loadtxt in numpy can read files with tabular data (any number of columns) and return the data in a two-dimensional array:

import numpy as np
# Read table of floats
data = np.loadtxt('xy.dat', dtype=np.float)
# Extract one-dim arrays from two-dim data
x = data[:,0]  # column with index 0
y = data[:,1]  # column with index 1

The present exercise asks you to implement a simplified version of loadtxt, but for later loading of a file with tabular data into an array you will certainly use loadtxt.

Exercise 15: Write function data to file

We want to dump \( x \) and \( f(x) \) values to a file, where the \( x \) values appear in the first column and the \( f(x) \) values appear in the second. Choose \( n \) equally spaced \( x \) values in the interval \( [a,b] \). Provide \( f \), \( a \), \( b \), \( n \), and the filename as input data on the command line.

Hint.

Filename: write_cml_function.

Exercise 16: Plot data from a file

The files density_water.dat and density_air.dat files in the folder src/plot contain data about the density of water and air (respectively) for different temperatures. The data files have some comment lines starting with # and some lines are blank. The rest of the lines contain density data: the temperature in the first column and the corresponding density in the second column. The goal of this exercise is to read the data in such a file and plot the density versus the temperature as distinct (small) circles for each data point. Let the program take the name of the data file as command-line argument. Apply the program to both files. Filename: read_density_data.

Exercise 17: Write table to file

Given a function of two parameters \( x \) and \( y \), we want to create a file with a table of function values. The left column of the table contains \( y \) values in decreasing order as we go down the rows, and the last row contains the \( x \) values in increasing order. That is, the first column and the last row act like numbers on an \( x \) and \( y \) axis in a coordinate system. The rest of the table cells contains function values corresponding to the \( x \) and \( y \) values for the respective rows and columns. For example, if the function formula is \( x+2y \), \( x \) runs from 0 to 2 in steps of 0.5, and \( y \) run from -1 to 2 in steps of 1, the table looks as follows:

    2    4  4.5    5  5.5    6
    1    2  2.5    3  3.5    4
    0    0  0.5    1  1.5    2
   -1   -2 -1.5   -1 -0.5    0

         0  0.5    1  1.5    2

The task is to write a function

def write_table_to_file(f, xmin, xmax, nx, ymin, ymax, ny,
                        width=10, decimals=None,
                        filename='table.dat'):

where f is the formula, given as a Python function; xmin, xmax, ymin, and ymax are the minimum and maximum \( x \) and \( y \) values; nx is the number of intervals in the \( x \) coordinates (the number of steps in \( x \) direction is then (xmax-xmin)/nx); ny is the number of intervals in the \( y \) coordinates; width is the width of each column in the table (a positive integer); decimals is the number of decimals used when writing out the numbers (None means no decimal specification), and filename is the name of the output file. For example, width=10 and decimals=1 gives the output format %10.1g, while width=5 and decimals=None implies %5g.

Here is a test function which you should use to verify the implementation:

def test_write_table_to_file():
    filename = 'tmp.dat'
    write_table_to_file(f=lambda x, y: x + 2*y,
                        xmin=0, xmax=2, nx=4,
                        ymin=-1, ymax=2, ny=3,
                        width=5, decimals=None,
                        filename=filename)
    # Load text in file and compare with expected results
    with open(filename, 'r') as infile:
        computed = infile.read()
    expected = """\ 
    2    4  4.5    5  5.5    6
    1    2  2.5    3  3.5    4
    0    0  0.5    1  1.5    2
   -1   -2 -1.5   -1 -0.5    0

         0  0.5    1  1.5    2"""
    assert computed == expected

Filename: write_table_to_file.

Exercise 18: Fit a polynomial to data points

The purpose of this exercise is to find a simple mathematical formula for how the density of water or air depends on the temperature. The idea is to load density and temperature data from file as explained in Exercise 16: Plot data from a file and then apply some NumPy utilities that can find a polynomial that approximates the density as a function of the temperature.

NumPy has a function polyfit(x, y, deg) for finding a best fit of a polynomial of degree deg to a set of data points given by the array arguments x and y. The polyfit function returns a list of the coefficients in the fitted polynomial, where the first element is the coefficient for the term with the highest degree, and the last element corresponds to the constant term. For example, given points in x and y, polyfit(x, y, 1) returns the coefficients a, b in a polynomial a*x + b that fits the data in the best way. (More precisely, a line \( y=ax+b \) is a best fit to the data points \( (x_i,y_i) \), \( i=0,\ldots,n-1 \) if \( a \) and \( b \) are chosen to make the sum of squared errors \( R=\sum_{j=0}^{n-1} (y_j - (ax_j + b))^2 \) as small as possible. This approach is known as least squares approximation to data and proves to be extremely useful throughout science and technology.)

NumPy also has a utility poly1d, which can take the tuple or list of coefficients calculated by, e.g., polyfit and return the polynomial as a Python function that can be evaluated. The following code snippet demonstrates the use of polyfit and poly1d:

coeff = polyfit(x, y, deg)
p = poly1d(coeff)
print p             # prints the polynomial expression
y_fitted = p(x)     # computes the polynomial at the x points
# Use red circles for data points and a blue line for the polyn.
plot(x, y, 'ro', x, y_fitted, 'b-',
     legend=('data', 'fitted polynomial of degree %d' % deg))

a) Write a function fit(x, y, deg) that creates a plot of data in x and y arrays along with polynomial approximations of degrees collected in the list deg as explained above.

b) We want to call fit to make a plot of the density of water versus temperature and another plot of the density of air versus temperature. In both calls, use deg=[1,2] such that we can compare linear and quadratic approximations to the data.

c) From a visual inspection of the plots, can you suggest simple mathematical formulas that relate the density of air to temperature and the density of water to temperature?

Filename: fit_density_data.

Exercise 19: Fit a polynomial to experimental data

Suppose we have measured the oscillation period \( T \) of a simple pendulum with a mass \( m \) at the end of a massless rod of length \( L \). We have varied \( L \) and recorded the corresponding \( T \) value. The measurements are found in a file src/plot/pendulum.dat. The first column in the file contains \( L \) values and the second column has the corresponding \( T \) values.

a) Plot \( L \) versus \( T \) using circles for the data points.

b) We shall assume that \( L \) as a function of \( T \) is a polynomial. Use the NumPy utilities polyfit and poly1d, as explained in Exercise 18: Fit a polynomial to data points, to fit polynomials of degree 1, 2, and 3 to the \( L \) and \( T \) data. Visualize the polynomial curves together with the experimental data. Which polynomial fits the measured data best?

Filename: fit_pendulum_data.

Exercise 20: Read acceleration data and find velocities

A file src/plot/acc.dat contains measurements \( a_0,a_1,\ldots,a_{n-1} \) of the acceleration of an object moving along a straight line. The measurement \( a_k \) is taken at time point \( t_k=k\Delta t \), where \( \Delta t \) is the time spacing between the measurements. The purpose of the exercise is to load the acceleration data into a program and compute the velocity \( v(t) \) of the object at some time \( t \).

In general, the acceleration \( a(t) \) is related to the velocity \( v(t) \) through \( v'(t)=a(t) \). This means that $$ \begin{equation} v(t) = v(0) + \int_0^t a(\tau) d\tau\tp \tag{22} \end{equation} $$ If \( a(t) \) is only known at some discrete, equally spaced points in time, \( a_0,\ldots,a_{n-1} \) (which is the case in this exercise), we must compute the integral in (22) numerically, for example by the Trapezoidal rule: $$ \begin{equation} v(t_k) \approx \Delta t\left( \frac{1}{2}a_0 + \frac{1}{2}a_{k} + \sum_{i=1}^{k-1} a_i\right),\quad 1 \leq k \leq n-1\tp \tag{23} \end{equation} $$ We assume \( v(0)=0 \) so that also \( v_0=0 \).

Read the values \( a_0,\ldots,a_{n-1} \) from file into an array, plot the acceleration versus time, and use (23) to compute one \( v(t_k) \) value, where \( \Delta t \) and \( k\geq 1 \) are specified on the command line. Filename: acc2vel_v1.

Exercise 21: Read acceleration data and plot velocities

The task in this exercise is the same as in Exercise 20: Read acceleration data and find velocities, except that we now want to compute \( v(t_k) \) for all time points \( t_k=k\Delta t \) and plot the velocity versus time. Now only \( \Delta t \) is given on the command line, and the \( a_0,\ldots, a_{n-1} \) values must be read from file as in Exercise 20: Read acceleration data and find velocities.

Hint.

Filename: acc2vel.

Exercise 22: Plot a trip's path and velocity from GPS coordinates

A GPS device measures your position at every \( s \) seconds. Imagine that the positions corresponding to a specific trip are stored as \( (x,y) \) coordinates in a file src/plot/pos.dat with an \( x \) and \( y \) number on each line, except for the first line, which contains the value of \( s \).

a) Plot the two-dimensional curve of corresponding to the data in the file.

Hint.

b) Plot the velocity in \( x \) direction versus time in one plot and the velocity in \( y \) direction versus time in another plot.

Hint.

Filename: position2velocity.

Exercise 23: Vectorize the Midpoint rule for integration

The Midpoint rule for approximating an integral can be expressed as $$ \begin{equation} \int_a^b f(x)dx \approx h\sum_{i=1}^n f(a - \frac{1}{2}h + ih), \tag{25} \end{equation} $$ where \( h=(b-a)/n \).

a) Write a function midpointint(f, a, b, n) to compute Midpoint rule. Use a plain Python for loop to implement the sum.

b) Make a vectorized implementation of the Midpoint rule where you compute the sum by Python's built-in function sum.

c) Make another vectorized implementation of the Midpoint rule where you compute the sum by the sum function in the numpy package.

d) Organize the three implementations above in a module file midpoint_vec.py. Equip the module with one test function for verifying the three implementations. Use the integral \( \int_2^4 2xdx = 12 \) as test case since the Midpoint rule will integrate such a linear integrand exactly.

e) Start IPython, import the functions from midpoint_vec.py, define some Python implementation of a mathematical function \( f(x) \) to integrate, and use the %timeit feature of IPython to measure the efficiency of the three alternative implementations.

Filename: midpoint_vec.

Remarks

The lesson learned from the experiments in e) is that numpy.sum is much more efficient than Python's built-in function sum. Vectorized implementations must always make use of numpy.sum to compute sums.

Exercise 24: Vectorize a function for computing the area of a polygon

The area of a polygon is given by a formula in the exercise "Compute the area of a polygon" in the document Functions and branching. Vectorize this formula such that there are no Python loops in the implementation. Make a test function that compares the scalar implementation in the referred exercise with the new vectorized implementation for some chosen polygons (the scalar version must then be available in a module so that the function can be imported).

Hint.

Filename: polygon_area_vec.

Exercise 25: Implement Lagrange's interpolation formula

Imagine we have \( n+1 \) measurements of some quantity \( y \) that depends on \( x \): \( (x_0,y_0), (x_1,y_1),\ldots,(x_{n}, y_n) \). We may think of \( y \) as a function of \( x \) and ask what \( y \) is at some arbitrary point \( x \) not coinciding with any of the points \( x_0,\ldots,x_n \). It is not clear how \( y \) varies between the measurement points, but we can make assumptions or models for this behavior. Such a problem is known as interpolation.

One way to solve the interpolation problem is to fit a continuous function that goes through all the \( n+1 \) points and then evaluate this function for any desired \( x \). A candidate for such a function is the polynomial of degree \( n \) that goes through all the points. It turns out that this polynomial can be written $$ \begin{equation} p_L(x) = \sum_{k=0}^n y_kL_k(x), \tag{26} \end{equation} $$ where $$ \begin{equation} L_k(x) = \prod_{i=0,i\neq k}^n\frac{x-x_i}{x_k -x_i} \tag{27} \thinspace . \end{equation} $$ The \( \prod \) notation corresponds to \( \sum \), but the terms are multiplied. For example, $$ \begin{equation*} \prod_{i=0,i\neq k}^n x_i = x_0x_1\cdots x_{k-1}x_{k+1}\cdots x_n\thinspace . \end{equation*} $$ The polynomial \( p_L(x) \) is known as Lagrange's interpolation formula, and the points \( (x_0,y_0),\ldots,(x_n,y_n) \) are called interpolation points.

a) Make functions p_L(x, xp, yp) and L_k(x, k, xp, yp) that evaluate \( p_L(x) \) and \( L_k(x) \) by (26) and (27), respectively, at the point x. The arrays xp and yp contain the \( x \) and \( y \) coordinates of the \( n+1 \) interpolation points, respectively. That is, xp holds \( x_0,\ldots,x_n \), and yp holds \( y_0,\ldots,y_n \).

b) To verify the program, we observe that \( L_k(x_k)=1 \) and that \( L_k(x_i)=0 \) for \( i\neq k \), implying that \( p_L(x_k)=y_k \). That is, the polynomial \( p_L \) goes through all the points \( (x_0,y_0),\ldots,(x_n,y_n) \). Write a function test_p_L(xp, yp) that computes \( |p_L(x_k)-y_k| \) at all the interpolation points \( (x_k,y_k) \) and checks that the value is approximately zero. Call test_p_L with xp and yp corresponding to \( 5 \) equally spaced points along the curve \( y=\sin(x) \) for \( x\in [0,\pi] \). Thereafter, evaluate \( p_L(x) \) for an \( x \) in the middle of two interpolation points and compare the value of \( p_L(x) \) with the exact one.

Filename: Lagrange_poly1.

Exercise 26: Plot Lagrange's interpolating polynomial

a) Write a function

def graph(f, n, xmin, xmax, resolution=1001):

for plotting \( p_L(x) \) in Exercise 25: Implement Lagrange's interpolation formula, based on interpolation points taken from some mathematical function \( f(x) \) represented by the argument f. The argument n denotes the number of interpolation points sampled from the \( f(x) \) function, and resolution is the number of points between xmin and xmax used to plot \( p_L(x) \). The \( x \) coordinates of the n interpolation points can be uniformly distributed between xmin and xmax. In the graph, the interpolation points \( (x_0,y_0),\ldots,(x_n,y_n) \) should be marked by small circles. Test the graph function by choosing 5 points in \( [0,\pi] \) and f as \( \sin x \).

b) Make a module Lagrange_poly2 containing the p_L, L_k, test_p_L, and graph functions. The call to test_p_L described in Exercise 25: Implement Lagrange's interpolation formula and the call to graph described above should appear in the module's test block.

Filename: Lagrange_poly2.

Exercise 27: Investigate the behavior of Lagrange's interpolating polynomials

Unfortunately, the polynomial \( p_L(x) \) defined and implemented in Exercise 25: Implement Lagrange's interpolation formula can exhibit some undesired oscillatory behavior that we shall explore graphically in this exercise. Call the graph function from Exercise 26: Plot Lagrange's interpolating polynomial with \( f(x)=|x| \), \( x\in [-2,2] \), for \( n=2,4,6,10 \). All the graphs of \( p_L(x) \) should appear in the same plot for comparison. In addition, make a new figure with results from calls to graph for \( n=13 \) and \( n=20 \). All the code necessary for solving this exercise should appear in some separate program file, which imports the Lagrange_poly2 module made in Exercise 26: Plot Lagrange's interpolating polynomial. Filename: Lagrange_poly2b.

Remarks

The purpose of the \( p_L(x) \) function is to compute \( (x,y) \) between some given (often measured) data points \( (x_0,y_0),\ldots,(x_n,y_n) \). We see from the graphs that for a small number of interpolation points, \( p_L(x) \) is quite close to the curve \( y=|x| \) we used to generate the data points, but as \( n \) increases, \( p_L(x) \) starts to oscillate, especially toward the end points \( (x_0,y_0) \) and \( (x_n,y_n) \). Much research has historically been focused on methods that do not result in such strange oscillations when fitting a polynomial to a set of points.

Exercise 28: Plot a wave packet

The function $$ \begin{equation} f(x, t) = e^{-(x-3t)^2}\sin\left( 3\pi (x-t)\right) \tag{28} \end{equation} $$ describes for a fixed value of \( t \) a wave localized in space. Make a program that visualizes this function as a function of \( x \) on the interval \( [-4,4] \) when \( t=0 \). Filename: plot_wavepacket.

Exercise 29: Judge a plot

Assume you have the following program for plotting a parabola:

import numpy as np
x = np.linspace(0, 2, 20)
y = x*(2 - x)
import matplotlib.pyplot as plt
plt.plot(x, y)
plt.show()

Then you switch to the function \( \cos (18\pi x) \) by altering the computation of y to y = cos(18*pi*x). Judge the resulting plot. Is it correct? Display the \( \cos (18\pi x) \) function with 1000 points in the same plot. Filename: judge_plot.

Exercise 30: Plot the viscosity of water

The viscosity of water, \( \mu \), varies with the temperature \( T \) (in Kelvin) according to $$ \begin{equation} \mu (T) = A\cdot 10^{B/(T-C)}, \tag{29} \end{equation} $$ where \( A=2.414\cdot 10^{-5}\hbox{ Pa s} \), \( B=247.8 \) K, and \( C=140 \) K. Plot \( \mu (T) \) for \( T \) between 0 and 100 degrees Celsius. Label the \( x \) axis with 'temperature (C)' and the \( y \) axis with 'viscosity (Pa s)'. Note that \( T \) in the formula for \( \mu \) must be in Kelvin. Filename: water_viscosity.

Exercise 31: Explore a complicated function graphically

The wave speed \( c \) of water surface waves depends on the length \( \lambda \) of the waves. The following formula relates \( c \) to \( \lambda \): $$ \begin{equation} c(\lambda) = \sqrt{{g\lambda\over 2\pi}\left( 1 + s{4\pi^2\over\rho g \lambda^2}\right) \tanh\left( {2\pi h\over\lambda}\right)}\tp \tag{30} \end{equation} $$ Here, \( g \) is the acceleration of gravity (\( 9.81\hbox{ m/s}^2 \)), \( s \) is the air-water surface tension (\( 7.9\cdot 10^{-2} \hbox{ N/m} \)) , \( \rho \) is the density of water (can be taken as \( 1000 \hbox{ kg/m}^3 \)), and \( h \) is the water depth. Let us fix \( h \) at 50 m. First make a plot of \( c(\lambda) \) (in m/s) for small \( \lambda \) (0.001 m to 0.1 m). Then make a plot \( c(\lambda) \) for larger \( \lambda \) (1 m to 2 km. Filename: water_wave_velocity.

Exercise 32: Plot Taylor polynomial approximations to \( \sin x \)

The sine function can be approximated by a polynomial according to the following formula: $$ \begin{equation} \sin x \approx S(x;n) = \sum_{j=0}^{n} (-1)^j {x^{2j+1}\over (2j+1)!}\tp \tag{31} \end{equation} $$ The expression \( (2j+1)! \) is the factorial (math.factorial can compute this quantity). The error in the approximation \( S(x;n) \) decreases as \( n \) increases and in the limit we have that \( \lim_{n\rightarrow\infty} S(x;n) = \sin x \). The purpose of this exercise is to visualize the quality of various approximations \( S(x;n) \) as \( n \) increases.

a) Write a Python function S(x, n) that computes \( S(x;n) \). Use a straightforward approach where you compute each term as it stands in the formula, i.e., \( (-1)^j x^{2j+1} \) divided by the factorial \( (2j+1)! \).

b) Plot \( \sin x \) on \( [0,4\pi] \) together with the approximations \( S(x;1) \), \( S(x;2) \), \( S(x;3) \), \( S(x;6) \), and \( S(x;12) \).

Filename: plot_Taylor_sin.

Exercise 33: Animate a wave packet

Display an animation of the function \( f(x,t) \) in Exercise 28: Plot a wave packet by plotting \( f \) as a function of \( x \) on \( [-6,6] \) for a set of \( t \) values in \( [-1,1] \). Also make an animated GIF file.

Hint.

Filename: plot_wavepacket_movie.

Exercise 34: Animate a smoothed Heaviside function

Visualize the smoothed Heaviside function \( H_{\epsilon}(x) \), defined in the exercise "Implement a smoothed Heaviside function" in the document Functions and branching, as an animation where \( \epsilon \) starts at 2 and then goes to zero. Filename: smoothed_Heaviside_movie.

Exercise 35: Animate two-scale temperature variations

We consider temperature oscillations in the ground as addressed in the section Example: Animating a function. Now we want to visualize daily and annual variations. Let \( A_1 \) be the amplitude of annual variations and \( A_2 \) the amplitude of the day/night variations. Let also \( P_1 = 365 \) days and \( P_2=24 \) h be the periods of the annual and the daily oscillations. The temperature at time \( t \) and depth \( z \) is then given by $$ \begin{equation} T(z,t) = T_0 + A_1e^{-a_1z}\sin (\omega_1 t - a_1z) + A_2e^{-a_2z}\sin (\omega_2 t - a_2z), \tag{32} \end{equation} $$ where $$ \begin{align*} \omega_1 &= {2\pi P_1},\\ \omega_2 &= {2\pi P_2},\\ a_1 &=\sqrt{\omega_1\over 2k},\\ a_2 &=\sqrt{\omega_2\over 2k}\tp\\ \end{align*} $$ Choose \( k = 10^{-6} \hbox{ m}^2/\hbox{s} \), \( A_1=15 \) C, \( A_2=7 \) C, and the resolution \( \Delta t \) as \( P_2/10 \). Modify the heatwave.py program in order to animate this new temperature function. Filename: heatwave2.

Remarks

We assume in this problem that the temperature \( T \) equals the reference temperature \( T_0 \) at \( t=0 \), resulting in a sine variation rather than the cosine variation in (18).

Exercise 36: Use non-uniformly distributed coordinates for visualization

Watching the animation in Exercise 35: Animate two-scale temperature variations reveals that there are rapid oscillations in a small layer close to \( z=0 \). The variations away from \( z=0 \) are much smaller in time and space. It would therefore be wise to use more \( z \) coordinates close to \( z=0 \) than for larger \( z \) values. Given a set \( x_0 < x_1 < \cdots < x_n \) of uniformly spaced coordinates in \( [a, b] \), we can compute new coordinates \( \bar x_i \), stretched toward \( x=a \), by the formula $$ \begin{equation*} \bar x_i = a + (b-a)\left({x_i-a\over b-a}\right)^s,\end{equation*} $$ for some \( s>1 \). In the present example, we can use this formula to stretch the \( z \) coordinates to the left.

a) Experiment with \( s\in [1.2,3] \) and few points (say 15) and visualize the curve as a line with circles at the points so that you can easily see the distribution of points toward the left end. Identify a suitable value of \( s \).

b) Run the animation with no circles and (say) 501 points with the found \( s \) value.

Filename: heatwave2a.

Exercise 37: Animate a sequence of approximations to \( \pi \)

The exercise named "Approximate $\pi$" in the document Functions and branching [4] outlines an idea for approximating \( \pi \) as the length of a polygon inside the circle. Wrap the code from that exercise in a function pi_approx(N), which returns the approximation to \( \pi \) using a polygon with \( N+1 \) equally distributed points. The task of the present exercise is to visually display the polygons as a movie, where each frame shows the polygon with \( N+1 \) points together with the circle and a title reflecting the corresponding error in the approximate value of \( \pi \). The whole movie arises from letting \( N \) run through \( 4,5,6,\ldots,K \), where \( K \) is some (large) prescribed value. Let there be a pause of 0.3 s between each frame in the movie. By playing the movie you will see how the polygons move closer and closer to the circle and how the approximation to \( \pi \) improves. Filename: pi_polygon_movie.

Exercise 38: Animate a planet's orbit

A planet's orbit around a star has the shape of an ellipse. The purpose of this exercise is to make an animation of the movement along the orbit. One should see a small disk, representing the planet, moving along an elliptic curve. An evolving solid line shows the development of the planet's orbit as the planet moves and the title displays the planet's instantaneous velocity magnitude. As a test, run the special case of a circle and verify that the magnitude of the velocity remains constant as the planet moves.

Hint 1.

Hint 2.

Filename: planet_orbit.

Exercise 39: Animate the evolution of Taylor polynomials

A general series approximation (to a function) can be written as $$ \begin{equation*} S(x; M, N)=\sum_{k=M}^{N} f_k(x)\tp \end{equation*} $$ For example, the Taylor polynomial of degree \( N \) for \( e^x \) equals \( S(x; 0, N) \) with \( f_k(x)=x^k/k! \). The purpose of the exercise is to make a movie of how \( S(x;M,N) \) develops and improves as an approximation as we add terms in the sum. That is, the frames in the movie correspond to plots of \( S(x; M, M) \), \( S(x; M, M+1) \), \( S(x; M,M+2) \), \( \ldots \), \( S(x; M, N) \).

a) Make a function

animate_series(fk, M, N, xmin, xmax, ymin, ymax, n, exact)

for creating such animations. The argument fk holds a Python function implementing the term \( f_k(x) \) in the sum, M and N are the summation limits, the next arguments are the minimum and maximum \( x \) and \( y \) values in the plot, n is the number of \( x \) points in the curves to be plotted, and exact holds the function that \( S(x) \) aims at approximating.

Hint.

b) Call the animate_series function for the Taylor series for \( \sin x \), where \( f_k(x) = (-1)^k {x^{2k+1}/(2k+1)!} \), and \( x\in [0,13\pi] \), \( M=0 \), \( N=40 \), \( y\in [-2,2] \).

c) Call the animate_series function for the Taylor series for \( e^{-x} \), where \( f_k(x)=(-x)^k/k! \), and \( x\in [0, 15] \), \( M=0 \), \( N=30 \), \( y\in [-0.5,1.4] \).

Filename: animate_Taylor_series.

Exercise 40: Plot the velocity profile for pipeflow

A fluid that flows through a (very long) pipe has zero velocity on the pipe wall and a maximum velocity along the centerline of the pipe. The velocity \( v \) varies through the pipe cross section according to the following formula: $$ \begin{equation} v(r) = \left({\beta\over 2\mu_0}\right)^{{1/ n}} {n \over n+1}\left( R^{1 + 1/n} - r^{1 + 1/n}\right), \tag{33} \end{equation} $$ where \( R \) is the radius of the pipe, \( \beta \) is the pressure gradient (the force that drives the flow through the pipe), \( \mu_0 \) is a viscosity coefficient (small for air, larger for water and even larger for toothpaste), \( n \) is a real number reflecting the viscous properties of the fluid (\( n=1 \) for water and air, \( n < 1 \) for many modern plastic materials), and \( r \) is a radial coordinate that measures the distance from the centerline (\( r=0 \) is the centerline, \( r=R \) is the pipe wall).

a) Make a Python function that evaluates \( v(r) \).

b) Plot \( v(r) \) as a function of \( r\in [0,R] \), with \( R=1 \), \( \beta =0.02 \), \( \mu_0 =0.02 \), and \( n=0.1 \).

c) Make an animation of how the \( v(r) \) curves varies as \( n \) goes from 1 and down to 0.01. Because the maximum value of \( v(r) \) decreases rapidly as \( n \) decreases, each curve can be normalized by its \( v(0) \) value such that the maximum value is always unity.

Filename: plot_velocity_pipeflow.

Exercise 41: Plot sum-of-sines approximations to a function

The function $$ \begin{equation*} S(t;n) = {4\over\pi}\sum_{i=1}^n {1\over 2i-1} \sin\left( {2(2i-1)\pi t\over T}\right)\tp \end{equation*} $$ is as approximation to $$ \begin{equation*} f(t) = \left\lbrace\begin{array}{ll} 1, & 0 < t < T/2,\\ 0, & t = T/2,\\ -1, & T/2 < t < T \end{array}\right. \end{equation*} $$ It can be shown that \( S(t;n)\rightarrow f(t) \) as \( n\rightarrow\infty \).

Plot \( S(t;1) \), \( S(t;3) \), \( S(t;20) \), \( S(t;200) \), and the exact \( f(t) \) function in the same plot. Use \( T=2\pi \). Filename: sinesum1_plot.

Exercise 42: Animate the evolution of a sum-of-sine approximation to a function

First perform Exercise 41: Plot sum-of-sines approximations to a function. A natural next step is to animate the evolution of \( S(t;n) \) as \( n \) increases. Create such an animation and observe how the discontinuity in \( f(t) \) is poorly approximated by \( S(t;n) \), even when \( n \) grows large (plot \( f(t) \) in each frame). This is a well-known deficiency, called Gibb's phenomenon, when approximating discontinuous functions by sine or cosine (Fourier) series. Filename: sinesum1_movie.

Exercise 43: Plot functions from the command line

For quickly getting a plot of a function \( f(x) \) for \( x\in [x_{\rm min},x_{\rm max}] \) it could be nice to a have a program that takes the minimum amount of information from the command line and produces a plot on the screen and saves the plot to a file tmp.png. The usage of the program goes as follows:

plotf.py "f(x)" xmin xmax

Plotting \( e^{-0.2x}\sin (2\pi x) \) for \( x\in [0,4\pi] \) is then specified as

plotf.py "exp(-0.2*x)*sin(2*pi*x)" 0 4*pi

Write the plotf.py program with as short code as possible (we leave it to Exercise 44: Improve command-line input to test for valid input).

Hint.

Filename: plotf.

Exercise 44: Improve command-line input

Equip the program from Exercise 43: Plot functions from the command line with tests on valid input on the command line. Also allow an optional fourth command-line argument for the number of points along the function curve. Set this number to 501 if it is not given. Filename: plotf2.

Exercise 45: Demonstrate energy concepts from physics

The vertical position \( y(t) \) of a ball thrown upward is given by \( y(t) = v_0t - \frac{1}{2}gt^2 \), where \( g \) is the acceleration of gravity and \( v_0 \) is the velocity at \( t=0 \). Two important physical quantities in this context are the potential energy, obtained by doing work against gravity, and the kinetic energy, arising from motion. The potential energy is defined as \( P=mgy \), where \( m \) is the mass of the ball. The kinetic energy is defined as \( K=\frac{1}{2}mv^2 \), where \( v \) is the velocity of the ball, related to \( y \) by \( v(t)=y'(t) \).

Make a program that can plot \( P(t) \) and \( K(t) \) in the same plot, along with their sum \( P+K \). Let \( t\in [0,2v_0/g] \). Read \( m \) and \( v_0 \) from the command line. Run the program with various choices of \( m \) and \( v_0 \) and observe that \( P+K \) is always constant in this motion. (In fact, it turns out that \( P+K \) is constant for a large class of motions, and this is a very important result in physics.) Filename: energy_physics.

Exercise 46: Plot a w-like function

Define mathematically a function that looks like the "w" character. Plot the function. Also write a formal test function that verifies the implementation. Filename: plot_w.

Exercise 47: Plot a piecewise constant function

Consider the piecewise constant function defined in The exercise named "Implement a piecewise constant function" in the document Functions and branching [4]. Make a Python function plot_piecewise(data, xmax) that draws a graph of the function, where data is the nested list explained in mentioned exercise and xmax is the maximum \( x \) coordinate. Use ideas from the section Piecewisely defined functions. Filename: plot_piecewise_constant.

Exercise 48: Vectorize a piecewise constant function

Consider the piecewise constant function defined in The exercise named "Implement a piecewise constant function" in the document Functions and branching [4]. Make a vectorized implementation piecewise_constant_vec(x, data, xmax) of such a function, where x is an array.

Hint.

Filename: piecewise_constant_vec.

Remarks

Plotting the array returned from piecewise_constant_vec faces the same problems as encountered in the section Piecewisely defined functions. It is better to make a custom plotting function that simply draws straight horizontal lines in each interval.

Exercise 49: Visualize approximations in the Midpoint integration rule

Consider the midpoint rule for numerical integration. Use Matplotlib to make an illustration of the midpoint rule as shown to the left in Figure 22.

The \( f(x) \) function used in Figure 22 is $$ \begin{equation*} f(x) = x(12-x) + \sin(\pi x),\quad x\in[0,10]\tp\end{equation*} $$

Hint.

Filename: viz_midpoint.

Exercise 50: Visualize approximations in the Trapezoidal integration rule

Redo Exercise 49: Visualize approximations in the Midpoint integration rule for the Trapezoidal rule to produce the graph shown to the right in Figure 22. Filename: viz_trapezoidal.

Exercise 51: Experience overflow in a function

We are give the mathematical function $$ \begin{equation*} v(x)={1-e^{x/\mu}\over 1-e^{1/\mu}},\end{equation*} $$ where \( \mu \) is a parameter.

a) Make a Python function v(x, mu=1E-6, exp=math.exp) for calculating the formula for \( v(x) \) using exp as a possibly user-given exponential function. Let the v function return the nominator and denominator in the formula as well as the fraction.

b) Call the v function for various x values between 0 and 1 in a for loop, let mu be 1E-3, and have an inner for loop over two different exp functions: math.exp and numpy.exp. The output will demonstrate how the denominator is subject to overflow and how difficult it is to calculate this function on a computer.

c) Plot \( v(x) \) for \( \mu =1,0.01, 0.001 \) on \( [0,1] \) using 10,000 points to see what the function looks like.

d) Convert x and eps to a higher precision representation of real numbers, with the aid of the NumPy type float96, before calling v:

import numpy
x = numpy.float96(x); mu = numpy.float96(e)

Repeat point b) with these type of variables and observe how much better results we get with float96 compared with the standard float value, which is float64 (the number reflects the number of bits in the machine's representation of a real number).

e) Call the v function with x and mu as float32 variables and report how the function now behaves.

Filename: boundary_layer_func1.

Remarks

When an object (ball, car, airplane) moves through the air, there is a very, very thin layer of air close to the object's surface where the air velocity varies dramatically, from the same value as the velocity of the object at the object's surface to zero a few centimeters away. This layer is called a boundary layer. The physics in the boundary layer is important for air resistance and cooling/heating of objects. The change in velocity in the boundary layer is quite abrupt and can be modeled by the functiion \( v(x) \), where \( x=1 \) is the object's surface, and \( x=0 \) is some distance away where one cannot notice any wind velocity \( v \) because of the passing object (\( v= 0 \)). The wind velocity coincides with the velocity of the object at \( x=1 \), here set to \( v=1 \). The parameter \( \mu \) is very small and related to the viscosity of air. With a small value of \( \mu \), it becomes difficult to calculate \( v(x) \) on a computer. The exercise demonstrates the difficulties and provides a remedy.

Exercise 52: Apply a function to a rank 2 array

Let \( A \) be the two-dimensional array $$ \begin{equation*} \left\lbrack\begin{array}{ccc} 0 & 2 & -1 \\ -1 & -1 & 0 \\ 0 & 5 & 0 \end{array}\right\rbrack \end{equation*} $$ Apply the function \( f \) from Exercise 5: Apply a function to a vector to each element in \( A \). Then calculate the result of the array expression A**3 + A*exp(A) + 1, and demonstrate that the end result of the two methods are the same. Filename: apply_arrayfunc.

Exercise 53: Explain why array computations fail

The following loop computes the array y from x:

>>> import numpy as np
>>> x = np.linspace(0, 1, 3)
>>> y = np.zeros(len(x))
>>> for i in range(len(x)):
...     y[i] = x[i] + 4

However, the alternative loop

>>> for xi, yi in zip(x, y):
...     yi = xi + 5

leaves y unchanged. Why? Explain in detail what happens in each pass of this loop and write down the contents of xi, yi, x, and y as the loop progresses. Filename: find_errors_arraycomp.

Exercise 54: Verify linear algebra results

When we want to verify that a mathematical result is true, we often generate matrices or vectors with random elements and show that the result holds for these "arbitrary" mathematical objects. As an example, consider testing that \( A+B=B+A \) for matrices \( A \) and \( B \):

def test_addition():
    n = 4  # matrix size
    A = matrix(random.rand(n, n))
    B = matrix(random.rand(n, n))

    tol = 1E-14
    result1 = A + B
    result2 = B + A
    assert abs(result1 - result2).max() < tol

Use this technique to write test functions for the following mathematical results:

1. \( (A + B)C = AC + BC \)
2. \( (AB)C = A(BC) \)
3. \( \text{rank}A=\text{rank}A^T \)
4. \( \det(AB) = \det A\det B \)
5. The eigenvalues if \( A \) equals the eigenvalues of \( A^T \) when \( A \) is square.

References

1. H. P. Langtangen. Random numbers and simple games, \emphhttp://hplgit.github.io/primer.html/doc/pub/random, http://hplgit.github.io/primer.html/doc/pub/random.
2. H. P. Langtangen. Python Scripting for Computational Science, 3rd edition, Texts in Computational Science and Engineering, Springer, 2009.
3. H. P. Langtangen. User input and error handling, \emphhttp://hplgit.github.io/primer.html/doc/pub/input, http://hplgit.github.io/primer.html/doc/pub/input.
4. H. P. Langtangen. Functions and branching, \emphhttp://hplgit.github.io/primer.html/doc/pub/funcif, http://hplgit.github.io/primer.html/doc/pub/funcif.

• ← Prev