.. !split Exercises ========= .. --- begin exercise --- .. _decay:exer:main2func: Exercise 1: Refactor a flat program in terms of a function ---------------------------------------------------------- For simple ODEs of the form .. math:: u' = f(t),\quad u(0)=I,\ t\in (0,T] we can find the solution by straightforward integration: .. math:: u(t) = \int_0^t f(\tau) d\tau{\thinspace .} To compute :math:`u(t)` for :math:`t\in [0,T]`, we introduce a uniform time mesh with points :math:`t_n=n\Delta t` and apply to Trapezoidal rule to approximate the integral. Suppose we have computed the numerical approximation :math:`u^n` to :math:`u(t_n)`. We have .. math:: u(t_{n+1}) = u(t_n) + \int_{t_n}^{t_{n+1}} f(\tau)d\tau{\thinspace .} Using the Trapezoidal rule we get .. math:: u^{n+1} = u^n + \frac{1}{2}\Delta t (f(t_n) + f(t_{n+1})){\thinspace .} The starting value is :math:`u^0=I`. A corresponding implementation for the case :math:`f(t)=2t+1` is given next. .. code-block:: python # f(t) is 2*t + 1 T = 2 from numpy import * dt = 0.2 # time step Nt = int(round(T/dt)) # no of mesh points u = zeros(Nt+1) t = linspace(0, T, Nt+1) for n in range(Nt-1): u[n+1] = u[n] + 0.5*dt*(2*t[n]+1 + 2*t[n+1]+1) This is a flat program. Refactor the program as a function ``solver(f, I, T, dt)``, where ``f`` is the Python implementation of the mathematical function :math:`f(t)` that is to be integrated. The return value of ``solver`` is the pair (``u``, ``t``) representing the solution values and the associated time mesh. Filename: ``integrate.py``. .. Closing remarks for this Exercise Remarks ~~~~~~~ Many prefer to do a first implementation of an algorithm as a flat program and hardcode formulas, here the :math:`f(t)`, into the algorithm. Unfortunately, this coding style makes it difficult to reuse a well-tested algorithm. When the flat program works, it is strongly recommended to *refactor* the code (i.e., rearrange the statements) such that general algorithms are encapsulated in Python functions. The function arguments should be chosen such that the function can be applied for a large class of problems, here all problems that can be expressed as :math:`u'=f(t)`,. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:plot:dtconst: Exercise 2: Compare methods for a given time mesh ------------------------------------------------- Make a program that imports the ``solver`` function from the ``decay_mod`` module and offers a function ``compare(dt, I, a)`` for comparing, in a plot, the methods corresponding to :math:`\theta=0,0.5,1` and the exact solution. This plot shows the accuracy of the methods for a given time mesh. Read input data for the problem from the command line using appropriate functions in the ``decay_mod`` module (the ``--dt`` option for giving several time step values can be reused: just use the first time step value for the computations). Filename: ``decay_compare_theta.py``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:doctest1: Problem 3: Write a doctest -------------------------- Type in the following program and equip the ``roots`` function with a doctest: .. 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) Make sure to test both real and complex roots. Write out numbers with 14 digits or less. Filename: ``doctest_roots.py``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:nosetest1: Problem 4: Write a nose test ---------------------------- Make a nose test for the ``roots`` function in :ref:`decay:exer:doctest1`. Filename: ``test_roots.py``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:module1: Problem 5: Make a module ------------------------ Let .. math:: q(t) = \frac{RAe^{at}}{R + A(e^{at} - 1)} {\thinspace .} Make a Python module ``q_module`` containing two functions ``q(t)`` and ``dqdt(t)`` for computing :math:`q(t)` and :math:`q'(t)`, respectively. Perform a ``from numpy import *`` in this module. Import ``q`` and ``dqdt`` in another file using the "star import" construction ``from q_module import *``. All objects available in this file is given by ``dir()``. Print ``dir()`` and ``len(dir())``. Then change the import of ``numpy`` in ``q_module.py`` to ``import numpy as np``. What is the effect of this import on the number of objects in ``dir()`` in a file that does ``from q_module import *``? Filename: ``q_module.py``. .. \frac{du}{dt}=au\left(1-\frac{u}{R}\right),\quad u(0)=A, .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:decay_class:exper: Exercise 6: Make use of a class implementation ---------------------------------------------- We want to solve the exponential decay problem :math:`u'=-au`, :math:`u(0)=I`, for several :math:`\Delta t` values and :math:`\theta=0,0.5,1`. For each :math:`\Delta t` value, we want to make a plot where the three solutions corresponding to :math:`\theta=0,0.5,1` appear along with the exact solution. Write a function ``experiment`` to accomplish this. The function should import the classes ``Problem``, ``Solver``, and ``Visualizer`` from the `decay_class `__ module and make use of these. A new command-line option ``--dt_values`` must be added to allow the user to specify the :math:`\Delta t` values on the command line (the options ``--dt`` and ``--theta`` implemented by the ``decay_class`` module have then no effect when running the ``experiment`` function). Note that the classes in the ``decay_class`` module should *not* be modified. Filename: ``decay_class_exper.py``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:decay_class2: Exercise 7: Generalize a class implementation --------------------------------------------- Consider the file `decay_class.py `__ where the exponential decay problem :math:`u'=-au`, :math:`u(0)=I`, is implemented via the classes ``Problem``, ``Solver``, and ``Visualizer``. Extend the classes to handle the more general problem .. math:: u'(t) = -a(t)u(t) + b(t),\quad u(0)=I,\ t\in (0,T], using the :math:`\theta`-rule for discretization. In the case with arbitrary functions :math:`a(t)` and :math:`b(t)` the problem class is no longer guaranteed to provide an exact solution. Let the ``u_exact`` in class ``Problem`` return ``None`` if the exact solution for the particular problem is not available. Modify classes ``Solver`` and ``Visualizer`` accordingly. Add test functions ``test_*()`` for the nose testing tool in the module. Also add a demo example where the environment suddenly changes (modeled as an abrupt change in the decay rate :math:`a`): .. math:: a(t) =\left\lbrace\begin{array}{ll} 1, & 0\leq t\leq t_p,\\ k, & t> t_p,\end{array}\right. where :math:`t_p` is the point of time the environment changes. Take :math:`t_p=1` and make plots that illustrate the effect of having :math:`k\gg 1` and :math:`k\ll 1`. Filename: ``decay_class2.py``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:decay_class3: Exercise 8: Generalize an advanced class implementation ------------------------------------------------------- Solve :ref:`decay:exer:decay_class2` by utilizing the class implementations in `decay_class_oo.py `__. Filename: ``decay_class3.py``. .. --- end exercise ---