Mixed Finite Element Methods for a Dirac Source: Divergence-Form Splitting and L^p Error Analysis

Abstract: For a mixed finite element method, a Dirac source is first a failure of duality, not of regularity: the conservation equation is tested against a Lebesgue space, and a Dirac measure lies in the dual of none. We therefore remove the measure from the conservation law by a divergence-form splitting. An explicit field whose divergence is the Dirac measure is subtracted from the physical flux, and the modified flux is taken as the mixed unknown, so that only the regular part of the load remains in the conservation equation. Equivalently, and independently of any discretization, the Dirac problem is rewritten as an elliptic equation whose data are in divergence form, generated by a field of L^p. The subtracted field depends on the location of the source alone and not on the coefficient. No coefficient-dependent singular solution and no discrete delta is needed, only the load vector of the RT_0-P_0 system changes, and the source may sit anywhere relative to the mesh: at a vertex, on a coefficient interface, or inside an element. Unless the splitting is matched to the operator at the source, the modified flux lies in L^p for every p<2 but not in L^2, so the flux error analysis has to leave the Hilbert scale. We prove a quasi-best approximation bound for the flux, and with it that on a quasi-uniform family the flux error is exactly of order h^2/p-1: the matching lower bound comes already from the single element carrying the source. Grading the mesh there restores first-order complexity, N^-1/2 in the number of elements, and the adaptive computations attain it. The scalar variable is limited only by piecewise constant approximation of the solution, which it attains. We also prove a residual norm equivalence in the Lebesgue scale, yielding a computable L^p estimator, reliable and locally efficient for the mixed flux together with a recovered potential.

• The paper introduces an exact divergence-form splitting that removes the Dirac measure from the conservation equation without regularizing the data, allowing mixed finite element methods to handle poles located at mesh vertices and coefficient interfaces.
• The analysis proves sharp quasi-uniform flux convergence of order h^(2/p−1) for 1<p<2, while graded meshes restore first-order complexity O(N^−1/2); scalar errors remain asymptotically first-order up to a logarithmic factor.
• The paper establishes an Lp residual estimator and numerical results confirm the predicted rates, while showing that divergent discrete L2 flux norms reflect the exact singularity rather than a failure of the method.

The failure is one of duality, not of regularity

The paper considers the diffusion problem $-(A\nabla u)=f_1+\delta_{x_0}$ in a bounded simply connected Lipschitz domain $\Omega\subset\mathbb{R}^2$ with Dirichlet data, where $A$ is symmetric uniformly elliptic. Its central diagnosis is that for a mixed method a Dirac source is not primarily an approximation problem but a duality problem (2608.17575). A mixed formulation tests the conservation equation against the scalar space, which by design is a Lebesgue space; since $\delta_{x_0}$ lies in the dual of no Lebesgue space, the pairing is undefined at the continuous level, before any mesh exists. The standard escapes—point evaluation on continuous spaces, discrete representers, or regularization—all act on the datum, and the representer route additionally requires $x_0$ to lie in the interior of an element, which conflicts with fitting the mesh to coefficient interfaces.

The divergence-form splitting

The remedy changes the unknown rather than the datum. An explicit field

$$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$

is subtracted from the physical flux, giving the modified flux $\sigma:=-A\nabla u-f_2$ with $\operatorname{div}\sigma=f_1$. The conservation law becomes an ordinary $L^p$ equation testable against the full Lebesgue space; the singular field survives only in the constitutive equation, where H\"older's inequality suffices. Eliminating $\sigma$ rewrites the problem as $\Omega\subset\mathbb{R}^2$0, exactly equivalent to the original problem with no parameter introduced and no datum perturbed, but with right-hand side in divergence form generated by an $\Omega\subset\mathbb{R}^2$1 field—the standard data class for second-order elliptic problems. The field $\Omega\subset\mathbb{R}^2$2 depends on $\Omega\subset\mathbb{R}^2$3 alone: not on $\Omega\subset\mathbb{R}^2$4, not on lower-order terms, and not on the discretization, so only load assembly changes and $\Omega\subset\mathbb{R}^2$5 may be a mesh vertex.

