Sharp superconvergence rate at Gauss-Lobatto points of the $C^0$-$Q^k$ finite element method for elliptic problems with Neumann boundary conditions and mixed second-order derivatives

For second-order elliptic boundary value problems on rectangular domains, the continuous ($C^0$-$Q^k$) finite element method exhibits function-value superconvergence of order $k+2$ at the tensor-product Gauss-Lobatto points (the quadrature nodes, which include the cell vertices) when the boundary conditions are of Dirichlet type; this has been proven for the exact bilinear form, for the interpolated-coefficient implementation, and for the $(k+1)$-point Gauss-Lobatto quadrature implementation. For Neumann-type boundary conditions, however, when the operator contains mixed second-order derivatives (a non-diagonal coefficient matrix), every existing proof stops at order $k+3/2$ in the discrete 2-norm over the Gauss-Lobatto points — a barrier that persists even for the exact finite element scheme without quadrature or coefficient interpolation. Numerical evidence is conflicting: elliptic tests on Cartesian grids show approximately $(k+2)$-th order accuracy, while spectral-element wave tests on smoothly perturbed grids exhibit roughly $k+5/3$. The problem asks for the sharp rate: prove $O(h^{k+2})$ under Neumann (or mixed Dirichlet-Neumann) boundary conditions with mixed second-order derivatives, or establish rigorously that the rate cannot reach $O(h^{k+2})$.

Background

The finite element method (FEM) approximates the solution of an elliptic boundary value problem by a piecewise-polynomial function. On a rectangular mesh of a two-dimensional domain, the $C^0$-$Q^k$ element uses continuous functions that are tensor products of degree-$k$ polynomials (the space $V^h$); $k\geq 2$ throughout. The $(k+1)$-point Gauss-Lobatto quadrature rule on a reference interval is the rule whose nodes include the two interval endpoints; its tensor product on each cell gives a set of grid points $Z_0$ (all Gauss-Lobatto points of all cells) which, on a uniform mesh, coincides with a finite difference grid. On rectangular meshes the $C^0$-$Q^k$ element is identical to what the spectral element method literature calls the $Q^k$ spectral element.

Superconvergence is the phenomenon that the FEM error at special points decays faster than the global rate: while $\|u-u_h\|_{L^2}=O(h^{k+1})$ in general, the function values at Gauss-Lobatto points converge at order $k+2$ in the discrete 2-norm $\|v\|_{2,Z_0}=(h^2\sum_{z\in Z_0}|v(z)|^2)^{1/2}$, and at order $O(h^{k+2}|\ln h|)$ in the maximum norm over $Z_0$ — classical results built on one-dimensional knot superconvergence, M-type projection (a tensor-product Gauss-Lobatto interpolation whose approximation properties underpin the analysis), and point-line-plane interpolation techniques.

A 2020 paper by Hao Li and Xiangxiong Zhang showed that this $(k+2)$-th order superconvergence for Dirichlet boundary conditions survives two practical simplifications: (i) replacing the variable coefficients by their piecewise $Q^k$ Lagrange interpolants at the Gauss-Lobatto points, which yields a finite-difference-type scheme requiring only grid-point values of the data (for $k=2$, a fourth-order scheme with quadratic-coefficient assembly); and (ii) evaluating all integrals by $(k+1)$-point Gauss-Lobatto quadrature — the cheapest and most common implementation, because the Lagrangian $Q^k$ basis functions are delta functions at those points. Statement (ii) was proven in their companion paper, whose key new tool is applying the Bramble-Hilbert lemma (a standard tool asserting that a functional annihilating polynomials of degree $k$ is bounded by the $k+1$-st derivative of its argument) after integration by parts so that quadrature errors cancel between neighboring cells; standard cell-wise quadrature estimates are insufficient.

For Neumann-type boundary conditions the picture changes. The duality/superconvergence argument needs test functions in the full space $V^h$ (not the subspace $V_0^h$ of functions vanishing on the boundary), and for boundary cells the M-type projection estimate carries a half-order, trace-type loss: schematically, $\int_\Omega a\,(u-u_p)_x v_y\,dxdy = O(h^{k+2-1/2})$ for $v\in V^h$ versus $O(h^{k+2})$ for $v\in V_0^h$, where $u_p$ is the M-type projection of the exact solution. This loss occurs only when the second-order operator contains mixed derivatives, i.e. the coefficient matrix $(a_{ij})$ is not diagonal: for diagonal coefficients the offending cross-term is absent and $(k+2)$-th order is provable for Neumann problems as well. Consequently, for a general elliptic operator with mixed second-order derivatives and Neumann (or mixed Dirichlet-Neumann) boundary conditions, the best proven rate at all Gauss-Lobatto points is $O(h^{k+3/2})$ in the discrete 2-norm — and this half-order barrier holds even for the full FEM with exact integration and exact coefficients, before any quadrature or interpolation simplification. Numerical experiments do not settle the question: the elliptic tests of Li and Zhang (with $a_{12}=a_{21}\neq 0$ and purely Neumann conditions on a Cartesian grid) show approximately fourth-order accuracy for $k=2$, while spectral-element wave-equation tests on smoothly perturbed grids — which induce mixed second-order derivatives after mapping to the reference mesh — show a rate of roughly $k+5/3$, strictly between $k+3/2$ and $k+2$.

Problem Statement

Let $\Omega$ be a rectangular domain in $\mathbb{R}^2$ with a rectangular mesh of mesh size $h$, and consider the second-order elliptic boundary value problem

