$$ \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\vexd}[1]{{v_{\small\mbox{e}, #1}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$ « Previous Next »

Exercises

Exercise 1: Refactor a flat program in terms of a function

For simple ODEs of the form $$ u' = f(t),\quad u(0)=I,\ t\in (0,T]$$ we can find the solution by straightforward integration: $$ u(t) = \int_0^t f(\tau) d\tau\tp$$ To compute \( u(t) \) for \( t\in [0,T] \), we introduce a uniform time mesh with points \( t_n=n\Delta t \) and apply to Trapezoidal rule to approximate the integral. Suppose we have computed the numerical approximation \( u^n \) to \( u(t_n) \). We have $$ u(t_{n+1}) = u(t_n) + \int_{t_n}^{t_{n+1}} f(\tau)d\tau\tp$$ Using the Trapezoidal rule we get $$ \begin{equation} u^{n+1} = u^n + \frac{1}{2}\Delta t (f(t_n) + f(t_{n+1}))\tp \end{equation} $$ The starting value is \( u^0=I \). A corresponding implementation for the case \( f(t)=2t+1 \) is given next. 

# 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 \( 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.

Remarks

Many prefer to do a first implementation of an algorithm as a flat program and hardcode formulas, here the \( 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 \( u'=f(t) \),.

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 \( \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.

Problem 3: Write a doctest

Type in the following program and equip the roots function with a doctest: 

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.

Problem 4: Write a nose test

Make a nose test for the roots function in Problem 3: Write a doctest. Filename: test_roots.py.

Problem 5: Make a module

Let $$ q(t) = \frac{RAe^{at}}{R + A(e^{at} - 1)} \tp $$ Make a Python module q_module containing two functions q(t) and dqdt(t) for computing \( q(t) \) and \( 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.

Exercise 6: Make use of a class implementation

We want to solve the exponential decay problem \( u'=-au \), \( u(0)=I \), for several \( \Delta t \) values and \( \theta=0,0.5,1 \). For each \( \Delta t \) value, we want to make a plot where the three solutions corresponding to \( \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 \( \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.

Exercise 7: Generalize a class implementation

Consider the file decay_class.py where the exponential decay problem \( u'=-au \), \( u(0)=I \), is implemented via the classes Problem, Solver, and Visualizer. Extend the classes to handle the more general problem $$ u'(t) = -a(t)u(t) + b(t),\quad u(0)=I,\ t\in (0,T],$$ using the \( \theta \)-rule for discretization.

In the case with arbitrary functions \( a(t) \) and \( 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 \( a \)): $$ a(t) =\left\lbrace\begin{array}{ll} 1, & 0\leq t\leq t_p,\\ k, & t> t_p,\end{array}\right. $$ where \( t_p \) is the point of time the environment changes. Take \( t_p=1 \) and make plots that illustrate the effect of having \( k\gg 1 \) and \( k\ll 1 \). Filename: decay_class2.py.

Exercise 8: Generalize an advanced class implementation

Solve Exercise 7: Generalize a class implementation by utilizing the class implementations in decay_class_oo.py. Filename: decay_class3.py.

« Previous Next »