# 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?

What you find is indeed correct. It is known that positivity is lost if very small time steps are chosen, see

https://doi.org/10.1515/cmam-2015-0018

This loss of positivity happens even for semi-discrete scheme.

The analysis for 1-d case is given in section 6 of this paper.

• thanks for the useful link
– mfnx
Mar 14, 2020 at 8:51
• @mfnx I did not try to give a detailed answer since it is not a simple answer, and the paper explains it well. If you think your question is answered by the analysis in the paper, you can close this by accepting my answer. Mar 15, 2020 at 4:53
• part of it is answered. I'd like to have an answer to the last questions
– mfnx
Mar 15, 2020 at 5:41