The Allen-Cahn equation with a time Caputo-Hadamard derivative: Mathematical and Numerical Analysis

1.   Introduction

Fractional calculus has been increasingly gaining attention from researchers due to its widespread applications in fields such as physics, chemistry and biology ^[1,2,3,4,5]. The most commonly used approaches entail the use of the Riemann-Liouville integral/derivative, Caputo derivative, Riesz derivative and fractional Laplacian. However, there is another type of fractional calculus that was discovered in 1892 but largely overlooked, and it is Hadamard calculus ^[6]. It was only in recent years that people realized its ability to provide more accurate descriptions of complex processes, such as Lomnitz's logarithmic creep law for special materials ^[7] and ultra-low diffusion processes ^[8]. As a result, it has gradually come into the spotlight ^[9,10].

The Allen-Cahn equation (AC equation) was originally proposed by Allen and Cahn in 1979 while studying the motion of phase boundaries in crystalline solids as a model for phase separation processes in binary alloys at a given temperature ^[11,12]. Over the years, it has become one of the most widely used phase field models for describing physical phenomena in materials science and fluid mechanics ^[13]. Early studies of the AC model primarily focused on the following integer-order partial differential equation:

where $\varepsilon > 0$ represents the interface width parameter and $F(u) = \frac{1}{4}(1-u^2)^2$ denotes a double-well potential ^[14,15,16]. The particle motion in model (1.1) follows Brownian motion, meaning that the mean square displacement (MSD) satisfies $\langle(x(t))^2\rangle\simeq t$), and the forces are spatially local, i.e., long-range interactions between particles are not considered. However, due to the heterogeneity of the medium, non-local interaction forces unavoidably exist in phase field models, which cannot be described in integer-order models. Therefore, an increasing number of researchers have started to focus on studying fractional AC equations:

where $_{C}{\rm D}_{0, t}^{\alpha}$ is the Caputo fractional derivative defined by

The fractional AC equation (1.2) describes subdiffusion phenomena in nature, characterized by the power-law growth of the MSD with time (i.e., $\langle(x(t))^2\rangle\simeq t^\alpha$).

However, an equally important and significant fact is that there are also many ultra-low diffusion behaviors in nature. This is because diffusing particles have a heavy-tailed waiting time distribution, which is slower than any power-law decay. Their MSDs exhibit logarithmic growth over time, i.e., $\langle(x(t))^2\rangle\simeq (\log t)^\alpha$. Describing these phenomena using Hadamard calculus would be more accurate. Therefore, we consider the following form of the Caputo-Hadamard-type time-fractional AC equation:

where $_{CH}{\rm{D}}_{a, t}^{\alpha} \, (0 < \alpha < 1)$ represents the Caputo-Hadamard fractional derivative defined in (2.6), and the domain $\Omega\subset\mathbb{R}^d\; (d\geq 1)$ is a bounded and convex polygon.

Although there have been many studies on the AC equation with a time Caputo derivative (see (1.2)), there seems to be no report on the time derivative in the Caputo-Hadamard sense. This paper will focus on the following two goals for this type of situation:

● Provide some theoretical results on the solution of problem (1.3), including the existence, uniqueness and regularity of the solution in certain spaces.

● Consider the numerical solution of problem (1.3), use the local discontinuous Galerkin (LDG) method for spatial approximation and discretize the time direction using a nonuniform difference formula to obtain the corresponding fully discrete scheme. The stability and convergence of this scheme are demonstrated through numerical examples.

The organization of this article is as follows. In Section 2, we provide some symbols and definitions. By utilizing the modified Laplace transform and its inverse form, the mild solution of (1.3) is derived. Section 3 mainly focuses on some theoretical analysis of (1.3). In Section 4, we present a nonuniform L1/LDG scheme for (1.3) with generalized alternating numerical fluxes. The stability analysis and error estimation of this scheme are investigated. Numerical examples are presented in Section 5. A brief summary is provided in the final section.

2.   Preliminaries

In this section, we review some common symbols and definitions and provide a representation for the solution of (1.3).

