.. !split Exercises (1) ====================== .. --- begin exercise --- .. _decay:exer:meshfunc: 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 :math:`f(t_0),\ldots,f(t_{N_t})`, where ``f`` is a Python function implementing a mathematical function ``f(t)`` and :math:`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. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **b)** Use ``mesh_function`` to compute the mesh function corresponding to .. math:: 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 :math:`t_n=n\Delta t` with :math:`\Delta t=0.1`. Plot the mesh function. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``mesh_function``. .. Closing remarks for this Exercise Remarks (1) ~~~~~~~~~~~~~~~~~~~~ In the section :ref:`decay:computing:error` 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 :math:`u^{\prime}=-au`, :math:`u(0)=1`, for :math:`t\in [0,4]`, where :math:`a=1` for :math:`t\in [0,3]` and :math:`a=3` for :math:`t\in [3,4]`. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:dudt: Problem 1.2: Differentiate a function ------------------------------------- .. index:: array arithmetics .. index:: array computing .. index:: vectorization Given a mesh function :math:`u^n` as an array ``u`` with :math:`u^n` values at mesh points :math:`t_n=n\Delta t`, the discrete derivative can be based on centered differences: .. _Eq:decay:exer:dudt:D2t: .. math:: \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: .. math:: d^0 = [D_t^+u]^n = \frac{u^{1}-u^{0}}{\Delta t}, and .. math:: 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 :math:`d^n` of the mesh function :math:`u^n`. The parameter ``dt`` reflects the mesh spacing :math:`\Delta t`. Write a corresponding test function ``test_differentiate()`` for verifying the implementation. .. --- begin hint in exercise --- **Hint.** The three differentiation formulas are exact for quadratic polynomials. Use this property to verify the program. .. --- end hint in exercise --- .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) **b)** A standard implementation of the formula :ref:`(58) ` is to have a loop over :math:`i`. For large :math:`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 .. math:: u^+ &= (u^2,u^3,\ldots,u^{N_t}), u^- &= (u^0,u^1,\ldots,u^{N_t-2}){\thinspace .} The formula :ref:`(58) ` can now be expressed as .. math:: (d^1,d^2,\ldots,d^{N_t-1}) = \frac{1}{2\Delta t}(u^+ - u^-){\thinspace .} The corresponding Python code reads .. code-block:: python 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``. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``differentiate``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:intdiv: Problem 1.3: Experiment with divisions -------------------------------------- Explain what happens in the following computations, where some are mathematically unexpected: .. code-block:: ipy >>> dt = 3 >>> T = 8 >>> Nt = T/dt >>> Nt 2 >>> theta = 1; a = 1 >>> (1 - (1-theta)*a*dt)/(1 + theta*dt*a) 0 .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``pyproblems``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:decay1err: Problem 1.4: Experiment with wrong computations ----------------------------------------------- Consider the ``solver`` function in the `decay_v1.py `__ file and the following call: .. code-block:: python u, t = solver(I=1, a=1, T=7, dt=2, theta=1) The output becomes .. code-block:: text 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 :math:`t=6` only even though we specified :math:`T=7`. .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``decay_v1_err``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:plot:error: Problem 1.5: Plot the error function ------------------------------------ Solve the problem :math:`u'=-au`, :math:`u(0)=I`, using the Forward Euler, Backward Euler, and Crank-Nicolson schemes. For each scheme, plot the error mesh function :math:`e^n = {u_{\small\mbox{e}}}(t_n)-u^n` for :math:`\Delta t=0.1, 0.05, 0.025`, where :math:`{u_{\small\mbox{e}}}` is the exact solution of the ODE and :math:`u^n` is the numerical solution at mesh point :math:`t_n`. .. --- begin hint in exercise --- **Hint.** Modify the `decay_plot_mpl.py `__ code. .. --- end hint in exercise --- .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``decay_plot_error``. .. --- end exercise --- .. --- begin exercise --- .. _decay:exer:inexact:output: 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 .. code-block:: text 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 .. code-block:: text 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 .. code-block:: text 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? .. removed !bsol ... !esol environment (because of the command-line option --without_solutions) Filename: ``decay_memsave_v2``. .. --- end exercise ---