$\alpha$-robust error analysis of two nonuniform schemes for Caputo-Hadamard fractional reaction sub-diffusion problems

1.   Introduction

In the past, researchers have conducted extensive studies in the fields of fractional calculus and fractional differential equations (FDEs)^[1,2,3,4]. To date, numerous scholars have extensively explored the realm of fractional integrals and derivatives in numerous forms, including Riemann-Liouville, Caputo, and Riesz integrals and derivatives, among others. However, there exists another type of fractional derivative that incorporates a logarithmic function in its definition, known as the Hadamard fractional derivative, which is defined as ^[2]

where $0 < \alpha < 1, \omega_{\beta}(t) = \frac{t^{\beta-1}}{\Gamma{(\beta)}}, \ \delta v(x, s) = \left(s\frac{\partial }{\partial s}\right)v(x, s),$ and $\Gamma(\cdot)$ representing the Gamma function.

Compared to the Riemann-Liouville and Caputo derivatives, the Caputo-Hadamard derivative, first introduced in 1892, more accurately captures some complex processes in practical applications. This includes Lomnitz logarithmic creep law ^[5,6] and ultra-slow mechanics ^[7,8], among others. In particular, the logarithmic increase in the mean square displacement of particles during ultra-slow diffusion has been demonstrated in ^[9,10,11]. As a result, the Hadamard fractional operator, whose kernel is a logarithmic function, has emerged as a natural model for ultra-slow diffusion processes and has garnered significant attention.

Let $\Omega = (x_l, x_r), \Lambda = (a, T)$ with $a > 0$, and we will pay attention to the numerical approximation for the following Caputo-Hadamard fractional reaction sub-diffusion equation:

Here, the real constant $\kappa\in \mathbb{R}$ is the reaction coefficient, and the source term $f(x, t)$ and the initial data $\varphi(x)$ are given functions.

Recent studies have explored various numerical methods for tackling Caputo-Hadamard FDEs, encompassing the L1 scheme ^[12], L1-2 scheme, and L2-1$\sigma$ scheme ^[13]. Moreover, Zhao et al. ^[14] introduced a spectral collocation method utilizing mapped Jacobi log orthogonal functions as basis functions, resulting in an efficient algorithm for solving Hadamard-type FDEs. Based on block-by-block approach, Ye et al. ^[15] proposed and analyzed a high order time stepping scheme having the convergence order more than three for the Caputo-Hadamard fractional differential equations.

Actually, the aforementioned studies mainly concentrate on the Caputo-Hadamard FDEs that possess smooth solutions. Recently, numerous efficient numerical methods have been developed for Caputo FDEs with weakly singular solutions, including the nonuniform L1 scheme ^[16,17], the nonuniform Alikhanov scheme ^[18,19,20], convolution quadrature method ^[21], and the spectral method ^[22]. Interested readers can also consult some recent references ^{[23,24,25,26,27]} for more numerical methods about FDE, such as the Alternating Direction Implicit method, extrapolation method, meshless method, and so on.

Nevertheless, numerical simulations for Caputo-Hadamard fractional reaction sub-diffusion equation (1.1)–(1.3) with weakly singular solutions remain relatively limited. For Eqs (1.1)–(1.3) without reaction term, Li et al. ^[28] proposed an L1 scheme on nonuniform meshes to approximate the time Caputo-Hadamard fractional derivative and employed the local discontinuous Galerkin method to approximate the spatial derivative. Later, the Alikhanov scheme with nonuniform time meshes for Caputo-Hadamard fractional sub-diffusion equations with an initial singularity was investigated in ^[29]. The stability and convergence of the resulting discrete scheme were analyzed, but the error bounds generally contain a constant factor $\Gamma(1-\alpha)$ or $1/(1-\alpha)$ which will blow up as $\alpha$ approaches $1^-$.

Zhang et al. ^[30] derived a novel $\alpha$-robust error analysis for convolution-type schemes with general nonuniform time step for Caputo fractional reaction sub-diffusion equations. By virtue of the ideas derived in ^[30], this paper will extend the nonuniform L1 and Alikhanov scheme presented in ^[28] and ^[29] to Caputo-Hadamard fractional sub-diffusion equations with reaction term, and then we will consider the $\alpha$-robust error analysis of the proposed schemes. This means that the derived error bounds will not contain any blowup factor and will remain valid as $\alpha\rightarrow 1^-$.