For any measurable subset $\omega$ of $\Omega$, let $(\cdot, \cdot)_\omega$ be the $L^2$ inner product on $\omega$ and $\|\cdot\|_{\omega}$ denote the $L^2(\omega)$ norm defined by $\|v\|_{\omega}^2 = (v, v)_\omega$. Here, we omit the subscript when $\omega = \Omega$. For each nonnegative integer $r$, $H^r(\omega)$ denotes the usual Sobolev space with its associated norm $\|\cdot\|_r$.

Suppose $\mathcal{S}$ is a real Banach space with norm $\|\cdot\|_{\mathcal{S}}$. $C\left([a, T]; \mathcal{S}\right)$ represents the space consisting of all continuous functions $v:[a, T]\rightarrow \mathcal{S}$, whose norm is defined as

To simplify the notation, in this paper we assume that $\varepsilon = 1$ in (1.3). Set $X = L^2(\Omega)$, $\mathcal{D} = H_0^1(\Omega)\cap H^2(\Omega)$ and $\mathcal{A} = \Delta$, with a homogeneous Dirichlet boundary condition. Then the operator $\mathcal{A}$ satisfies the following resolvent estimate ^[17]:

where $\|\cdot\|_{X\rightarrow X}$ is the operator norm on the space $X$, $\mathcal{I}$ is the identity operator and $\Sigma_\theta: = \{z\in \mathbb{C}\backslash\{0\}:|\arg (z)|\leq\theta\}$. One can thus get from (2.2) that

Definition 2.1 (^[18]). The Hadamard fractional integral of a given function $f(t)$ with order $\alpha > 0$ is defined as

where $\Gamma(\cdot)$ is the usual Gamma function.

Definition 2.2 (^[18]). The Hadamard fractional derivative of a given function $f(t)$ with order $\alpha (n-1 < \alpha < n\in\mathbb{Z}^+)$ is defined as

where $\delta^n g(x) = \left(x{\rm d}/{\rm d}x\right)^n g(x) = \delta (\delta^{n-1}g(x))$.

Definition 2.3 (^[8,19]). The Caputo-Hadamard fractional derivative of a given function $f(t)$ with order $\alpha (n-1 < \alpha < n\in\mathbb{Z}^+)$ is defined as

For more information on Definitions 2.1–2.3, one may refer to ^{[20,21,22,23]}.

Definition 2.4 (^[21]). The modified Laplace transform of a given function $f(t)$ with $t\in[a, +\infty)\, (a > 0)$ is defined by

The inverse modified Laplace transform is given by

Definition 2.5 (^[21]). For the given functions $f(t)$ and $g(t)$ defined on $[a, +\infty)\, (a > 0)$, the convolution is defined by

From the definition of operator $\delta$, it is easy to see that

Let $w = u-u_a$. Then, one can get from (1.3) that $w$ satisfies the following equations:

By virtue of the modified Laplace transform, one has

which further yields

According to the inverse modified Laplace transform and the convolution rule, one gets

where the operators $\mathcal{E}(t), \, \mathcal{F}(t):X\rightarrow X$ are defined as

For fixed $\vartheta > 0$ and $\theta\in(\frac{\pi}{2}, \pi)$, $\Gamma_{\theta, \vartheta}$ is given by

in which Im $s$ increases.

As a consequence, the mild solution of (1.3) can be derived, that is,

Here and hereafter $u(t) = u({\bf x}, t)$ with ${\bf x}\in\Omega$.

Now we study the properties of the operators $\mathcal{F}(t)$ and $\mathcal{E}(t)$.

Lemma 2.1. The operators $\mathcal{E}(t)$ and $\mathcal{F}(t)$ defined in (2.14) satisfy the following properties:

(ⅰ) For $t\in[a, T]$, $\mathcal{F}(t):X\rightarrow \mathcal{D}$ is continuous and $\mathcal{A}\mathcal{F}(a) = 0$.

(ⅱ) For $t\in(a, T]$, it holds that

(ⅲ) For $t\in(a, T]$, it holds that

(ⅳ) For $t\in(a, T]$, $\mathcal{E}(t):X\rightarrow \mathcal{D}$ is continuous.

Proof. (ⅰ) For $t\in[a, T]$, Sakamoto and Yamamoto have shown in [24, Theorem 2.1] that $\mathcal{A}\mathcal{F}(t) = \mathcal{F}(t)\mathcal{A}:X\rightarrow X$ is continuous. Thus, $\mathcal{F}(t):X\rightarrow \mathcal{D}$ is continuous with respect to $t\in[a, T]$. Letting $f(u) = 0$ in (2.16), and taking the limit as $t\rightarrow a$, one can deduce $\mathcal{A}\mathcal{F}(a) = 0$.

