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:
\begin{equation} \frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2} \end{equation}
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:
\begin{equation} \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 \end{equation}
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:
\begin{equation} \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 \end{equation}
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]$):
\begin{equation} u_1^1 = -\frac{\left(\frac{1}{6} - \alpha \theta \right) \hat{u}_0}{\frac{2}{3} + 2\theta \alpha} \end{equation}
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:
- The boundary value at node $0$ differs from the initial solution at node $0$.
- 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:
- Is this a known problem? I guess it is, but I've had trouble finding info on it.
- 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$:
\begin{equation} \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}} \end{equation}
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?
