Homework 2

Due: Thursday, October 11th in class

Monte Carlo Sampling

Note: Read the lecture notes here before doing this assignment.

Description: In this assignment you will consider a stochastic ODE problem and study the performance of the standard Monte Carlo sampling technique (MC) for solving the problem.

For a vector-valued stochastic process \({\bf u}(t,y) = [u_1(t,y),u_2(t,y)]^{\top}\), consider the stochastic ODE problem:

\begin{eqnarray*} \left\{ \begin{array}{l} {\bf u}_t(t,y) = A(y) \, {\bf u}(t,y), \qquad t \in [0,T]\\ \\ {\bf u}(0,y) = [1, 0]^{\top}, \hskip 2.2cm t=0 \end{array} \right. \end{eqnarray*}

where \(A(y)\) is a stochastic matrix given in terms of one uniform random variable \(y\):

$$ A(y) = \left[ \begin{array}{c c} 0 & \, \, 1+y\\ -1-y & \, \, 0 \end{array} \right], \qquad y \sim U[0,1]. $$

Consider the following stochastic quantity of interest (QoI):

$$Q(y) = u_1(T, y).$$

We want to approximate the expected value of \(Q\) (i.e. \({\mathbb E}[Q]\)) by MC.

The stochastic problem can indeed be solved analytically. The analytical solution is given by \(u_1 = \sin (t (1+y))\) and \(u_2 = \cos (t (1+y))\). Hence, the exact expected value of the QoI is

$${\mathbb E}[Q] = (\cos(T)-\cos(2T))/T.$$

This value will be used to measure the total error. Throughout this homework we take \(T=1\).

Tasks:

  1. In the lecture notes, we have seen the cost-error optimization strategy used to obtain the optimal number of samples \(M_{opt}\) and the optimal time step size \(h_{opt}\), while keeping the order \(q\) of the discretization scheme fixed. Perform a similar strategy and find the optimal \(M_{opt}\) and the optimal \(q_{opt}\), while keeping \(h\) fixed.

  2. Download the MATLAB files mc_conv.m and mc.m and ode_taylor.m from the subdirectory hw2 in the public Bitbucket repository https://bitbucket.org/motamed/uq2017 .

    The main code the you are supposed to modify is mc_conv.m. It consists of 3 major parts: 1) Input parameters; 2) Computations; and 3) Plots. It calls the other two codes. The deterministic solver, i.e. the \(q\) -th order accurate Taylor’s method, is implemented in ode_taylor.m. For the following taks, you will need to modify only the main code. Do not modify the other two files, but you are welcome to take a look and try to figure out what they do.

    Let the order of accuracy of deterministic solver be fixed (see the values of \(q\) below). Consider a set of tolerances \(\varepsilon_{TOL} = 20 \times 10^{-4}, 19 \times 10^{-4}, 18\times 10^{-4}, \dotsc, 1 \times 10^{-4}\). In the main code this can be done for instance by writing Eps=(20:-1:1)*(1e-4)}.

    Generate two figures:

    • Figure 1. Let \(q=2\). Plot the total error in MC (y-axis) versus tolerances \(\varepsilon_{TOL}\) (x-asis). Use loglog to get a logarithmic scale on both axes. The plot should have labels on each axis. Use large fonts. Also plot the line \(\varepsilon_{TOL}\) versus \(\varepsilon_{TOL}\) in the same figure window. From your figure you should be able to see if the total error generated by the the MC method remain below \(\varepsilon_{TOL}\) or not. If not, why? Can you explain?
    • Figure 2. Plot the total computational cost of MC (y-axis) versus tolerances \(\varepsilon_{TOL}\) (x-asis) for three different values of \(q=1,2,3\). Plot all three curves in the same figure wondow. Use different markers/colors for the errors generated by different orders \(q\). The plot should have labels on each axis and a legend to describe each curve. Use loglog to get a logarithmic scale on both axes. Use large fonts. Comment and discuss on the convergence rate of the computational cost with respect to tolerance \(\varepsilon_{TOL}\). Compare it with the theoretical results in the lecture notes.

  1. Extra 10 bonus points. If in Figure 1 you observe that the total MC error exceeds the tolerance at some points, can you propose a strategy (in words) to fix or improve the problem? In other words what to do to force the MC error to remain below the tolerance?