A second–order approximation scheme for Caputo–Hadamard derivative and its application in fractional Allen–Cahn equation

1.   Introduction

Fractional calculus, an extension of classical calculus, offers a more precise and flexible tool for modeling complex phenomena by introducing fractional derivatives and integrals ^[1,2,3]. It demonstrates unique advantages in capturing memory effects and dealing with processes involving long-range dependencies, which are often difficult to accurately describe within the traditional framework of integer-order calculus. Hadamard's fractional calculus, proposed by Hadamard in 1892, shares similarities in form with the commonly used Riemann–Liouville integral/derivative and Caputo derivative, yet it has received relatively less attention. However, Hadamard's calculus can offer more accurate descriptions of phenomena such as Lomnitz's logarithmic creep law for specific substances ^[4] and ultra-slow diffusion processes ^[5]. For further applications of Hadamard's calculus, readers may refer to the literature ^{[6,7,8,9,10,11]}.

It is widely acknowledged that determining exact solutions for fractional differential equations, irrespective of the type of fractional derivative involved, presents a significant challenge. Consequently, there has been an increased emphasis on researching numerical methods to address these solutions. In the context of Caputo–Hadamard fractional differential equations, scholars have devised several effective numerical approaches for their resolution. For instance, Li et al. employed the finite difference method as discussed in ^[12], while the local discontinuous Galerkin finite element method was explored in ^[13]. Zhao et al. introduced a spectral collocation method in ^[14]. Ou et al. examined the regularity and logarithmic decay of solutions to the Caputo–Hadamard fractional diffusion-wave equation in ^[15] and resolved it using the difference method. They later introduced a fitted scheme in ^[16], which demonstrated a convergence order of $\min\{2r\alpha, 2\}$, where $r$ is the grading mesh parameter. Fan et al. proposed a series of approximation formulas for the Caputo–Hadamard fractional derivative in ^[5], including the L1-2 and L2-$1_\sigma$ formulas with the $(3-\alpha)$th order of convergence for $\alpha \in (0, 1)$, and the H2N2 formula also exhibiting the $(3-\alpha)$th order of convergence for $\alpha \in (1, 2)$. On the basis of the L1 and L2-$1_\sigma$ formulas approximating the Caputo derivative, Wang et al. ^[17,18] extended them to the discrete Caputo–Hadamard derivative, also known as the L1 and L2-$1_\sigma$ formulas in the logarithmic sense. It is worth mentioning that the main difference between the Caputo–Hadamard derivative and the Caputo derivative lies in its integral kernel, which is composed of a logarithmic function. Therefore, if a logarithmic transformation is applied to the variable in the Caputo–Hadamard derivative, it can be converted into the corresponding Caputo analog ^[19,20,21]. Consequently, the approximation schemes in the aforementioned literature were derived under this idea. To provide approximation formulas specifically for the Caputo–Hadamard derivative, Wang et al. recently proposed the general L1 formula ^[22] and the L2-$1_\sigma$ formula ^[23], which, although similar in form to those in the logarithmic sense, differ in the selection of grid points. Mustapha ^[24] first introduced the time-stepping L1 scheme for discretizing the Riemann–Liouville fractional derivative, achieving second-order accuracy. Liao et al. proposed the L1⁺ scheme for the Caputo derivative in ^[25], and subsequently established its discrete gradient structure in ^[26]. Since then, scholars have successfully applied this formula to the numerical solutions of many models and have gradually perfected its theoretical framework ^[27,28,29]. However, research on the Caputo–Hadamard fractional derivative has not yet been addressed. Therefore, the first task of this paper is to generalize it by constructing the L1⁺ formula for the Caputo–Hadamard fractional derivative.

The Allen–Cahn equation is a significant nonlinear partial differential equation that plays a crucial role in explaining phase transition phenomena and addressing interfacial diffusion issues. This equation was initially proposed by Allen and Cahn ^[30], and it has been extensively studied in the field of nonlinear dynamics, demonstrating its broad application value in various domains such as materials science ^[31] and grain growth ^[32]. In an effort to better understand and describe the anomalous diffusion transport behavior in heterogeneous porous materials, scholars have introduced the time-fractional Allen–Cahn equation. Liu et al. ^[33] considered a finite difference approximation for the time-fractional Allen–Cahn equation. Liao et al. developed a series of efficient numerical schemes for solving the time-fractional Allen–Cahn equation, such as the backward Euler and stabilized semi-implicit scheme ^[34], the Crank-Nicolson-type scheme with variable steps ^[35], the variable-step L1 scheme preserving a compatible energy law ^[36], and the nonuniform L2-$1_\sigma$ scheme ^[37]. Huang and Stynes ^[38] employed the Alikhanov L2-$1_\sigma$ scheme and a standard finite element method (FEM) to solve the time-fractional Allen–Cahn equation. Fan and Li numerically solved the Allen–Cahn equation with different time derivatives ^[39]. In this paper, we consider the following time-fractional Allen–Cahn equation:

