Figure 22: Visualization of numerical integration rules, with the Midpoint rule to the left and the Trapezoidal rule to the right. The filled areas illustrate the deviations in the approximation of the area under the curve.
This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.
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. You may adapt the example in the section Using lists for collecting function data.
Filename: fill_lists.
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.
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.
Make a plot of the function in Exercise 1: Fill lists with function values
for \( x\in [-4,4] \).
Filename: plot_Gaussian.
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.
Suppose x and t are two arrays of the same length, entering a
vectorized expression
y = cos(sin(x)) + exp(1/t)
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.
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.
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.
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.
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. You need a different vector of \( t \) coordinates for each curve.
Filename: plot_ball2.
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.
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.
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.
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
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.
Read the file line by line, split each line
into words, convert to float, and append to x and y.
The computations with y are simpler if the list is converted
to an array.
Filename: read_2columns.
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
loadtxt, but for later loading of a file with tabular data into
an array you will certainly use loadtxt.
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.
You may use the StringFunction tool (see the sections ref{sec:input:StringFunction} and
Vectorization of StringFunction objects) to turn the textual expression
for \( f \) into a Python function.
Filename: write_cml_function.
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.
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'):
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
write_table_to_file.
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.
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.
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.
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.
Repeated use of (23) for all \( k \) values
is very inefficient. A more efficient formula arises if we add the area of
a new trapezoid to the previous
$$
\begin{equation}
v(t_{k}) = v(t_{k-1}) + \int\limits_{t_{k-1}}^{t_k} a(\tau) d\tau
\approx v(t_{k-1}) + \Delta t \frac{1}{2}(a_{k-1} + a_{k}),
\tag{24}
\end{equation}
$$
for \( k=1,2,\ldots,n-1 \), while
\( v_0=0 \).
Use this formula to fill an array v with velocity values.
Filename: acc2vel.
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.
Load \( s \) into a float variable
and then the \( x \) and \( y \) numbers into two arrays. Draw a straight line
between the points, i.e., plot the \( y \) coordinates versus the \( x \)
coordinates.
b) Plot the velocity in \( x \) direction versus time in one plot and the velocity in \( y \) direction versus time in another plot.
Hint.
If \( x(t) \) and \( y(t) \) are the coordinates of the positions as a function
of time, we have that the velocity in \( x \) direction is \( v_x(t)=dx/dt \),
and the velocity in \( y \) direction is \( v_y=dy/dt \).
Since \( x \) and \( y \) are only known for some discrete times,
\( t_k=ks \), \( k=0,\ldots,n-1 \), we must use numerical differentiation. A
simple (forward) formula is
$$
\begin{equation*} v_x(t_k) \approx {x(t_{k+1})-x(t_k)\over s},\quad
v_y(t_k) \approx {y(t_{k+1})-y(t_k)\over s},\quad k=0,\ldots,n-2\tp \end{equation*}
$$
Compute arrays vx and vy with velocities based on the
formulas above for
\( v_x(t_k) \) and \( v_y(t_k) \), \( k=0,\ldots,n-2 \).
Filename: position2velocity.
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.
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.
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.
Observe that the formula \( x_1y_2+x_2y_3 + \cdots + x_{n-1}y_n =
\sum_{i=0}^{n-1}x_iy_{i+1} \) is the dot product of two vectors,
x[:-1] and y[1:],
which can be computed as numpy.dot(x[:-1], y[1:]).
Filename: polygon_area_vec.
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.
a) Write a function
def graph(f, n, xmin, xmax, resolution=1001):
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.
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.
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.
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.
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()
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.
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.
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.
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.
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. A suitable resolution can be 1000 intervals (1001 points) along the \( x \) axis, 60 intervals (61 points) in time, and 6 frames per second in the animated GIF file. Use the recipe in the section Making animations and remember to remove the family of old plot files in the beginning of the program.
Filename: plot_wavepacket_movie.
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.
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.
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).
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.
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.
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. The points \( (x,y) \) along the ellipse are given by the expressions $$ \begin{equation*} x=a\cos (\omega t),\quad y = b\sin (\omega t),\end{equation*} $$ where \( a \) is the semi-major axis of the ellipse, \( b \) is the semi-minor axis, \( \omega \) is an angular velocity of the planet around the star, and \( t \) denotes time. One complete orbit corresponds to \( t\in [0, 2\pi/\omega] \). Let us discretize time into time points \( t_k = k\Delta t \), where \( \Delta t = 2\pi/(\omega n) \). Each frame in the movie corresponds to \( (x,y) \) points along the curve with \( t \) values \( t_0,t_1,\ldots,t_i \), \( i \) representing the frame number (\( i=1,\ldots,n \)).
Hint 2. The velocity vector is $$ \begin{equation*} (\frac{dx}{dt}, \frac{dy}{dt}) = (-\omega a\sin(\omega t), \omega b\cos(\omega t)),\end{equation*} $$ and the magnitude of this vector becomes \( \omega\sqrt{a^2\sin^2(\omega t) + b^2\cos^2(\omega t)} \).
Filename: planet_orbit.
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)
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.
Here is some more information on how to write the animate_series
function. The function must accumulate the \( f_k(x) \) terms in a
variable \( s \), and for each \( k \) value, \( s \) is plotted against \( x \)
together with a curve reflecting the exact function. Each plot must
be saved in a file, say with names tmp_0000.png, tmp_0001.png, and
so on (these filenames can be generated by tmp_%04d.png, using an
appropriate counter). Use the movie function to combine all the plot
files into a movie in a desired movie format.
In the beginning of the animate_series function, it is necessary to
remove all old plot files of the form tmp_*.png. This can be done by
the glob module and the os.remove function as exemplified in
the section Making animations.
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.
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.
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.
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.
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
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.
Make \( x \) coordinates from the second and third command-line arguments and
then use eval
on the first argument.
Filename: plotf.
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.
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.
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.
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.
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.
You can use ideas from the Nv1 function in the section Vectorization of a hat function.
However, since the number of intervals is not known, it is necessary to
store the various intervals and conditions in lists.
Filename: piecewise_constant_vec.
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.
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.
Figure 22: Visualization of numerical integration rules, with the Midpoint rule to the left and the Trapezoidal rule to the right. The filled areas illustrate the deviations in the approximation of the area under the curve.
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.
Look up the documentation of the Matplotlib function fill_between
and use this function to create the filled areas between \( f(x) \)
and the approximating rectangles.
Note that the fill_between requires the two curves to have the same
number of points. For accurate visualization of \( f(x) \) you need quite
many \( x \) coordinates, and the rectangular approximation to \( f(x) \) must
be drawn using the same set of \( x \) coordinates.
Filename: viz_midpoint.
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.
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)
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.
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.
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.
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
>>> for xi, yi in zip(x, y):
... yi = xi + 5
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.
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
verify_linalg.