Exercises

Exercise 1.1: Define a mesh function and visualize it

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

\[\begin{split}f(t) = \left\lbrace \begin{array}{ll} e^{-t},& 0\leq t\leq 3,\\ e^{-3t}, & 3 < t\leq 4 \end{array}\right.\end{split}\]

Choose a mesh \(t_n=n\Delta t\) with \(\Delta t=0.1\). Plot the mesh function.

Filename: mesh_function.

Remarks

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]\).

Problem 1.2: Differentiate a function

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:

\[\tag{58} d^n = [D_{2t}u]^n = \frac{u^{n+1}-u^{n-1}}{2\Delta t},\quad n=1,\ldots,N_t-1{\thinspace .}\]

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}{\thinspace .}\]

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.

Hint. The three differentiation formulas are exact for quadratic polynomials. Use this property to verify the program.

b) A standard implementation of the formula (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{split}u^+ &= (u^2,u^3,\ldots,u^{N_t}), u^- &= (u^0,u^1,\ldots,u^{N_t-2}){\thinspace .}\end{split}\]

The formula (58) can now be expressed as

\[(d^1,d^2,\ldots,d^{N_t-1}) = \frac{1}{2\Delta t}(u^+ - u^-){\thinspace .}\]

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.

Problem 1.3: Experiment with divisions

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.

Problem 1.4: 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 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.

Problem 1.5: 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 mesh function \(e^n = {u_{\small\mbox{e}}}(t_n)-u^n\) for \(\Delta t=0.1, 0.05, 0.025\), where \({u_{\small\mbox{e}}}\) is the exact solution of the ODE and \(u^n\) is the numerical solution at mesh point \(t_n\).

Hint. Modify the decay_plot_mpl.py code.

Filename: decay_plot_error.

Problem 1.6: 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

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.