For simple ODEs of the form
we can find the solution by straightforward integration:
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
Using the Trapezoidal rule we get
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
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
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\)):
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.