🔢 mathematics

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

This paper introduces a divergence-form splitting technique that removes Dirac measures from the conservation law of mixed finite element methods, enabling $L^p$ error analysis and optimal adaptive convergence for the flux without requiring singular solutions or discrete delta functions, regardless of the source's position relative to the mesh.

Original authors: Yueyao Wu, Shun Zhang

Published 2026-08-19

📖 7 min read🧠 Deep dive

CC BY 4.0

Original authors: Yueyao Wu, Shun Zhang

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the world of mathematical modeling, scientists often try to predict how things like heat, electricity, or fluid pressure move through a material. To do this, they break the material down into a grid of tiny shapes, like a mosaic, and solve equations for each piece. This approach, known as the finite element method, is a workhorse of modern engineering. However, it hits a wall when the source of the energy is not spread out, but concentrated at a single, infinitesimal point. Imagine trying to measure the temperature of a room where the heat comes from a pinprick of fire. In the language of mathematics, this is called a Dirac source. The problem is that the standard tools used to solve these equations are built on the assumption that the source can be averaged out over a small area. When the source is a single point, those tools break down because the math requires a value at that exact spot, but the equations themselves do not allow for such a sharp, singular point to exist within the system's rules.

For decades, researchers have tried to fix this by smoothing out the point source or by forcing the grid to align perfectly with the point. These workarounds often require the mesh to be specially designed so that the point sits safely inside a single tile, or they involve replacing the point with a small, artificial blob of data. While these methods can produce results, they are essentially patching over a fundamental flaw in how the problem is set up. They treat the difficulty as a problem of approximation, assuming that if we just make the grid fine enough or the blob small enough, the answer will emerge. But the real issue is deeper: the mathematical framework itself is not designed to handle a point source without breaking the rules of duality, a core principle that ensures the equations make sense together.

A team of mathematicians at the City University of Hong Kong has proposed a different way to look at the problem. Instead of trying to force the point source to fit into the existing grid or smoothing it out, they decided to move the point source out of the main equation entirely. They realized that the difficulty arises because the conservation law—the rule that says what goes in must come out—is being tested against a space that cannot handle a single point. Their solution was to split the problem into two parts. They introduced a specific, known field that carries the weight of the point source and subtracted it from the physical flow. By doing this, they removed the singular point from the main conservation equation, leaving behind a standard, well-behaved problem that any grid can handle. The point source is no longer a mysterious, unmanageable spike in the equation; it has been transformed into a known, explicit field that can be calculated directly.

This change is profound because it means the method no longer depends on where the point sits relative to the grid. In previous approaches, if the point happened to land on a corner where several grid tiles met, the calculation would fail or become ambiguous. With this new splitting technique, the point can sit anywhere—inside a tile, on an edge, or right at a corner—and the math works the same way. The researchers did not need to change the grid, the elements, or the boundary conditions. They simply changed the load vector, the part of the calculation that represents the input force. This makes the method incredibly robust and flexible, allowing it to be applied to complex materials with changing properties without needing to rebuild the entire mathematical structure for every new scenario.

However, removing the point source from the main equation came with a trade-off. While the conservation law became clean and standard, the modified flow field that the researchers were now solving for lost some of its smoothness. In the language of mathematics, this field no longer belongs to the "Hilbert scale," a category of functions that are very smooth and easy to work with using standard tools. Instead, it belongs to a slightly rougher category called the Lebesgue scale. This might sound like a minor technicality, but it is significant. It means that the usual error analysis tools, which rely on the smoothness of the Hilbert scale, could no longer be used. The researchers had to develop a completely new set of tools to measure how close their solution was to the true answer, working in this rougher, less forgiving mathematical space.

Despite this hurdle, the team proved that their method works with high precision. They showed that on a standard, uniform grid, the error in the flow field decreases at a specific rate that depends on the roughness of the field. More importantly, they demonstrated that by grading the mesh—making the grid tiles smaller and smaller as they get closer to the point source—they could recover the best possible rate of convergence. This means that with the right mesh refinement, the method becomes as efficient as the best possible approximation theory allows. They also developed a way to estimate the error after the calculation is done, creating a reliable indicator that tells engineers exactly where the solution needs more refinement, without needing complex extra calculations.

