Discrete Normal Derivative Analysis

• Discrete normal derivative is a boundary-sensitive operator that approximates the outward flux in discrete models, central to finite element and variational formulations.
• It enables finite element realizations where standard and variational approaches yield precise L2 error estimates, especially in polygonal domains affected by corner singularities.
• Additional constructions include directional stencils, spectral analogues on lattices, and fractal methods, broadening its applicability to nonlocal and graph-based settings.

Searching arXiv for papers on discrete normal derivatives and closely related formulations. A discrete normal derivative is a boundary-oriented discrete operator that approximates, represents, or generalizes the outward normal derivative in a discrete model. In the finite element analysis of elliptic PDEs, the term usually refers either to the classical boundary flux obtained from the discrete gradient on boundary elements or to a variational trace quantity defined by a discrete Green identity; for the Poisson problem on polygonal domains, both notions admit precise $L^2(\Gamma)$ error estimates and are strongly influenced by corner singularities and mesh grading (Pfefferer et al., 2018). In broader usage, closely related constructions occur as directional derivative stencils on grids, graph- and vertex-based fluxes on fractals, and nonlocal spectral analogues on lattices, so the term denotes a class of boundary-sensitive discrete operators rather than a unique universal formula (Cao et al., 2016).

1. Continuous antecedent and weak boundary-flux interpretation

For the model problem

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$

with $\Omega\subset \mathbb R^2$ polygonal and $f\in L^2(\Omega)$, the weak formulation is

$$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$

If $u\in H^2(\Omega)$, the outward normal derivative is the classical boundary flux

$$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$

For weak solutions on polygons, the normal derivative is understood through Green’s identity: $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$ In the weighted regularity regime used for polygonal domains, this identity yields $\partial_n u\in L^2(\Gamma)$ (Pfefferer et al., 2018).

The regularity of $\partial_n u$ is controlled by corner singularities. If the corners are $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$0 with interior angles $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$1, the singular exponents are

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$2

These exponents limit the global regularity of both $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$3 and $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$4. In weighted Sobolev spaces $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$5, one obtains precise corner-dependent regularity, and pointwise normal derivatives may vanish at convex corners and may blow up at non-convex ones (Pfefferer et al., 2018).

A complementary boundary-kernel viewpoint is furnished by Green functions. For a smooth bounded domain $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$6, with Green function $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$7 and Poisson kernel

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$8

the normal derivative of a function $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$9 satisfying $\Omega\subset \mathbb R^2$0 on $\Omega\subset \mathbb R^2$1 obeys

$\Omega\subset \mathbb R^2$2

This identifies the normal derivative as a boundary functional generated by an interior operator and a boundary kernel (Martin, 2012).

2. Finite element realizations

For a conforming triangulation $\Omega\subset \mathbb R^2$3, the standard linear finite element spaces are

$\Omega\subset \mathbb R^2$4

and the discrete state $\Omega\subset \mathbb R^2$5 solves

$\Omega\subset \mathbb R^2$6

Within this setting, two distinct discrete normal derivatives are used (Pfefferer et al., 2018).

Variant	Definition	Representation
Standard discrete normal derivative	$\Omega\subset \mathbb R^2$7 on a boundary edge $\Omega\subset \mathbb R^2$8 of element $\Omega\subset \mathbb R^2$9	elementwise constant on boundary edges
Variational discrete normal derivative	$f\in L^2(\Omega)$0 for all $f\in L^2(\Omega)$1	$f\in L^2(\Omega)$2

The standard derivative is the natural elementwise flux. Because $f\in L^2(\Omega)$3 is piecewise linear, $f\in L^2(\Omega)$4 is piecewise constant, so the boundary flux on each boundary edge is immediately available from the adjacent element.

The variational discrete normal derivative is the $f\in L^2(\Omega)$5-Riesz representation of the discrete boundary flux functional

$f\in L^2(\Omega)$6

Its coefficients are obtained from a boundary mass matrix system: with $f\in L^2(\Omega)$7 the boundary nodal basis,

$f\in L^2(\Omega)$8