$$-\nabla\cdot(a\,\nabla u)+\sum_i b_i\,\partial_{x_i}u+cu=f\quad\text{in }\Omega,$$

where the coefficients are smooth in the general operator form — $a_{ij}, b_i, c \in W^{k+2,\infty}(\Omega)$ with the coefficient matrix $a=(a_{ij})$ uniformly positive definite, and with mixed second-order derivatives present (i.e. $a_{12}=a_{21}\neq 0$ on a set of positive measure) — together with homogeneous Neumann boundary conditions $(a\nabla u)\cdot n=0$ on all of $\partial\Omega$, or mixed homogeneous Dirichlet-Neumann conditions. For a purely Neumann problem, assume the usual compatibility condition on $f$ and normalize $u$ (and $u_h$) to have zero mean over $\Omega$. Let $V^h$ be the $C^0$-$Q^k$ finite element space with $k\geq 2$, let $Z_0$ be the set of all tensor-product $(k+1)$-point Gauss-Lobatto points of all mesh cells, and let $\|v\|_{2,Z_0}=(h^2\sum_{z\in Z_0}|v(z)|^2)^{1/2}$ be the discrete 2-norm. Let $u_h\in V^h$ be the finite element solution — for the exact bilinear form, and for the two standard practical implementations: all integrals computed by $(k+1)$-point Gauss-Lobatto quadrature, and coefficients replaced by their piecewise $Q^k$ Lagrange interpolants at the Gauss-Lobatto points. Assume the exact solution satisfies $u\in H^{k+2}(\Omega)$ together with the elliptic regularity that the standard superconvergence duality argument requires.

Determine the sharp convergence rate of the function values at the Gauss-Lobatto points under Neumann-type boundary conditions with mixed second-order derivatives. Specifically: is

$$\|u-u_h\|_{2,Z_0}=O(h^{k+2}),$$

with the corresponding $O(h^{k+2}|\ln h|)$ bound in the maximum norm over $Z_0$, as it is for Dirichlet boundary conditions — or does the discrete 2-norm rate provably fall short of $O(h^{k+2})$? That is, is the $(k+2)$-th order discrete 2-norm bound provable for the full finite element scheme with exact integration and exact coefficients, under Neumann-type boundary conditions (pure Neumann or mixed Dirichlet-Neumann), with the assumptions stated — the case that remains open, whose resolution determines the two practical implementations — or is there a construction of coefficient fields, right-hand sides $f$ and mesh sequences for which the discrete 2-norm error at $Z_0$ provably decays strictly slower than $h^{k+2}$, showing the sharp rate to be strictly below $k+2$ (whether it is $k+3/2$ or an intermediate rate such as the $\approx k+5/3$ suggested by some spectral-element tests)? The question is open even for the full finite element scheme with exact integration and exact coefficients; resolving that case resolves the barrier shared by all three implementations.

References

1. Hao Li and Xiangxiong Zhang, Superconvergence of $C^0$-$Q^k$ Finite Element Method for Elliptic Equations with Approximated Coefficients, Journal of Scientific Computing 82, Article 1 (2020). DOI: 10.1007/s10915-019-01102-1, https://doi.org/10.1007/s10915-019-01102-1 (preprint: arXiv:1902.00945, https://arxiv.org/abs/1902.00945)
2. Hao Li and Xiangxiong Zhang, Superconvergence of High Order Finite Difference Schemes Based on Variational Formulation for Elliptic Equations, Journal of Scientific Computing 82, Article 36 (2020). DOI: 10.1007/s10915-020-01144-w, https://doi.org/10.1007/s10915-020-01144-w (preprint: arXiv:1904.01179, https://arxiv.org/abs/1904.01179)
3. Hao Li, Daniel Appelö and Xiangxiong Zhang, Accuracy of Spectral Element Method for Wave, Parabolic, and Schrödinger Equations, SIAM Journal on Numerical Analysis 60 (2022) 339-363. DOI: 10.1137/21m1401760, https://doi.org/10.1137/21m1401760 (preprint: arXiv:2103.00400, https://arxiv.org/abs/2103.00400)
4. Hao Li and Xiangxiong Zhang, On the monotonicity and discrete maximum principle of the finite difference implementation of $C^0$-$Q^2$ finite element method, Numerische Mathematik 145 (2020) 437-472. DOI: 10.1007/s00211-020-01110-6, https://doi.org/10.1007/s00211-020-01110-6 (preprint: arXiv:1905.06422, https://arxiv.org/abs/1905.06422)
5. Logan J. Cross and Xiangxiong Zhang, On the Monotonicity of $Q^2$ Spectral Element Method for Laplacian on Quasi-Uniform Rectangular Meshes, Communications in Computational Physics 35 (2024) 160-180. DOI: 10.4208/cicp.oa-2023-0206, https://doi.org/10.4208/cicp.oa-2023-0206
6. Logan J. Cross and Xiangxiong Zhang, On the Monotonicity of $Q^3$ Spectral Element Method for Laplacian, Annals of Applied Mathematics 40 (2024) 161-190. DOI: 10.4208/aam.oa-2024-0007, https://doi.org/10.4208/aam.oa-2024-0007
7. Xiangxiong Zhang, Recent Progress on $Q^k$ Spectral Element Method: Accuracy, Monotonicity and Applications, in Spectral and High-Order Methods for Partial Differential Equations ICOSAHOM 2023, Lecture Notes in Computational Science and Engineering, Springer (2025) 113-138. DOI: 10.1007/978-3-031-76988-7_5, https://doi.org/10.1007/978-3-031-76988-7_5