The researchers tested their theory with several numerical experiments, including cases where the point source sits on a sharp interface between two different materials and cases where the material properties change direction. In every scenario, the method held up. They found that the flow field they calculated converged to the true answer at the predicted rate, even though the field itself was too rough for standard methods to handle. They also confirmed that the error estimator they created was accurate, reliably pointing out the areas where the grid needed to be finer. One striking result was that while the flow field converged beautifully in their new mathematical framework, it would have appeared to diverge or fail if analyzed using the old, standard tools. This proved that the difficulty was not with the method itself, but with the mathematical lens through which it was being viewed.

The work also clarified what happens when the point source is perfectly matched to the material properties. In some special cases, the singularity can be canceled out completely, leaving a smooth solution. But in the more general, unmatched cases that are common in real-world applications, the singularity remains. The researchers showed that their method handles these unmatched cases naturally, without needing to know the exact shape of the singularity in advance. This is a crucial advantage because, in many practical problems, the exact nature of the singularity is unknown or too complex to calculate. By using a field that depends only on the location of the point and not on the material properties, the method remains simple and universal.

Ultimately, this research provides a new way to think about problems involving point sources. It shifts the focus from trying to approximate the impossible to reformulating the problem so that the impossible part is handled explicitly and removed from the main calculation. The result is a method that is not only mathematically sound but also practically robust, capable of handling complex geometries and material interfaces without the need for special mesh alignment. By proving that the error can be controlled and estimated even in this rougher mathematical landscape, the authors have opened the door for more accurate simulations of phenomena ranging from heat transfer in composite materials to the behavior of electrical fields around sharp defects. The paper stands as a demonstration that sometimes, the best way to solve a difficult problem is not to work harder at the approximation, but to change the question itself.

Technical Summary: Mixed Finite Element Methods for a Dirac Source

Problem Statement
The paper addresses the numerical solution of the diffusion equation with a Dirac source term in a bounded, simply connected Lipschitz domain $\Omega \subset \mathbb{R}^2$ :
$-\nabla \cdot (A\nabla u) = f_1 + \delta_{x_0} \quad \text{in } \Omega, \quad u = g \quad \text{on } \partial\Omega.$
The authors identify that for mixed finite element methods (MFEM), the primary difficulty with a Dirac source is not merely the low regularity of the solution (which typically lies in $W^{1,p}$ for $p<2$ rather than $H^1$ ), but a fundamental failure of duality. Standard mixed formulations test the conservation equation against a scalar space designed to be a Lebesgue space ( $L^q$ ). However, the Dirac measure $\delta_{x_0}$ does not belong to the dual of any Lebesgue space; consequently, the pairing $\langle \delta_{x_0}, v \rangle$ is undefined for general test functions $v \in L^q$ . This failure occurs at the continuous level, independent of discretization, and cannot be repaired by mesh refinement or element choice alone.

Methodology: Divergence-Form Splitting
To resolve the duality failure without altering the mesh or the approximation space, the authors introduce a divergence-form splitting. They subtract an explicit vector field $f_2$ from the physical flux such that $\nabla \cdot f_2 = \delta_{x_0}$ . Specifically, they define:
$f_2(x) := \frac{1}{2\pi} \frac{x - x_0}{|x - x_0|^2},$
which satisfies $\nabla \cdot f_2 = \delta_{x_0}$ in the sense of distributions and belongs to $L^p(\Omega)^2$ for all $1 < p < 2$ .

The mixed unknown is redefined as the modified flux:
$\sigma := -A\nabla u - f_2.$
This transformation rewrites the problem as a first-order system:
$A^{-1}\sigma + \nabla u = -A^{-1}f_2, \quad \nabla \cdot \sigma = f_1.$
In this new formulation:

The conservation law $\nabla \cdot \sigma = f_1$ involves only the regular data $f_1$ , allowing it to be tested against the full Lebesgue space as required by MFEM.
The singularity is confined to the constitutive equation, where it pairs with the flux test function. Since $f_2 \in L^p$ and test functions $\tau \in L^{p'}$ (where $1/p + 1/p' = 1$ ), Hölder's inequality ensures the pairing is well-defined.
The field $f_2$ depends only on the source location $x_0$ , not on the coefficient $A$ or the mesh. Thus, no coefficient-dependent singular solutions or discrete delta approximations are required. The method works regardless of whether $x_0$ lies in the interior of an element or on a vertex.

