Structure-Preserving Numerical Schemes for a Time-Fractional Porous Medium Equation with Nonlocal Pressure

Introduction

Nonlocal porous medium equations combine a density-dependent diffusivity with interactions acting over long spatial scales. In the fractional-pressure model, these interactions are encoded with the pressure given by the fractional gradient $\nabla(-\Delta)^{-\sigma}u$. In addition, replacing the classical time derivative by a Caputo derivative introduces memory into the evolution. The resulting equation combines nonlinear degeneracy with nonlocality in both space and time, making the preservation of positivity, mass, and stability particularly relevant to its numerical approximation. There are several difficulties in constructing such an approximation. The diffusivity vanishes where the density is zero, while the pressure gradient depends on the density in the whole space. Moreover, the time derivative depends on the previous history of the solution. Hence, the discretizations of the flux and of the memory term have to be considered together when proving stability and convergence.

In this paper, we introduce and analyze numerical approximations of

where $m>1$, $0<\sigma<1$, and $0<\alpha<1$. We take non-negative initial data $u_0\in L^1(\mathbb{R}^{d})\cap H^r(\mathbb{R}^{d})$, with $r\ge d/2+3$, and take the regularity assumption $u\in C([0,T];L^1(\mathbb{R}^{d})\cap H^r(\mathbb{R}^{d}))$. Existence of the solutions of such regularity was shown for $\alpha=1$ and $m=1$ by Choi and Yeong [6]. Also, when additional spatial and temporal regularity is required, we impose it where needed for consistency and convergence. In particular, the convergence estimates below are conditional on the stated consistency bounds. We do not establish existence or uniqueness for the continuous problem in this generality.

For sufficiently regular functions, the Caputo derivative can be written as

On the other hand, the pressure is formally $p=(-\Delta)^{-\sigma}u$ and since only its gradient enters (1), we define the spatial interaction directly by the Fourier multiplier

For smooth, rapidly decaying functions, the corresponding singular-integral representation is

where the normalization constants $c_{d,\sigma}>0$ are known and symmetric limits are used when needed. Defining the gradient directly also avoids the low-frequency normalization of the pressure in one dimension when $\sigma\ge1/2$. This operator was first introduced by Horváth [12]. We write

For the classical time derivative and $m=2$, Caffarelli and Vázquez [3] established the existence of bounded weak solutions with finite propagation speed for compactly supported data. Stan, del Teso, and Vázquez [17] studied the more general mobility $u^{m-1}$ and the dependence of propagation on $m$. The derivation of this equation for the classical time derivative from physical principles was provided by the first author [16]. In one dimension, the use of integrated cumulative variables $v(x,t)=\int_{-\infty}^{x}u(y,t)\,dy$ transforms the density equation to a nondivergence form equation that involves a fractional Laplacian. This formulation has the advantage of providing access to comparison arguments. In addition, a related viscosity theory is developed in [5]. In the context of the mathematical analysis of dislocation density, this equation with $m = 2$ and $\sigma= 1/2$ was analyzed by Biler, Karch, and Monneau [2]. For our time-fractional equation, the same integration in space gives, under sufficient regularity and vanishing flux at $-\infty$,

Notice that $v$ is bounded and non-decreasing when $u$ is non-negative and integrable. This gives a useful class of functions on which to study the comparison and stability of the numerical method. It also allows us to recover the density approximation by discrete differentiation. The relation between these two formulations is used throughout our one-dimensional analysis.

The combination of fractional pressure and a fractional time derivative was studied by Allen, Caffarelli, and Vasseur [1]. For $m = 2$ and $0 < \sigma< 1/2$, they proved the existence and Hölder regularity of weak solutions under suitable decay assumptions on the initial data and forcing. Their results provide an analytical foundation for the time-fractional model, but do not supply a general uniqueness theorem or the strong regularity assumed in our error analysis. Viscosity methods for parabolic problems with a Caputo derivative have also been developed by Topp and Yangari [20].

The closest numerical predecessor is the work of del Teso and Jakobsen [10] for the classical-time equation. They constructed a monotone upwind approximation of the one-dimensional cumulative problem and recovered a conservative density scheme by numerical differentiation. Their convergence theory includes nonsmooth and measure-valued initial data. The cumulative formulation, the choice of the upwind derivative, and the recovery of the density are therefore already available in the classical-time setting. In the present work, we study how this construction can be combined with the L1 history term, both for implicit and explicit time stepping. This requires comparison estimates that retain the contribution of all previous time levels.

A different approach is taken by Carrillo, Fronzoni, and Süli [4], who analyze a finite element method for a parabolically regularized model on bounded domains in dimensions two and three, with pressure defined through a spectral fractional Neumann Laplacian. For the distinct class of equations $\partial_t u - \mathcal{L}[\varphi(u)] = f$, robust numerical methods are developed by del Teso, Endal, and Jakobsen [9]. In that class, a linear nonlocal operator acts on a nonlinear function of the density. In (1), the nonlinearity instead enters through the product of the mobility and the pressure gradient. Thus, although positivity and conservation are common numerical objectives, spatial operators have different structures. The analysis in the following uses the conservative flux in the density variable and the comparison property in the cumulative variable. Our time discretization uses the L1 formula, whose stability and accuracy for time-fractional diffusion were studied, for example, by Lin and Xu [15].