The key error identity is

$f\in L^2(\Omega)$9

This identity is the starting point for the $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$0 analysis of the variational flux and explains why $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$1 is especially natural in variational boundary formulations (Pfefferer et al., 2018).

3. Mesh dependence and error estimates

On general shape-regular quasi-uniform meshes, the convergence of discrete normal derivatives is limited by corner regularity. If $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$2, then for both $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$3 and $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$4 the known quasi-uniform rates are $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$5 up to logarithms when $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$6, and otherwise $$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$7 (Pfefferer et al., 2018).

The principal refinement strategy in the finite element analysis is boundary concentration. Writing

$$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$8

the mesh grading condition is

$$(\nabla u,\nabla v)_{L^2(\Omega)}=(f,v)_{L^2(\Omega)} \qquad \forall v\in H_0^1(\Omega).$$9

Thus, elements touching the boundary have diameter $u\in H^2(\Omega)$0, while elements at distance $u\in H^2(\Omega)$1 from the boundary have size $u\in H^2(\Omega)$2. The total number of elements is $u\in H^2(\Omega)$3, only a logarithmic factor above the $u\in H^2(\Omega)$4 complexity of quasi-uniform meshes (Pfefferer et al., 2018).

The analysis uses weighted errors involving

$u\in H^2(\Omega)$5

together with weighted regularity, dyadic decomposition by boundary distance, localized energy estimates, and trace inequalities. For the standard boundary flux, a local estimate on a boundary edge $u\in H^2(\Omega)$6 with adjacent element $u\in H^2(\Omega)$7 is

$u\in H^2(\Omega)$8

The resulting global rates on boundary-concentrated meshes are substantially sharper. In convex polygons with $u\in H^2(\Omega)$9, the standard derivative satisfies

$$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$0

while the variational derivative satisfies

$$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$1

In non-convex polygons, the best-case order drops to $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$2 in the corresponding estimates. In particular, if $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$3, then boundary concentration yields $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$4 up to logarithms, whereas quasi-uniform meshes yield only $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$5 up to logarithms. This is the sense in which boundary-concentrated meshes double the order in $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$6 for normal-derivative approximation (Pfefferer et al., 2018).

A common source of confusion is that this grading is boundary-concentrated, not corner-concentrated. The refinement law is governed by distance to $$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$7, not only by distance to re-entrant corners (Pfefferer et al., 2018).

4. Boundary kernels, perturbations, and control problems

Normal derivatives enter directly into boundary representation formulas. For the Laplacian, the Poisson kernel is

$$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$8

and the solution of the homogeneous Dirichlet problem is represented by

$$\partial_n u=\nabla u\cdot n \quad \text{on }\Gamma.$$9

This identifies the normal derivative of the Green function as a boundary density reproducing interior harmonic values (Martin, 2012).

The same boundary-kernel viewpoint persists under operator perturbation. For the Schrödinger operator $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$0, the outward normal derivative of the Green function satisfies

$$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$1

For the Laplace–Beltrami operator $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$2 with $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$3, the first variation contains a bulk term involving $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$4 and local multiplicative corrections proportional to $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$5 (Martin, 2012). This operator-theoretic structure explains why variational discrete normal derivatives often appear as natural boundary unknowns.

A direct application is Dirichlet boundary control. With desired state $$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$6, the unconstrained problem is

$$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$7

subject to

$$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$8

The continuous optimality system includes

$$(\partial_n u,w)_{L^2(\Gamma)} = (\nabla u,\nabla w)_{L^2(\Omega)}-(f,w)_{L^2(\Omega)} \qquad \forall w\in H^1(\Omega).$$9

so the control is $\partial_n u\in L^2(\Gamma)$0. In the discrete system,

$\partial_n u\in L^2(\Gamma)$1

hence $\partial_n u\in L^2(\Gamma)$2 in $\partial_n u\in L^2(\Gamma)$3. The control error is then limited by the discrete variational normal derivative, and on boundary-concentrated meshes one obtains

$\partial_n u\in L^2(\Gamma)$4

