Exercises and Problems

Exercise 1: Derive schemes for Newton’s law of cooling

Show in detail how we can apply the ideas of the Forward Euler, Backward Euler, Crank-Nicolson, and \(\theta\)-rule discretizations to derive explicit computational formulas for new temperature values in Newton’s law of cooling (see the section Newton’s law of cooling):

(1)\[ {dT\over dt} = -k(T-T_s),\quad T(0)=T_0{\thinspace .}\]

Here, \(T\) is the temperature of the body, \(T_s\) is the temperature of the surroundings, \(t\) is time, \(k\) is the heat transfer coefficient, and \(T_0\) is the initial temperature of the body.

Filename: schemes_cooling.pdf.

Exercise 2: Implement schemes for Newton’s law of cooling

Formulate a \(\theta\)-rule for the three schemes in Exercise 1: Derive schemes for Newton’s law of cooling such that you can get the three schemes from a single formula by varying the \(\theta\) parameter. Implement the \(\theta\) scheme in a function cooling(T0, k, T_s, t_end, dt, theta=0.5), where T0 is the initial temperature, k is the heat transfer coefficient, T_s is the temperature of the surroundings, t_end is the end time of the simulation, dt is the time step, and theta corresponds to \(\theta\). The cooling function should return the temperature as an array T of values at the mesh points and the time mesh t. Construct verification examples to check that the implementation works.

Hint. For verification, try to find an exact solution of the discrete equations. A trick is to introduce \(u=T-T_s\), observe that \(u^{n}=(T_0-T_s)A^n\) for some amplification factor \(A\), and then express this formula in terms of \(T^n\).

Filename: cooling.py.

Exercise 3: Find time of murder from body temperature

A detective measures the temperature of a dead body to be 26.7 C at 2 pm. One hour later the temperature is 25.8 C. The question is when death occurred.

Assume that Newton’s law of cooling (9.5) is an appropriate mathematical model for the evolution of the temperature in the body. First, determine \(k\) in (9.5) by formulating a Forward Euler approximation with one time steep from time 2 am to time 3 am, where knowing the two temperatures allows for finding \(k\). Assume the temperature in the air to be 20 C. Thereafter, simulate the temperature evolution from the time of murder, taken as \(t=0\), when \(T=37\hbox{ C}\), until the temperature reaches 25.8 C. The corresponding time allows for answering when death occurred. Filename: detective.py.

Exercise 4: Experiment with integer division

Explain what happens in the following computations, where some are mathematically unexpected:

>>> dt = 3
>>> T = 8
>>> Nt = T/dt
>>> Nt
2
>>> theta = 1; a = 1
>>> (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
0

Filename: pyproblems.txt.

Exercise 5: Experiment with wrong computations

Consider the solver function in the decay_v1.py file and the following call:

u, t = solver(I=1, a=1, T=7, dt=2, theta=1)

The output becomes

t= 0.000 u=1
t= 2.000 u=0
t= 4.000 u=0
t= 6.000 u=0

Print out the result of all intermediate computations and use type(v) to see the object type of the result stored in v. Examine the intermediate calculations and explain why u is wrong and why we compute up to \(t=6\) only even though we specified \(T=7\). Filename: decay_v1_err.py.

Exercise 6: Plot the error function

Solve the problem \(u'=-au\), \(u(0)=I\), using the Forward Euler, Backward Euler, and Crank-Nicolson schemes. For each scheme, plot the error function \(e^n = {u_{\small\mbox{e}}}(t_n)-u^n\) for \(\Delta t\), \(\frac{1}{4}\Delta t\), and \(\frac{1}{8}\Delta t\), where \({u_{\small\mbox{e}}}\) is the exact solution of the ODE and \(u^n\) is the numerical solution at mesh point \(t_n\). Filename: decay_plot_error.py.

Exercise 7: 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.

Exercise 8: Change formatting of numbers and debug

The decay_memsave.py program writes the time values and solution values to a file which looks like

0.0000000000000000E+00  1.0000000000000000E+00
2.0000000000000001E-01  8.3333333333333337E-01
4.0000000000000002E-01  6.9444444444444453E-01
6.0000000000000009E-01  5.7870370370370383E-01
8.0000000000000004E-01  4.8225308641975323E-01
1.0000000000000000E+00  4.0187757201646102E-01
1.2000000000000000E+00  3.3489797668038418E-01
1.3999999999999999E+00  2.7908164723365347E-01

Modify the file output such that it looks like

0.000  1.00000
0.200  0.83333
0.400  0.69444
0.600  0.57870
0.800  0.48225
1.000  0.40188
1.200  0.33490
1.400  0.27908

Run the modified program

Terminal> python decay_memsave_v2.py --T 10 --theta 1 \
          --dt 0.2 --makeplot

The program just prints Bug in the implementation! and does not show the plot. What went wrong? Filename: decay_memsave_v2.py.

Problem 9: 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 10: Write a nose test

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

Problem 11: Make a module

Let

\[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 \(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 12: 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 13: 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 exact_solution 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\)):

\[\begin{split}a(t) =\left\lbrace\begin{array}{ll} 1, & 0\leq t\leq t_p,\\ k, & t> t_p,\end{array}\right.\end{split}\]

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 14: Generalize an advanced class implementation

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