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.


1 Answer 1


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.

  • $\begingroup$ 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 $\endgroup$ Commented Nov 29, 2013 at 8:40
  • $\begingroup$ (ex: $u(x,y,t)=x^{3}+y^{3}+t$) I hit an accuracy barrier and my errors stagnate around $10^{-7}$. $\endgroup$ Commented Nov 29, 2013 at 8:41
  • $\begingroup$ 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). $\endgroup$ Commented Nov 29, 2013 at 8:46
  • $\begingroup$ 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. $\endgroup$ Commented Dec 1, 2013 at 23:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.