.. !split


.. admonition:: Goal

   This document illustrates *best practice* for developing scientific software
   in an efficient and reliable way. Not only will the outlined techniques
   save a lost of human time, but they will also help assure reproducible science
   and higher quality of computational investigations. Key questions to
   be answered are
   
    * How should I organize a program?
   
    * How can I efficiently and safely provide input data and run my code?
   
    * How can I verify that the implementation is correct?
   
    * How should I reliably work with files and documents?
   
    * How should I conduct large numerical experiments?




<font color="red">[<b>hpl 1</b>: <em>Need to cover
functions, classes, modules,
cml, GUI, hand calc, trivial problems, exact num sol, MMS, qualitative results,
Git, bitbucket/github, scripting, report generation</em>]</font>

Sample problem and code
=======================

This first introduction to good programming habits in scientific
computing will make use of a very simple mathematical problem to keep
the mathematical details at the lowest possible level while
introducing a series of computer science concepts. The simplicity of
the mathematical problem obviously prevents us from treating several
techniques that are only meaningful for complex scientific software.

Mathematical problem
--------------------

We consider a linear ordinary differential equation with variable
coefficients:

.. _Eq:softeng1:problem:ode:

.. math::
   :label: softeng1:problem:ode
        
        u'(t) = -a(t)u(t) + b(t),\quad u(0)=U_0,\quad t\in (0,T]{\thinspace .}
        
        