In addition, a relevant issue in its analysis is the regularity of the solution near the initial time. The usual $O(\tau^{2-\alpha})$ estimate requires time regularity, which does not need to hold even for smooth initial data. Jin, Lazarov, and Zhou [13] analyzed the L1 method for linear subdiffusion with smooth and nonsmooth initial data. Stynes, O’Riordan, and Gracia [19] studied the effect of the initial singularity on finite difference approximations on uniform and graded meshes. Kopteva [14] provided a framework for L1 error estimates in the maximum and $L^2$ norms. These works explain why requiring the consistency error to tend to zero uniformly over all time levels can be too restrictive. They also provide the background for the weighted temporal estimate used here.

More recently, del Teso and Płociniczak [11] combined L1 time stepping with monotone spatial discretizations for diffusion equations nonlocal in space and time, obtaining convergence and regularity results in a viscosity framework. Their spatial operator has a prescribed space-time dependent diffusivity. In (3), the coefficient $\lvert\partial_x v\rvert^{m-1}$ depends on the unknown solution and can vanish. For this reason, the monotonicity of the complete nonlinear spatial discretization has to be checked before applying a temporal comparison argument.

The spatial construction used here is based on the fractional powers of the discrete Laplacian studied by Ciaurri et al. [7, 8]. These works also provide the kernel formulas and consistency estimates used in our one-dimensional analysis. For the semigroup pressure-gradient weights, taking the negative discrete divergence gives the fractional power $(-\Delta_h)^{1-\sigma}$. The resulting operator has a difference representation with non-negative weights. In several dimensions, the relevant sign condition concerns the sum over the spatial directions. The individual directional coefficients need not be non-negative. We verify the summed sign condition and the required cancellation and summability properties in the appendix.

The present work combines a conservative upwind flux with implicit and explicit L1 time discretizations. We first study the implicit density scheme in an arbitrary space dimension. Under the stated assumptions on the weights and summability of the fluxes, we establish structural estimates for the implicit scheme, including non-negativity, mass conservation, maximum-norm and energy bounds, and existence of discrete solutions. These results use the upwind form of the flux, the summed sign condition on the spatial weights, and the positivity of the L1 history coefficients. The implicit estimates do not require a CFL restriction.

In one dimension, the cumulative formulation yields comparison and contraction estimates. For $m \ge2$, these estimates give a uniform convergence of the cumulative approximation and a convergence of the density against compactly supported smooth test functions, provided that the assumed spatial and weighted temporal consistency bounds hold. More precisely, the comparison estimate for the nonlinear resolvent gives a scalar L1 inequality for the maximum-norm error. A weighted discrete comparison estimate then controls the accumulation of the consistency errors. This yields the convergence result in Theorem 4.

The temporal analysis allows for the bounds $\lVert\partial_t v(t)\rVert_\infty\le C_v t^{\alpha-1}$ and $\lVert\partial_{tt}v(t)\rVert_\infty\le C_v t^{\alpha-2}$ for $t > 0$. Under these bounds and the stated spatial assumptions, the implicit error estimate takes the form of

Here, $h$ and $\tau$ are the spatial and temporal mesh sizes, $t_n=n\tau$, $V^n$ is the discrete cumulative solution, and $v^n$ denotes the sampled exact solution. The quantity $\eta_A(h)$ is the consistency error of the fractional discrete Laplacian, and $C$ is independent of $h$ and $\tau$. The estimate is uniform over all time levels, including the first one. It follows from a time-weighted residual bound and does not require the unweighted first-step residual to tend to zero. The stated regularity of the cumulative solution remains an assumption on the continuous problem.

We also derive a sufficient CFL condition for the explicit scheme throughout $0<\sigma<1$. The restriction involves $\tau^\alpha$ and reflects the spatial scaling on either side of $\sigma=1/2$, with a logarithmic correction at this value. Under this condition, comparison and contraction give an analogous error estimate, with an additional temporal modulus for the continuous spatial operator. This modulus accounts for evaluating the numerical flux at the previous time level; its required properties are stated in Theorem 5. Thus, the analysis extends the cumulative upwind approach to Caputo time stepping while retaining explicit assumptions on the regularity of the solution being approximated.

The paper is organized as follows. § 2 recalls a continuous maximum principle. § 3 introduces the schemes and develops their structural properties, consistency, and one-dimensional convergence analysis. § 4 presents numerical experiments. Appendix § 6 describes the discretizations of the pressure-gradient and the connection with the discrete fractional Laplacian.

Preliminaries

then

Proof. Since $u(\cdot,t)$ has a global maximum at $x_0$, we have

Expanding $\mathcal{A}u$, we obtain

If $u(x_0,t)=0$, then the non-negativity of $u$ and the maximality of $u(x_0,t)$ imply that $u(\cdot,t)\equiv0$, and the conclusion is immediate. Hence, assume that $u(x_0,t)>0$. Then

Therefore, the first term vanishes at $x_0$. For the second term from the identity

and by the singular-integral representation of the fractional Laplacian, we have

Since $x_0$ is a global maximum, it follows that

and hence

Since $u(x_0,t)^{m-1}\geq0$, we conclude that

which proves the claim.