How to refine $h$ and $\Delta t$ for convergence tests on evolution PDE

Question

Setting

I am solving for $u(x,y,t)$ the wave equation $\partial_{tt} u - \partial_{xx} u = f$ on $(x,y)\in\Omega=[0,1]\times \mathbb{R}$ by splitting it into an equivalent first order system: $$\partial_t \begin{pmatrix} u\\ \phi \end{pmatrix}+ \begin{pmatrix}0&b\cdot\nabla\\ b\cdot\nabla&0 \end{pmatrix}\begin{pmatrix} u\\ \phi \end{pmatrix}=0 ,\qquad b=\begin{pmatrix} 1\\ 0\end{pmatrix}$$

and applying a Crank-Nicolson method to this system (which is equivalent to the Newmark scheme for the original wave equation). I use FEM (fenicsx software library) with polynomials of various degrees on a fully-structured triangle-element mesh for the space discretization.

Expected convergence rates (from FEM theory):

The error $u-u_h$ should converge of rate 2 in the $L^2$ norm, when using polynomials of order 1 and the Crank-Nicolson scheme, i.e. $$\epsilon = \Vert u-u_h\Vert_{L^2(\Omega)} \leq C( h^2 + \Delta t^2)$$

When i instead use polynomials of order 2 i expect to get $$\epsilon = \Vert u-u_h\Vert_{L^2(\Omega)} \leq C( h^3 + \Delta t^2)$$

Now i would like to measure these rates in a numerical experiment. In particular, I would like to measure $h^3$, be refining $\Delta t$ faster (so that I dont have to switch my time-discretization scheme).

Computing the convergence rate $r$ using the numerical solution

Suppose we have a stationary problem, then (C independent of $h$ or $i$) $$\epsilon_i = Ch_i^r \implies \left(\frac{h_{i+1}}{h_i}\right)^r=\frac{\epsilon_{i+1}}{\epsilon_i} \implies r=\frac{ \log (\epsilon_{i+1}/\epsilon_{i})}{\log(h_{i+1}/h_i)}$$

However my wave problem has 2 discretization parameter: $h$ and $\Delta t$. So i cant find a formula for the rates $r_h$ or $r_t$ using the above approach. But i can refine $\Delta t$ so fast (that it "vanishes" asymptotically), that I actually observe $r_h=3$ and then the above formula for $r$ is correct.

Lets fix $h_{i+1}=\frac{1}{2}h_i \implies \epsilon_{i+1}=\frac{1}{8}\epsilon_i$.

Does this mean, the refinement factor for the time scheme to reveal this fast FEM convergence is $$\frac{\Delta t_{i+1}}{\Delta t_i}=\frac {1}{\sqrt 8}$$

so that we end up with $$\frac{\epsilon_{i+1}}{\epsilon_i} \leq \frac{C}{C}\frac{h_{i+1}^3 + \Delta t_{i+1}^2}{h_i^3 + \Delta t_i^2} =\frac{(h_i\frac{1}{2})^3 + (\Delta t_i \frac{1}{\sqrt 8})^2}{h_i^3 + \Delta t_i^2} = \frac{1}{8}\frac{h_i^3 + \Delta t_i^2}{h_i^3 + \Delta t_i^2}=\frac{1}{8}?$$

Is this derivation correct?

A numerical experiment reveals, that I can not observe the fast rate $r_h=3$ with the $\sqrt 8$ time-refinement.

I would like to understand if

• either my theory is fine and my code is wrong
• or my theory is wrong

Answer 1

Here are a few solutions that you could explore to determine the orders in space and time.

1) Separate study of time error

You can use a given spatial mesh, and perform multiple simulations with finer and finer time steps (take them typically uniformly distributed in log space, i.e. $\Delta t_{i+1} = k \Delta t_i, i=2..N$ with $k\in [0,1]$ and a chosen value for $\Delta t_0$. Take the solution with the smallest time step as the reference solution. You can compute the error as:

$$\epsilon_i = ||u_i(t_f) - u_N(t_f)||_{L^p((\Omega)}$$

with $p=1$ or 2, i.e. a norm of the error in space at the final time (only useful is you final solution is still in sufficiently unsteady situation).

$$\epsilon_i = ||u_i(t_f) - u_N(t_f)||_{L^p((\Omega \times \mathbb{R})}$$ i.e. a norm of the error in space and time across the whole transient. This latter form gives in my experience smoother curves, but the order observed order is the same. Its disadvantage is that it requires some kind of decimation (easy if $k=1/2$ for instance) or time interpolation to be computed.

Note that I am not sure that this process can be performed conversely for a fixed time mesh a varying space mesh refinement. Indeed the error constant in the temporel error term most likely will depend on the semi-discrete problem you seek to integrate, which will hence be modified as you refine the space mesh.aonversely, you could do the same thing with a fixed time mesh and a varying space mesh. I have never tried this personally, and most likely it would also fail with some time schemes, which might suffer from e.g. the increased stiffness of the problem as the mesh is refined.

2) Space study with negligible time error

For that study, you refine more and more the spatial mesh, and for every mesh, you integrate the corresponding semi-discrete system in time with a "sufficiently" small time step (to be manually determined), or at best with an adaptive scheme and tight integration tolerances so that you mimic an exact temporal integration. Note that in the case of a linear system, you can get an exact temporal integration with an exponential integrator.

Then, simply compute the error in space at the final time (or in space-time across the whole transient) with respect to the solution with the finest mesh. The observed order will be the order of the space scheme.

3) Combined space-time convergence study

Given that you know the expected accuracy of your space and time schemes, you know that the error asymptotically behaves as $\epsilon = O(\Delta t^p) + O(\Delta x^q)$, with $(p,q) \in \mathbb{N}^2$. If you gradually refine $\Delta x$ and, for each of its value, you ensure $\Delta t^p = a \Delta x^q$, with a chosen constant (use the same for all meshes), then both error terms behave as $(\Delta x^q)$. If you do not see the expected convergence with order $q$, then something is most likely wrong, either in the theoretical determination of the spatial order $q$ or that of the time order $p$, or with the implementation. Or you may also be in the non-asymptotic regime, where the error may behave differently than or vanishing space and time steps.

In your case, this would indeed lead to $\Delta t_{i+1} = \frac{1}{\sqrt{8}} \Delta t_i$, hence I would suggest trying both previous points to assess whether the error is coming from the space or time discretisation. Maybe also try computing the space-time error also, it may help recover the expected orders.

4) Method of manufactured solution

You can also look up this particular method, which involves adding a source term in your equations, determined so that the solution is driven to a chosen analytical profile. This method may be simple to implement and test.