Significance of p-convergence studies

Question

Consider a method (e.g., FEM) with variable approximation order $p$. Now, we know that the optimal order of convergence is given by $$e = C h^{p+1},$$ where $h$ denotes the mesh size and a constant $C$ that does not depend on $h$. As a result, we can study the experimental order convergence (or $h$-convergence) on a series of grids $h, \frac{h}{2}, \ldots$ because you can directly compare two errors $e_1$ and $e_2$ on different grids.

On the other hand, I've always been puzzled by the way people perform and evaluate $p$-convergence studies in this context. Here, we have a fixed grid with mesh size $h$, but different approximation orders $p=1,2,3,\ldots$. My problem now is that the constant $C$ does not depend on $h$, but it does depend on $p$ (otherwise, we would see that the lines in a log-log-plot of the $h$-convergence results for different values of $p$ would intersect in a single point - which they, in my experience, don't do). Now, what is the significance of a $p$-convergence study having realized that we cannot really compare two errors, $e_1$ and $e_2$, using the above formula for the error, given that they have been obtained using different polynomial degrees?

Answer 1

It isn't so much that we want to compare the $p$-refinement and $h$-refinement errors directly, instead we want to compare the convergence properties (e.g. speed) of each refinement strategy. This requires more knowledge of the constant in the apriori error estimate. We'll illustrate by looking at the apriori error estimate of a discontinuous Galerkin spectral method (which can be interpreted as a high-order FEM).

It is true that the constant $C$ in the inequality $$ e \leq Ch^{p+1}, $$ depends on the approximation order $p$. Thus, it is not immediately clear what to expect when we increase the order $p$. However, the error inequality above hides a lot of information about the convergence behaviour in that constant. In $p$-refinement it becomes important how the $C$ depends on $p$.

In a 1D spectral element method the error estimate on a given element $k$ has a form (derived in a similar fashion as in Canuto et. al.)$$ e \leq \mathcal{C}h^{\min(m,p+1)}p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)}, $$ where $\mathcal{C}$ is independent of $h$ and $p$, $q^k$ is the true solution on element $k$, and $H^m$ is the $m$ order Sobolev space. In we apply the numerical method to a problem with a "smooth" solution it will generally be true that $p+1\leq m$. Now we can more readily see the difference in $p$ versus $h$ refinement.

If we fix $p$ then the terms $$ p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)} $$ are a constant that we absorb into $\mathcal{C}$ to obtain the $C$ in the first error estimate. Now we can perform a standard $h$-refinement argument and see that the convergence will be similar to that of a finite difference scheme, i.e., some polynomial order on $h$.

However, if we fix the grid $h$ then the term $h^{\min(m,p+1)}$ is constant and the error estimate looks like $$ e \leq \hat{\mathcal{C}}p^{\frac{3}{2}-m}\|q^k\|_{H^m(\Omega)},$$ where $\hat{\mathcal{C}} = \mathcal{C}h^{\min(m,p+1)}$. Now we immediately see that the rate of convergence changes in the $p$-refinement study. Provided $m>\frac{3}{2}$ the error decays as a power of $p$. For "very smooth" solutions, e.g. infinitely differentiable, the decay in the error is faster than any polynomial order. Sometimes this is referred to as exponential convergence because if you look at a semilog plot of the error against $p$ you see a straight line.

As a quick aside we note that in spectral methods the convergence rate depends on the smoothness of the solution. Heuristically, we can interpret this as "the larger $m$ is in the error estimate, the smoother the solution," because the true solution is in a higher order Sobolev space, i.e., the true solution is square integrable up to the $m^{\textrm{th}}$ order derivative. The convergence of a Sobolev norm is often viewed in terms of the decay of Fourier coefficients of the solution $q^k$ (particularly for fractional derivatives). This is one reason the family of methods are referred to as spectral methods.

Answer 2

You can conduct your $p$ convergence experiments for different $h$. Since the "constant" $C$ depends on $p$, but not on $h$, this allows you to separate the effects of the $h^{p+1}$ and the $C(p)$.

The estimate of the form $C(p) h^{p+1}$ is justified by the fact that the $C(p)$ is typically polynomial in $p$, and thus for small $h$ and large $p$, the estimate is dominated by the exponential term.

Answer 3

This might be answering a different question that what you asked, but I think a chart like that is more about calling out to a reader whether p-adaption or h-adaption is more effective for converging the underlying model. Although the x-axis is often taken to be mesh granularity (h) or number of unknowns, if you normalize it to something that quantifies "work" (e.g. RAM consumption, flops required, elapsed wallclock), then you could gather at a glance whether p-refinement or h-refinement is the better strategy to hit a given error bound with minimum "work". Certainly this is problem dependent (in particular, whether or not the underlying fields are smooth or singular). But a couple of graphs like that ("work" vs error) could immediately demonstrate e.g. how coarse-h high-p is a better strategy on smooth problems, but fine-h low-p is better on say a singularity-dominated re-entrant corner. (Granted, on real applications that mix those two regimes, the answer itself is mixed - hence the motivation towards adaptive refinements)