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):
$$ \begin{equation} {dT\over dt} = -k(T-T_s),\quad T(0)=T_0\tp \tag{74} \end{equation} $$ 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.
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.
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 (74) is an
appropriate mathematical model for the evolution of the temperature in
the body. First, determine \( k \) in (74) 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.
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.
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.
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 = \uex(t_n)-u^n \) for \( \Delta t \), \( \frac{1}{4}\Delta t \), and
\( \frac{1}{8}\Delta t \), where \( \uex \) is the exact solution of the ODE and
\( u^n \) is the numerical solution at mesh point \( t_n \).
Filename: decay_plot_error.py.
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.
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.
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 9: 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 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 \)):
$$ 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 13: Generalize a class implementation by utilizing the
class implementations in
decay_class_oo.py.
Filename: decay_class3.py.