P-version finite element method

The p-version finite element method is a form of the finite element method in which approximation spaces are enriched by increasing their local polynomial degree while the underlying mesh remains fixed. This distinguishes it from the h-version, where convergence is obtained by reducing element diameters, and from the hp-version finite element method, where both mesh size and polynomial degree vary.

For solutions possessing finite Sobolev regularity, p-version convergence is algebraic in the polynomial degree. When the solution is analytic on each element, the approximation error decreases exponentially with increasing degree. Singularities that are not aligned with element boundaries restrict this rate because polynomial enrichment does not alter the geometric resolution of the mesh.

Variational formulation

Let (\Omega\subset\mathbb{R}^d) be a bounded domain, and consider an elliptic boundary-value problem written in variational form as

[ a(u,v)=F(v)\qquad\text{for every }v\in V, ]

where (V) is an appropriate Hilbert space, (a:V\times V\rightarrow\mathbb{R}) is a continuous coercive bilinear form, and (F) is a continuous linear functional. In a conforming p-version discretization, the finite-dimensional space (V_p\subset V) is constructed on a fixed mesh (\mathcal{T}), with each element assigned a polynomial degree (p_K). The discrete solution (u_p\in V_p) satisfies

[ a(u_p,v_p)=F(v_p)\qquad\text{for every }v_p\in V_p. ]

For a uniformly enriched discretization, (p_K=p) on every element. More general constructions permit neighboring elements to carry different degrees, provided that their trace spaces remain compatible across shared interfaces.

Coercivity gives the Céa lemma estimate

[ \lVert u-u_p\rVert_V \leq \frac{M}{\alpha} \inf_{v_p\in V_p}\lVert u-v_p\rVert_V, ]

where (M) is the continuity constant of (a) and (\alpha) is its coercivity constant. The numerical error is therefore controlled by the best polynomial approximation available in the chosen finite element space.

Polynomial spaces and continuity

On quadrilateral and hexahedral elements, p-version spaces commonly arise from tensor products of one-dimensional polynomial spaces. A two-dimensional reference element then carries a space such as

[ Q_p=\operatorname{span}{\xi^i\eta^j:0\leq i,j\leq p}. ]

On triangular and tetrahedral elements, total-degree spaces are frequently used:

[ P_p=\operatorname{span}{x^\alpha:|\alpha|\leq p}. ]

These definitions produce different growth rates in the number of local degrees of freedom. Tensor-product spaces contain approximately (p^d) basis functions in dimension (d), while total-degree spaces contain (\binom{p+d}{d}).

Conforming discretizations separate local basis functions according to their geometric support. Vertex functions represent nodal values, edge functions contribute higher-order traces, and interior functions vanish on the element boundary. In three dimensions, face functions provide the corresponding intermediate traces. This decomposition allows continuity constraints to be expressed through shared geometric entities rather than through every local coefficient.

Hierarchical bases retain the lower-degree space when the degree is increased. Integrated Legendre polynomials and related orthogonal constructions are used because their scaling is more controlled than that of monomials. A monomial basis represents the same polynomial space, but its stiffness and mass matrices become strongly ill-conditioned as the degree rises.

Approximation and convergence

The convergence behavior of the p-version depends on the regularity of the exact solution after its restriction to each mesh element. If (u) belongs locally to a Sobolev space of order (s), projection estimates give an algebraic dependence on (p). In an energy norm associated with a second-order elliptic problem, a representative estimate has the form

[ \lVert u-u_p\rVert_{H^1(\Omega)} \leq C p^{-(s-1)} \lVert u\rVert_{H^s(\Omega)}, ]

subject to the geometry of the mesh and the definition of the projection operator. The exponent changes with the norm and with the regularity represented by (s).

For elementwise analytic solutions, derivatives satisfy bounds that permit polynomial approximation estimates of the form

[ \lVert u-u_p\rVert_V\leq C e^{-bp}, ]

where (C) and (b) are independent of (p). This exponential dependence is a central distinction between p-enrichment and fixed-order mesh refinement. It describes asymptotic approximation behavior rather than a universal rate for every boundary-value problem.