(ⅱ) Notice that $\mathcal{A}(z^\alpha-\mathcal{A})^{-1} = -\mathcal{I}+z^\alpha(z^\alpha-\mathcal{A})^{-1}$, and choose $\vartheta = (\log t-\log a)^{-1}$ in the contour $\Gamma_{\theta, \vartheta}$; then, for any nonnegative integers $k$, one has

So we get (ⅱ).

(ⅲ) Because $\mathcal{E}(t) = \delta \mathcal{F}(t)$, (ⅲ) follows from (2.17).

(ⅳ) Using the equivalent norm $\|v\|_{\mathcal{D}}\sim \|v\|_X+\|\mathcal{A}v\|_X, \; \forall v\in \mathcal{D}$, the conclusion in (iv) can be obtained.

3.   Regularity Analysis

In order to provide a theoretical basis for the numerical analysis that follows, we will consider the existence, uniqueness and regularity of solutions to (1.3) in the present section.

Theorem 3.1. Suppose $f:\mathbb{R}\rightarrow \mathbb{R}$ is a Lipschitz continuous function, i.e.,

Assume that $u_a\in H_0^1(\Omega)\cap H^2(\Omega)$. Then, (1.3) has a unique solution $u$ that satisfies

where $C$ is a positive constant.

Proof. First, we prove the existence of a unique solution to (1.3). For any fixed $\lambda > 0$, we denote by $C\left([a, T]; X\right)_\lambda$ the weighted norm space of function $v\in C\left([a, T]; X\right)$, equipped with the norm

Let $\mathscr{M}: C\left([a, T]; X\right)_\lambda\rightarrow C\left([a, T]; X\right)_\lambda$ be a nonlinear map defined as

Then, for any $v_1(t), v_2(t)\in C\left([a, T]; X\right)_\lambda$, by virtue of Lemma 2.1 and the Lipschitz continuity of $f$, we have

Thus, by choosing a sufficiently large $\lambda$, the following inequality holds:

which means that $\mathscr{M}$ is a contractive mapping on the space $C\left([a, T]; X\right)_\lambda$. Then, based on the contraction mapping principle and the equivalence of spaces $C\left([a, T]; X\right)_\lambda$ and $C\left([a, T]; X\right)$, we obtain that (2.16) has a unique fixed point $u\in C\left([a, T]; X\right)$.

We now prove that there exists a positive constant $C$ such that

When $s = t$, (3.6) obviously holds. We focus on proving the case of $a\leq s < t\leq T$. The same can be obtained for the other case. According to (2.16), one has

For the first term, we apply Lemma 2.1–(ⅱ) and the Minkowski inequality to get

For the second term, we can obtain from Lemma 2.1–(ⅱ) that

Similarly, the third term can be bounded as

Denoting

and substituting (3.8)–(3.10) into (3.7) yield

Hence, by choosing a sufficiently large $\lambda$, we can achieve the desired result.

Applying the operator $\mathcal{A}$ to both sides of (2.16), and noting that

we get

By directly applying Lemma 2.1–(ⅰ), we can obtain $\Phi_1(t)\in C\left([a, T]; X\right)$ and

In view of Lemma 2.1–(ⅲ), one has

which implies that $\Phi_2(t)$ is continuous at $t = a$. Meanwhile, by using Lemma 2.1–(ⅳ), we know that $\Phi_2(t)$ is continuous for $t\in (a, T]$. Therefore, $\Phi_2(t)\in C\left([a, T]; X\right)$. Combining the previous three estimates, we can see that

This, together with (1.3), also show that $_{CH}{\rm D}_{a, t}^{\alpha} u\in C\left([a, T]; L^2(\Omega)\right)$.

Finally, the term $\|\delta u(t)\|$ has not been estimated yet. By differentiating (2.16) with respect to variable $t$, we obtain

Multiplying both sides of this equation by $e^{-\lambda(\log t-\log a)}(\log t-\log a)^{1-\alpha}$, and by using Lemma 2.1, one gets

By taking maximum of the left-hand side over $t\in[a, T]$ and choosing a sufficiently large $\lambda$, we obtain