where $\Omega\subset\mathbb{R}^2$ represents a bounded rectangular domain with the boundary $\partial\Omega$. The nonlinear term $f(u) = u^3-u$ with $u_a(\mathbf x)$ is the initial function, and $_{CH}{\rm{D}}^{\alpha}_{a, t}$ is the Caputo–Hadamard fractional derivative of order $\alpha$, which is defined as ^[40]

with $\omega_{\beta}(t) = \frac{t^{\beta-1}}{\Gamma(\beta)}$ and $\delta^n v(s) = \left(s\frac{\rm d}{{\rm d}s}\right)^nv(s) = \delta(\delta^{n-1}v(s))$, and $\forall n\in\mathbb{Z}^+$.

The nonconforming FEM boasts several advantages, including the ability to handle complex geometries and the flexibility of adaptive mesh refinement. It is well-known that the nonconforming Wilson element exhibits superior convergence properties compared with the conforming bilinear element. However, the Wilson element only converges on rectangular and parallelogram meshes. To enhance its applicability, modifications to the Wilson element have been attempted, as seen in ^[41], among others. Specifically, Jiang and Cheng introduced a quasi-Wilson element in ^[42], whose shape functions are independent of the element's geometric shape and converge on arbitrary quadrilateral meshes. Therefore, the second task of this paper is to apply the nonconforming quasi-Wilson FEM to the numerical solution of the time-fractional Allen–Cahn equation (1.1). It is widely recognized that solutions to time-fractional partial differential equations may exhibit a weak singularity at the initial time ^[43,44]. Additionally, it is crucial to note that the solution may display anisotropic behavior in space, meaning that the solution changes significantly only in certain directions while varying slowly in others. In such cases, using anisotropic meshes can reduce computational effort and clearly reflect the anisotropic characteristics of the solution. Consequently, we will employ an anisotropic nonconforming quasi-Wilson FEM for the spatial discretization of Eq. (1.1), and for the temporal direction, we will use the L1⁺ formula approximation. We will analyze the error in the $L^2$-norm, the superclose property in the $H^1$-norm, and the global superconvergence in the $H^1$-norm. Finally, we provide numerical examples to validate the theoretical results and simulate solutions to the Allen–Cahn equation with different time derivatives, observing the impact of the time derivative on the solution's behavior.

The structure of this paper is as follows. In Section Section 2, we propose a nonuniform L1⁺ formula for approximating the Caputo–Hadamard derivative and conduct a detailed analysis of its truncation error. Through a specific constructive example, we verify the convergence and numerical accuracy of this formula. In Section Subsection 3, we combine the spatially nonconforming finite element method with the L1⁺ formula on a nonuniform time mesh to obtain the fully–discrete scheme of Eq. (1.1). We analyze the $L^2$-norm error and the $H^1$-norm superclose estimate of this scheme, and further give the global superconvergence estimate in the $H^1$-norm. In Section Subsection 4, we provide some numerical examples to test the correctness of the theoretical analysis in Section Subsection 3. The last section is the conclusion.

2.   Nonuniform L1⁺ approximation for Caputo–Hadamard fractional derivative

Throughout the paper, $C$ denotes a generic constant that can take different values in different places but is independent of the mesh parameters. The notation $A\approx B$ means that a constant $C$ exists such that $A = CB$. For a given finite time $T$, let $\mathcal{I}_t = \left\{I_n = (t_{n-1}, t_n)\right\}_{n = 1}^N$ be a partition of $(a, T)$, where $a = t_0 < t_1 < \cdots < t_N = T$, $t_n = a+(T-a)(n/N)^r$, and the grading mesh parameter $r\geq 1$. From ^[22], we know that

2.1.   Nonuniform L1⁺ discretization

