# Measure the convergence rate of a discretization of a wave equation

I'm currently trying to approximate the following type of wave equation (in weak formulation):

Let $\Omega \subset \mathbb{R}^d$ ($d=2$) be some polygonal domain. We seek a function $u \in L^2\left(0, T; H^1_0(\Omega)\right)$ with $u' \in L^2\left(0, T; L^2(\Omega)\right)$, $u'' \in L^2\left(0, T; H^{-1}(\Omega)\right)$ and

$$$$\langle u''(t), v \rangle_{H^{-1}(\Omega) \times H^1_0(\Omega)} + \langle \nabla u(t), \nabla v\rangle_{L^2(\Omega)} = \langle f(t), v\rangle_{L^2(\Omega)},$$$$ for almost every $t \in (0,T)$ and all $v \in H^1_0(\Omega)$ and with initial conditions $$$$u(0) = u_0, \ \ u'(0)=v_0,$$$$ with data $f \in L^2\left(0, T; L^2(\Omega)\right)$, $u_0 \in H^1_0(\Omega)$ and $v_0 \in L^2(\Omega)$.

Discretization: First, we just discretize the problem in space by using a finite element space $V_h$ of continuous, piecewise-linear polynomials and then use an (explicit) centred difference quotient of second order for the time derivative "$u''(t)$" with time step size $\Delta t = \frac{T}{M}$, $M \in \mathbb{N}$, and end up with a discrete scheme:

$$\label{leapfrogscheme} \left\langle \frac{\left( u_h^{n-1} - 2 u_h^{n} + u_h^{n+1} \right)}{(\Delta t)^2}, v_h \right\rangle_{L^2} + \langle \nabla u_h^{n}, \nabla v_h \rangle_{L^2} = \langle f^n, v_h \rangle_{L2},$$

for all $v_h \in V_h$ and with given start values $u_h^{0}$ and $u_h^{1}$, where the approximation of the Galerkin solution $u_h$ at time $t_n = n \Delta t$ is denoted by $u_h^n$ and we define $f^n := f(t_n)$.

As start values, we take for instance the $L^2$-projection (for all $v_h \in V_h$): \begin{align} \langle u_h^0, v_h \rangle_{L^2} &:= \langle u_0 , v_h \rangle_{L^2}, \\ \langle u_h^1, v_h \rangle_{L^2} &:= \left\langle u_0 + (\Delta t) v_0 + \frac{(\Delta t)^2}{2}\left( f^{0} + \Delta u_h^0 \right), v_h \right\rangle_{L^2}, \end{align} where we used the second order Taylor expansion \begin{align} u_h(\Delta t) \approx u_h(0) + (\Delta t) u_h'(0) + \frac{(\Delta t)^2}{2} u_h''(0) \end{align} in time $t=0$ to define $u_h^{1}$.

Error analysis: Here comes the part i'm having some trouble with, because i haven't found any good literature that covers this topic.
I basically want to figure out, what kind of convergence rate in a norm i can optimally expect (given the solution is smooth enough) and which discrete norm i should take to verify that my computed solutions resolve this rate.

In the paper Garth A. Baker - "Error Estimates for Finite Element Methods for Second Order Hyperbolic Equations", i found that, if additionally $u \in L^\infty\left(0, T; H^2(\Omega)\right)$, $u' \in L^2\left(0, T; H^2(\Omega)\right)$ and the third and fourth time derivative of $u$ exists and is in $L^2\left(0, T; H^2(\Omega)\right)$, then

$$\max_{n} \| u(t_n) - u_h^{n} \|_{L^2(\Omega)} \leq C \, ( h^2 + \Delta t^2 ),$$ where $h>0$ is the spatial mesh size of an uniform mesh that defines the finite element space $V_h$, i.e. all elements $T$ are of size $\mathrm{diam}(T) \approx h$.

For stability reasons, one has to take a sufficient small time step length $\Delta t \approx h_{\min}$ (CFL-condition). So, if we replace $\Delta t$ by $h$ and rewrite (in dimension d=2) $h = N^{-\frac{1}{2}}$, where $N$ denotes the degrees of freedom (here: number of inner nodes of the finite element mesh), then one should hopefully see an optimal error of $$\max_{n} \| u(t_n) - u_h^{n} \|_{L^2(\Omega)} \leq C \, N^{-1}.$$ So, this is basically an error analysis in the $L^\infty\left(0, T; L^2(\Omega)\right)$-norm.

Is this - so far - how people would measure the convergence of such a hyperbolic problem?

But now it gets a little bit more complicated by considering domains with reentrant corners, so that we can only expect the solution to be in $u(t) \in H^{1+\delta}(\Omega)$, for some $0 < \delta < 1$. This would damp the $L^2$-error of $u(t) - u_h(t)$ for a fix $t \in [0,T]$ down to $\mathcal{O}\left(h^{1+\delta}\right)$ instead of $\mathcal{O}\left(h^{2}\right)$ (i.e. down to $\mathcal{O}(N^{-\frac{1+\delta}{2}}$) if we would use an uniform mesh of size $h$).

To get around this drawback, i would like to use local refined meshes to still get the optimal $L^2$-error of $\mathcal{O}\left(N^{-1}\right)$ for a fix $t$. The problem now is that i still have to ensure stability by taking $\Delta t \approx h_{\min}$ - this depends on the smallest element diameter, which decreases drastically for locally refined meshes. So, the computational work we put into the time discretization isn't really comparable anymore with our number of interior nodes.

How would someone measure the convergence for this kind of discretization? I would also be grateful for some literature in general that considers the error analysis of the wave equation by using finite elements.

## 1 Answer

You are asking very complicated questions for which there are likely no answers that can be rigorously proven.

If you go back for a second and ask the same question for the solution of the Laplace equation, then you are asking whether, given that $u\in H^{1+\delta}$, you can find an algorithm that constructs a mesh with $N$ cells so that $\|u-u_h\|_{L_2} = {\cal O}(N^{-\frac{1+\delta}{2}})$ in 2d. I don't know whether that's actually true for the $L_2$ norm, but we know since the work of Binev, Dahmen, DeVore that there are algorithms to construct meshes so that $\|u-u_h\|_{H^1} = {\cal O}(N^{-1})$ if the solution is in $H^2$. That's already a very deep result that was only proven in the 2000s, if I recall correctly. Definitely not trivial.

Conversely, it will not be very difficult to show experimentally that a mesh created by a less sophisticated algorithm will not satisfy this guarantee. As a consequence, the concrete choice of the adaptive mesh matters -- a lot!

Now you go further and ask whether this all can be translated into similar results for the wave equation. I bet that this is not possible with today's mathematical strategies. The solution of the wave equation is less regular, because the wave equation is hyperbolic. So there is less to work with in order to prove the embeddings of the solution spaces of the wave equation into the approximation spaces used in the theorem mentioned above. You then also want to couple the mesh size and the time step size, which further complicates things. I suspect that nobody in the world today can prove results for this. Or, at the very least, the half dozen people who could do so, already know all of the things I mentioned above.