Throughout this paper, we employ $C$ to represent a generic constant that is independent of the mesh, and it may take different values at different places. Additionally, $C$ exhibits $\alpha$-robustness, meaning it is influenced by $\alpha$, yet as $\alpha$ approaches 1, the value of $C$ remains finite, avoiding any potential explosion. The proposed method in this paper is analyzed under the following regularity assumptions: for all $(x, t)\in \Omega\times \Lambda$, it holds

with the regularity parameter $\sigma\in(0, 1)\cup(1, 2)$, and $\delta^{l}u(x, t) = (t\frac{\partial}{\partial t})^l u(x, t)$.

The remainder of this paper is organized as follows. In Section 2, we describe the detailed construction of the general convolution-type scheme. After that, the abstract result for graded mesh is applied to two typical numerical schemes, i.e., the widely used L1 scheme and Alikhanov's scheme. In Sections 3 and 4, we give a rigorous analysis of the stability and convergence of the L1 scheme and Alikhanov's scheme, and derive $\alpha$-robust error estimates under specific regularity conditions imposed on the exact solution. In Section 5, some numerical examples are provided to support the theoretical statement. Some concluding remarks are given in the final section.

2.   Preliminaries

To develop a finite difference scheme for solving (1.1)–(1.3), we first divide the spatial interval $[x_l, x_r]$ into $M$ subintervals with grid size $h = \frac{x_r-x_l}{M}$. Set discrete grid $\Omega_h = \{x_i|0 \leq i \leq M\}$ with $x_i = x_l+ih.$ For any grid function $w = \{w_i|w_i = w(x_i), \ x_i\in \Omega_h\}$, let $\delta_x^2 w_i$ be the standard second-order approximation of $\partial_x^2 w(x_i)$, i.e.,

Denote the space of grid functions $\mathcal{W} = \{w|w_0 = w_M = 0\}.$ For any $w, v\in \mathcal{W}$, the discrete $L^2$ inner product and the associate $L^2$ norm are given as

Denote

for any grid functions $w, v\in \mathcal{W}$, and it holds

We now proceed to the discretization of time. First, we partition the interval $[a, T]$ arbitrarily with $a = t_0 < t_1 < \cdots < t_{k-1} < t_k < \cdots < t_N = T,$ and set

Based on this partition, we derive L1 and Alikhanov's scheme, which are given in (3.1) and (4.1), respectively. Second, we further study and discuss these two formulas in the following divisions on the interval $[a, T]$:

and, correspondingly,

Define $v^n: = v(x, t_n), t_{n-\theta}: = \theta t_{n-1}+(1-\theta)t_{n}$, and $v^{n-\theta}: = \theta v^{n-1}+(1-\theta)v^{n}$ with an offset parameter $\theta\in [0, 1)$. Then, the Caputo-Hadamard derivative operator in the problem (1.1) can be approximated as a convolution as detailed in the succeeding article ^[28]:

where the difference operator is $\nabla_{\tau}v^{k} = v^{k}-v^{k-1}$ for $k\geq 1$. To conduct our error analysis, we require the following three assumptions.

A1. The discrete kernels are monotone, meaning that

A2. There exists a constant $\pi_{A} > 0$ such that the discrete kernels satisfy a lower bound

A3. There exists a constant $\rho > 0$ for which the step size ratios $\rho_k$ satisfy

For the graded mesh, it can be checked that the step size ratio is $\rho_{k} < 1$. In fact, by using the mean value theorem, we have

with $\eta_1\in (k-1, k), \; \; \eta_2\in(k, k+1).$

Next, to derive the global consistency error, we introduce the discrete complementary convolution (DCC) kernels ^[16]

which are specifically chosen to enforce the identity

Futher, a discrete fractional Gronwall inequality is given as follows.

Lemma 2.1. Assume that A1-A3 hold, $\theta\in[0, 1)$, and the nonnegative sequences $(g^{n})_{n = 1}^{N},$ $(\lambda^{l})_{l = 0}^{N-1}, (v^{k})_{k = 0}^{N}$ satisfy

If a constant $\Pi$ satisfies $\Pi\geq\sum_{l = 0}^{N-1}\lambda_{l}$ and the maximum step size satisfies

then it holds that

