.. !split

Exercises
=========

.. --- begin exercise ---

.. _softeng1:exer:derivative:

Problem 1: Make a tool for differentiating curves
-------------------------------------------------

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

.. math::
        
        d_0 = \frac{y_1-y_0}{\Delta x},
        

.. math::
          
        d_i = \frac{y_{i+1}-y_{i-1}}{2\Delta x},\quad i=1,\ldots,n-1,
        

.. math::
          
        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 :math:`d_i`, :math:`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 :math:`d_i` are exact for points on
   a straight line, as all the :math:`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 :math:`d_i`, :math:`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 :math:`y=\sin(\pi x)` and :math:`n+1`
points in :math:`[0,1]`. Apply ``differentiate`` four times and plot the
resulting curve and the exact :math:`y=\sin\pi x` for :math:`n=6, 11, 21, 41`.

.. Using a 2nd-order backward formula at x=1 does not improve the

.. results much, one gets large errors at the end points.

Filename: ``curvediff``.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:integral:flat:

Problem 2: Make solid software for the Trapezoidal rule
-------------------------------------------------------

An integral

.. math::
         \int_a^b f(x)dx 

can be numerically approximated by the Trapezoidal rule,

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

where :math:`x_i` is a set of uniformly spaced points in :math:`[a,b]`:

.. math::
         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 :math:`\int_0^\pi \sin^2x\, dx`:

.. code-block:: python

        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 :math:`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 :math:`\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:

.. code-block:: text

        Terminal> py.test -s -v trapezoidal

.. --- begin hint in exercise ---

**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 :math:`f(x)`.

.. --- end hint in exercise ---

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

.. code-block:: text

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

Here, :math:`a=0`, :math:`b=\pi`, and :math:`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.

.. --- begin hint in exercise ---

**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:

.. code-block:: python

        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

.. code-block:: python

        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 :math:`a` and :math:`b` to be ``math`` expressions like ``pi/4``
and ``exp(4)``, do

.. code-block:: text

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

.. --- end hint in exercise ---

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

Filename: ``trapezoidal``.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:integral:flat2:

Problem 3: Implement classes for the Trapezoidal rule
-----------------------------------------------------

We consider the same problem setting as in :ref:`softeng1:exer:integral:flat`. 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 :ref:`softeng1:exer:integral:flat`.
Filename: ``trapezoidal_class``.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:doctest1:

Problem 4: Write a doctest and a test function
----------------------------------------------

Type in the following program:

.. code-block:: python

        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``.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:tol:

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 :math:`10^k`, :math:`k=-5,-9,\ldots,20` and ``a=1.0/49*b*49``.
Filename: ``tolerance``.

.. Closing remarks for this Problem

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
:math:`10^6` m, ``tol`` must be about :math:`10^{-9}`, while ``a`` and ``b`` around unity can
have ``tol`` of size :math:`10^{-15}`.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:class:dts:

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 :ref:`softeng1:prog:se:class`.
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``.

.. --- end exercise ---

.. --- begin exercise ---

.. _softeng1:exer:logistic:

Exercise 7: Make solid software for a difference equation
---------------------------------------------------------

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

.. _Eq:softeng1:exer:logistic:eq:

.. math::
   :label: softeng1:exer:logistic:eq
        
        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, :math:`n` is a counter in time, :math:`\Delta t` is time between time levels
:math:`n` and :math:`n+1` (assumed constant), :math:`r` is a net reproduction rate
for the specie,
and :math:`M^n` is the upper limit of the population that the environment can
sustain at time level :math:`n`.
Filename: ``logistic``.

.. --- end exercise ---