All of this completes the proof.

Lemma 3.1. Let $Dv(t) = (\log t-\log a)\delta v(t)$. Then, the following relation holds:

Proof. Recalling Definition 2.4, the following relation holds:

Lemma 3.2 (Gronwall inequality ^[25]). Suppose that $f(t)$ and $g(t)$ are nonnegative integrable functions on $[a, b]$. If there exists a nonnegative constant $C_1$ such that

then

In particular, if $g(t)$ is non-decreasing, then

Based on the above lemmas, we next show that the solution $u({\bf x}, t)$ of (1.3) satisfies higher regularity.

Theorem 3.2. Assume that $f$ satisfies the condition in Theorem 3.1 and $u_a\in H_0^1(\Omega)\cap H^4(\Omega)$. Then, for $t\in(a, T]$, it holds that

and

where $C$ is a positive constant.

Proof. Step 1. By applying the operator $D$ on both sides of (2.16) and utilizing Lemma 3.1, one gets

Applying the Laplace operator $\mathcal{A}$ to both sides of (3.17) further implies the following

Recalling the Sobolev embedding formula $\|u\|_{L^\infty(\Omega)}+\|\nabla u\|_{L^4(\Omega)}\leq C\|u\|_2$, and by using the fact that $u\in C\left([a, T]; H_0^1(\Omega)\cap H^2(\Omega)\right)$ from (3.2), one has

Likewise,

Combining (3.18)–(3.20) and Lemma 2.1, one gets

By virtue of Lemma 3.2, we obtain

which confirms the case of $l = 1$ in (3.15).

Step 2. In view of the definition of the Caputo-Hadamard fractional derivative in (2.6), we obtain

where the second inequality is derived by (3.22). As a result, we get from (1.3) and (3.19) that

By virtue of (1.3), we know that $_{CH}{\rm D}_{a, t}^{\alpha} u|_{\partial\Omega} = \mathcal{A} u|_{\partial\Omega} = 0$ for $t\in[a, T]$. Noting that $\Omega$ is a bounded convex polygonal domain, we can therefore prove that (3.16) holds according to (3.23) and (3.24).

Step 3. Now, we demonstrate the case $l = 2$ in (3.15). Apply operator $\delta$ on both sides of (3.17) to obtain that

Then, applying the Laplace operator $\mathcal{A}$ on both sides of (3.25) further leads to

From (3.20) and (3.22), one can deduce that

Again use the embedding theorem $\|u\|_{L^\infty(\Omega)}+\|\nabla u\|_{L^4(\Omega)}\leq C\|u\|_2$ to show that

Similar to the proof of (3.20), one can get that

Therefore, by the assumption $u_a\in H^4(\Omega)$, (3.26)–(3.29) and Lemma 2.1, one can derive that

This, together with Lemma 3.2, leads to the desired result.

4.   Nonuniform L1/LDG Discretization of the Time-Fractional AC Equation

As shown in Theorems 3.1 and 3.2, the solution $u({\bf x}, t)$ of problem (1.3) may behave as weakly regular at the starting time $t = a$. Thus, we utilize the L1 scheme on nonuniform meshes (see ^[21] for more information about this scheme) to discretize the time Caputo-Hadamard derivative and by using the LDG method in space. Without loss of generality, suppose that the bounded domain $\Omega = (-1, 1)^d$ in (1.3) and $f(u)$ satisfies

where $C$ is a positive constant.

In the next analysis, we will consider the cases $d = 1$ and $2$. The more general case $d > 2$, which can also be obtained by changing the tensor product structure of the mesh, is omitted here.

4.1.   Nonuniform L1 Approximation in the Logarithmic Sense

For $r\geq 1$, denote $t_{n} = a\left(T/a\right)^{(n/M)^{r}}$, where $n = 0, 1, \cdots, M$, $M\in \mathbb{N}$. We divide the interval $[a, T]$ into a grading mesh in the logarithmic sense, that is, $\log a = \log t_{0} < \log t_{1} < \cdots < \log t_{n-1} < \log t_{n} < \cdots < \log t_{M} = \log T$ with

Let $\tau_n = \log t_n-\log t_{n-1}$, $n = 1, \ldots, M$ be the time mesh sizes.

