In lectures we saw that we can solve the heat equation with periodic boundary conditions:

$$u_t = u_{\theta\theta}, \qquad u(0,\theta) = u_0(\theta), \qquad u \hbox{ is periodic}$$where initial condition is also periodic. This proceeded by either using either `expm`

, `ode45`

, or discretizing the forward Euler method in time

or the backward Euler method in time

$$\mathbf w_{k+1} = (I-Δt D_n^2)^{-1}\mathbf w_k$$In this lecture, we will investigate the heat equation with Dirichlet boundary conditions:

$$u_t = u_{xx}, \qquad u(0,x) = u_0(x), u(t,0) = u(t,1) =0$$where $u_0(x)$ also satisfies the Dircihlet conditions: $u_0(0) = u_0(1) = 0$.

`expm`

and `ode45`

¶Define the $n-1 \times n+1$ matrix

$$\tilde D_n^2 \triangleq D_{n-1} D_n = {1 \over h^2} \begin{pmatrix} 1 & -2 & 1 \cr & 1 & -2 & 1\cr &&\ddots & \ddots & \ddots \cr &&&1 & -2 & 1\end{pmatrix}$$**Question 1(a)** Construct a function `D2(h,n)`

that returns the $n-1 \times n+1$ matrix $\tilde D_n^2$, without checking the lecture notes.

In [ ]:

```
```

**Question 1(b)** Why is it not the case that

when we have Dirichlet boundary conditions? Try evaluating `expm`

on `D2(0.1,10)`

to confirm your explanation.

In [ ]:

```
```

**Question 1(c)** Why can't we solve the equation with `ode45`

?

We will now implement the forward Euler method. However, we need to modify it to incorporate boundary conditions.

**Question 2(a)** Why can't you implement the forward Euler scheme

**Question 2(b)** Construct a forward Euler scheme that only modifies the interior points, that is the `2:end-1`

entries of the vector, for $u_0(x) = x(1-x)\exp(x)$. Plot the solution at time $t=1$, where $n=99$, $\Delta x = 1/(n+1)$ and $\Delta t = \Delta x^2/2$.

In [ ]:

```
using PyPlot
n=99
Δx=1/(n+1)
Δt=Δx^2/2
x=linspace(0.,1.,n+1) # Grid points
K=20000 # number of time steps to get to T=1
u0=x->x.*(1-x).*exp(x) # initial condition
w=Array(Vector{Float64},K+1)
w[1]=u0(x)
##TODO: Fill w with each time step from Forward Euler
plot(w[K+1])
```

**Question 2(c)** What goes wrong if $\Delta t = \Delta x^2$?

In [ ]:

```
```

We will now show that backward Euler easily allows for boundary conditions. But we need to put in some thought into how boundary conditions are incorporated:

**Question 3(a)** Why does the following backward Euler construction break down:
$$\mathbf w_{k+1} = (I-Δt \tilde D_n^2)^{-1}\mathbf w_k$$

**Question 3(b)** How would you incorporate the boundary conditions in the linear solve? Repeat **2(b)** with the backward Euler method.

In [ ]:

```
```

**Question 3(c)** Do you see the same issues as in forward Euler in **2(c)** when using $\Delta t =\Delta x^2$ in backward Euler?

In [ ]:

```
```