a)
Write a function mesh_function(f, t) that returns an array with
mesh point values \( f(t_0),\ldots,f(t_{N_t}) \), where f is a Python
function implementing a mathematical function f(t) and \( t_0,\ldots,t_{N_t} \)
are mesh points stored in the array t. Use a loop over the mesh
points and compute one mesh function value at the time.
b)
Use mesh_function to compute the mesh function corresponding to
$$
f(t) = \left\lbrace
\begin{array}{ll}
e^{-t},& 0\leq t\leq 3,\\
e^{-3t}, & 3 < t\leq 4
\end{array}\right.
$$
Choose a mesh \( t_n=n\Delta t \) with \( \Delta t=0.1 \).
Plot the mesh function.
Filename: mesh_function.
In the section Computing the numerical error as a mesh function we show how easy it is to
compute a mesh function by array arithmetics (or array computing).
Using this technique, one could simply implement mesh_function(f,t)
as return f(t). However, f(t) will not work if there are
if tests involving t inside f as is the case in b). Typically,
if t < 3 must have t < 3 as a boolean expression, but if t is
array, t < 3, is an array of boolean values, which is not legal
as a boolean expression in an if test.
Computing one element
at a time as suggested in a) is a way of out of this problem.
We also remark that the function in b) is the solution of \( u^{\prime}=-au \), \( u(0)=1 \), for \( t\in [0,4] \), where \( a=1 \) for \( t\in [0,3] \) and \( a=3 \) for \( t\in [3,4] \).
Given a mesh function \( u^n \) as an array u with \( u^n \) values at mesh
points \( t_n=n\Delta t \), the discrete derivative can be based on
centered differences:
$$
\begin{equation}
d^n = [D_{2t}u]^n =
\frac{u^{n+1}-u^{n-1}}{2\Delta t},\quad n=1,\ldots,N_t-1\tp
\tag{1.58}
\end{equation}
$$
At the end points we use forward and backward differences:
$$ d^0 = [D_t^+u]^n = \frac{u^{1}-u^{0}}{\Delta t},$$
and
$$ d^{N_t} = [D_t^-u]^n = \frac{u^{N_t}-u^{N_t-1}}{\Delta t}\tp$$
a)
Write a function
differentiate(u, dt) that returns the discrete derivative \( d^n \) of the
mesh function \( u^n \). The parameter dt reflects the
mesh spacing \( \Delta t \). Write a corresponding test function
test_differentiate() for verifying the implementation.
The three differentiation formulas are exact for quadratic polynomials. Use this property to verify the program.
b) A standard implementation of the formula (1.58) is to have a loop over \( i \). For large \( N_t \), such loop may run slowly in Python. A technique for speeding up the computations, called vectorization or array computing, replaces the loop by array operations. To see how this can be done in the present mathematical problem, we define two arrays $$ \begin{align*} u^+ &= (u^2,u^3,\ldots,u^{N_t}), u^- &= (u^0,u^1,\ldots,u^{N_t-2})\tp \end{align*} $$ The formula (1.58) can now be expressed as $$ (d^1,d^2,\ldots,d^{N_t-1}) = \frac{1}{2\Delta t}(u^+ - u^-)\tp$$ The corresponding Python code reads
d[1:-1] = (u[2:] - u[0:-2])/(2*dt)
# or
d[1:N_t] = (u[2:N_t+1] - u[0:N_t-1])/(2*dt)
Recall that an array slice u[1:-1] contains the elements in u starting
with index 1 and going all indices up to, but not including, the last one
(-1).
Use the ideas above to implement a vectorized version of the
differentiate function without loops. Make a corresponding
test function that compares the result with that of
differentiate.
Filename: differentiate.
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.
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 some variable 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.
Solve the problem \( u'=-au \), \( u(0)=I \), using the Forward Euler, Backward Euler, and Crank-Nicolson schemes. For each scheme, plot the error mesh function \( e^n = \uex(t_n)-u^n \) for \( \Delta t=0.1, 0.05, 0.025 \), where \( \uex \) is the exact solution of the ODE and \( u^n \) is the numerical solution at mesh point \( t_n \).
Modify the decay_plot_mpl.py code.
Filename: decay_plot_error.
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
If you have just modified the formatting of numbers in the file, running the modified program
Terminal> python decay_memsave_v2.py --T 10 --theta 1 \
--dt 0.2 --makeplot
leads to printing of the message Bug in the implementation! in the
terminal window. Why?
Filename: decay_memsave_v2.