The nonuniform L1 approximation in the logarithmic sense for the Caputo-Hadamard derivative ^[21] at $t = t_n$ is defined as

where the discrete coefficients and the local truncation error are given, respectively, by

and

Denote $a^{(n)}_{n-k} = b_{n, n-k+1}/\Gamma(2-\alpha)$, $k = 1, 2, \cdots, n$, and

Letting $\omega_\beta(t) = t^{\beta-1}/\Gamma(\beta)$, we use (4.2) to obtain

Similar to [26, Lemma 2.1–(ⅱ)], one can prove that

In view of the integral mean-value theorem, one has

For simplicity, we denote $u^n = u({\bf x}, t_n)$; then, the nonuniform L1 approximation scheme given in (4.2) can be rewritten as

Lemma 4.1. ^[21] Let the function $u({\bf x}, t)$ satisfy that $\left|\delta^{l} u(\cdot, t)\right|\leq C\left(1+(\log t-\log a)^{\alpha-l}\right)$ for $l = 0, 1, 2$ and all $t\in(a, T]$. Then, it holds that

Lemma 4.2. Assume that $u(\cdot, t)\in C^2(a, T]$ and $\left|\delta^{l} u(\cdot, t)\right|\leq C\left(1+(\log t-\log a)^{\alpha-l}\right)$ for $l = 0, 1, 2$ and all $t\in(a, T]$. Then for $n = 1, 2, \cdots, M$, the following inequality holds:

Proof. The proof of this lemma is analoguous to that of (3.12) in ^[26], so is omitted here.

Lemma 4.3 (Discrete Gronwall inequality). Let $(\lambda_l)^{M-1}_{l = 0}$ be a nonnegative sequence and there exist a constant $\lambda$ independent of time steps such that $\sum_{l = 0}^{M-1} \lambda_l\leq \lambda$. Assume that the sequences $\{\phi^n\}_{n = 1}^M$ and $\{\psi^n\}_{n = 1}^M$ are nonnegative, and that the grid function $\{v^n\}_{n = 1}^M$ satisfies

If the maximum time-step $\tau_M\leq (2\Gamma(2-\alpha)\lambda)^{-1/\alpha}$, the following holds:

where $E_{\alpha, 1}(z)$ is the well-known Mittag-Leffler function.

Proof. Denote

If $v^n\leq \Psi^\ast: = \sqrt{\Gamma(1-\alpha)} \max_{1\leq k \leq n}\{(\log t_k-\log a)^{\alpha/2}\psi^k\}$, then (4.9) is directly obtained from $\mathcal{E}_\alpha^n\geq 2$. For the alternative case $v^n > \Psi^\ast$, we have $v^n > \sqrt{\Gamma(1-\alpha)}\{(\log t_n-\log a)^{\alpha/2}\psi^k\}$, and the inequality (4.8) can be rewritten as

Using [27, Lemma 3.6] with

we get from (4.4) that

The proof is completed.

Remark 4.1. The conclusion in Lemma 4.3 provides the theoretical support for the numerical approach to the Caputo-Hadamard fractional differential equation. The results are almost identical to the usual nonuniform L1 formula (for Caputo fractional derivative, see [28, Theorem 2.3] for details).

4.2.   Notations and Projections of the LDG Method

Let us denote by $\Omega_h = \{K\}$ a shape-regular subdivision of $\Omega$, and set $\partial\Omega_h = \{\partial K: K\in\Omega_h\}$. Suppose that the "left" and "right" elements $K_L$ and $K_R$ share a face $e$, and $\varphi$ is a function defined on $K_L$ and $K_R$, but may be discontinuous on $e$. Then, we use $\varphi_L$ and $\varphi_R$ to denote the traces of $e$ from the left and right direction, respectively. The finite element space associated with the mesh $\Omega_h$ is of the form

where $\mathcal{Q}^{k}(K)$ is a tensor product space defined over $K$ with maximal $k$-th polynomial. When d = 1, $\mathcal{Q}^{k}(K) = \mathcal{P}^k(K)$.

