🔢 mathematics

Positive generalized convolution quadrature for variable-order advection-diffusion: discontinuous Galerkin approximation and compression bounds

This paper presents a positive, stable, and first-order convergent discontinuous Galerkin method for variable-order advection-diffusion equations that utilizes backward-Euler generalized convolution quadrature with composite Gaussian compression to efficiently handle time-dependent fractional orders on arbitrary grids.

Original authors: Lijie Liu, Leilei Wei

Published 2026-10-05

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Lijie Liu, Leilei Wei

Original paper licensed under CC BY 4.0 (https://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

Many natural processes do not move in a straight, predictable line. Instead, they carry a memory of their past, where the speed and direction of change depend on how the system behaved moments ago, hours ago, or even days ago. This is the realm of fractional evolution equations, mathematical tools used to describe phenomena like the slow spread of pollutants through porous rock or the way certain materials relax after being stretched. In these systems, the "memory" is not fixed; it can change over time, becoming stronger or weaker as conditions shift. To model this accurately, scientists use a specific type of mathematical operator that weighs the entire history of a process, but with a twist: the weight given to past events changes based on the current moment. This creates a complex calculation where every new step requires looking back at everything that happened before, making the math incredibly heavy and slow to solve on a computer.

Researchers Lijie Liu and Leilei Wei have developed a new way to solve these difficult equations that is both fast and reliable. They focused on a specific type of problem involving the movement of fluids and heat, where the memory of the system changes over time. Their approach combines a method for handling the history of the system with a technique for breaking space into small, manageable pieces. The core of their innovation is a way to compress the massive amount of historical data needed for each calculation without losing accuracy. Instead of storing every single past value, which would quickly overwhelm a computer's memory, they use a clever mathematical shortcut that summarizes the distant past into a few key numbers while keeping the most recent changes exact. This allows them to run simulations on complex, irregular grids that would have been impossible to handle with previous methods.

The team tested their method on a variety of scenarios, including situations where the rules of the system change abruptly or where the solution starts with a sharp, difficult-to-capture spike. They found that their approach remains stable and accurate even when the time steps between calculations vary wildly, a situation that often causes other methods to fail. In one test, they used a grid where the time steps jumped by a factor of nearly nine thousand between consecutive moments, yet their method held steady. They also compared their technique to a widely used alternative known as the L1 formula. While the L1 method was slightly more accurate for very smooth, predictable problems, the new method proved superior in its ability to handle irregular time steps and sudden changes in the system's behavior without crashing.

A critical part of their work is the ability to prove, while the computer is running, how close the result is to the true answer. They created a system that generates a "certificate" of accuracy, a number that tells the user exactly how much error might be present due to the compression of the history. In their tests, this certificate was remarkably tight, often showing that the difference between their compressed result and a highly precise reference was less than one part in ten billion. This means that scientists can trust the results of these complex simulations without needing to run them multiple times with different settings to check for consistency. The method works by ensuring that the mathematical structure of the solution remains positive and stable, preventing the small errors that usually accumulate in such calculations from growing out of control.

The researchers also explored how their method performs when the solution is not perfectly smooth, which is common in real-world physics where materials might have cracks or sudden changes in density. They found that by adjusting the density of their calculation points in time, they could capture these sharp changes effectively. When they used a standard, evenly spaced grid, the accuracy dropped significantly for these difficult cases. However, by clustering more calculation points near the start of the process where the changes are most rapid, they recovered the high level of accuracy. This flexibility suggests that their method can be adapted to a wide range of physical problems, from modeling the flow of blood in arteries to simulating the diffusion of heat in materials with varying properties.

Ultimately, this work provides a robust framework for simulating systems with changing memory laws. It removes the strict requirement that time steps must be uniform or that the system's behavior must change slowly. By separating the memory states from the changing rules of the system, the researchers created a tool that is both computationally efficient and mathematically sound. The ability to handle arbitrary time steps and sudden changes in the order of the memory operator opens the door to more realistic simulations of complex physical phenomena. The method does not just offer a faster way to calculate; it offers a way to calculate things that were previously too unstable to model with confidence, providing a new level of certainty for scientists studying the memory-laden processes of our world.

Technical Summary: Positive Generalized Convolution Quadrature for Variable-Order Advection-Diffusion

Problem Statement
The paper addresses the numerical solution of variable-order advection-diffusion equations involving a Caputo derivative where the fractional order $\alpha(t)$ is evaluated at the current observation time. The governing equation is defined on a bounded domain $\Omega \subset \mathbb{R}^d$ ( $d=2,3$ ) with homogeneous Dirichlet boundary conditions. A key challenge in this formulation is that changing the order $\alpha(t)$ reweights the entire history of the solution, complicating the construction of efficient and stable numerical schemes. Standard convolution quadrature (CQ) methods often struggle with non-uniform time steps and variable orders, while existing fast history evaluation techniques may not preserve the necessary positivity and stability properties required for variable-order operators.

Methodology
The authors propose a fully discrete scheme combining a Discontinuous Galerkin (DG) spatial discretization with a specialized temporal discretization:

Spatial Discretization: The method employs the Symmetric Interior Penalty Galerkin (SIPG) method for diffusion and an upwind flux for convection. The penalty parameter is chosen conservatively to ensure spatial accretivity, independent of mesh size or time-step ratios.
Temporal Discretization (Generalized CQ): The time derivative is approximated using a Backward-Euler Generalized Convolution Quadrature (CQ). The operator is represented via a diffusive identity involving an integral of discrete resolvent products.
History Compression: To reduce the computational cost of the non-local history term, the authors apply composite Gauss-Legendre quadrature to the resolvent products.
- Auxiliary States: The history is compressed into a set of auxiliary states that depend on the time grid but are independent of the fractional order $\alpha(t)$ .
- Exact Diagonal Correction: Crucially, the scheme retains the exact coefficient for the current time increment ( $\tau_n^{-\alpha_n}$ ) rather than approximating it via quadrature. This ensures the newest memory contribution is preserved exactly, confining the compression error to older increments.
- Positivity Preservation: The integrand's complete monotonicity ensures that the Gaussian quadrature underestimates the integral, and the diagonal correction remains non-negative. This preserves the "positive memory" structure essential for stability.
Error Certification: The paper introduces computable scalar states (certificates) that provide upper bounds for the Gaussian quadrature remainder and the omitted tails of the history integral. These bounds, combined with the algebraic residual, allow for a rigorous estimate of the discrepancy between the compressed scheme and the exact-quadrature scheme.

Key Contributions
The paper makes three primary theoretical and practical contributions:

Stability without Restrictions: The authors prove $L^2$ -norm stability for the scheme on arbitrary positive time grids. Unlike many existing methods, this stability holds without restrictions on the ratio of adjacent time steps or the variation of the fractional order $\alpha(t)$ . The proof relies on a Hilbert-space comparison argument utilizing the positive, increasing nature of the increment weights within each time row.
Convergence Analysis:
- For smooth solutions, the method achieves first-order temporal convergence ( $O(\tau^*)$ ) and optimal spatial convergence ( $O(h^{p+1})$ ), where $\tau^*$ is the maximum time step and $p$ is the polynomial degree.
- The analysis explicitly handles nonzero initial slopes and weak initial singularities, the latter of which is addressed numerically using graded time grids.
Computable Compression Bounds: The paper derives explicit, computable bounds for the error introduced by history compression. These bounds depend on the solution's increments and derivative states, allowing the user to certify the accuracy of the compressed solution relative to the exact-quadrature scheme without needing to know the exact solution.

Numerical Results
The proposed method is validated through a series of numerical experiments in two dimensions:

Temporal Accuracy: Tests with smooth and weakly singular solutions confirm the theoretical first-order convergence rate. The method is compared against the non-uniform L1 formula; while L1 shows slightly higher accuracy for smooth data (due to higher approximation order), the proposed CQ method demonstrates robust first-order behavior and superior stability on highly irregular grids.
Spatial Accuracy: Experiments on perturbed triangular meshes confirm optimal spatial convergence rates ( $p+1$ in $L^2$ and $p$ in energy norms).
Compression and Stability: The method successfully handles highly irregular time steps (adjacent ratios up to ~8800) and rapidly oscillating or jumping fractional orders. The computed error bounds consistently overestimate the actual difference between the compressed and reference solutions, validating the certification mechanism.
Efficiency: The compression reduces the storage and computational cost from $O(N^2)$ to $O(NQ)$ , where $Q$ is the number of quadrature states, though the paper notes that for high-accuracy requirements, $Q$ may need to be large.

Significance and Claims
The paper claims that its construction successfully decouples the dependence on the fractional order from the memory states, allowing for efficient reweighting of history without reconstructing past values. The preservation of the positive row structure through exact diagonal correction is identified as the key mechanism enabling stability on arbitrary grids.

The authors modestly frame their work as a foundation for variable-order problems, noting that while the compression analysis applies to arbitrary discrete sequences, the total PDE error estimates rely on specific regularity assumptions. They acknowledge that the current first-order CQ construction involves a trade-off in accuracy compared to higher-order L1 methods for smooth data, but emphasize the robustness of the stability properties. The paper concludes that while extensions to higher-order CQ and sharper nonsmooth analysis remain open, the proposed framework provides a stable, certifiable, and efficient approach for variable-order advection-diffusion problems with arbitrary time stepping.

Technical Summary: Positive Generalized Convolution Quadrature for Variable-Order Advection-Diffusion

More like this