How to measure the error of Finite Element approximation in satisfying the PDE?

Question

In Galerkin methods, we seldom can measure the accuracy of an approximation by tracking the value of the residual. For example, take the wave equation:

$ u_{tt} = u_{xx}$, $(x,t) \in (0,L) \times (0,T)$ with approximation $u(x,t) \approx u_N(x,t) = \sum_{i=1}^N\phi_i(x)q_i(t)$

The residual for this approximation is defined as follows:

$R(x,t,N) = \sum_{i=1}^N\phi_i''(x)q_i(t) - \phi_i(x)q_i''(t)$

for which we expect $R(x,t) = 0$ everywhere for $(x,t) \in (0,L) \times (0,T)$. The error metric in this case will be

$ \epsilon(N) = \int_0^T \int_0^L R(x,t,N)^2 dx dt $

However, in the Finite Element method, developing the above integral, the term $\phi_i''(x)$ persists. This term is not well defined in the FEM as we have $\phi_i(x)$ continuous and differentiable once. Is there a way to estimate the residual without necessitating $\phi_i''(x)$?

My Work

I understand that in the FEM, we are mostly concerned with the error metric $||u-u_N||$ rather than the error in residual. The other solution will be to compute $u_\infty$ (ie, $N$ very large, assuming convergence occurs) and use the error metric $||u_\infty-u_N||$ instead. However, in my case, as I increase $N$, I do not necessarily keep my initial conditions constant. My understanding is the $u_\infty$ is computed for given initial and boundary conditions while, for myself, only the boundary conditions are given and the initial conditions are found such that the solution is periodic $u_N(x,t) = u(x,t+T)$. So the initial conditions are not given per-se. I would still prefer if there existed a way to compute the residual.

Answer 1

The residual error estimater can help in your case.

Let $u_h\in V_h$ be the FE solution given at a time $t$ and assume you know $\frac{\partial^2 u_h}{\partial t^2}$ at the same time $t$ then you can define the residual at an element $T$ as:

\begin{align} \Delta u_h|_T - \frac{\partial^2 u_h}{\partial t^2}|_T \end{align}

As you mentioned $\Delta u_h$ is not well defined since $\nabla u_h$ can have jumps at the edges of our elements in the mesh. But what we can still do is measure the violation of the continuity of $\nabla u_h$ at an edge $E=T_1\cap T_2$ by the normal-jump:

\begin{align} \left[\left|\frac{\partial u_h}{\partial n} \right|\right]:= \frac{\partial u_h}{\partial n_1}|_{T_1} + \frac{\partial u_h}{\partial n_2}|_{T_2} \end{align}

By that we obtain the residual error estimator as: \begin{align} \eta^{res}(u_h,\frac{\partial^2 u_h}{\partial t^2})^2 := \sum_{T\in\mathcal{T}}\left[ h_T^2 ||\Delta u_h-\frac{\partial^2 u_h}{\partial t^2}||_{L_2(T)}^2 + \sum_{E:\text{ edge in } \partial T} h_E || \left[\left|\frac{\partial u_h}{\partial n} \right|\right] ||_{L_2(E)}^2\right] \end{align} where $h_T$ is the maximum distance between two points in element $T$ (longest edge for triangles) and $h_E$ respectively the longest distance between two points on the edge $E$ (edge lenght for triangles). To see that the residual error estimator is a reasonable choice to recover the unknown error one can prove the following:

\begin{align} c||u-u_h||_{V(\Omega)} \leq \eta^{res}(u_h,\frac{\partial^2 u_h}{\partial t^2}) \leq C||u-u_h||_{V(\Omega)} \end{align} where $c,C>0$ are positive constants and $u \in V(\Omega)$ is the exact solution of $\Delta u=\frac{\partial^2 u_h}{\partial t^2}$ at time $t$ in the Domain $\Omega$. We usually call an error estimator reliable if the left inequality holds and efficient if the right inequality holds. Since the residual error estimator is efficient and reliable, you have some guarantee, that it is at most a constant factor apart from your error $||u_h-u||_{V(\Omega)}$ and if your error blows up or goes to zero it will do the same.