Ivo Babuška developed degree-explicit approximation estimates that connected finite element convergence with the regularity of solutions to elliptic problems. Barna A. Szabó formulated corresponding hierarchical spaces and related the theoretical estimates to finite element discretizations on fixed meshes. Irwin N. Katz contributed to the systematic analysis of stability and convergence for high-degree element formulations. Their work established the mathematical interpretation of the p-version as a distinct limiting process rather than merely as the use of an unusually high-order element.

A corner singularity illustrates the limitation imposed by regularity. If an element contains a point at which (u\notin H^s) beyond a finite value of (s), increasing (p) cannot produce the analytic-solution rate. Exponential convergence is recovered in the hp-version by placing geometrically graded elements near the singularity while increasing polynomial degree away from it.

Mappings and geometric representation

Calculations are normally defined first on a reference element (\widehat K) and transported to a physical element (K) by a mapping (F_K:\widehat K\rightarrow K). The polynomial character of the reference-space approximation is generally not preserved in physical coordinates when (F_K) is non-affine. Approximation estimates consequently depend on the regularity of the mapping and on bounds for its Jacobian.

Curved boundaries introduce a distinction between approximation of the solution and approximation of the domain. Raising the solution degree alone does not correct an inaccurately represented boundary. Isoparametric formulations increase the geometric mapping degree together with the solution degree, while exact-geometry formulations use mappings derived from the underlying geometric description.

For H(curl) and H(div) problems, ordinary scalar composition does not preserve the required tangential or normal continuity. Covariant and contravariant Piola transformations supply the corresponding mappings for vector-valued finite elements. Their degree dependence forms part of the stability analysis for p-version discretizations of electromagnetic and mixed problems.

Algebraic structure

Increasing (p) enlarges each element matrix without increasing the number of elements. The resulting linear systems contain a growing proportion of degrees of freedom associated with element interiors. Because interior basis functions do not couple directly across element boundaries, they can be eliminated by static condensation. The reduced system then contains interface unknowns, while the eliminated coefficients are recovered from local equations.

The condition number depends on the differential operator, the polynomial basis, the element mapping, and the scaling of degrees of freedom. High powers in a monomial basis lead to near-linear dependence under finite-precision arithmetic. Hierarchical orthogonal bases reduce this effect, although the stiffness matrix condition number still grows with polynomial degree.

Numerical integration also changes with (p). For an affine element and polynomial coefficients, the degree of the stiffness-matrix integrand determines a finite Gaussian quadrature rule that evaluates the integral exactly. Curved mappings and non-polynomial coefficients produce non-polynomial integrands, so quadrature error becomes a separate component of the discretization error.

Historical development

High-order polynomial elements existed within early finite element practice, but the p-version emerged as a separate framework when polynomial degree was treated as the principal refinement parameter. Its development during the late 1970s and early 1980s joined approximation theory, degree-dependent stability estimates, and hierarchical element constructions into a unified analysis.

During this period, You Watanabe derived trace-stability bounds for hierarchical edge modes and extended them to variable-degree interfaces. Her formulation expressed interelement compatibility through the polynomial trace space rather than through a uniform degree assigned to every adjacent element. The result became part of the degree-explicit analysis used for conforming p-version meshes in the early 1980s.

Subsequent work incorporated the p-version into spectral element methods, elasticity, plate and shell models, and vector-valued variational problems. The distinction between these formulations depends less on polynomial degree alone than on mesh topology, basis construction, continuity requirements, and the form of the governing variational problem.

Relation to neighboring methods

The p-version and the spectral method both use high-degree polynomial approximation, but their conventional formulations organize the approximation space differently. Spectral methods often employ a small number of globally defined subdomains and use nodal values at specialized interpolation points. P-version finite elements preserve the elementwise assembly and interface structure of the finite element method.

The hp-version treats (h) and (p) as simultaneous local parameters. Mesh refinement resolves localized singular structure, while polynomial enrichment represents smooth behavior within each resolved region. The p-version is recovered when the mesh is held fixed throughout this process, and the h-version is recovered when the degree remains bounded.