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.
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) \),.
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.
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.
Make a nose test for the roots function in Problem 3: Write a doctest.
Filename: test_roots.py.
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.
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.
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.
Solve Exercise 7: Generalize a class implementation by utilizing the
class implementations in
decay_class_oo.py.
Filename: decay_class3.py.