where $E_{\alpha}(\cdot) = \sum\limits_{j = 0}^{\infty}\frac{(\cdot)^{j}}{\Gamma{(j\alpha+1)}}$ represents the special Mittag-Leffler function.

The proof of this lemma is similar to that in ^[19].

Lemma 2.2. ^[31] For $\nu > 0, \beta > 0, b > a > 0,$ there holds

Lemma 2.3. ^[19] If $g$ is monotone increasing and $h$ is monotone decreasing on the interval $[a, b]$, and if both functions are integrable, then

In the following, we derive two important lemmas which are useful in the convergence analysis later on.

Lemma 2.4. Assume that A1-A2 hold, then we have

Proof. Denote $h(t) = \omega_{1+\alpha}(\log t-\log a)$, then it holds $\delta h(t) = \omega_{\alpha}(\log t-\log a).$ In one side, it follows from the definition of the Caputo-Hadamard derivative that

It is easy to check that $\omega_{1-\alpha}(\log t_{j}-\tau)$ is monotone increasing and $\omega_{\alpha}(\tau-\log a)$ is monotone decreasing on the interval $[\log t_{k-1}, \log t_{k}]$. Thus, by Lemma 2.3 and A2, we obtain

In another side, by the definition of the Caputo-Hadamard derivative and Lemma 2.2, we get again that

Using (2.9), (2.10), and (2.5), we have

□

Lemma 2.5. Assume that A1–A2 hold, and for any positive sequence $(\upsilon^{k})_{k = 1}^{n}$, one has

Proof. It follows from A2 that

Thus, using (2.12) and (2.5), we have

□

3.   $\alpha$-robust error analysis of the L1 scheme

We are now in a position to consider the $\alpha$-robust error estimate of the L1 scheme. Following ^[28], the L1 approximation to the Caputo-Hadamard derivative $^{CH}D_{a, t}^{\alpha} v^{n}$ is given by

where the discrete coefficients $A_{n-k}^{(n)}$ are defined by

It is easy to see that Assumption A2 holds with $\pi_{A} = 1$. Using the integral mean-value theorem, one has

with $\xi_{k+1}\in(t_{k}, t_{k+1}), \; \; \xi_{k}\in(t_{k-1}, t_{k})$. Thus Assumption A1 holds.

Let $u_i^n$ be the discrete approximation of solution $u(x_{i}, t_{n})$ for $x_i\in \Omega_h, 0\leq n\leq N.$ The fully discrete scheme for problem (1.1)–(1.3) is given as

where $f_{i}^{n} = f(x_{i}, t_{n})$.

We aim to demonstrate the stability and convergence of the scheme (3.4)–(3.6). The stability is established in the following theorem.

Theorem 3.1. Assume that the assumptions (1.4), (1.5), and A3 hold. Let $\kappa_{+}: = \max\{-\kappa, 0\}$, and $v_i^n$ be the solutions of the following difference equation:

If the maximum step size is $\tau\leq \mathfrak{B}_1$ with

we have

Proof. After taking the inner product with $2v^{n}$ on both sides of (3.7), we obtain

By employing the definition of the discrete Caputo-Hadamard derivative given in (3.1) and utilizing the monotone property given in (3.3), it can be inferred that

It follows that

Noting (2.1), and using the Schwarz inequality, we obtain

By applying Lemma 2.1 to the above inequality, we obtain the desired estimate (3.11), thereby completing the proof. □

We now consider the convergence of the scheme. To this end, let $u(x, t)$ be the exact solution of (1.1)–(1.3), and let $\{u_{i}^{n}|0\leq i\leq M, 0\leq n\leq N\}$ be the solution of problem (3.4)–(3.6). Set

It is straightforward to obtain the following error equation:

where

After conducting stability analysis in Theorem 3.1, if the maximum step size satisfies (3.10), one has, with $\bar{C} = 2E_{\alpha}(4\kappa_{+}\max\{1, \rho\} (\log t_{n})^{\alpha})$,

Under the assumption of spatial regularity given by (1.4), the Taylor expansion provides a direct demonstration that $\|(R_{s})^{n}\|\leq Ch^2.$ When combined with Lemma 2.4, we can further deduce that

Now, we only need to estimate the term $\sum_{j = 1}^{n}P_{n-j}^{(n)}\|(R_{t})^{n}\|$, which will be achieved through the following lemmas.