Case A ($d = 1$): For an arbitrary element $K: = I_j = (x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}})$ with $j = 1, 2, \cdots, N$, we denote $x_j = (x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2$, $h_j = x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ and $h = \max_{1\leq j\leq N} h_j$. Obviously, $x_{\frac{1}{2}} = -1$ and $x_{N+\frac{1}{2}} = 1$. Let $\Omega_h$ be a quasi-uniform mesh, that is, there exists a fixed positive constant $\rho$ independent of $h$ such that $\rho h\leq h_j\leq h$ for $j = 1, 2, \cdots, N$ as $h\rightarrow 0$.

Let $\mathscr{P}_{h}:L^2(\Omega)\rightarrow V_h$ represent the standard $L^2$ projection, defined as

The Gauss-Radau projections $\mathscr{P}_{h}^\pm: H^{1}(\Omega)\rightarrow V_{h}$ are given by ^[29]

and

Case B ($d = 2$): For an arbitrary rectangular element $K: = K_{ij} = I_i\times J_j = (x_{i-\frac{1}{2}}-x_{i+\frac{1}{2}})\times(y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}})$, we denote $h_i^x = x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $h_j^y = y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$. Analogous to the one-dimensional case, $h_{ij} = \max\{h_i^x, h_j^y\}$ and $h = \max_{K_{ij}\in \Omega_h}h_{ij}$ are well defined. We also list the projections that will be used ^[30].

$-$ The projection $\Pi_h^-:H^1(\Omega)\rightarrow {\bf V}_h$ for scalar functions is defined as

where $\mathscr{P}_{h, x}^-$ and $\mathscr{P}_{h, y}^-$ represent the one-dimensional projection $\mathscr{P}_{h}^{-}$ given in (4.14) on a two-dimensional rectangular element $K_{ij}$.

$-$ Let $\mathscr{P}_{h, x}$ and $\mathscr{P}_{h, y}$ be the standard $L^2$ projections in the $x$ and $y$ directions, respectively.

$-$ The projection ${\bf \Pi_h^+} = \mathscr{P}_{h, x}^+\otimes\mathscr{P}_{h, y}:\, [H^1(\Omega)]^2\rightarrow {\bf V}_h$ for vector-valued functions is defined as

where $\vec{{\bf n}}$ denotes the outward unit normal vector.

As shown in [31, Lemma 2.4], the projections mentioned above satisfy the following approximation properties:

where $\mathbb{Q}_h = \mathscr{P}_{h}^\pm$, $\Pi_h^-$, or ${\bf\Pi_h^+}$. Moreover, the projection $\Pi_h^-$ also has the following superconvergence property (see [31, Lemma 3.7]):

The "hat" term here is the numerical flux, which will be given later.

4.3.   Numerical Analysis

Rewrite (1.3) into the following equivalent first-order system:

Then the weak form of (4.16) at $t_n$ can be formulated as follows:

in which $v$ and ${\bf w}$ are test functions. The fully discrete nonuniform L1/LDG scheme is defined as follows: find $(u_h^n, {\bf p}_h^n)\in V_h\times{\bf V}_h$ such that

hold for any $(v_h, {\bf w}_h)\in V_h\times{\bf V}_h$. The alternating numerical fluxes are chosen, namely,

or

It is now time to present the stability and error estimate for the scheme (4.20) in the $L^2$-norm.

Theorem 4.1. (Stability) Assume that $u_h^n$ and ${\bf p}_h^n$ $(n = 1, 2, \cdots, M)$ are the LDG solutions of (4.20) with numerical flux (4.21). Then, it holds that

Proof. Taking $(v_h, {\bf w}_h) = (u_h^n, {\bf p}_h^n)$ in (4.20), and by summing over all $K$, one has

Adding them together and using integration by parts, one gets

which means that

Notice that

This, together with (4.24), yields

Therefore, utilizing Lemma 4.3 with $v^n = \|u_h^n\|$ and $\phi^n = \psi^n = 0$, one has

provided that the maximum time step $\tau_M\leq (4\Gamma(2-\alpha))^{-1/\alpha}$. The proof is completed.

Theorem 4.2. (Error estimate) Let $u({\bf x}, t_n)$ be the exact solution of problem (1.3), which satisfies that $\left|\delta^{l} u(\cdot, t)\right|\leq C\left(1+(\log t-\log a)^{\alpha-l}\right)$ for $l = 0, 1, 2$ and all $t\in(a, T]$. $u_h^n$ and ${\bf p}_h^n$ $(n = 1, 2, \cdots, M)$ are the LDG solutions of (4.20), with numerical flux given by (4.21). Suppose that $f(u)$ satisfies the condition (4.1). Then, there exists a positive constant $C$ independent of $M$ and $h$ such that

