Now we would like to consider the time dependent heat equation \begin{equation*} u_t = \kappa \nabla \cdot \nabla u \end{equation*} in 1 spatial dimension. Note that \(\kappa\) is the diffusivity (the rate of diffusion) so in terms of physical problems, if \(\kappa\) is small then the diffusion occurs slowly and if \(\kappa\) is large then the diffusion occurs quickly.

In 1 spatial dimension, the heat equation is simply \begin{equation*} u_t = \kappa u_{xx} \end{equation*} and we can approximate the derivatives with an Euler-type approximation of the time and a central difference in space: \begin{equation*} \frac{U_i^{n+1} - U_i^n}{\Delta t} = \kappa \left( \frac{U_{i+1}^n - 2U_i^n + U_{i+1}^n}{\Delta x^2} \right). \end{equation*}

Here we are taking \(U_i^n \approx u(t_n,x_i)\) (superscripts represent time step and subscripts represent spatial steps). Rearranging we see that \begin{gather*} U_i^{n+1} = U_i^n + \frac{\kappa \Delta t}{\Delta x^2} \left( U_{i+1}^n - 2U_i^n + U_{i+1}^n \right). \end{gather*}

Implement <<Unresolved xref, reference "eqn_heat_euler"; check spelling or use "provisional" attribute>> in MATLAB to approximate the solution to the following problem: \begin{equation*} \text{ Solve: } u_t = 0.5u_{xx} \text{ with } x \in (0,1), \, u(0,x) = \sin(2 \pi x), \, u(t,0) = 0, \, \text{ and } \, u(t,1) = 0. \end{equation*}

You may have noticed in the previous problem that you will have terribly unstable solutions for certain choices of \(\Delta x\) and \(\Delta t\). Set \(\kappa = 1\) in the previous problem and experiment with choices for \(\Delta x\) and \(\Delta t\) to find where <<Unresolved xref, reference "eqn_heat_euler"; check spelling or use "provisional" attribute>> gives a stable numerical solution to the heat equation. For each choice of \(\Delta x\) and \(\Delta t\) report the value of \(\frac{\Delta t}{\Delta x^2}\).

The instabilities of the heat equation with and Euler-type time discretization and a central differencing scheme is maddening. Thankfully, we can avoid this issue almost entirely by considering an implicit scheme called the Crank-Nicolson method. In this method we approximate the temporal derivative with an Euler-type approximation, but we approximate the spatial derivative as the average of the central difference at the old time step and the central difference at the new time step. That is: \begin{equation*} \frac{U_i^{n+1} - U_i^n}{\Delta t} = \frac{1}{2} \left[\kappa \left( \frac{U_{i+1}^n - 2U_i^n + U_{i+1}^n}{\Delta x^2}\right) +\kappa \left(\frac{U_{i+1}^{n+1} - 2U_i^{n+1} + U_{i+1}^{n+1}}{\Delta x^2} \right) \right]. \end{equation*}

Letting \(r = \kappa \Delta t / (2\Delta x^2)\) we can rearrange to get \begin{equation*} -r U_{i-1}^{n+1} + (1+2r) U_{i}^{n+1} - r U_{i+1}^{n+1} = r U_{i-1}^{n} + (1-2r) U_{i}^{n} + r U_{i+1}^{n}. \end{equation*}

This can now be viewed as a system of equation that can be solved for all \(U_j^{n+1}\) all at once with one linear solve. Write MATLAB code to solve the heat equation from the previous problems with the Crank-Nicolson method.

In this problem we will solve a more realistic 1D heat equation. We will allow the diffusivity to change spatially, so \(\kappa = \kappa(x)\) and we want to solve \begin{equation*} u_t = \left( \kappa(x) u'(x) \right)' \end{equation*} on \(x \in (0,1)\) with Dirichlet boundary conditions \(u(t,0) = u(t,1) = 0\) and initial condition \(u(0,x) = \sin(2 \pi x)\). In this problem we will take \(\kappa(x)\) to be the parabola \(\kappa(x)= x(1-x)\). We start by doing some calculus to rewrite the differential equation: \begin{equation*} u_t = \kappa(x) u''(x) + \kappa'(x) u'(x). \end{equation*}

Your job is:

1. Write an explicit scheme to solve this problem by using centered differences for the spatial derivatives and an Euler-type discretization for the temporal derivative. Write a clear and thorough explanation for how you are doing the discretization as well as a discussion for the errors that are being made with each discretization.

2. Write a MATLAB script to find an approximate solution to this problem. Upload this script, but be sure that it is well-commented.

3. Write a clear and thorough discussion about how your will choose \(\Delta x\) and \(\Delta t\) to give stable solutions to this equation.

4. Graphically compare your solution to this problem with a heat equation where \(\kappa\) is taken to be the average diffusivity: \(\kappa = 0.5.\) How does the changing diffusivity change the shape of the solution?

Write code to solve the 2D Dirichlet problem \begin{equation*} \text{ Solve: } u_t = \kappa \nabla \cdot \nabla u \text{ in } x \in \Omega \text{ with } u(0,x,y) = 2000 (x-0.5)(y-0.5)\text{ exp } \left( -\frac{(x-0.5)^2 + (y-0.5)^2}{0.02} \right) \end{equation*} subject to homogeneous Dirichlet boundary conditions.

For simplicity we suggest that you take \(\Delta x = \Delta y\). You should also be very careful of the stability conditions for the heat equation.