Lemma 3.1. [29, Lemma 3.1] Suppose that $f(x)$ has a continuous $\delta$-derivative of $n + 1$ order in some field of point $x_0$. A Taylor-like formula with integral remainder is given by

Lemma 3.2. For any function $u\in C^3((a, T])$, the local truncation errors $(R_{t})_{i}^{k}$ satisfy

where

Lemma 3.3. Assume that A3 holds, and $u(x, t)$ satisfies the regularity assumption (1.5), then

where $C = (1+\rho)\Gamma{(2-\alpha)}$.

The proofs of Lemmas 3.2 and 3.3 are left to Appendix 7.1 and 7.2 for brevity.

Theorem 3.2. Let $u(x, t)$ be the exact solution of (1.1)–(1.3), and let $\{u_{i}^{n}|0\leq i\leq M, 0\leq n\leq N\}$ be the solution of (3.4)–(3.6). Suppose the regularity assumptions (1.4)–(1.5) hold. Set

If the maximum step size is $\tau\leq \mathfrak{B}_1$, then the discrete solution is convergent in the $L^2$-norm with $\bar{C} = 2E_{\alpha}(4\kappa_{+}\max\{1, \rho\}\left(\log t_n\right)^{\alpha})$ such that

where

In particular, if graded mesh is used, then it holds

with

Proof. By combining (3.18), (3.19) and Lemma 3.3, we arrive at (3.22). We now proceed to examine the global approximation error on the graded mesh. It is easy to verify $\tau_{k}\leq r \log\frac{T}{a}N^{-r}k^{r-1}$ as follows:

Note that $\sigma > 0$, it follows from (3.23) and (3.25) that

Taking into account that

(3.26) can be simplified to

For $\sigma r\geq2-\alpha$, we have from (3.27) that

For $\sigma r < 2-\alpha$, it follows from (3.27) that

Note that

then we get from (3.29) that

By utilizing (3.28), (3.30), and (3.23), we achieve the desired global approximation error on the graded mesh. □

4.   $\alpha$-robust error analysis of the Alikhanov's scheme

In this section, we delve into the $\alpha$-robust error estimate pertaining to the Alikhanov scheme. Referring to ^[28], the Alikhanov scheme for the Caputo-Hadamard derivative is expressed as follows:

with $\theta = \frac{\alpha}{2}$. Here the discrete convolution kernel $A_{n-k}^{(n)}$ is defined as: $A_{0}^{(n)} = a_{0}^{(n)}$ if $n = 1$ and, for $n\geq2$,

Moreover, the discrete coefficients $a_{n-k}^{(n)}$ and $b_{n-k}^{(n)}$ are determined by:

where

As proved in ^[29], one can verify that Assumptions A1 and A2 hold with $\pi_{A} = 11/4$.

Next, the difference scheme is established. Considering (1.1) at mesh point $(x_{i}, t_{n-\theta})$, the fully discrete scheme for problem (1.1) is given as

We proceed to assess the stability and convergence of Alikhanov's scheme.

Theorem 4.1. Assume that the assumptions (1.5) and A3 hold. Let $\kappa_{+}: = \max\{-\kappa, 0\}$ and $v_i^n$ be the solution of the following differential equation:

If the maximum step size is $\tau\leq \mathfrak{B}_2$ with

we have

Proof. Taking the inner products with $2v_{i}^{n-\theta}$ on both sides of (4.5), we obtain

Utilizing (3.12), we can further derive

Taking into account (2.1) and applying the Schwarz inequality, we arrive at

By applying Lemma 2.1 to the aforementioned inequality, we successfully derive the desired estimate (4.9), thus, conclusively completing the proof. □

We now consider the convergence of the scheme. To this end, let $u(x, t)$ be the exact solution of (1.1)–(1.3), and let $\{u_{i}^{n-\theta}|0\leq i\leq M, 0\leq n\leq N\}$ be the solution of problem (4.2)–(4.4). Set

It is straightforward to obtain the following error equation:

where $R_{i}^{n} = (R_{t})_{i}^{n}+(R_{s})_{i}^{n}+(R_{l})_{i}^{n}$ with

By the stability analysis in Theorem 4.1, if the maximum step size fulfills the condition stated in (4.8), one can derive the following inequality:

