Exercises

Problem 1: Make a tool for differentiating curves

Suppose we have a curve specified through a set of discrete coordinates \((x_i,y_i)\), \(i=0,\ldots,n\), where the \(x_i\) values are uniformly distributed with spacing \(\Delta x\): \(x_i=\Delta x\). The derivative of this curve, defined as a new curve with points \((x_i, d_i)\), can be computed via finite differences:

\[d_0 = \frac{y_1-y_0}{\Delta x},\]

\[d_i = \frac{y_{i+1}-y_{i-1}}{2\Delta x},\quad i=1,\ldots,n-1,\]

\[d_n = \frac{y_n-y_{n-1}}{\Delta x}{\thinspace .}\]

a) Write a function differentiate(x, y) for differentiating a curve with coordinates in the arrays x and y, using the formulas above. The function should return the coordinate arrays of the resulting differentiated curve.

b) Since the formulas for differentiation used here are only approximate, with unknown approximation errors, it is challenging to construct test cases. Here are three approaches, which should be implemented in three separate test functions.

1. Consider a curve with three points and compute \(d_i\), \(i=0,1,2\), by hand. Make a test that compares the hand-calculated results with those from the function in a).
2. The formulas for \(d_i\) are exact for points on a straight line, as all the \(d_i\) values are then the same, equal to the slope of the line. A test can check this property.
3. For point lying on a parabola, the values for \(d_i\), \(i=1,\ldots,n-1\), should equal the exact derivative of the parabola. Make a test based on this property.

c) Start with a curve corresponding to \(y=\sin(\pi x)\) and \(n+1\) points in \([0,1]\). Apply differentiate four times and plot the resulting curve and the exact \(y=\sin\pi x\) for \(n=6, 11, 21, 41\).

Filename: curvediff.

Problem 2: Make solid software for the Trapezoidal rule

An integral

\[\int_a^b f(x)dx\]

can be numerically approximated by the Trapezoidal rule,

\[\int_a^b f(x)dx \approx \frac{h}{2}(f(a) + f(b)) + h\sum_{i=1}^{n-1} f(x_i),\]

where \(x_i\) is a set of uniformly spaced points in \([a,b]\):

\[h = \frac{b-a}{n},\quad x_i=a + ih,\ i=1,\ldots,n-1{\thinspace .}\]

Somebody has used this rule to compute the integral \(\int_0^\pi \sin^2x\, dx\):

from math import pi, sin
np = 20
h = pi/np
I = 0
for k in range(1, np):
    I += sin(k*h)**2
print I

a) The “flat” implementation above suffers from serious flaws:

1. A general numerical algorithm (the Trapezoidal rule) is implemented in a specialized form where the formula for \(f\) is inserted directly into the code for the general integration formula.
2. A general numerical algorithm is not encapsulated as a general function, with appropriate parameters, which can be reused across a wide range of applications.
3. The lazy programmer dropped the first terms in the general formula since \(\sin(0)=\sin(\pi)=0\).
4. The sloppy programmer used np (number of points?) as variable for n in the formula and a counter k instead of i. Such small deviations from the mathematical notation are completely unnecessary. The closer the code and the mathematics can get, the easier it is to spot errors in formulas.

Write a function trapezoidal that fixes these flaws. Place the function in a module trapezoidal.

b) Write a test function test_trapezoidal. Call the test function explicitly to check that it works. Remove the call and run pytest on the module:

Terminal> py.test -s -v trapezoidal

Hint. Note that even if you know the value of the integral, you do not know the error in the approximation produced by the Trapezoidal rule. However, the Trapezoidal rule will integrate linear functions exactly (i.e., to machine precision). Base a test function on a linear \(f(x)\).

c) Add functionality such that we can compute \(\int_a^b f(x)dx\) by providing \(f\), \(a\), \(b\), and \(n\) as positional command-line arguments to the module file:

Terminal> python trapezoidal.py 'sin(x)**2' 0 pi 20

Here, \(a=0\), \(b=\pi\), and \(n=20\).

Note that the trapezoidal.py file must still be a valid module file, so the interpretation of command-line data and computation of the integral must be performed from calls in a test block.

Hint. To translate a string formula on the command line, like sin(x)**2, into a Python function, you can wrap a function declaration around the formula and run exec on the string to turn it into live Python code:

import math, sys
formula = sys.argv[1]
f_code = """
def f(x):
    return %s
""" % formula
exec(code, math.__dict__)

The result is the same as if we had hardcoded

from math import *

def f(x):
    return sin(x)**2

in the program. Note that exec needs the namespace math.__dict__, i.e., all names in the math module, such that it understands sin and other mathematical functions. Similarly, to allow \(a\) and \(b\) to be math expressions like pi/4 and exp(4), do

a = eval(sys.argv[2], math.__dict__)
b = eval(sys.argv[2], math.__dict__)

d) Write a test function for verifying the implementation of data reading from the command line.

Filename: trapezoidal.

Problem 3: Implement classes for the Trapezoidal rule

We consider the same problem setting as in Problem 2: Make solid software for the Trapezoidal rule. Make a module with a class Problem representing the mathematical problem to be solved and a class Solver representing the solution method. The rest of the functionality of the module, including test functions and reading data from the command line, should be as in Problem 2: Make solid software for the Trapezoidal rule. Filename: trapezoidal_class.

Problem 4: Write a doctest and a test function

Type in the following program:

import sys
# This sqrt(x) returns real if x>0 and complex if x<0
from numpy.lib.scimath import sqrt

def roots(a, b, c):
    """
    Return the roots of the quadratic polynomial
    p(x) = a*x**2 + b*x + c.

    The roots are real or complex objects.
    """
    q = b**2 - 4*a*c
    r1 = (-b + sqrt(q))/(2*a)
    r2 = (-b - sqrt(q))/(2*a)
    return r1, r2

a, b, c = [float(arg) for arg in sys.argv[1:]]
print roots(a, b, c)

a) Equip the roots function with a doctest. Make sure to test both real and complex roots. Write out numbers in the doctest with 14 digits or less.

b) Make a test function for the roots function. Perform the same mathematical tests as in a), but with different programming technology.

Filename: test_roots.

Problem 5: Experiment with tolerances in comparisons

When we replace a comparison a == b, where a and/or b are float objects, by a comparison with tolerance, abs(a-b) < tol, the appropriate size of tol depends on the size a and b. Investigate how the size of abs(a-b) varies when b takes on values \(10^k\), \(k=-5,-9,\ldots,20\) and a=1.0/49*b*49. Filename: tolerance.

Remarks

You will experience that if a and b are large, as they can be in geophysical applications where lengths measured in meters can be of size \(10^6\) m, tol must be about \(10^{-9}\), while a and b around unity can have tol of size \(10^{-15}\).

Exercise 6: Make use of a class implementation

Implement the experiment_compare_dt function from decay.py using class Problem and class Solver from the section Classes for problem and solution method. The parameters I, a, T, the scheme name, and a series of dt values should be read from the command line. Filename: experiment_compare_dt_class.

Exercise 7: Make solid software for a difference equation

We have the following evolutionary difference equation for the number of individuals \(u^n\) of a certain specie at time \(n\Delta t\):

(1)\[ u^{n+1} = u^n + \Delta t r u^n\left(1 - \frac{u^n}{M^n}\right), \quad u^0=U_0{\thinspace .} \\]

Here, \(n\) is a counter in time, \(\Delta t\) is time between time levels \(n\) and \(n+1\) (assumed constant), \(r\) is a net reproduction rate for the specie, and \(M^n\) is the upper limit of the population that the environment can sustain at time level \(n\). Filename: logistic.