Let $\Pi_{\log}v(t)$ represent the linear interpolation function of $v(t)$ in the interval $[t', t'']$, i.e.,

Writing $v^n = v(t_n)$, we then present the L1⁺ discretization for the Caputo–Hadamard fractional derivative at time $t = t_n$,

where

For simplicity, we write

Set $\rho_n = \frac{\log t_{n-1}-\log t_{n-2}}{\log t_{n}-\log t_{n-1}}$ for $n = 2, 3, \cdots, N$. Our grid partition implies that $0 < \rho_n\le1$.

Lemma 1. Regarding $b^n_{n-k}, 1\le k\le n\le N$. We have the following properties:

(ⅰ) Assume $\alpha\ge (5-\sqrt{17})/2\approx0.43845$, if $\rho_n\ge \frac{2(1-\alpha)}{\alpha(3-\alpha)}$ for $1\leq n\leq N$, we have

(ⅱ) $b_{n-k}^n\ge \frac{1}{(2-\alpha)(\log t_k-\log t_{k-1})} {\int}_{t_{k-1}}^{t_k}\omega_{1-\alpha}(\log t_n-\log s)\frac{{\rm d}s}{s}$, $1\le k\le n\le N$.

Proof. Part (ⅰ) can be obtained along the same lines as in Section 5 of ^[29].

Next, we prove that Part (ⅱ) holds. According to the definition of $b_{n-k}^n\ (1\leq k\leq n)$, when $k = n$, the result can be obtained directly. When $1\le k\le n-1$, a $\xi_n\in(t_{n-1}, t_n)$ exists such that

The proof has been completed. □

2.2.   Truncation error analysis of L1⁺ scheme

In this subsection, we analyze the truncation error of the L1⁺ formula. For $t\in (a, T]$, suppose $u(t)$ satisfies the initial weak singularity condition

Assume that

For $1\le n\le N$ and $t\in I_n$, let $\psi(t) = v(t)-\Pi_{\log}v(t)$. We can then obtain the following from (2.8)

For $2\le n\le N$ and $t\in I_n$, according to the mean value theorem, $\xi_n\in I_n$ exists such that

For $n = 1$, since $v(t)$ satisfies the condition (2.7), $\delta v(t)$ may blow up as $t\rightarrow a^+$. Therefore, it needs to be estimated separately. Utilizing (2.4), we have

which further yields

For $1\leq n\leq N$, let

Our next task is to prove that

From (2.3), it follows that for $r(1+\sigma)\le 2$, we have

Therefore, to prove that (2.14) is true, it suffices to show that

Next, we will prove that (2.15) holds for different values of $n$ through several lemmas.

Lemma 2. For $r\ge 1$, when $n = 1, 2, 3$, we have $|R_1^n|\le C\Phi$.

Proof. When $n = 1$, by applying inequality (2.12), we can deduce that

For the case of $n = 2$, using $\log {t_2-\log t_1}\approx\log t_1-\log a$ along with the inequality (2.11), we can similarly infer that

When $n = 3$, the proof process is analogous to the case of $n = 2$, ultimately yielding

where $\xi_1\in(a, t_1)$ and $\xi_2\in (t_1, t_2)$. By combining the estimated results of these three scenarios and employing (2.15), the proof of this lemma can be completed. □

Lemma 3. For $r\ge 1$, when $n\ge 4$, $|R_1^n|\le C\Phi$ is true.

Proof. Firstly, we decompose $R_1^n$ into $R_1^n: = R_{11}^n+R_{12}^n+R_{13}^n$, where

with $n_0 = [n/2]$. We then prove one by one that for $i = 1, 2, 3$, $|R_{1i}^n|\le C\Phi$ holds.

Step 1: From (2.16), we can deduce that a $\xi_n\in(t_{n-1}, t_n)$ exists such that

where for the last inequality, we have utilized $\log t_{n-1}-\log a\approx\log t_{n}-\log a$, $(\log t_{n}-\log t_1)^{-\alpha-1} = \left(\log t_{n}-\log a-(\log t_1-\log a)\right)^{-\alpha-1}\le C(\log t_{n}-\log a)^{-\alpha-1}$, and $\log t_{n-1}-\log t_{n_0} = C(N^{-r}((n-1)^r-[n/2]^r))\ge C(\log t_n-\log a)$.

Since

we can obtain

Furthermore, according to $(\log t_n-\log a)^{1+\sigma} \approx (n/N)^{r(1+\sigma)}$ and $(\log t_n-\log a)^{-2/r}\approx n^{-2}N^{2}$, we can see that $|R_{11}^n|\leq C\Phi$.

Step 2: Using (2.1)–(2.3), one can derive $\log t_{n_0}-\log a\approx(\log t_n-\log a)$ for $n\ge1$, and $\log t_n-\log t_{n_0+1}\ge C(\log t_n-\log a)$ for $n\ge 4$. From (2.10), we obtain

which, combined with (2.15), implies that $|R_{12}^n|\leq C\Phi$.

Step 3: Let

and

From (2.18), we can derive $R_{13}^n = \sum_{k = n_0+1}^{n-1}r^k$. Applying the weak singularity condition (2.7), it follows that for $j = k, \, k+1$, a $\xi_j\in (t_{j-1}, t_j)$ exists such that

Combining (2.20)–(2.22), we deduce that

which further yields

Next, we aim to estimate $\sum_{k = n_0+1}^{n-1}\widetilde{r}^k$. To achieve this, we initially divide it into

where

To estimate $\sum_{k = n_0+1}^{n-1}\widetilde{r}_1^k$, we first perform a substitution for $s$ with $s = s + \log t_{k+1} - \log t_k$ in (2.25). Then, $\widetilde{r}_1^k$ can be rewritten as

where

Since $\widetilde{r}_{11}^{n-1} = 0$, we only need to estimate $\widetilde{r}_{11}^k$ for $n_0+1\le k\le n-2$. From (2.7), we have

Then we can bound the first term as follows:

For the second term $\widetilde{r}_{12}^k$, we can show that

Substituting (2.28) and (2.29) into (2.27), we know from (2.15) that

Finally, the estimation of $\sum_{k = n_0+1}^{n-1}\widetilde{r}_2^k$ remains unaddressed. Applying (2.19) and noticing that $\log t_{k+1}-\log t_{k}-(\log t_{k}-\log t_{k-1})\le CN^{-r}n^{r-2}\le CN^{-2}(\log t_n-\log a)^{1-2/r}$ yields

Combining (2.24), (2.30), and (2.31), we can obtain

which, together with (2.23) and the triangle inequality, leads to

By synthesizing the conclusions of Steps 1–3, we have completed the proof of the lemma. □

Theorem 1. When $r\ge 1$, for $n = 1, 2, \cdots, N$, it holds that

Proof. The conclusion of the theorem can be directly derived from Lemmas Lemma 2 and Lemma 3. □

Next, we present an example to validate the error and convergence order of the proposed nonuniform L1⁺ formula (2.5).

Example 1. We solve the following fractional differential equation numerically:

where the exact solution is $y(t) = (\log t)^{\alpha}+(\log t)^3$, and the function $f(t)$ can be directly computed.

It is clear that when $\sigma = \alpha$, $y(t)$ satisfies the initial weak singularity condition (2.7). Let $Y^n$ be an approximation of $y(t)$ at $t = t_n$, in which case, we have the following approximation scheme:

Analogous to the analysis process in Theorem 6.2 of ^[29], it can be shown that the maximum norm error

Figure 1 illustrates the maximum norm errors and numerical accuracy for various values of $\alpha$ when the grading mesh parameter $r = 1$. Clearly, the numerical accuracy is $\mathcal{O}(N^{-\alpha})$. When $r = \frac{2}{\alpha}$, the numerical results are shown in Figure 2, in which second-order convergence can be observed. These results align with the conclusions of Theorem 1.

3.   The fully–discrete numerical scheme for the problem (1.1)

In the spatial domain, let $(\cdot, \cdot)$ denote the $L^2$-inner product on $\Omega$, and let $\|\cdot\|$ represent the corresponding $L^2$-norm. For a non-negative integer $p$, we define $H^p(\Omega)$ as the Sobolev space with the norm $\|\cdot\|_p$ and the seminorm $|\cdot|_p$.

3.1.   Anisotropic nonconforming quasi-Wilson FEM

We utilize $\Gamma_h$ to stand for a rectangular partition of the domain $\Omega$, with $\overline{\Omega}$ being the union of all elements $e$ in $\Gamma_h$, i.e., $\overline{\Omega} = \cup_{e\in \Gamma_h} e$. For any element $e\in\Gamma_h$, let $O_e = (x_e, y_e)$ represent its center. The lengths of the two edges of $e$ in the $x$ and $y$ directions are denoted $2h_{x, e}$ and $2h_{y, e}$, respectively. Here, $\Gamma_h$ is a family of anisotropic grids without any regularity or quasi-uniform assumptions. In other words, $\Gamma_h$ is not required to satisfy $\frac{h_{x, e}}{h_{y, e}}\leq C$ or $\frac{h_e}{h_{y, e}}\leq C$, where $h_e = \max\{h_{x, e}, h_{y, e}\}$, and $h = \max_{e\in\Gamma_h}\{h_e\}$.

Assume that $\hat{e} = [-1, -1]^2$ is the reference element on the $\hat{x}-\hat{y}$ plane. Its four vertices are $\hat{a}_1 = (-1, -1)$, $\hat{a}_2 = (1, -1)$, $\hat{a}_3 = (1, 1)$, and $\hat{a}_4 = (-1, 1)$. Similar to the work in ^[45,46], we define the nonconforming quasi-Wilson finite element $\{\hat{e}, \hat{P}, \hat{\Sigma}\}$ as follows:

where $N_i(\hat{x}, \hat{y}) = \frac{1}{4}(1+\hat{x}_i\hat{x})(1+\hat{y}_i\hat{y})$, with $(\hat{x}_1, \hat{x}_2, \hat{x}_3, \hat{x}_4) = (-1, 1, 1, -1)$ and $(\hat{y}_1, \hat{y}_2, \hat{y}_3, \hat{y}_4) = (-1, -1, 1, 1)$; $\hat{\Psi}(s) = \frac{1}{2}(s^2-1)-\frac{5}{12}(s^4-1)$. Also, $\hat{v}_i = \hat{v}(\hat{A}_i)\, (i = 1, 2, 3, 4)$, $\hat{v}_5 = \frac{1}{|\hat{e}|} {\int}_{\hat{e}}\frac{\partial^2 \hat{v}}{\partial \hat{x}^2}{\rm d}\hat{x}{\rm d}\hat{y}$, and $\hat{v}_6 = \frac{1}{|\hat{e}|} {\int}_{\hat{e}}\frac{\partial^2 \hat{v}}{\partial \hat{y}^2}{\rm d}\hat{x}{\rm d}\hat{y}$. For any $\hat{v}\in\hat{P}$, the associated interpolation can be uniquely expressed as $\hat{I}_h\hat{v}(\hat{x}, \hat{y}) = \sum_{i = 1}^4N_i(\hat{x}, \hat{y})\hat{v}_i+\hat{\Psi}(\hat{x})\hat{v}_5+\hat{\Psi}(\hat{y})\hat{v}_6$. The affine mapping $F_e:\hat{e}\rightarrow e$ is defined by $x = x_e+h_{x, e}\hat{x}$ and $y = y_e+h_{y, e}\hat{y}$.

We define the finite element space by

For any given $v_h = \bar{v}_h+\tilde{v}_h$, with $\bar{v}_h$ and $\tilde{v}_h$ denoting the conforming and nonconforming components of $V_h$ respectively. As can be obtained from ^[47], we have

The Ritz projection, denoted $R_h: H^1_0(\Omega)\rightarrow V_h$, is defined in such a way that for any $v_h\in V_h$, the equation

holds. Let $I_h: H^1(\Omega)\rightarrow V_h$ be an interpolation operator satisfying $I_h|_e = I_e$ and $I_ev = \hat{I}_h\hat{v}\circ F_e^{-1}$. The $L^2$-projection operator $P_h: L^2(\Omega)\rightarrow V_h$ and the discrete Laplacian operator $\Delta_h: V_h\rightarrow V_h$ are defined as

respectively. According to ^[48], it is known that

where $\vec{n}$ is the unit outward normal vector on $\partial e$.

Lemma 4 (^[49]). For $v \in H_0^1(\Omega) \cap H^2(\Omega)$, it is true that

If $v \in H_0^1(\Omega) \cap H^3(\Omega)$, one can derive

If $v \in H_0^1(\Omega) \cap H^4(\Omega)$, one has

Integrating the first equation of (1.1) over the time interval $I_n$, we obtain

For simplicity, let $\overline{g^n} = \frac{1}{\log t_n-\log t_{n-1}} {\int}_{t_{n-1}}^{t_n}g(t)\frac{{\rm d}t}{t}$. Then the equation above can be rewritten as

where $u^n = u(\mathbf{x}, t_n)$ for $0\leq n\leq N$. For the nonlinear term $({u^n})^3$, we employ the Newton linearization method (see, e.g., ^[38]), namely $({u^n})^3 \approx (u^{n-1})^3 + 3(u^{n-1})^2(u^n - u^{n-1})$, to derive the time semi-discrete scheme

where $H(u^n) = (u^{n-1})^3+\frac{3}{2}(u^{n-1})^2(u^n-u^{n-1})$, and $u^{n-\frac{1}{2}} = (u^n+u^{n-1})/2$.

On the basis of (3.1) and the equation above, we can now describe the fully–discrete numerical scheme to find $\{U_h^n\}_{n = 0}^N\in V_h$ such that

where $U_h^n = U_h(t_n)$ and $U_h^{n-1/2} = \frac{U_h^n+U_h^{n-1}}{2}$.

3.2.   Global superconvergence analysis

In this subsection, we discuss the convergence of the scheme (3.7) under the $L^2$-norm and the superclose convergence under the $H^1$-norm. Subsequently, on the basis of the interpolation postprocessing technique defined in ^[50], we analyze the global superconvergence of this scheme under the $H^1$-norm. Suppose that the solution $u$ in (1.1) satisfies the following regularity condition:

Lemma 5. Assume

Under the condition (3.8), a constant $C$ exists such that

Proof. From ^[29], it follows that $H(u^n)$ provides a second-order approximation for the nonlinear term $\frac{1}{\log t_n - \log t_{n-1}} {\int}_{t_{n-1}}^{t_n} \left(u(t) - u^3(t) \right) \frac{\mathrm{d}t}{t}.$

On the basis of (2.10)–(2.12), we readily derive

Combining this result with Lemmas 3.1 and 3.2 in ^[23], we obtain the desired conclusion. □

Theorem 2. Suppose $u$ and $\{U_h^n\}_{n = 1}^N$ are the solutions of equations (1.1) and (3.7), respectively. Assume that $u$ satisfies the conditions in (3.8) and $_{CH}{\rm D}_{a, t}^{\alpha} u \in L^\infty(a, T;H^2(\Omega))$. Then the following estimates hold

Proof. To begin, let $u^n-U_h^n = R_hu^n-U_h^n+u^n-R_hu^n: = \xi^n+\eta^n.$ By applying (3.5)-(3.7), we derive that

Notice that

From (3.9), we then get

where $\Phi^{n} = \Phi_1^{n}(\xi^{n-1}+\eta^{n-1})+\Phi_2^{n}(\xi^{n}+\eta^{n})$, $\Phi_1^{n} = \big((u^{n-1})^2+(U_h^{n-1})^2+u^{n-1}U_h^{n-1}+\frac{3}{2}(u^{n}-u^{n-1})(u^{n-1}+U_h^{n-1})-\frac{3}{2}(U_h^{n-1})^2$, and $\Phi_2^{n} = \frac{3}{2}(U_h^{n-1})^2$.

Taking $v_h = 2\xi^{n-\frac{1}{2}}$ in (3.10), we obtain

Applying Lemma A.1 in Appendix, one can get

Substituting (3.12) into (3.11) and noting that

we have

Utilizing the Sobolev embedding theorem, we find that a constant $C$ exists such that $\|u^n\|_{L^\infty(\Omega)} \leq C$. Following a similar argument to that in [51, Theorem 3.2], it can be shown that $U_h^n$ for $1 \leq n \leq N$ satisfies $\|U_h^n\| \leq C$. Using (3.2) together with the definition of $\eta^n$ and an application of the Cauchy–Schwarz inequality yields

In view of (3.5), we have

where $u^{n-\frac{1}{2}} = \frac{u^n+u^{n-1}}{2}$.

Through the three estimates above and again using the Cauchy–Schwarz inequality, we obtain

By applying Lemma 5, the following holds:

Substituting (3.17) into (3.16) indicates that

Using the discrete Gronwall inequality (see Theorem B.1) and (A.4), we obtain

Thus, the first assertion of the theorem follows by applying (3.3) and the triangle inequality.

Taking the test function $v_h = -2\Delta_h\xi^{n-\frac{1}{2}}$ in (3.10) and using (3.12), we have

By employing (3.2) and leveraging the projection property (3.3), we derive

and

Thus, it follows from (3.20)–(3.23) that

By repeating the discussion process of (3.19), we can conclude that

This, together with (3.4), leads to the following result:

The second assertion in the theorem has been proven completed. □

In order to examine the superconvergence of the fully–discrete scheme (3.7), a new mesh family $\hat{\Gamma}_h = \{\hat{e}\}$ is introduced. Each new element $\hat{e}$ in $\hat{\Gamma}_h$ is formed by combining four adjacent elements from the original mesh $\Gamma_h$. Let $a_i \, (i = 1, \cdots, 9)$ be the nine vertices of the four small elements. For any element $\hat{e}$ in $\hat{\Gamma}_h$, an interpolation postprocessing operator $\hat{I}_h$ is defined as stated in ^[50]:

where $Q_{22}(\hat{e}) = \text{span}\{1, x, y, xy, x^2y, xy^2, x^2, y^2, x^2y^2\}$.

Theorem 3. Assuming that the conditions of Theorem 2 are satisfied, for $1\le n\le N$, we have

Proof. According to the conclusion of Theorem 2 and using (3.26), we deduce that

This completes the proof. □

Remark 1. Although, theoretically, in order to achieve second-order accuracy in time, the scheme (3.7) requires the fractional derivative order $\alpha$ to satisfy the condition $0.43845\approx(5-\sqrt{17})/2 < \alpha < 1$ (see Theorems 2 and 3), our numerical experiments have shown that this condition is not absolutely necessary. In fact, for all $0 < \alpha < 1$, as long as $r\ge 2/\sigma$, second-order convergence can be achieved. For details, see Example Example 2.

4.   Numerical examples

In this section, we aim to validate the effectiveness and accuracy of the scheme (3.7) for solving the Allen–Cahn equation with a Caputo–Hadamard time derivative through numerical experiments. Furthermore, we conduct numerical simulations of the Allen–Cahn equation using various time derivatives (including the integer order derivative, the Caputo derivative, and the Caputo–Hadamard derivative) to examine the influence of the time derivative on the solution's behavior.

Example 2. We study the following Allen–Cahn equation:

where $\Omega = (0, 1)\times(0, 1)$. The function $g(\mathbf{x}, t)$ is selected to fulfill the exact solution $u(\mathbf{x}, t) = (\log t)^{\alpha} (1 - x)(1 - e^{-x})y(1 - y)$.

It is easy to observe that for Eq. (4.1), the solution exhibits initial weak singularity in time and anisotropic behavior in space. We take the regularization parameter $\sigma = \alpha$ from (2.7). According to Theorem 2 and 3, to achieve optimal second-order accuracy in time, the grading mesh parameter $r$ should satisfy $r \geq \frac{2}{\alpha}$. In Table 1, we present the errors and temporal convergence rates for different values of $\alpha$ when $r = \frac{2}{\alpha}$, demonstrating that the proposed scheme achieves second-order convergence. Table 2 shows the numerical results in the spatial direction, which are also consistent with the theoretical analysis.

Example 3. (Mean curvature flow problem ^[52]) We study the following equation with Neumann boundary conditions:

where $\Omega = (0, 1)\times (0, 1)$ and $\epsilon = 0.01$. If $\alpha = 1$, the operator $\partial^\alpha_tu(\mathbf x, t) = \frac{\partial u(\mathbf x, t)}{\partial t}$. For $\alpha\in(0, 1)$, $\partial^\alpha_tu(\mathbf x, t)$ can be defined using either the Caputo fractional derivative or the Caputo–Hadamard fractional derivative.

For the case of $\alpha = 1$, Eq. (4.2) represents the classical Allen–Cahn equation. Figure 3 depicts the changing process of the numerical solution. Clearly, the circle is shrinking as time passes, which is in line with the numerical calculation results in Choi et al. ^[52].

If the time derivative is the Caputo fractional derivative and $\alpha = 0.4$ and $0.8$, the numerical results at $t = 0$, $2$, and $10$ are shown in Figure 4. As time increases, the circle shrinks. Moreover, the larger the value of $\alpha$, the faster it shrinks. However, it shrinks more slowly than the integer-order model.

If the time derivative is the Caputo–Hadamard fractional derivative and $\alpha = 0.4$ and $0.8$, the numerical results at $t = 1$, $100$, and $300$ are shown in Figure 5. The circle also shrinks as time goes on. The larger the value of $\alpha$, the faster it shrinks. But it shrinks more slowly than the Caputo time derivative model.

In conclusion, the shrinking speeds under different derivatives are as follows: integer-order Allen–Cahn equation $>$ Caputo fractional Allen–Cahn equation $>$ Caputo–Hadamard fractional Allen–Cahn equation. In the case of the same fractional derivative, the larger the fractional derivative $\alpha$, the faster the shrinking. This is consistent with the conclusion obtained by Fan and Li ^[39].

5.   Concluding remarks

The primary contributions of this paper are as follows.

$(1)$ We introduce a high-order approximation scheme for approximating the Caputo–Hadamard fractional derivative, known as the nonuniform L1⁺ formula. When the grading mesh parameter satisfies certain conditions, this method can achieve second-order accuracy.

$(2)$ We explore a numerical algorithm for the Allen–Cahn equation with time as the Caputo–Hadamard fractional derivative. The algorithm employs the nonuniform L1⁺ formula for temporal discretization and uses the nonconforming quasi-Wilson FEM for spatial approximation. We analyze the error of the algorithm in the $L^2$-norm and the superclose error in the $H^1$-norm. Furthermore, by applying interpolation postprocessing techniques, we demonstrate global superconvergence results in the $H^1$-norm.

$(3)$ Numerical experiments validate the theoretical convergence rates. Furthermore, we numerically solve the Allen–Cahn equation with different types of time derivatives, including integer-order derivatives, Caputo derivatives, and Caputo–Hadamard derivatives, to observe their diffusion processes. The results indicate that the relationship of diffusion speed is as follows: integer-order Allen–Cahn equation $>$ Caputo fractional Allen–Cahn equation $>$ Caputo–Hadamard fractional Allen–Cahn equation.

Author contributions

Luhan Sun: Writing-original draft; Zhen Wang: Writing-review & editing, Methodology; Yabing Wei: Validation.

Use of AI tools declaration

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

Acknowledgments

The research of Zhen Wang is supported by the National Natural Science Foundation of China (No. 12101266). The research of Yabing Wei is supported by the Natural Science Foundation of Jiangsu Province (No. BK20240832).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Properties of the discrete coefficients in (2.5)

Lemma A.1. For $\alpha \in \left(\frac{5-\sqrt{17}}{2}, 1\right)$, the operator ${_{CH}{\delta}}_{a, t}^{\alpha}$ defined in (2.6) satisfies

Proof. Following the proof idea of Corollary 1 in ^[53], the corresponding result can be obtained. □

Lemma A.2. Define the complementary discrete kernels $Q_{n-k}^n$ ($1 \le k \le n$) as follows:

For $\alpha \in \left(\frac{5-\sqrt{17}}{2}, 1\right)$, the following properties hold:

(ⅰ)

and

(ⅱ) Consider a function $v: [a, T] \rightarrow \mathbb{R}$ that is continuous and piecewise $C^1$. If $\delta u$ is monotonically decreasing and non-negative, then the following holds:

(ⅲ) For any real number $\mu > 0$, we have

Here, $E_{\alpha} = \sum_{k = 0}^{\infty}\frac{z^k}{\Gamma(k\alpha+1)}$ represents the Mittag–Leffler function.

Proof. (ⅰ) Applying Lemma 1 (ⅱ), (A.2) is proven. Since $Q_{n-k}^n b_0^k \leq \sum_{j = k}^{n} Q_{n-j}^n b_{j-k}^j = 1$, then $Q_{n-k}^n \leq \frac{1}{b_0^k} \leq \Gamma(3-\alpha) (\log t_n - \log a)^\alpha$. Because $\omega_{1-\alpha}(\log t_j - \log a)$ is monotonically decreasing, we have $\sum_{k = 1}^{n} Q_{n-k}^n \omega_{1-\alpha}(\log t_n - \log a) \leq \sum_{k = 1}^{n} Q_{n-k}^n \omega_{1-\alpha}(\log t_k - \log a) \leq (2-\alpha)\sum_{k = 1}^{n} Q_{n-k}^n b_{k-1}^k = 2-\alpha$, yielding (A.3). Similar to the proof of Lemma 3.5 in ^[18], it can be seen that (A.4) holds.

(ⅱ) Using Lemma 1 (ⅱ) and Lemma A.2, we obtain

Thus,

(ⅲ) The proof process closely mirrors that of Lemma 3.4 in ^[23], which is omitted here for brevity. □

B.   Gronwall inequality

Theorem B.1. Let $\{\lambda_k\}_{k = 0}^{N-1}$ be a non-negative sequence, and a constant $\Lambda$ exists such that $\sum_{k = 0}^{N-1}\lambda_k \leq \Lambda$. Let $\{\xi^n, \eta^n : 1\leq n\leq N\}$ be non-negative sequences, and the grid function $\{v^n: n\geq 0\}$ satisfies the condition

If $\log t_N - \log t_{N - 1} \leq (2(2-\alpha)\Gamma(3-\alpha)\Lambda)^{-1/\alpha}$, then the following conclusion holds

Proof. The proof of the lemma follows a similar approach to that outlined in Lemma 3.6 of ^[18] and is thus omitted here. □