Numerical recovery of the diffusion coefficient in diffusion equations from terminal measurement

Introduction

Inverse problems in partial differential equations (PDEs) represent a cornerstone of modern computational science. One of the most challenging tasks in this field is the recovery of a space-dependent diffusion coefficient from a single terminal measurement. This problem, often encountered in hydrology, heat transfer, and anomalous diffusion modeling, is inherently ill-posed and nonlinear. In a recent study published in Numerische Mathematik, researchers Bangti Jin, Xiliang Lu, Qimeng Quan, and Zhi Zhou present a rigorous stability and error analysis for this problem, covering both standard and fractional diffusion models.

The Inverse Problem

Consider a system where a substance spreads through a medium. The process is governed by a diffusion or subdiffusion equation, where the diffusion coefficient $q(x)$ varies across the spatial domain $\Omega$. In many real-world applications, continuous monitoring is impossible. Instead, we only have access to the state of the system at a final time $T$, denoted as $u(x, T)$. The objective is to reconstruct $q(x)$ from this single snapshot.

This task is complicated by three primary factors:

1. Ill-posedness: Small perturbations or noise in the terminal data can lead to disproportionately large errors in the estimated coefficient.
2. Nonlinearity: The mapping from the coefficient to the terminal state is non-linear, requiring sophisticated optimization techniques.
3. Fractional Dynamics: In subdiffusion models, which utilize fractional time derivatives to account for "memory" effects in complex media, standard parabolic analysis fails.

Methodological Framework

The authors employ an output least-squares formulation to tackle the inversion. To ensure stability, they utilize Tikhonov regularization with an $H^1(\Omega)$-seminorm penalty. This approach penalizes excessive oscillations in the reconstructed coefficient, favoring physically plausible, smoother solutions.

The numerical implementation involves:

• Spatial Discretization: A Galerkin Finite Element Method (FEM) using piecewise linear elements.
• Temporal Discretization: A Backward Euler Convolution Quadrature (BECQ) scheme, which is particularly effective for handling the non-local nature of fractional derivatives.

Theoretical Contributions

The core contribution of the paper is the establishment of a Hölder-type conditional stability estimate. The researchers proved that for a sufficiently large terminal time $T$, the difference between two potential coefficients can be bounded by the difference in their corresponding terminal states. This provides a mathematical guarantee that the reconstruction is stable under controlled conditions.

Furthermore, the paper provides the first rigorous $L^2$ error bounds for the numerical approximation. Unlike many previous studies that rely on "source-type conditions"—abstract mathematical assumptions that are difficult to verify in practice—this analysis provides concrete convergence rates based on observable parameters like noise level ($\delta$), mesh size ($h$), and regularization parameter ($\gamma$).

Numerical Validation

The theoretical findings were validated through synthetic experiments in 1D and 2D domains.

• 1D Case: Recovery of $q^\dagger(x) = 2 + \frac{1}{2}\sin(2\pi x)$ on $(0, 1)$.
• 2D Case: Recovery of $q^\dagger(x, y) = 1 + \frac{1}{2}y(1-y)\sin(\pi x)$ on a square domain.

The experiments confirmed that the stability of the recovery significantly improves as the terminal time $T$ increases. This "large $T$ advantage" allows the system's evolution to sufficiently imprint the characteristics of the diffusion coefficient onto the final state. Additionally, the empirical convergence rates often exceeded the theoretical $O(\delta^{1/4})$ prediction, suggesting that the derived bounds are robust.

Practical Implications

For practitioners, the paper offers specific algorithmic guidelines. To achieve optimal results, the regularization parameter should scale as $\gamma \approx \delta^2$ and the mesh size as $h \approx \sqrt{\delta}$. These heuristics provide a clear roadmap for setting up numerical solvers for diffusion-based inverse problems.

Conclusion and Future Directions

While the paper marks a significant advancement, it also highlights remaining challenges. The current stability results require $T$ to be large; when $T$ is very small, the information content in the terminal data is insufficient for accurate recovery. Future research may focus on extending these results to non-convex domains, discontinuous coefficients, or more complex boundary conditions. Nonetheless, this work provides a much-needed theoretical foundation for the reliable numerical recovery of diffusion parameters from terminal data.