# Backward Euler time step for finite elements

For the backward Euler discretization in time: $$\left( \frac{u^{(k)}-u^{(k-1)}}{\Delta t}, v\right) + a(u^{(k)},v) = \ell(v)$$

where $a(\cdot,\cdot)$ is the bilinear operator associated with the discretization in space. I’m considering linear elliptic problems in 2D. So I rearrange the above equation: $$\left(\frac{1}{\Delta t}u^{(k)},v\right) + a(u^{(k)},v) = \ell(v) + \left(\frac{1}{\Delta t}u^{(k-1)},v\right)$$ since $u^{(k-1)}$ is a known quantity from the previous time step. For non-time dependent problems, the convergence in space is typically on the order of $\mathcal{O}(h^{p+1})$, with $p$ being the basis degree. Now if I have a time-dependent problem as given above, how small does my time step have to be in order to test for convergence in space? If I take a large time step, the method is stable but yields error that is large, on the order of $10^{3}$.

Taking my time step extremely small yields proper convergence rates, but I don’t think this tells me anything since for small $\Delta t$, the time-dependent terms will dominate and I’m essentially solving $$\left(\frac{1}{\Delta t}u^{(k)},v\right) = \left(\frac{1}{\Delta t}u^{(k-1)},v\right)$$ and if the time step is small, we can most likely expect $u^{(k)} \approx u^{(k-1)}$, so solving the system will yield the “correct solution”.

I’m just wondering if my implementation error of Backwards euler is incorrect. I know typically we solve the linear system: $$(\mathbf{M} + \Delta t\mathbf{K})u^{(k)} = \mathbf{M}u^{(k-1)} + \Delta t\mathbf{f}$$

with $\mathbf{M}$ and $\mathbf{K}$ being the mass and stiffness matrices and $\mathbf{f}$ being the source term. But is the implementation I gave above incorrect? I originally implemented the first way I listed, then tried the second for comparison but they yield the same result in both cases, so I'm guessing that's not the issue.

Remember that the backward Euler is first-order in time so your error will look something like: \begin{equation*} \|u-u_{h,\Delta t}\| = c_1 h^{p+1} + c_2 \Delta t + \cdots \end{equation*} If you want to test for spatial accuracy then $c_2 \Delta t$ needs to be an order smaller than $c_1 h^{p+1}$ on the finest grid and vice versa. If you obtain the theoretical orders for a manufactured solution it confirms that the implementation is correct.
• I feel a bit silly for asking, but on a 2D uniform mesh of regular triangles, $h$ would be the shortest edge length, correct? I’m a bit confused how to translate the theory to implementation. Let $h_{\text{min}}$ be the minimum edge length on the finest grid. Then on each grid, I’m taking $\Delta t = h_{\text{min}}^{p+1}$ with 10 time steps and not observing the theoretical orders. By an order smaller, should I just divide by 10? Using $\Delta t = h_{\text{min}}^{p+1}$ for manufactured polynomial solutions, I get close to the correct convergence rates, but when I should have zero error – Justin Dong Nov 29 '13 at 8:40
• (ex: $u(x,y,t)=x^{3}+y^{3}+t$) I hit an accuracy barrier and my errors stagnate around $10^{-7}$. – Justin Dong Nov 29 '13 at 8:41
• Dividing by 10 seems to yield better/correct convergence rates but the accuracy barrier around $10^{-7} - 10^{-9}$ depending on $p$ is still observed (my code doesn’t exhibit this barrier in the steady state case). – Justin Dong Nov 29 '13 at 8:46
• You forget about the presence of the factors $c_1, c_2$. I don't know what values they have, but they definitely have physical units, so you can't just assume that they can be set to one. – Wolfgang Bangerth Dec 1 '13 at 23:40