under the regularity assumptions stated for the theorem (Pfefferer et al., 2018).

5. Grid-based directional constructions

On structured grids, a discrete normal derivative is often formed from discrete partial derivatives and a normal vector. One high-order construction is the Vandermonde-based discrete differential operator. Using $\partial_n u\in L^2(\Gamma)$5 samples $\partial_n u\in L^2(\Gamma)$6, one solves

$\partial_n u\in L^2(\Gamma)$7

where $\partial_n u\in L^2(\Gamma)$8 is the Vandermonde matrix, and sets $\partial_n u\in L^2(\Gamma)$9. The theorem gives

$\partial_n u$0

In two dimensions, the estimated gradient components $\partial_n u$1 and $\partial_n u$2 are extracted from the coefficient vector, and the discrete normal derivative is then

$\partial_n u$3

This construction computes all derivatives up to the chosen order simultaneously and makes the normal derivative a post-processing of the discrete gradient (Wang et al., 13 Jul 2025).

A scale-space formulation uses Gaussian derivative operators. For a unit normal $\partial_n u$4,

$\partial_n u$5

The analysis of hybrid discretizations distinguishes normalized sampled Gaussian plus central differences and integrated Gaussian plus central differences from direct sampled or integrated Gaussian derivatives and from the genuinely discrete Gaussian $\partial_n u$6. The reported conclusions are that direct sampled or integrated derivative kernels match continuous spatial spread substantially better than the hybrid methods, while genuinely discrete derivatives behave better than hybrids; hybrid methods remain attractive because one smoothing stage can be shared across multiple derivative orders (Lindeberg, 2024).

A separate directional framework is Tao General Difference. Its first-order finite-window operator is

$\partial_n u$7

and in multiple dimensions the directional operator takes the convolution form

$\partial_n u$8

Choosing $\partial_n u$9 to be the unit normal yields a discrete normal derivative. The orthogonal construction replaces the derivative kernel in the normal direction by smoothing kernels in tangential directions, for example

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$00

which is particularly suited to axis-aligned boundaries on Cartesian grids (Tao et al., 2023).

6. Fractal and lattice analogues, and conceptual distinctions

On p.c.f. fractals, discrete normal derivatives appear as one component of Strichartz’s derivatives at vertices. For a boundary vertex $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$01,

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$02

and the case $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$03 is, up to normalization, the normal derivative. Here $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$04 and $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$05 are left eigenvectors and eigenvalues of the local harmonic extension matrix $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$06. At junction vertices, the branchwise normal derivatives satisfy the compatibility condition

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$07

For functions in $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$08, normal derivatives are uniformly bounded on all vertices. If the normal derivative vanishes at a fixed vertex, then the normal derivatives at neighboring vertices decay to zero with explicit rates governed by the comparison among $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$09 and $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$10, with the three cases

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$11

according to whether $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$12 is greater than, equal to, or less than $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$13 (Cao et al., 2016).

On $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$14, a different, nonlocal notion arises from the derivative at $$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$15 of the fractional discrete Laplacian. In one dimension,

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$16

and more generally

$$-\Delta u=f \quad \text{in }\Omega,\qquad u=0 \quad \text{on }\Gamma=\partial\Omega,$$17

The source explicitly interprets this as a discrete logarithmic Laplacian and a nonlocal boundary/normal derivative-type operator associated with the discrete Laplacian and its fractional powers (Li et al., 1 Mar 2026). This usage is not equivalent to the finite element boundary flux, but it extends the same boundary-oriented vocabulary into a spectral and nonlocal setting.

The main conceptual distinction is therefore structural. In finite elements, the discrete normal derivative is a boundary flux tied to Green’s identity and trace spaces. On structured grids, it is typically a directional derivative projected onto a normal vector. On p.c.f. fractals, it is a renormalized limit of graph-level boundary combinations. On lattices with fractional operators, it can become a nonlocal spectral derivative-type operator. These objects are analogous because each one encodes outward flux, directional variation, or boundary response in a discrete medium, but they are not interchangeable.