Key Contributions and Analysis

Regularity and Error Analysis in $L^p$ Scales:
Because the modified flux $\sigma$ generally retains a singularity of order $|x-x_0|^{-1}$ (unless the coefficient $A$ is a scalar multiple of the identity near $x_0$ ), $\sigma$ belongs to $L^p(\Omega)^2$ for $1 < p < 2$ but not to $L^2(\Omega)^2$ . Consequently, the error analysis must be conducted below the Hilbert scale ( $p=2$ ).
- Flux Error: The authors prove a quasi-best approximation bound for the flux. On quasi-uniform meshes, the flux error converges at the rate $O(h^{2/p - 1})$ . This rate is shown to be sharp and is dictated by the local approximation of the singularity on the element containing the pole.
- Scalar Variable Error: Using a discrete $L^{p'}$ inf–sup condition and a duality argument, the error in the scalar variable $u$ is shown to be asymptotically dominated by the best piecewise constant approximation of $u$ . On quasi-uniform meshes, this yields an $L^2$ error of order $O(h |\log h|^{1/2})$ .
Construction of Comparison Fluxes:
Standard canonical interpolation (Raviart–Thomas) may fail for the modified flux if the singularity's angular component causes non-integrable normal moments on edges passing through $x_0$ . The authors construct an equilibrated comparison flux via local flux balance on a patch surrounding the pole. This ensures the existence of a valid comparison field for the error analysis even when canonical degrees of freedom are undefined.
Mesh Grading and Complexity:
The sub-optimal rate $h^{2/p - 1}$ on quasi-uniform meshes is a property of the singularity, not the method. By employing graded meshes that refine geometrically toward $x_0$ , the authors demonstrate that the flux error can be restored to first-order complexity, $O(N^{-1/2})$ in terms of the number of elements $N$ .
A Posteriori Error Estimation:
The paper establishes a residual norm equivalence in the Lebesgue scale, analogous to first-order system least-squares methods. This leads to a computable $L^p$ estimator for the augmented pair (modified flux, recovered potential). The estimator is shown to be reliable and locally efficient without requiring bubble functions or jump terms.

Results and Numerical Experiments
Numerical experiments using the lowest-order Raviart–Thomas ( $RT_0$ ) and piecewise constant ( $P_0$ ) elements confirm the theoretical findings:

Adaptive Runs: Adaptive refinement driven by the $L^p$ estimator achieves the optimal convergence rate of $O(N^{-1/2})$ for the flux, effectively removing the penalty of the singularity.
Uniform Runs: On uniform meshes, the flux error exhibits the predicted rate $h^{2/p - 1}$ , which varies with $p$ (e.g., $h^{0.667}$ for $p=1.2$ ).
Scalar Convergence: The scalar variable converges at order $O(h)$ (ignoring logarithmic factors) on graded meshes, independent of $p$ , confirming that the duality argument successfully recovers the optimal rate for the potential.
$L^2$ Divergence: The experiments verify that while the $L^p$ error decays, the $L^2$ norm of the discrete flux diverges logarithmically, confirming that the problem cannot be analyzed within the standard Hilbert framework.

Significance
The paper claims that the divergence-form splitting is a robust, mesh-independent remedy for the duality failure inherent in mixed methods with Dirac sources. By shifting the singularity from the conservation law to the constitutive law, the method allows standard mixed formulations to be applied to Dirac problems without modifying the discrete system matrix or introducing artificial discrete deltas. The analysis acknowledges that the cost is a shift to $L^p$ error analysis ( $p<2$ ), but demonstrates that this cost is affordable and yields optimal convergence rates when combined with appropriate mesh grading. The approach is presented as applicable to a broad class of discretizations (conforming, DG, finite volume) by treating the Dirac source as a divergence-form load generated by an $L^p$ field.

Technical Summary: Mixed Finite Element Methods for a Dirac Source

More like this