This problem is numerically solved by the so-called :math:`\theta`-rule,
which is a convenient way to merge different formulas for the
well-known Forward Euler, Backward Euler, and Crank-Nicolson
(midpoint/central) schemes. We introduce a uniform time mesh
:math:`t_n=n\Delta t`, :math:`n=0,1,\ldots,N_t`, and seek :math:`u(t)` at the mesh
points. The numerical approximation to :math:`u(t_n)` is denoted
:math:`u^n`. Since we will use the symbol :math:`u` both for the exact analytical
solution of :eq:`softeng1:problem:ode` and for the numerical
approximation, we sometimes introduce :math:`{u_{\small\mbox{e}}}(t)` to help distinguish
the two types of solutions (i.e., subscript e for "exact") [#e-subscript]_.

.. [#e-subscript] In the literature, it is more common to put a subscript
   (like :math:`u_\Delta` or :math:`u_h`)
   on the numerical solution to distinguish it from the exact solution.
   However, we will use the variable ``u`` in the code for the numerical
   approximation to be computed, and therefore adjust the mathematical
   notation to convenient conventions in the code such that we can have
   as close correspondence as possible between the implementation and
   the mathematics.

The :math:`\theta`-rule leads to an explicit updating formula for :math:`u^{n+1}`,
given :math:`u^n`:

<font color="red">[<b>hpl 2</b>: <em>not ready. Need more here!!!!!!!</em>]</font>

Flat program implementation
---------------------------

Function implementation
-----------------------

``from module import *``

Prefixing functions by module names
-----------------------------------

``import module as m``

User interfaces
===============

It is good programming practice to let programs read input from the
user rather than require the user to edit the source code when trying
out new values of input parameters. One reason is that any edit of the
code has a danger of introducing bugs. Another reason is that it is
easier and less manual work to supply data to a program instead of
editing the program code. A third reason is that a program that reads
input can easily be run by another program, and in this way we can
automate a large number of runs in scientific investigations.


.. admonition:: Tip

   We shall make it a habit to equip any implementation of a
   numerical solver with an appropriate user interface before testing out
   the code.




Reading input data can be done in many ways. We have to decide on
desired *user interface*, i.e., how we want to operate the program
when providing input, and then use appropriate tools to implement
the user interface. There are four basic types of user interface
of relevance to our programs, listed here with increasing complexity
of the implementation:

1. Questions and answers in the terminal window

2. Command-line arguments

3. Reading data from file

4. Graphical user interfaces

Although conceptually simple, alternative 1 involves more typing than
the other alternatives and is therefore abandoned. Below, we shall
address alternative 2 and 4, which are most appropriate for the
present problem.

[[[
.. _decay:commandline:

Creating command-line interfaces
--------------------------------

.. index:: user interfaces to programs

.. index:: command-line arguments

Reading input from the command line is a simple and flexible way of interacting
with the user. Python stores all the command-line arguments in
the list ``sys.argv``, and there are, in principle, two ways of programming with
command-line arguments in Python:

 * Decide upon a sequence of parameters on the command line and read
   their values directly from the ``sys.argv[1:]`` list (``sys.argv[0]`` is
   the just program name).

 * Use option-value pairs (``--option value``) on
   the command line to override default values of input parameters,
   and utilize the ``argparse.ArgumentParser`` tool to interact with
   the command line.

Both strategies will be illustrated next.

Reading a sequence of command-line arguments
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. index:: reading the command line

The `decay_plot_mpl.py <http://tinyurl.com/jvzzcfn/decay/decay_plot_mpl.py>`_
program needs the following input data: :math:`I`, :math:`a`, :math:`T`, an option to
turn the plot on or off (``makeplot``), and a list of :math:`\Delta t` values.

The simplest way of reading this input from the command line is to say
that the first four command-line arguments correspond to the first
four points in the list above, in that order, and that the rest of the
command-line arguments are the :math:`\Delta t` values.  The input given for
``makeplot`` can be a string among ``'on'``, ``'off'``, ``'True'``, and
``'False'``. The code for reading this input is most conveniently put in
a function:

.. code-block:: python

        import sys
        
        def read_command_line():
            if len(sys.argv) < 6:
                print 'Usage: %s I a T on/off dt1 dt2 dt3 ...' % \ 
                      sys.argv[0]; sys.exit(1)  # abort
        
            I = float(sys.argv[1])
            a = float(sys.argv[2])
            T = float(sys.argv[3])
            makeplot = sys.argv[4] in ('on', 'True')
            dt_values = [float(arg) for arg in sys.argv[5:]]
        
            return I, a, T, makeplot, dt_values

.. index:: list comprehension

.. index:: sys.argv

One should note the following about the constructions in the program above:

  * Everything on the command line ends up in a *string* in
    the list ``sys.argv``. Explicit conversion to, e.g., a ``float`` object is
    required if the string as a number we want to compute with.

  * The value of ``makeplot`` is determined from a boolean expression,
    which becomes ``True`` if the command-line argument is either ``'on'`` or
    ``'True'``, and ``False`` otherwise.

  * It is easy to build the list of :math:`\Delta t` values: we simply run through
    the rest of the list, ``sys.argv[5:]``, convert each command-line argument
    to ``float``, and collect these ``float`` objects in a list, using the
    compact and convenient *list comprehension* syntax in Python.

The loops over :math:`\theta` and :math:`\Delta t` values can be coded in a ``main`` function:

.. code-block:: python

        def main():
            I, a, T, makeplot, dt_values = read_command_line()
            for theta in 0, 0.5, 1:
                for dt in dt_values:
                    E = explore(I, a, T, dt, theta, makeplot)
                    print '%3.1f %6.2f: %12.3E' % (theta, dt, E)

The complete program can be found in `decay_cml.py <http://tinyurl.com/jvzzcfn/decay/decay_cml.py>`_.

Working with an argument parser
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. index:: argparse (Python module)

.. index:: ArgumentParser (Python class)

.. index:: option-value pairs (command line)

.. index:: command-line options and values

.. index:: reading the command line

Python's ``ArgumentParser`` tool in the ``argparse`` module makes it easy
to create a professional command-line interface to any program. The
documentation of `ArgumentParser <http://docs.python.org/library/argparse.html>`_ demonstrates its
versatile applications, so we shall here just list an example
containing basic features.  On the command line we want to specify
option-value pairs for :math:`I`, :math:`a`, and :math:`T`, e.g., ``--a 3.5 --I 2 --T
2``. Including ``--makeplot`` turns the plot on and excluding this option
turns the plot off.  The :math:`\Delta t` values can be given as ``--dt 1 0.5
0.25 0.1 0.01``.  Each parameter must have a sensible default value so
that we specify the option on the command line only when the default
value is not suitable.

We introduce a function for defining the mentioned command-line options:

.. code-block:: python

        def define_command_line_options():
            import argparse
            parser = argparse.ArgumentParser()
            parser.add_argument('--I', '--initial_condition', type=float,
                                default=1.0, help='initial condition, u(0)',
                                metavar='I')
            parser.add_argument('--a', type=float,
                                default=1.0, help='coefficient in ODE',
                                metavar='a')
            parser.add_argument('--T', '--stop_time', type=float,
                                default=1.0, help='end time of simulation',
                                metavar='T')
            parser.add_argument('--makeplot', action='store_true',
                                help='display plot or not')
            parser.add_argument('--dt', '--time_step_values', type=float,
                                default=[1.0], help='time step values',
                                metavar='dt', nargs='+', dest='dt_values')
            return parser

Each command-line option is defined through the ``parser.add_argument``
method. Alternative options, like the short ``--I`` and the more
explaining version ``--initial_condition`` can be defined. Other arguments
are ``type`` for the Python object type, a default value, and a help
string, which gets printed if the command-line argument ``-h`` or ``--help`` is
included. The ``metavar`` argument specifies the value associated with
the option when the help string is printed. For example, the option for
:math:`I` has this help output:

.. code-block:: text

        Terminal> python decay_argparse.py -h
          ...
          --I I, --initial_condition I
                                initial condition, u(0)
          ...

The structure of this output is

.. code-block:: text

          --I metavar, --initial_condition metavar
                                help-string

The ``--makeplot`` option is a pure flag without any value, implying a
true value if the flag is present and otherwise a false value. The
``action='store_true'`` makes an option for such a flag.

Finally, the ``--dt`` option demonstrates how to allow for more than one
value (separated by blanks) through the ``nargs='+'`` keyword argument.
After the command line is parsed, we get an object where the values of
the options are stored as attributes. The attribute name is specified
by the ``dist`` keyword argument, which for the ``--dt`` option is
``dt_values``. Without the ``dest`` argument, the value of an option ``--opt``
is stored as the attribute ``opt``.

The code below demonstrates how to read the command line and extract
the values for each option:

.. code-block:: python

        def read_command_line():
            parser = define_command_line_options()
            args = parser.parse_args()
            print 'I={}, a={}, T={}, makeplot={}, dt_values={}'.format(
                args.I, args.a, args.T, args.makeplot, args.dt_values)
            return args.I, args.a, args.T, args.makeplot, args.dt_values

The ``main`` function remains the same as in the ``decay_cml.py`` code based
on reading from ``sys.argv`` directly. A complete program featuring the
demo above of ``ArgumentParser`` appears in the file `decay_argparse.py <http://tinyurl.com/jvzzcfn/decay/decay_argparse.py>`_.

Creating a graphical web user interface
---------------------------------------

The Python package `Parampool <https://github.com/hplgit/parampool>`_
can be used to automatically generate a web-based *graphical user interface*
(GUI) for our simulation program. Although the programming technique
dramatically simplifies the efforts to create a GUI, the forthcoming
material on equipping our ``decay_mod`` module with a GUI is quite technical
and of significantly less importance than knowing how to make
a command-line interface (the section :ref:`decay:commandline`).
There is no danger in jumping right to the section :ref:`decay:convergence:rate`.

Making a compute function
~~~~~~~~~~~~~~~~~~~~~~~~~

The first step is to identify a function
that performs the computations and that takes the necessary input
variables as arguments. This is called the *compute function* in
Parampool terminology. We may start with a copy of the basic file
`decay_plot_mpl.py <http://tinyurl.com/jvzzcfn/decay/decay_plot_mpl.py>`_,
which has a ``main`` function displayed in
the section :ref:`decay:plotting` for carrying out simulations and plotting
for a series of :math:`\Delta t` values. Now we want to control and view the same
experiments from a web GUI.

To tell Parampool what type of input data we have,
we assign default values of the right type to all arguments in the
main function and call it ``main_GUI``:

.. code-block:: python

        def main_GUI(I=1.0, a=.2, T=4.0,
                 dt_values=[1.25, 0.75, 0.5, 0.1],
                 theta_values=[0, 0.5, 1]):

The compute function must return the HTML code we want for displaying
the result in a web page. Here we want to show plots of the numerical
and exact solution for different methods and :math:`\Delta t` values.
The plots can be organized in a table with :math:`\theta` (methods) varying
through the columns and :math:`\Delta t` varying through the rows.
Assume now that a new version of the ``explore`` function
not only returns the error ``E`` but also HTML code containing the
plot. Then we can write the ``main_GUI`` function as

.. code-block:: python

        def main_GUI(I=1.0, a=.2, T=4.0,
                 dt_values=[1.25, 0.75, 0.5, 0.1],
                 theta_values=[0, 0.5, 1]):
            # Build HTML code for web page. Arrange plots in columns
            # corresponding to the theta values, with dt down the rows
            theta2name = {0: 'FE', 1: 'BE', 0.5: 'CN'}
            html_text = '<table>\n'
            for dt in dt_values:
                html_text += '<tr>\n'
                for theta in theta_values:
                    E, html = explore(I, a, T, dt, theta, makeplot=True)
                    html_text += """
        <td>
        <center><b>%s, dt=%g, error: %s</b></center><br>
        %s
        </td>
        """ % (theta2name[theta], dt, E, html)
                html_text += '</tr>\n'
            html_text += '</table>\n'
            return html_text

Rather than creating plot files and showing the plot on the screen,
the new version of the ``explore`` function makes a string with the PNG code of
the plot and embeds that string in HTML code. This action is
conveniently performed by Parampool's ``save_png_to_str`` function:

.. code-block:: python

        import matplotlib.pyplot as plt
        ...
        # plot
        plt.plot(t, u, r-')
        plt.xlabel('t')
        plt.ylabel('u')
        ...
        from parampool.utils import save_png_to_str
        html_text = save_png_to_str(plt, plotwidth=400)

Note that we now write ``plt.plot``, ``plt.xlabel``, etc.
The ``html_text`` string is long and contains all the characters that
build up the PNG file of the current plot. The new ``explore``
function can make use of the above code snippet and return
``html_text`` along with ``E``.

Generating the user interface
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The web GUI is automatically generated by
the following code, placed in a file `decay_GUI_generate.py <http://tinyurl.com/jvzzcfn/decay/decay_GUI_generate.py>`_

.. code-block:: python

        from parampool.generator.flask import generate
        from decay_GUI import main
        generate(main,
                 output_controller='decay_GUI_controller.py',
                 output_template='decay_GUI_view.py',
                 output_model='decay_GUI_model.py')

Running the ``decay_GUI_generate.py`` program results in three new
files whose names are specified in the call to ``generate``:

 1. ``decay_GUI_model.py`` defines HTML widgets to be used to set
    input data in the web interface,

 2. ``templates/decay_GUI_views.py`` defines the layout of the web page,

 3. ``decay_GUI_controller.py`` runs the web application.

We only need to run the last program, and there is no need to look into
these files.

Running the web application
~~~~~~~~~~~~~~~~~~~~~~~~~~~

The web GUI is started by

.. code-block:: text

        Terminal> python decay_GUI_controller.py

Open a web browser at the location ``127.0.0.1:5000``. Input fields for
``I``, ``a``, ``T``, ``dt_values``, and ``theta_values`` are presented.  Setting
the latter two to ``[1.25, 0.5]`` and ``[1, 0.5]``, respectively, and
pressing *Compute* results in four plots, see Figure
:ref:`decay:fig:GUI`. With the techniques demonstrated here, one can
easily create a tailored web GUI for a particular type of application
and use it to interactively explore physical and numerical effects.

.. _decay:fig:GUI:

.. figure:: decay_GUI.png
   :width: 800

   *Automatically generated graphical web interface*

Verification
============

Comparison with hand calculations
---------------------------------

One of the simplest and most powerful methods for verifying numerical
codes is to perform some steps of the algorithm by hand and compare the
results with those produced by the code.
In the present case, we may choose some test problem and run three
steps by hand. Picking :math:`a(t)=t^2`... <font color="red">[<b>hpl 3</b>: <em>Not ready. Time-dep :math:`a`?</em>]</font>

Test function
-------------

.. Note: the admon function needs a raw string to handle \theta (\t, \b, etc)


.. admonition:: Caution: choice of parameter values

   For the choice of values of parameters in verification tests one should
   stay away from integers, especially 0 and 1, as these can
   simplify formulas too much for test purposes. For example, with
   :math:`\theta =1` the nominator in the formula for :math:`u^n` will be the same for
   all :math:`a` and :math:`\Delta t` values. One should therefore choose more
   "arbitrary" values, say :math:`\theta =0.8` and :math:`I=0.1`.




Comparison with an exact discrete solution
------------------------------------------

Sometimes it is possible to find a closed-form
*exact discrete solution* that fulfills the discrete finite
difference equations. The implementation can then be verified against
the exact discrete solution. This is usually the best technique for
verification.

.. Not so limited, will later guess that linear functions and MMS can

.. be used in the discrete eqs as well!

Define

.. math::
         A = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a \Delta t}{\thinspace .}  

Manual computations with the :math:`\theta`-rule results in

.. math::
        
        u^0 &= I,\\ 
        u^1 &= Au^0 = AI,\\ 
        u^2 &= Au^1 = A^2I,\\ 
        &\vdots\\ 
        u^n &= A^nu^{n-1} = A^nI {\thinspace .}
        

We have then established the exact discrete solution as

.. _Eq:decay:un:exact:

.. math::
   :label: decay:un:exact
        
        u^n = IA^n
        
        {\thinspace .}
        


.. admonition:: Caution

   One should be conscious about the different meanings of the notation
   on the left- and right-hand side
   of :eq:`decay:un:exact`: on the left, :math:`n` in :math:`u^n`
   is a superscript reflecting a counter
   of mesh points (:math:`t_n`), while on the right, :math:`n`
   is the power in the exponentiation :math:`A^n`.




Comparison of the exact discrete solution and the computed
solution is done in the following function:

.. code-block:: python

        def verify_exact_discrete_solution():
        
            def exact_discrete_solution(n, I, a, theta, dt):
                A = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
                return I*A**n
        
            theta = 0.8; a = 2; I = 0.1; dt = 0.8
            Nt = int(8/dt)  # no of steps
            u, t = solver(I=I, a=a, T=Nt*dt, dt=dt, theta=theta)
            u_de = array([exact_discrete_solution(n, I, a, theta, dt)
                          for n in range(Nt+1)])
            difference = abs(u_de - u).max()  # max deviation
            tol = 1E-15  # tolerance for comparing floats
            success = difference <= tol
            return success

The complete program is found in the file `decay_verf2.py <http://tinyurl.com/jvzzcfn/decay/decay_verf2.py>`_ (``verf2`` is a short name for "verification,
version 2").


.. admonition:: Local functions

   One can define a function inside another function, here called
   a *local function* (also known as *closure*) inside a *parent function*.
   A local function is invisible outside the parent function.
   A convenient property is that any local function has access to all
   variables defined in the parent function, also if we send the
   local function to some other function as argument (!).
   In the present example, it means that the local function
   ``exact_discrete_solution`` does not need its five arguments as the
   values can alternatively be accessed through the local variables defined
   in the parent function ``verify_exact_discrete_solution``. We can send
   such an ``exact_discrete_solution`` without arguments to any other
   function and ``exact_discrete_solution`` will still have access to
   ``n``, ``I``, ``a``, and so forth defined in its parent function.




.. _decay:convergence:rate:

Computing convergence rates
---------------------------

.. index:: convergence rate

We expect that the error :math:`E` in the numerical solution is
reduced if the mesh size :math:`\Delta t` is decreased. More specifically,
many numerical methods obey a power-law relation between :math:`E` and
:math:`\Delta t`:

.. _Eq:decay:E:dt:

.. math::
   :label: decay:E:dt
        
        E = C\Delta t^r,
        
        

where :math:`C` and :math:`r` are (usually unknown) constants independent of :math:`\Delta t`.
The formula :eq:`decay:E:dt` is viewed as an asymptotic model valid for
sufficiently small :math:`\Delta t`. How small is normally hard to estimate
without doing numerical estimations of :math:`r`.

The parameter :math:`r` is known as the *convergence rate*. For example,
if the convergence rate is 2, halving :math:`\Delta t` reduces the error by
a factor of 4. Diminishing :math:`\Delta t` then has a greater impact on
the error compared with methods that have :math:`r=1`. For a given value of :math:`r`,
we refer to the method as of :math:`r`-th order. First- and second-order
methods are most common in scientific computing.

Estimating :math:`r`
~~~~~~~~~~~~~~~~~~~~

There are two alternative ways of estimating :math:`C` and :math:`r` based on a set of
:math:`m` simulations with corresponding pairs :math:`(\Delta t_i, E_i)`, :math:`i=0,\ldots,m-1`,
and :math:`\Delta t_{i} < \Delta t_{i-1}` (i.e., decreasing cell size).

 1. Take the logarithm of :eq:`decay:E:dt`, :math:`\ln E = r\ln \Delta t + \ln C`,
    and fit a straight line to the data points :math:`(\Delta t_i, E_i)`,
    :math:`i=0,\ldots,m-1`.

 2. Consider two consecutive experiments, :math:`(\Delta t_i, E_i)` and
    :math:`(\Delta t_{i-1}, E_{i-1})`. Dividing the equation
    :math:`E_{i-1}=C\Delta t_{i-1}^r` by :math:`E_{i}=C\Delta t_{i}^r` and solving
    for :math:`r` yields

.. _Eq:decay:conv:rate:

.. math::
   :label: decay:conv:rate
        
        r_{i-1} = \frac{\ln (E_{i-1}/E_i)}{\ln (\Delta t_{i-1}/\Delta t_i)}
        
        

for :math:`i=1,\ldots,m-1`.

The disadvantage of method 1 is that :eq:`decay:E:dt` might not be valid
for the coarsest meshes (largest :math:`\Delta t` values). Fitting a line
to all the data points is then misleading.  Method 2 computes
convergence rates for pairs of experiments and allows us to see
if the sequence :math:`r_i` converges to some value as :math:`i\rightarrow m-2`.
The final :math:`r_{m-2}` can then be taken as the convergence rate.
If the coarsest meshes have a differing rate, the corresponding
time steps are probably too large for :eq:`decay:E:dt` to be valid.
That is, those time steps lie outside the asymptotic range of
:math:`\Delta t` values where the error behaves like :eq:`decay:E:dt`.

Implementation
~~~~~~~~~~~~~~

It is straightforward to extend the ``main`` function in the program
``decay_argparse.py`` with statements for computing :math:`r_0, r_1, \ldots, r_{m-2}`
from :eq:`decay:E:dt`:

.. code-block:: python

        from math import log
        
        def main():
            I, a, T, makeplot, dt_values = read_command_line()
            r = {}  # estimated convergence rates
            for theta in 0, 0.5, 1:
                E_values = []
                for dt in dt_values:
                    E = explore(I, a, T, dt, theta, makeplot=False)
                    E_values.append(E)
        
                # Compute convergence rates
                m = len(dt_values)
                r[theta] = [log(E_values[i-1]/E_values[i])/
                            log(dt_values[i-1]/dt_values[i])
                            for i in range(1, m, 1)]
        
            for theta in r:
                print '\nPairwise convergence rates for theta=%g:' % theta
                print ' '.join(['%.2f' % r_ for r_ in r[theta]])
            return r

The program containing this ``main`` function is called `decay_convrate.py <http://tinyurl.com/jvzzcfn/decay/decay_convrate.py>`_.

.. index:: dictionary

The ``r`` object is a *dictionary of lists*. The keys in this
dictionary are the :math:`\theta` values. For example,
``r[1]`` holds the list of the :math:`r_i` values corresponding to
:math:`\theta=1`. In the loop ``for theta in r``, the loop variable ``theta``
takes on the values of the keys in the dictionary ``r`` (in an
undetermined ordering). We could simply do a ``print r[theta]``
inside the loop, but this would typically yield output of
the convergence rates with 16 decimals:

.. code-block:: text

        [1.331919482274763, 1.1488178494691532, ...]

Instead, we format each number with 2 decimals, using a list
comprehension to turn the list of numbers, ``r[theta]``, into
a list of formatted strings. Then we join these strings
with a space in between to get a sequence of rates on one line
in the terminal window. More generally, ``d.join(list)`` joins the
strings in the list ``list`` to one string, with ``d``
as delimiter between ``list[0]``, ``list[1]``, etc.

Here is an example on the outcome of the convergence rate computations:

.. code-block:: text

        Terminal> python decay_convrate.py --dt 0.5 0.25 0.1 0.05 0.025 0.01
        ...
        Pairwise convergence rates for theta=0:
        1.33 1.15 1.07 1.03 1.02
        
        Pairwise convergence rates for theta=0.5:
        2.14 2.07 2.03 2.01 2.01
        
        Pairwise convergence rates for theta=1:
        0.98 0.99 0.99 1.00 1.00

The Forward and Backward Euler methods seem to have an :math:`r` value which
stabilizes at 1, while the Crank-Nicolson seems to be a second-order
method with :math:`r=2`.

.. index:: verification

Very often, we have some theory that predicts what :math:`r` is for a numerical
method. Various theoretical error measures for the :math:`\theta`-rule point to
:math:`r=2` for :math:`\theta =0.5` and :math:`r=1` otherwise. The computed estimates of :math:`r` are
in very good agreement with these theoretical values.


.. admonition:: Why convergence rates are important

   The strong practical application of computing convergence rates is for
   verification: wrong convergence rates point to errors in the code, and
   correct convergence rates brings evidence that the implementation is
   correct. Experience shows that bugs in the code easily destroy the
   expected convergence rate.




Debugging via convergence rates
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Let us experiment with bugs and see the implication on the convergence
rate. We may, for instance, forget to multiply by ``a`` in the denominator
in the updating formula for ``u[n+1]``:

.. code-block:: python

        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt)*u[n]

Running the same ``decay_convrate.py`` command as above gives the expected
convergence rates (!). Why? The reason is that we just specified
the :math:`\Delta t` values are relied on default values for other
parameters. The default value of :math:`a` is 1. Forgetting the factor
``a`` has then no effect. This example shows how important it is to
avoid parameters that are 1 or 0 when verifying implementations.
Running the code ``decay_v0.py`` with :math:`a=2.1` and :math:`I=0.1` yields

.. code-block:: text

        Terminal> python decay_convrate.py --a 2.1 --I 0.1  \ 
                  --dt 0.5 0.25 0.1 0.05 0.025 0.01
        ...
        Pairwise convergence rates for theta=0:
        1.49 1.18 1.07 1.04 1.02
        
        Pairwise convergence rates for theta=0.5:
        -1.42 -0.22 -0.07 -0.03 -0.01
        
        Pairwise convergence rates for theta=1:
        0.21 0.12 0.06 0.03 0.01

This time we see that the expected convergence rates for the Crank-Nicolson and
Backward Euler methods are not obtained, while :math:`r=1` for the Forward Euler
method. The reason for correct rate in the latter case is that :math:`\theta=0`
and the wrong ``theta*dt`` term in the denominator vanishes anyway.

The error

.. code-block:: python

        u[n+1] = ((1-theta)*a*dt)/(1 + theta*dt*a)*u[n]

manifests itself through wrong rates :math:`r\approx 0` for all three methods.
About the same results arise from an erroneous initial condition, ``u[0] = 1``,
or wrong loop limits, ``range(1,Nt)``. It seems that in this simple
problem, most bugs we can think of are detected by the convergence rate
test, provided the values of the input data do not hide the bug.

A ``verify_convergence_rate`` function could compute the dictionary of
list via ``main`` and check if the final rate estimates (:math:`r_{m-2}`)
are sufficiently close to the expected ones. A tolerance of 0.1
seems appropriate, given the uncertainty in estimating :math:`r`:

.. code-block:: python

        def verify_convergence_rate():
            r = main()
            tol = 0.1
            expected_rates = {0: 1, 1: 1, 0.5: 2}
            for theta in r:
                r_final = r[theta][-1]
                diff = abs(expected_rates[theta] - r_final)
                if diff > tol:
                    return False
            return True  # all tests passed

We remark that ``r[theta]`` is a list and the last element in any list
can be extracted by the index ``-1``.