with $\tilde{C} = E_{\alpha}\left(11\kappa_{+}\max\{1, \rho\}(\log t_{n})^{\alpha}\right).$

Now, our sole task is to estimate the terms $\sum_{j = 1}^{n}P_{n-j}^{(n)}\|(R_{l})^{n}\|$ and $\sum_{j = 1}^{n}P_{n-j}^{(n)}\|(R_{t})^{n}\|$, which we will accomplish by utilizing the lemmas outlined below.

Lemma 4.1. Assume that A3 holds, and $u(x, t)$ satisfies the regularity assumption (1.5), then

Proof. Set $v(t) = (\delta_{x}^2-\kappa)u(x_{i}, t)$. Using the Taylor-like formula given in Lemma 3.1, we derive

By combining (4.18) and (4.19), we obtain

This leads to the following expression for $(R_{l})_{i}^{j}$:

For $j = 1$, utilizing the regularity assumption (1.5), we deduce from (4.20) that

Analogously, for $2\leq j\leq N$, we have

It follows from (2.4), A1, and A2 that

By integrating the information from (4.21) to (4.23) along with Lemma 2.4, we arrive at the conclusion stated in (4.17). □

Lemma 4.2. Assume that A3 holds and $u(x, t)$ satisfies the regularity assumption (1.5), then

Proof. It follows from ^[29] that the local truncation error defined in (4.13) can be bounded by

where

By the regularity assumption (1.5), it is easy to obtain

and

Similar to (7.8), we have

Applying Lemma 2.5 by taking $\upsilon_{k} = \mathcal{G}^{k}$ in (2.11) to (4.29), we get

According to (4.26), (4.27), and (4.28), we have

and

Therefore, combining (4.30), (4.31), and (3.21) we obtain the desired result. □

Theorem 4.2. Let $u(x, t)$ be the exact solution of (1.1)–(1.3), and let $\{u_{i}^{n}|0\leq i\leq M, 0\leq n\leq N\}$ be the solution of (4.2)–(4.4). Suppose the regularity assumptions (1.4)–(1.5) hold. Let

If the maximum time step is $\tau\leq \mathfrak{B}_2,$ then the discrete solution is convergent in the $L^2$-norm with $\tilde{C} = E_{\alpha}\left(11\kappa_{+}\max\{1, \rho\}(\log t_{n})^{\alpha}\right)$ such that

where

In particular, if graded mesh is used, then it holds that

Proof. By combining (3.19), (4.16), Lemma 4.1 and Lemma 4.2, we arrive at (4.33). We now proceed to examine the global approximation error on the graded mesh. Due to (2.2) and (3.25), we can estimate the righthand side of (4.33) term by term as follows. It easy to check that

Further, it holds that

For the third term in the righthand side of (4.33), we have

If $r\sigma\geq 2$, then $\max\limits_{2\leq k\leq n}k^{r\sigma-2} = n^{r\sigma-2}\leq N^{r\sigma-2}.$ If $r\sigma < 2$, then $\max\limits_{2\leq k\leq n}k^{r\sigma-2} = 2^{r\sigma-2} < 1.$ Thus, we obtain

Moreover, we have

Combining (4.34) and (4.39), we obtain

If $r\sigma\geq3-\alpha$, we have from (4.40) that

If $r\sigma < 3-\alpha$, we get from (4.40) that

Note that

then we have from (4.42) that

By combining (4.36), (4.37), (4.38), (4.41), and (4.44), we achieve the desired global approximation error (4.35) on the graded mesh. □

5.   Numerical expertiments

In this section, we present several numerical examples to validate the theoretical result stated in Theorem 3.2 and Theorem 4.2. In our computations, the spatial domain $\Omega = [0, \pi]$ is uniformly divided into $M$ parts, and the time interval $\Lambda$ is divided into $N$ subintervals using graded meshes $t_k = a\left(\frac{T}{a}\right)^{(k/N)^{r}}$. The grading constant $r\geq 1$ controls the extent to which the time levels are concentrated near $t = a$. As $r$ increases, the initial step sizes become smaller than the later ones, which can be visually observed in Figure 1.

Example 5.1. Consider the problem (1.1)–(1.3) with

It can be verified that the corresponding exact solution is $u(x, t) = (\sin x) (\log{t})^{\alpha}$ and $\sigma = \alpha$.