The paper contrasts this with the ideal splitting using $\Omega\subset\mathbb{R}^2$6 from a singular solution of the operator itself, which returns $\Omega\subset\mathbb{R}^2$7 to $\Omega\subset\mathbb{R}^2$8 when available (e.g., $\Omega\subset\mathbb{R}^2$9 a constant scalar matrix near $A$0). Matching is defined as square-integrability of the modified flux near the pole; for a discontinuous coefficient it is a condition coupling traces from all sides, and matching the singular potential is shown to differ from matching the singular flux. In every unmatched configuration analyzed—coefficient jumps or anisotropy at the pole—the mismatch retains a homogeneous $A$1 part, so that $A$2 for every $A$3 but $A$4. The entire error analysis must therefore leave the Hilbert scale.

Flux error analysis below the Hilbert exponent

The analytic input is Dur\'an's $A$5 stability estimate for the solenoidal projection onto the discrete kernel, stated as Hypothesis (SP) or its two-term form (SP$A$6), in which the degenerating $A$7-dependence multiplies only a data-oscillation term. Under (SP) the paper proves a quasi-best approximation bound over the equilibrated affine subspace $A$8 with no data-oscillation term; under (SP$A$9) the leading constant does not degenerate as $\delta_{x_0}$0. A notable limitation is conceded here explicitly: Dur\'an's result requires a smooth domain, quasi-uniform meshes, and a Lipschitz scalar coefficient, and its stream-function reduction is two-dimensional; outside that setting the stability estimate is kept as a hypothesis.

On quasi-uniform fitted meshes the flux error is exactly of order $\delta_{x_0}$1. The upper bound combines the patch construction with elementwise interpolation; the matching lower bound comes from the single element carrying the pole via a scaling argument requiring only that the leading $\delta_{x_0}$2 part of $\delta_{x_0}$3 be nonzero there. This establishes that the slow rate is a property of the singularity, not of the method. Two further technical contributions are needed: first, an example of a piecewise anisotropic coefficient (two sectors meeting at the pole) for which the canonical Raviart–Thomas normal moments of $\delta_{x_0}$4 are non-integrable on edges through the pole, so the degrees of freedom themselves are undefined—an equilibrated comparison flux is instead constructed by local flux balance on the patch, solvable because the splitting has removed the measure from the divergence theorem; second, grading past the barrier. Balancing local indicators yields graded meshes with exponent $\delta_{x_0}$5 restoring first-order complexity, $\delta_{x_0}$6, versus $\delta_{x_0}$7 on quasi-uniform families ($\delta_{x_0}$8 at $\delta_{x_0}$9).

The scalar variable

Two routes transfer the flux estimates to the potential. The unconditional route uses a discrete $x_0$0 inf–sup condition obtained by lifting the divergence directly in the graph-norm dual, avoiding the inverse inequality whose factor $x_0$1 would destroy convergence entirely. Under dual regularity (piecewise $x_0$2 of the dual problem), a Douglas–Roberts duality argument gains the factor $x_0$3, converting the flux rate into a first-order bound on $x_0$4 independent of $x_0$5. Combined with the Pythagorean identity splitting $x_0$6 into best approximation plus the projected error, the total scalar error is asymptotically the best $x_0$7 approximation of $x_0$8, of order $x_0$9 under a nonvanishing logarithmic leading term. The logarithm is accumulated over all dyadic annuli between $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$0 and a fixed radius—a property of the singularity and of $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$1, not of the method—and grading removes it.

A posteriori estimation in the Lebesgue scale

The a posteriori theory rests on a residual norm equivalence in $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$2: for $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$3 and $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$4,

$$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$5

