# FEM diffusion: inaccurate results small time steps

I wrote some FEM code and found some strange results when using a very small time step. So, I decided to analyze the discrete equations.

Consider the following linear diffusion problem in 1 dimension:

$$$$\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2}$$$$

where k is a positive constant.

I will solve this problem on a 1D mesh of 2 linear elements, that is: 3 nodes located at $$x_0$$, $$x_1$$ and $$x_2$$ where $$x_{i+1} = x_i + h$$ and $$h$$ is the spacing (a positive constant). The first and last node constitute the Dirichlet boundary (so, there is no Neumann boundary). The initial and boundary conditions read:

\begin{align} u(t = 0, x) &= u_i \\\ u(t, x_0) &= \hat{u}_0 \\\ u(t, x_2) &= \hat{u}_2 \end{align}

Using the classic Galerkin method (with shape functions $$(1 - \xi, \xi)$$), discretization of the spatial domain yields:

$$$$\frac{h}{6} \begin{pmatrix} 2 & 1 & 0 \\\ 1 & 4 & 1 \\\ 0 & 1 & 2 \end{pmatrix}\dot{\bf{u}} + \frac{k}{h} \begin{pmatrix} 1 & -1 & 0 \\\ -1 & 2 & -1 \\\ 0 & -1 & 1 \end{pmatrix} \bf{u} = 0$$$$

where $$\bf{u}$$ is the nodal vector $$(u_0, u_1, u_2)^T$$ and $$\dot{\bf{u}}$$ is the derivative of $$\bf{u}$$ with respect to time.

Using the $$\theta$$-method for the discretization of the time domain, we get:

$$$$\left[ \frac{h}{6 \Delta t} \begin{pmatrix} 2 & 1 & 0 \\\ 1 & 4 & 1 \\\ 0 & 1 & 2 \end{pmatrix} + \theta\frac{k}{h} \begin{pmatrix} 1 & -1 & 0 \\\ -1 & 2 & -1 \\\ 0 & -1 & 1 \end{pmatrix} \right] {\bf{u}^{n+1}}\\\ = \\\ \left[ \frac{h}{6 \Delta t} \begin{pmatrix} 2 & 1 & 0 \\\ 1 & 4 & 1 \\\ 0 & 1 & 2 \end{pmatrix} - (1 - \theta)\frac{k}{h} \begin{pmatrix} 1 & -1 & 0 \\\ -1 & 2 & -1 \\\ 0 & -1 & 1 \end{pmatrix} \right] \bf{u}^n$$$$

where $$\Delta t$$ is the time step and super-indexes refer to the time, such that $$t_{n+1} = t_n + \Delta t$$.

Finally, defining $$\alpha = \frac{k \Delta t}{h^2}$$ and imposing the boundary conditions, we get:

\begin{align} u_0^{n+1} &= \hat{u}_0 \\\ u_1^{n+1} &= \frac{1}{\frac{2}{3} + 2\theta \alpha} \left[ \left(\frac{1}{6} - \alpha \theta \right)\left(u_0^n - \hat{u}_0\right) + \left(\frac{1}{6} - \alpha \theta \right)\left(u_2^n - \hat{u}_2\right) + \\\ \alpha \left( u_0^n + u_2^n \right) + \left(\frac{2}{3} - 2\alpha (1 - \theta)\right)u_1^n \right]\\\ u_2^{n+1} &= \hat{u}_2 \end{align}

where sub-indexes refer to the nodes.

For the sake of simplicity, let's assume that:

\begin{align} u_i &= 0 \\\ \hat{u}_0 &> 0 \\\ \hat{u}_2 &= 0 \end{align}

Now, the solution at time $$t_{n+1}$$ is given by:

\begin{align} u_0^{n+1} &= \hat{u}_0 \\\ u_1^{n+1} &= \frac{1}{\frac{2}{3} + 2\theta \alpha} \left[ \left(\frac{1}{6} - \alpha \theta \right)\left(u_0^n - \hat{u}_0\right) + \alpha u_0^n + \left(\frac{2}{3} - 2\alpha (1 - \theta)\right)u_1^n \right]\\\ u_2^{n+1} &= \hat{u}_2 \end{align}

Now, under those conditions, we know that the value of $$u_1$$, at any time, lies in the interval $$\left[ 0, \frac{\hat{u}_0}{2} \right]$$, and in particular, when $$\Delta t$$ tends to $$\infty$$, it is easy to show that $$u_1$$ tends to $$\frac{\hat{u}_0}{2}$$.

However, we can see that the solution $$u^{n+1}_1$$ can take negative values if $$\alpha$$ is small enough (which, for instance, can be achieved by choosing a small enough time step $$\Delta t$$). This can be easily seen when calculating the solution at node 1 at time $$t_1$$ (recall that the initial solution was chosen to be 0 $$\forall x \in [0,2]$$):

$$$$u_1^1 = -\frac{\left(\frac{1}{6} - \alpha \theta \right) \hat{u}_0}{\frac{2}{3} + 2\theta \alpha}$$$$

If we choose a time step such that $$\alpha < \frac{1}{6\theta}$$, then $$u^1_1$$ becomes negative, which is clearly an inaccurate result.

To summarize, this problem comes from the difference $$u_0^n - \hat{u}_0$$ and may arise in the following situations:

1. The boundary value at node $$0$$ differs from the initial solution at node $$0$$.
2. The boundary value at node $$0$$ at time $$t_n$$ differs from the boundary value at node $$0$$ at time $$t_{n+1}$$.

Both points are actually equivalent. So, here are 2 questions:

1. Is this a known problem? I guess it is, but I've had trouble finding info on it.
2. If it is not a known problem, where did I go wrong? Did I apply the boundary values wrongly?

This second question made me calculate the solution $$u^{n+1}_1$$ when $$\alpha$$ tends to $$0$$:

$$$$\lim_{\alpha \to 0} u^{n+1}_1 = \frac{\frac{1}{6}\left(u_0^n - \hat{u}_0\right) + \frac{2}{3}u^n_1}{\frac{2}{3}}$$$$

The limit should be equal to $$u^n_1$$, which is only true if $$u_0^n = \hat{u}_0$$. So, here goes another question:

When applying a boundary value on node $$0$$ at time $$t_{n+1}$$, do I need to set $$u_0^n$$ to that same value as well? This would eliminate the problem of having negative values, but am I still solving the same problem here?