Since the spacial accuracy is standard for the second order central difference scheme, we only explore the convergence rate of the time stepping scheme. To this end, we calculate the $L^2$ errors between the exact and numerical solutions

In Tables 1–3, we list the temporal $L^2$ errors for the L1 scheme by taking fixed $M$ and increasing $N$ for different $\alpha$. Based on the obtained numerical errors, we further estimate the order of temporal convergence using the formula as

Results are also listed in Tables 1–3. As can be observed, the temporal convergence order for the L1 scheme is close to $\min\{r\sigma, 2-\alpha\}$ for all cases, aligning with the theoretical analysis presented in Theorem 3.2. In Tables 4–7, we list the temporal $L^2$ errors for Alikhanov's scheme as a function of $N$ for different $\alpha$. Also shown are the corresponding decay rates based on graded meshes about Alikhanov's scheme. From Tables 4–7, it is observed that the convergence rate is close to $\min\{r\sigma, 2\}$ in time.

Example 5.2. Consider the problem (1.1) with

Tables 8 and 9 display the temporal $L^2$ errors by taking fixed $M$ and increasing $N$ for the L1 scheme and Alikhanov's scheme with $\kappa = -1$, respectively. Additionally, the maximum time step size $\tau$, as well as the conditions $\mathfrak{B}_1$ and $\mathfrak{B}_2$ defined in (3.10) and (4.8), respectively, are also presented. The convergent results displayed show that the rate of convergence in time is in agreement with Theorem 3.2 and Theorem 4.2. It is worth noting that when $\alpha$ approaches $1^-$, the convergence order in the time direction is also in accordance with the theoretical analysis.

6.   Concluding remarks

In this paper, we have proposed L1 and Alikhanov schemes with nonuniform time steps for solving Caputo-Hadamard fractional reaction sub-diffusion equations. We conduct a rigorous analysis of the stability and convergence of these two schemes, and further derive $\alpha$-robust error estimates under specific regularity conditions imposed on the exact solution. The derivation of these regularity assumptions is currently under active investigation and will be the subject of our forthcoming research.

7.   Appendix

7.1.   The proof of Lemma 3.2

Proof. Let $v(t) = u(x_i, t)$, for $1\leq k\leq n$, and using the Taylor-like formula given in Lemma 3.1, we derive

Thus the truncation error at time $t = t_n$ is given as

Exchanging the order of integration, we have

For the sake of simplicity, we denote $\varpi_n(y) = \omega_{2-\alpha}(\log t_n-\log y)$. Note that the linear logarithmic interpolation function of $\varpi_n(y)$ with respect to the nodes $t_{k-1}$ and $t_k$ is given by

Let

then (7.1) can be rewritten as follows:

Similar to the proof of Theorem 3.1 in ^[29], $\widetilde{L}_{\log, 1}^{k}\varpi_n(y)$ can be rewritten as

where $\xi_k(y, s) = \max\{\log y-\log s, 0\}-(\log y-\log t_{k-1})(\log t_k-\log s)/\tau_k$ such that

For $k = n$, since the function $\varpi_n(y)$ is decreasing with respect to $y$ and $\delta^2\varpi_n(y) < 0$, one has

Thus,

For $1\leq k\leq n-1$, using (7.3) and (7.4), one has

Then we obtain for $1\leq k\leq n-1$ that

Combining (7.2), (7.5), and (7.7), the proof is complete. □

7.2.   The proof of Lemma 3.3

Proof. After multiplying the inequality (3.20) by $P_{n-k}^{(n)}$ and summing the index $k$ from 1 to $n$, it is possible to switch the order of summation and apply the definition (2.5) of $P_{n-k}^{(n)}$ to obtain

Applying Lemma 2.5 by taking $\upsilon_{k} = \mathcal{G}^{k}$ in (2.11) to (7.8), one has

On one hand, it follows from (3.2) that

On the other hand, the assumption of regularity (1.5) leads to the conclusions that

and

for $2\leq k\leq n.$ By combining (7.9), (7.10), and (7.11) with (7.12), we obtain the desired result. □

Use of AI tools declaration

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

Acknowledgments

This work is supported by the Natural Science Foundation of Fujian Province of China (2022J01338, 2024J01119) and Fujian Alliance of Mathematics of China (2024SXLMMS03).

Conflict of interest

The authors declare that they have no competing interests.