Problem

Solve the heat equation on the unit square (square.m) with homogeneous Dirichlet boundary conditions, $a=1$, initial condition $u_0=0$, $T=\pi/2$, and the right-hand side given by

$$f(x,t) = x_1 (1-x_1) x_2 (1-x_2) \cos t + 2 ( x_1 (1-x_1) + x_2 (1-x_2) ) \sin t,$$ (5)

using a time step of size $k = T/10$.

Compare your computed solution $U(x,t)$ with the exact solution $u(x,t)$. How large is the error in the maximum norm at the final time? The maximum norm of the error is given by

$$\Vert e(T)\Vert _{\infty} = \Vert U(T) - u(T)\Vert _{\infty} = \max_{x\in\Omega} \vert U(x,T) - u(x,T)\vert.$$ (6)

How large is the error when you refine the mesh (square_refined.m)? How large is the error when you divide the time step with a factor $2$?

To help you solve the problem, we will guide you through the implementation, step by step. Don't look at these instructions if you want to solve the problem in your own way. (You could of course glance occasionally at the instructions.)

1. We start by implementing the time-stepping loop. Write a loop that takes $10$ time steps and updates a variable called time, starting at time = 0. Note that you should name the time variable time instead of t, since t is used for the mesh.
2. Next, we add two vectors U0 and U1, representing the solution at time $t=t_{n-1}$ and time $t=t_n$ respectively. Initialize both vectors to zero (of the correct size), and update the vector U0 to be equal to U1 at the beginning of every time step. Remember to load the mesh at the beginning of the program.
3. We now need to assemble the matrix A for the left-hand side and the vector b for the right-hand side of the variational formulation. The matrix is not time-dependent, so it is enough if we assemble it one time, before the time-stepping begins. Since the vector depends on time ($f$ is time-dependent), we need to assemble it in each step of the loop.

   Note that the old value $U_{n-1}$ (given by U0) is present in the right-hand side of the variational formulation (2). We thus need to supply this value to the assembler. This is done by using the extra argument W of the program AssembleVector.m. In your program Heat.m, you then obtain the value of $U_{n-1}$ in the variable w.

4. Now solve the linear system to compute the new value of U1 in each step of the loop. Compute the exact solution at final time and compute the error.

Check your answer: You should obtain the following errors in the maximum norm for the two different meshes and the different values of $k$:

	$k = T/10$	$k = T/20$
grid.m	0.000998	0.000876
grid_refined.m	0.000445	0.000318