Proof. Let us first denote

Here, the projectors are selected as

Subtracting (4.20) from (4.19) and summing over all $K$ yield the following error equations:

Taking $(v_h, {\bf w}_h) = (Pe_u^n, \Pi e_{\bf p}^n)$ in (4.29), and by noticing the error decomposition (4.27), we obtain

where $\Upsilon^{n} = {_{CH}{\rm D}_{a, t}^{\alpha}}u^n- \Lambda_{\log}^{\alpha} u^n$. By virtue of (4.21) and the projection properties (4.12)–(4.15), we have

Then, from the Cauchy-Schwarz inequality and superconvergence property given by (4.17) that

On the other hand, for the nonlinear term in (4.32), we can obtain

where $\xi = \theta u^n+(1-\theta)Pu^n, \, \theta\in [0, 1]$. Employing the Cauchy-Schwarz inequality and interpolation property (4.16), we have

Notice that $f(u)-f(v) = f'(u)(u-v)-(u-v)^3+3u(u-v)^2.$ Hence, we derive from (4.1) that

Substituting (4.34)–(4.35) into (4.33) and applying (4.32), we have

This, combined with (4.25), further results in

As a consequence, according to Lemma 4.3 with $v^n = \|Pe_u^n\|$, $\phi^n = 2\|\Upsilon^n\|$, $\psi^n = \sqrt{2C}h^{k+1}$, $\lambda_0 = 2C$ and $\lambda_j = 0$ for $j = 1, 2, \cdots, M-1$, as long as $\tau_M\leq\left(4C\Gamma(2-\alpha)\right)^{-1/\alpha}$, we will obtain

Then, Lemma 4.2 leads to

By combining the above estimate with the triangle inequality, the desired result can be obtained.

5.   Numerical Experiments

The main purpose of this section is to give a numerical example to demonstrate the validity of the proposed scheme (4.20).

Example 5.1.

where $\Omega = (-1, 1)\times(-1, 1)$ and the source term $g(x, y, t)$ is chosen such that the exact solution of the problem is $u(x, y, t) = \big((\log t)^\alpha+(\log t)^2\big)(x+1)^2(x-1)^2(y+1)^2(y-1)^2$.

We apply the nonuniform L1/LDG scheme (4.20) to solve problem (5.1). Table 1 gives the $L^2$-errors and convergence orders versus $M$ for different values of $\alpha\, (\alpha = 0.4, 0.6, 0.8)$ and grading parameter $r \, (r = 1, \frac{2-\alpha}{2\alpha}, \frac{2-\alpha}{\alpha})$ when taking $t = 2$ and $M = N_x = N_y$, from which it is obvious that the convergence order in time is $\min\{2-\alpha, r\alpha\}$. To investigate the spatial convergence order, we utilize (4.20) to solve (5.1) by using both linear and quadratic finite element approximations, respectively. The $L^2$ errors and convergence order are listed in Table 2. The results show that the spatial convergence orders for the $L^2$-norm are close to $(k+1)$.

Comparisons between the numerical solution and the exact solution are depicted in Figures 1–3, and it can be seen that the numerical solution is in good agreement with the exact solution. The numerical solution surfaces for different times $t\, (t = 1.2, 1.4, 1.6, 1.8)$ and $\alpha \, (\alpha = 0.1, 0.5, 0.9)$ are shown in Figures 4–7. We can observe that the diffusion behavior of $u_h$ increases with time, and the maximum peak always appears in the center of the region. But if $\alpha$ is smaller, the diffusion process changes more slowly.

6.   Summary

The article first investigates the existence, uniqueness and regularity of solutions to (1.3). Then, a nonuniform L1/LDG scheme is constructed, and its stability and convergence are proven. Finally, the theoretical analysis is validated through numerical examples. In future work, we will focus on showcasing the physical properties of this numerical scheme and explore the implications of different definitions of $\alpha$-order fractional derivatives in the original problem. Additionally, we will examine which definition yields better results in terms of effectiveness.

Acknowledgments

The work was supported by the National Natural Science Foundation of China (No. 12101266).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.