with the $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$6 well-posedness of the Dirichlet problem replacing the Hilbert-space coercivity unavailable for $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$7. Filling the flux slot with $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$8 yields a computable estimator combining the constitutive residual against a recovered potential $$f_2(x):=\frac{1}{2\pi}\frac{x-x_0}{|x-x_0|^2},\qquad \operatorname{div}f_2=\delta_{x_0},$$9 with the data oscillation $\sigma:=-A\nabla u-f_2$0—no bubble functions, jump terms, or flux reconstruction. It is reliable and locally efficient for the augmented pair $\sigma:=-A\nabla u-f_2$1; reliability for the flux alone follows, while efficiency for the flux alone would require control of the recovery error, which the equivalence does not supply. The paper also notes that reliability and efficiency alone say nothing about an adaptive loop driven by the estimator without a priori control of the solve-and-recover pair.

Numerical verification

Three problems with $\sigma:=-A\nabla u-f_2$2 are computed with the $\sigma:=-A\nabla u-f_2$3–$\sigma:=-A\nabla u-f_2$4 pair at $\sigma:=-A\nabla u-f_2$5: a coefficient jump across an interface through the pole, a constant anisotropic matrix, and the sector coefficient for which the canonical interpolant does not exist. On adaptive meshes all reported orders lie between $\sigma:=-A\nabla u-f_2$6 and $\sigma:=-A\nabla u-f_2$7 in $\sigma:=-A\nabla u-f_2$8, attaining the first-order complexity predicted by the graded analysis, with effectivity indices between $\sigma:=-A\nabla u-f_2$9 and $\operatorname{div}\sigma=f_1$0 stable in $\operatorname{div}\sigma=f_1$1. On uniform refinement the measured flux order agrees with $\operatorname{div}\sigma=f_1$2 to three decimals at $\operatorname{div}\sigma=f_1$3 and $\operatorname{div}\sigma=f_1$4 (e.g., $\operatorname{div}\sigma=f_1$5 vs. $\operatorname{div}\sigma=f_1$6), and $\operatorname{div}\sigma=f_1$7 against $\operatorname{div}\sigma=f_1$8 at $\operatorname{div}\sigma=f_1$9 with a level-by-level rise consistent with the far-field correction scaling $L^p$0. The scalar order is $L^p$1 and $L^p$2-independent, exhibiting the refined duality estimate over both the naive bound (which predicts no convergence) and the basic estimate. Strikingly, the same computation produces a divergent $L^p$3 growing like $L^p$4—matching the exact-flux law—while the $L^p$5 error decays at order $L^p$6: no $L^p$7 error analysis can describe these computations, and the divergence reflects the norm, not a defect of the method.

Limitations and open questions

Several assumptions bear directly on the results. Hypothesis (SP) is unverified on polygons, adaptive meshes, and discontinuous or anisotropic coefficients. Assumption (A3)—the $L^p$8 well-posedness used for reliability—is invoked at $L^p$9 for the computed coefficients without being established there. Dual regularity fails on nonconvex polygons. Efficiency for the flux alone, plain convergence of the adaptive loop built on this estimator, extension beyond two dimensions (where Dur\'an's technique does not apply), the full operator with lower-order terms, and method-specific analyses for DG, hybridizable, finite volume, and nonconforming discretizations applied to the divergence-form equation are all left open, several deferred to separate papers.

Conclusion

The paper replaces a measure-valued load by an exact divergence-form lifting that depends on the pole alone, returning the Dirac problem to a standard data class while paying the price of analyzing the modified flux in $\sigma$0, $\sigma$1. Within that scale it delivers sharp a priori rates on quasi-uniform and graded meshes, first-order scalar accuracy asymptotically equal to best approximation, and a reliable, locally efficient residual estimator, with numerical results confirming each predicted rate—including the exact exponent $\sigma$2 and the logarithmic divergence of the discrete flux in $\sigma$3.