High-Order Convolution Quadrature Scheme for the Subdiffusion Equation with Nonsmooth Data

1. Introduction

This paper is concerned with the development and analysis of high-order time-stepping discrete methods for the following initial-boundary value problem of the subdiffusion equation [1,2]:

$$\begin{array}{cc}\hfill {}_0{}^{C}D_t^{\alpha }u(x,t)+Au(x,t)& =f(x,t)\mathrm{in}\Omega \times (0,T],\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill u(x,t)& =0\mathrm{on}\partial \Omega \times (0,T],\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill u(x,0)& =v\left(x\right)\mathrm{in}\Omega ,\hfill \end{array}$$

(1c)

where $\Omega \subset {\mathbb{R}}^d$ ($d\in \{1,2,3\}$) is a convex polygonal or polyhedral domain with boundary $\partial \Omega$, and $T>0$ is a fixed number. Here, $f:\overline{\Omega }\times [0,T]\to \mathbb{R}$ is a given source function, and $v\in L^2\left(\Omega \right)$ is a given initial data. The operator $A=-\Delta$ is the negative Laplacian defined on the domain $D\left(A\right)=H^2\left(\Omega \right)\cap H_0^1\left(\Omega \right)\subset L^2\left(\Omega \right)$. ${}_0{}^{C}D_t^{\alpha }u$ is the left-sided Caputo fractional derivative of order $\alpha \in (0,1)$, defined by [9]

$${}_0{}^{C}D_t^{\alpha }u(x,t)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^t{(t-s)}^{-\alpha }{\displaystyle \frac{\partial u(x,s)}{\partial s}}ds,$$

(2)

where $\Gamma (·)$ denotes the Euler Gamma function. The subdiffusion equation is a valuable tool for describing slow transport processes in complex systems, including subsurface environments, turbulent plasmas, and heterogeneous porous media [1]. To capture such anomalous behavior accurately, high-order time discretization schemes are often essential, especially when the data are nonsmooth [2].

Finite difference methods for the Caputo fractional derivative mainly fall into two classes: L1-type and Grünwald-Letnikov-type schemes. The former handles variable coefficients and nonsmooth data more flexibly. Within the L1 framework, the standard scheme achieves $(2-\alpha )$ order for smooth solutions, whereas the L1-2 formula raises this to $(3-\alpha )$ via quadratic interpolation [4]. However, these polynomial-based schemes demand high regularity. As noted in [3], solutions typically exhibit weak singularities near $t=0$. Without strict compatibility conditions, this singular behavior severely degrades the convergence rate of high-order methods down to first order.

The Laplace (or Fourier) transform of fractional calculus yields $s^{\alpha }$ (or ${\left(i\omega \right)}^{\alpha }$) [9], motivating the approximation of fractional derivatives directly from the frequency domain. Generating functions facilitate this approach by mapping discrete operators into the Laplace domain. Recent developments heavily utilize the Riemann zeta and polylogarithm functions. Dimitrov et al. [7] constructed a $(2-\alpha )$-order approximation using the polylogarithm, and subsequently [8] derived an equivalent-order scheme via asymptotic expansions involving $\zeta \left(\alpha \right)$. These approaches endow the discrete weights with explicit integral representations. Such exactness significantly simplifies the subsequent error analysis.

We extend the zeta/polylogarithm framework of Dimitrov et al. [7,8] by constructing a family of generating functions of order $k-\alpha$ ($k=3,4,5,6,7$). This leads to our new time-stepping methods, the Zeta schemes. We then incorporate the starting correction technique of [12] directly into the discrete convolution. This specific modification alters both the generating function and the Laplace-domain integration contour, compensating for the initial singular behavior. To further suppress truncation errors stemming from the dominant singularity, we introduce an additional correction factor. A linear system involving the Riemann zeta function determines the necessary correction coefficients to cancel the leading errors. By employing contour integral techniques, we prove that the corrected schemes uniformly recover the optimal $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ rate for any fixed $t>0$. Remarkably, this optimal convergence holds for both homogeneous and inhomogeneous problems, even when dealing with incompatible initial data $v\in L^2\left(\Omega \right)$ and nonsmooth sources.

The remainder of this paper is organized as follows. Section 2 and Section 3 detail the $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ time discretization. Here, we introduce the starting correction terms and compute the coefficients needed for the high-order generating function. Section 4 and Section 5 contain the core theoretical analysis: deriving uniform error estimates for the homogeneous and inhomogeneous cases, respectively, alongside the high-order correction factor. Section 6 presents numerical experiments that validate our theoretical bounds. Section 7 gives a brief conclusion.

2. High-Order CQ schemes

2.1. Pre Knowledge

Let $\tau =T/N$ and $u_n=u\left(t_n\right)$ be the value of the function $u\left(t\right)$ at the point $t_n=n\tau$. The L1 approximation of the Caputo derivative is constructed by dividing the interval $[0,t]$ into subintervals of equal length $\tau$ and approximating the first derivative on each subinterval using a second-order central difference approximation.

For any sequence ${\left\{g^n\right\}}_{n=0}^{\infty }\in {\ell }^2\left(L^2\left(\Omega \right)\right)$, let $\tilde{g}\left(\xi \right)$ denote the generating function of the sequence defined by

$$\tilde{g}\left(\xi \right)={\displaystyle \sum _{n=0}^{\infty }}{g}^n{\xi }^n\mathrm{for}\xi \in \mathbb{D}.$$

(3)

The L1 approximation formula is given by

$${}_0{}^{C}D_t^{\alpha }u\left(t_n\right)={\displaystyle \frac{{\tau }^{-\alpha }}{\Gamma (2-\alpha )}}{\displaystyle \sum _{j=0}^n}{b}_j{u}^{n-j}+\mathcal{O}\left({\tau }^{2-\alpha }\right),$$

(4)

where $b_0=1$, $b_n^{\left(\alpha \right)}={(n-1)}^{1-\alpha }-n^{1-\alpha }$, and $b_j={(j+1)}^{1-\alpha }-2j^{1-\alpha }+{(j-1)}^{1-\alpha }$ ($j=2,\dots ,n-1$). The popular L1 scheme [3] is associated with the generating symbol

$${\displaystyle \frac{{(1-\zeta )}^2}{\zeta \Gamma (2-\alpha )}}{\mathrm{Li}}_{\alpha -1}\left(\zeta \right)={\displaystyle \sum _{j=0}^{\infty }}{b}_j{\zeta }^j,$$

(5)

where

$${\mathrm{Li}}_p\left(z\right)={\displaystyle \sum _{j=1}^{\infty }}{\displaystyle \frac{z^j}{j^p}}=z+{\displaystyle \frac{z^2}{2^p}}+\dots +{\displaystyle \frac{z^n}{n^p}}+\dots$$

(6)

is the polylogarithmic function, which is well defined for $\left|z\right|<1$ and can be analytically continued to the split complex plane $\mathbb{C}\setminus [1,\infty )$.

The CQ generated by the k-th order BDF [10,11] is given by formal power series expansions. To establish our high-order CQ approximations, we rely on the following asymptotic summation behavior and singularity expansions.

Lemma 1. 

[13] If $u\left(t\right)\in C^k[0,T]$, then

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\tau }^{\alpha }}}{\displaystyle \sum _{j=1}^{n-1}}{\displaystyle \frac{u(t-j\tau )}{j^{1+\alpha }}}& =\Gamma (-\alpha ){}_{0}D_t^{\alpha }u\left(t\right)+{\displaystyle \sum _{m=0}^k}{\displaystyle \frac{{(-1)}^m}{m!}}\zeta (1+\alpha -m)u^{\left(m\right)}\left(t\right){\tau }^{m-\alpha }\hfill \\ \hfill & -{\displaystyle \sum _{m=0}^k}{\displaystyle \frac{B_{2m}}{\left(2m\right)!}}\left({\displaystyle \sum _{l=0}^{2m-1}}\left({\displaystyle \genfrac{}{}{0pt}{}{2m-1}{l}}\right){\displaystyle \frac{\Gamma (1+\alpha +l)}{\Gamma (1+\alpha )}}{\displaystyle \frac{u^{(2m-1-l)}\left(0\right)}{t^{1+\alpha +l}}}\right){\tau }^{2m}.\hfill \end{array}$$

(7)

In particular, if $u^{\left(m\right)}\left(0\right)=0$ for $m=0,1,\dots ,k$, then

$${\displaystyle \frac{1}{{\tau }^{\alpha }}}{\displaystyle \sum _{j=1}^{n-1}}{\displaystyle \frac{u(t-j\tau )}{j^{1+\alpha }}}=\Gamma (-\alpha ){}_{0}D_t^{\alpha }u\left(t\right)+{\displaystyle \sum _{m=0}^k}{\displaystyle \frac{{(-1)}^m}{m!}}\zeta (1+\alpha -m)u^{\left(m\right)}\left(t\right){\tau }^{m-\alpha }.$$

(8)

Lemma 2. 

[24] For $p\ne 1,2,\dots$, the function ${\mathrm{Li}}_p\left(e^{-z}\right)$ satisfies the singular expansion

$${\mathrm{Li}}_p\left(e^{-z}\right)\sim \Gamma (1-p)z^{p-1}+{\displaystyle \sum _{l=0}^{\infty }}{(-1)}^l\zeta (p-l){\displaystyle \frac{z^l}{l!}}\mathrm{as}z\to 0,$$

(9)

where ζ is the Riemann zeta function.

Proof. 

To establish the absolute convergence of the series in a neighborhood of $z=0$, we employ the functional equation for the Riemann zeta function:

$$\zeta (p-l)={\displaystyle \frac{2}{{\left(2\pi \right)}^{1-p+l}}}cos\left({\displaystyle \frac{\pi (p-l)}{2}}\right)\Gamma (1-p+l)\zeta (1-p+l).$$

(10)

As $l\to \infty$, the argument $1-p+l\to +\infty$, which implies $\zeta (1-p+l)\to 1$. Furthermore, by Stirling’s formula (or the asymptotic ratio of Gamma functions), we have ${\displaystyle \frac{\Gamma (1-p+l)}{l!}}\sim l^{-p}$ as $l\to \infty$. Consequently, the general term of the series satisfies

$$\left|\zeta (p-l){\displaystyle \frac{z^l}{l!}}\right|\le C{\displaystyle \frac{l^{-p}}{{\left(2\pi \right)}^l}}{\left|z\right|}^l,$$

(11)

for some constant $C>0$. By the ratio test, this power series has a radius of convergence $R=2\pi$. Since we are concerned with the asymptotic behavior as $z\to 0$, the series converges absolutely in this domain.    □

In particular, taking $p=\alpha -1$ in Lemma 2, we obtain the specific singular expansion

$${\mathrm{Li}}_{\alpha -1}\left(e^{-z}\right)\sim \Gamma (2-\alpha )z^{\alpha -2}+{\displaystyle \sum _{l=0}^{\infty }}{(-1)}^l\zeta (\alpha -1-l){\displaystyle \frac{z^l}{l!}}\mathrm{as}z\to 0.$$

(12)

2.2. High-Order Time Discretization Schemes

Inspired by the corrected scheme in [13] and drawing on the idea of polynomial interpolation, we construct new corrected schemes for both smooth and nonsmooth cases based on Lemma 1. The key idea is to regard the convolution weights of the discrete scheme as the coefficients of a power series $\tilde{b}\left(\xi \right)$, which is understood in the sense of formal power series for the purpose of asymptotic expansion.

Theorem 1. 

Let $0<\alpha <1$ and let $k\ge 2$ be a given integer. Assume that the function $u\in C^k[0,T]$ and its derivatives up to order k vanish at $t=0$. Construct the generating function

$$\Gamma (-\alpha )\tilde{b}\left(\xi \right)=\Gamma (-\alpha ){\displaystyle \sum _{j=0}^{\infty }}{b}_j{\xi }^j={\displaystyle \sum _{j=0}^{\infty }}{\displaystyle \frac{{\xi }^j}{j^{1+\alpha }}}+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{\xi }^j={Li}_{1+\alpha }\left(\xi \right)+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{\xi }^j,$$

(13)

where the coefficients ${\left\{c_j\right\}}_{j=0}^{k-1}$ are uniquely determined by the equations

$$\zeta (\alpha +1-m)+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{j}^m=0,m=0,1,\dots ,k-1.$$

(14)

Then the convolution scheme whose weight coefficients $b_j$ are derived by $\tilde{b}\left(\xi \right)$,

$${\displaystyle \frac{1}{{\tau }^{\alpha }}}{\displaystyle \sum _{j=0}^{\infty }}{b}_{j}u(t_n-j\tau )={}_{0}D_t^{\alpha }u\left(t_n\right)+\mathcal{O}\left({\tau }^{k-\alpha }\right),$$

(15)

has $k-\alpha$ order accuracy.

Proof. 

We carry out the analysis in the frequency domain. Define the error symbol

$$E(z,\tau ):={\tau }^{-\alpha }\Gamma (-\alpha ){\displaystyle \sum _{j=0}^{\infty }}{b}_j{e}^{-jz\tau }-\Gamma (-\alpha )z^{\alpha },$$

(16)

where z lies on a suitable contour in the left half-plane. It suffices to show that $E(z,\tau )=\mathcal{O}\left({\tau }^{k-\alpha }\right)$ as $\tau \to 0$, uniformly on the contour. The desired time-domain estimate then follows by applying the inverse Laplace transform together with standard Paley–Wiener type arguments.

Since $\Gamma (-\alpha ){\sum }_{j=0}^{\infty }{b}_j{e}^{-jz\tau }=F\left(e^{-z\tau }\right)$, we have

$$E(z,\tau )={\tau }^{-\alpha }F\left(e^{-z\tau }\right)-\Gamma (-\alpha )z^{\alpha }.$$

(17)

Substituting the explicit form of F, we obtain

$$F\left(e^{-z\tau }\right)={Li}_{1+\alpha }\left(e^{-z\tau }\right)+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{e}^{-jz\tau }.$$

(18)

Using the expansion

$${Li}_{1+\alpha }\left(e^{-z}\right)=\Gamma (-\alpha )z^{\alpha }+{\displaystyle \sum _{m=0}^{\infty }}{(-1)}^m{\displaystyle \frac{\zeta (\alpha +1-m)}{m!}}{z}^m$$

(19)

and

$$e^{-jz}={\displaystyle \sum _{m=0}^{\infty }}{\displaystyle \frac{{(-j)}^m}{m!}}{z}^m,$$

(20)

we have

$${\tau }^{-\alpha }F\left(e^{-z\tau }\right)=\Gamma (-\alpha )z^{\alpha }+{\displaystyle \sum _{m=0}^{\infty }}\left({(-1)}^m{\displaystyle \frac{\zeta (\alpha +1-m)}{m!}}+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{\displaystyle \frac{{(-j)}^m}{m!}}\right)z^m{\tau }^{m-\alpha }.$$

(21)

Truncating high order terms for accuracy yields

$${\tau }^{-\alpha }F\left(e^{-z\tau }\right)=\Gamma (-\alpha )z^{\alpha }+{\displaystyle \sum _{m=0}^{k-1}}\left({(-1)}^m{\displaystyle \frac{\zeta (\alpha +1-m)}{m!}}+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{\displaystyle \frac{{(-j)}^m}{m!}}\right)z^m{\tau }^{m-\alpha }+\mathcal{O}\left({\tau }^{k-\alpha }\right).$$

(22)

Thus,

$$E(z,\tau )={\displaystyle \sum _{m=0}^{k-1}}\left({(-1)}^m{\displaystyle \frac{\zeta (\alpha +1-m)}{m!}}+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{\displaystyle \frac{{(-j)}^m}{m!}}\right)z^m{\tau }^{m-\alpha }+\mathcal{O}\left({\tau }^{k-\alpha }\right).$$

(23)

To achieve $E(z,\tau )=\mathcal{O}\left({\tau }^{k-\alpha }\right)$, it is necessary and sufficient to let the first several terms of the above expansion vanish. Factoring out the common ${(-1)}^m$ multiplier from the summation, the condition reduces to:

$$\zeta (\alpha +1-m)+{\displaystyle \sum _{j=0}^{k-1}}{c}_j{j}^m=0,m=0,1,\dots ,k-1.$$

(24)

Finally, applying the inverse Laplace transform gives the time-domain error estimate

$${\displaystyle \frac{1}{{\tau }^{\alpha }}}{\displaystyle \sum _{j=0}^{\infty }}{b}_{j}u(t_n-j\tau )={}_{0}D_t^{\alpha }u\left(t_n\right)+\mathcal{O}\left({\tau }^{k-\alpha }\right).$$

(25)

This completes the proof.    □

Since $u\left(0\right)=0$, we extend the solution by zero for negative times to derive the numerical scheme. We have

$${\displaystyle \sum _{j=0}^{\infty }}{b}_{j}u(t_n-j\tau )={\displaystyle \sum _{j=0}^n}{b}_{j}u(t_n-j\tau ),$$

(26)

i.e.,

$${\displaystyle \frac{1}{{\tau }^{\alpha }}}{\displaystyle \sum _{j=0}^n}{b}_{j}u(t_n-j\tau )={}_{0}D_t^{\alpha }u\left(t_n\right)+\mathcal{O}\left({\tau }^{k-\alpha }\right).$$

(27)

We now discuss special cases corresponding to indices $k=2,3,4,5,6,7$. Solving the equations

$$\left\{\begin{array}{c}\zeta (\alpha +1)+c_0+c_1=0,\hfill \\ \zeta \left(\alpha \right)+c_1=0\hfill \end{array}\right.$$

(28)

by taking $k=2$, we get $c_0=\zeta \left(\alpha \right)-\zeta (1+\alpha )$ and $c_1=-\zeta \left(\alpha \right)$, yielding the weight coefficients

$$\begin{array}{cc}\hfill b_0^{\left(\alpha \right)}& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}(\zeta \left(\alpha \right)-\zeta (1+\alpha )),b_1^{\left(\alpha \right)}={\displaystyle \frac{1}{\Gamma (-\alpha )}}(1-\zeta \left(\alpha \right)),b_j^{\left(\alpha \right)}={\displaystyle \frac{1}{\Gamma (-\alpha )j^{1+\alpha }}}(2\le j\le n).\hfill \end{array}$$

(29)

By means of the generating function, we obtain the same result as in [13].

If we take $k=3$, the system of equations is rewritten as

$$\left\{\begin{array}{c}\zeta (\alpha +1)+c_0+c_1+c_2=0,\hfill \\ \zeta \left(\alpha \right)+c_1+2c_2=0,\hfill \\ \zeta (\alpha -1)+c_1+4c_2=0.\hfill \end{array}\right.$$

(30)

Solving the equations yields $c_0={\displaystyle \frac{1}{2}}(-\zeta (-1+\alpha )+3\zeta \left(\alpha \right)-2\zeta (1+\alpha ))$, $c_1=\zeta (-1+\alpha )-2\zeta \left(\alpha \right)$, and $c_2=-{\displaystyle \frac{1}{2}}\zeta (-1+\alpha )+{\displaystyle \frac{1}{2}}\zeta \left(\alpha \right)$, which imply

$$\begin{array}{cc}\hfill b_0& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{2}}(-\zeta (-1+\alpha )+3\zeta \left(\alpha \right)-2\zeta (1+\alpha ))\right],\hfill \\ \hfill b_1& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[1+\zeta (-1+\alpha )-2\zeta \left(\alpha \right)\right],\hfill \\ \hfill b_2& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{2^{1+\alpha }}}-{\displaystyle \frac{1}{2}}\zeta (-1+\alpha )+{\displaystyle \frac{1}{2}}\zeta \left(\alpha \right)\right],\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{j^{1+\alpha }\Gamma (-\alpha )}},j=3,4,\dots \hfill \end{array}$$

(31)

If we take $k=4$ and solve the equations, namely,

$$\left\{\begin{array}{c}\zeta (\alpha +1)+c_0+c_1+c_2+c_3=0,\hfill \\ \zeta \left(\alpha \right)+c_1+2c_2+3c_3=0,\hfill \\ \zeta (\alpha -1)+c_1+4c_2+9c_3=0,\hfill \\ \zeta (\alpha -2)+c_1+8c_2+27c_3=0,\hfill \end{array}\right.$$

(32)

we get

$$\begin{array}{cc}\hfill c_0& ={\displaystyle \frac{1}{6}}\zeta (-2+\alpha )-\zeta (-1+\alpha )+{\displaystyle \frac{11}{6}}\zeta \left(\alpha \right)-\zeta (1+\alpha ),\hfill \\ \hfill c_1& ={\displaystyle \frac{1}{2}}\left(-\zeta (-2+\alpha )+5\zeta (-1+\alpha )-6\zeta \left(\alpha \right)\right),\hfill \\ \hfill c_2& ={\displaystyle \frac{1}{2}}\zeta (-2+\alpha )-2\zeta (-1+\alpha )+{\displaystyle \frac{3}{2}}\zeta \left(\alpha \right),\hfill \\ \hfill c_3& =-{\displaystyle \frac{1}{6}}\zeta (-2+\alpha )+{\displaystyle \frac{1}{2}}\zeta (-1+\alpha )-{\displaystyle \frac{1}{3}}\zeta \left(\alpha \right),\hfill \end{array}$$

(33)

which imply the weight coefficients

$$\begin{array}{cc}\hfill b_0& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{6}}\zeta (-2+\alpha )-\zeta (-1+\alpha )+{\displaystyle \frac{11}{6}}\zeta \left(\alpha \right)-\zeta (1+\alpha )\right],\hfill \\ \hfill b_1& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[1+{\displaystyle \frac{1}{2}}\left(-\zeta (-2+\alpha )+5\zeta (-1+\alpha )-6\zeta \left(\alpha \right)\right)\right],\hfill \\ \hfill b_2& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{2^{1+\alpha }}}+{\displaystyle \frac{1}{2}}\zeta (-2+\alpha )-2\zeta (-1+\alpha )+{\displaystyle \frac{3}{2}}\zeta \left(\alpha \right)\right],\hfill \\ \hfill b_3& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{3^{1+\alpha }}}-{\displaystyle \frac{1}{6}}\zeta (-2+\alpha )+{\displaystyle \frac{1}{2}}\zeta (-1+\alpha )-{\displaystyle \frac{1}{3}}\zeta \left(\alpha \right)\right],\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{j^{1+\alpha }}},(j\ge 4).\hfill \end{array}$$

(34)

To remove the stringent requirement that $u^{\left(k\right)}\left(0\right)=0$ for all $k\le 4$, we introduce a corrected scheme that is applicable to smooth solutions without such boundary restrictions. Designed along the local singularity matching principle of Jin et al. [12], this strategy applies specific algebraic weights directly to the right-hand side. Denoting the spatial operator by $A=-{\Delta }_h$ and setting $V^n=u_h\left(t_n\right)-u_h\left(0\right)$, the corrected scheme for the homogeneous equation reads

$${\tau }^{-\alpha }{\displaystyle \sum _{j=1}^n}{b}_{n-j}{V}^j+A{V}^n=-\left(1+{\overline{c}}_n^{\left(4\right)}\right)Av,1\le n\le 3,$$

(35)

where the correction coefficients are

$${\overline{c}}_1^{\left(4\right)}={\displaystyle \frac{31}{24}},{\overline{c}}_2^{\left(4\right)}=-{\displaystyle \frac{7}{6}},{\overline{c}}_3^{\left(4\right)}={\displaystyle \frac{3}{8}},$$

(36)

and ${\overline{c}}_n^{\left(4\right)}=0$ for $n\ge 4$.

For the inhomogeneous equation the correction involves the source function $g_h\left(t\right)=f_h\left(t\right)-A{u}_h\left(0\right)$ and its derivatives at $t=0$:

$${\tau }^{-\alpha }{\displaystyle \sum _{j=1}^n}{b}_{n-j}{V}^j+A{V}^n=g_h\left(t_n\right)+{\overline{c}}_n^{\left(4\right)}{g}_h\left(0\right)+{\displaystyle \sum _{\ell =1}^3}{\tau }^{\ell }{d}_{\ell ,n}^{\left(4\right)}{g}_h^{(\ell )}\left(0\right),1\le n\le 3.$$

(37)

For $n\ge 4$ the right-hand side reduces to $g_h\left(t_n\right)$. Matching the Taylor expansion of $g_h\left(t\right)$ at $t_n$ with the homogeneous correction factors gives the relation

$$d_{\ell ,n}^{\left(4\right)}={\displaystyle \frac{n^{\ell }}{\ell !}}{\overline{c}}_n^{\left(4\right)}.$$

(38)

If we take $k=5$, we have

$$\begin{array}{cc}\hfill b_0& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{24}}\left[-\zeta (-3+\alpha )+10\zeta (-2+\alpha )-35\zeta (-1+\alpha )+50\zeta \left(\alpha \right)-24\zeta (1+\alpha )\right],\hfill \\ \hfill b_1& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[1+{\displaystyle \frac{1}{6}}\left(\zeta (-3+\alpha )-9\zeta (-2+\alpha )+26\zeta (-1+\alpha )-24\zeta \left(\alpha \right)\right)\right],\hfill \\ \hfill b_2& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{2^{1+\alpha }}}+{\displaystyle \frac{1}{4}}\left(-\zeta (-3+\alpha )+8\zeta (-2+\alpha )-19\zeta (-1+\alpha )+12\zeta \left(\alpha \right)\right)\right],\hfill \\ \hfill b_3& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{3^{1+\alpha }}}+{\displaystyle \frac{1}{6}}\left(\zeta (-3+\alpha )-7\zeta (-2+\alpha )+14\zeta (-1+\alpha )-8\zeta \left(\alpha \right)\right)\right],\hfill \\ \hfill b_4& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left[{\displaystyle \frac{1}{4^{1+\alpha }}}+{\displaystyle \frac{1}{24}}\left(-\zeta (-3+\alpha )+6\zeta (-2+\alpha )-11\zeta (-1+\alpha )+6\zeta \left(\alpha \right)\right)\right],\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{j^{1+\alpha }}},(j\ge 5).\hfill \end{array}$$

(39)

If we take $k=6$, we get

$$\begin{array}{cc}\hfill b_0& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{120}}[\zeta (-4+\alpha )-15\zeta (-3+\alpha )+85\zeta (-2+\alpha )-225\zeta (-1+\alpha )\hfill \\ \hfill & +274\zeta \left(\alpha \right)-120\zeta (1+\alpha )],\hfill \\ \hfill b_1& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[1+{\displaystyle \frac{1}{24}}(-\zeta (-4+\alpha )+14\zeta (-3+\alpha )-71\zeta (-2+\alpha )+154\zeta (-1+\alpha )\hfill \\ \hfill & -120\zeta \left(\alpha \right)\left)\right],\hfill \\ \hfill b_2& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{2^{1+\alpha }}}+{\displaystyle \frac{1}{12}}\zeta (-4+\alpha )-{\displaystyle \frac{13}{12}}\zeta (-3+\alpha )+{\displaystyle \frac{59}{12}}\zeta (-2+\alpha )\hfill \\ \hfill & -{\displaystyle \frac{107}{12}}\zeta (-1+\alpha )+5\zeta \left(\alpha \right)],\hfill \\ \hfill b_3& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{3^{1+\alpha }}}-{\displaystyle \frac{1}{12}}\zeta (-4+\alpha )+\zeta (-3+\alpha )-{\displaystyle \frac{49}{12}}\zeta (-2+\alpha )\hfill \\ \hfill & +{\displaystyle \frac{13}{2}}\zeta (-1+\alpha )-{\displaystyle \frac{10}{3}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_4& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{4^{1+\alpha }}}+{\displaystyle \frac{1}{24}}\zeta (-4+\alpha )-{\displaystyle \frac{11}{24}}\zeta (-3+\alpha )+{\displaystyle \frac{41}{24}}\zeta (-2+\alpha )\hfill \\ \hfill & -{\displaystyle \frac{61}{24}}\zeta (-1+\alpha )+{\displaystyle \frac{5}{4}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_5& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{5^{1+\alpha }}}-{\displaystyle \frac{1}{120}}\zeta (-4+\alpha )+{\displaystyle \frac{1}{12}}\zeta (-3+\alpha )-{\displaystyle \frac{7}{24}}\zeta (-2+\alpha )\hfill \\ \hfill & +{\displaystyle \frac{5}{12}}\zeta (-1+\alpha )-{\displaystyle \frac{1}{5}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{j^{1+\alpha }}},(j\ge 6).\hfill \end{array}$$

(40)

If we take $k=7$, we get

$$\begin{array}{cc}\hfill b_0& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{720}}[-\zeta (-5+\alpha )+21\zeta (-4+\alpha )-175\zeta (-3+\alpha )+735\zeta (-2+\alpha )\hfill \\ \hfill & -1624\zeta (-1+\alpha )+1764\zeta \left(\alpha \right)-720\zeta (1+\alpha )],\hfill \\ \hfill b_1& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[1+{\displaystyle \frac{1}{120}}(\zeta (-5+\alpha )-20\zeta (-4+\alpha )+155\zeta (-3+\alpha )-580\zeta (-2+\alpha )\hfill \\ \hfill & +1044\zeta (-1+\alpha )-720\zeta \left(\alpha \right)\left)\right],\hfill \\ \hfill b_2& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{2^{1+\alpha }}}-{\displaystyle \frac{1}{48}}\zeta (-5+\alpha )+{\displaystyle \frac{19}{48}}\zeta (-4+\alpha )-{\displaystyle \frac{137}{48}}\zeta (-3+\alpha )\hfill \\ \hfill & +{\displaystyle \frac{461}{48}}\zeta (-2+\alpha )-{\displaystyle \frac{117}{8}}\zeta (-1+\alpha )+{\displaystyle \frac{15}{2}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_3& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{3^{1+\alpha }}}+{\displaystyle \frac{1}{36}}\zeta (-5+\alpha )-{\displaystyle \frac{1}{2}}\zeta (-4+\alpha )+{\displaystyle \frac{121}{36}}\zeta (-3+\alpha )\hfill \\ \hfill & -{\displaystyle \frac{31}{3}}\zeta (-2+\alpha )+{\displaystyle \frac{127}{9}}\zeta (-1+\alpha )-{\displaystyle \frac{20}{3}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_4& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{4^{1+\alpha }}}-{\displaystyle \frac{1}{48}}\zeta (-5+\alpha )+{\displaystyle \frac{17}{48}}\zeta (-4+\alpha )-{\displaystyle \frac{107}{48}}\zeta (-3+\alpha )\hfill \\ \hfill & +{\displaystyle \frac{307}{48}}\zeta (-2+\alpha )-{\displaystyle \frac{33}{4}}\zeta (-1+\alpha )+{\displaystyle \frac{15}{4}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_5& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{5^{1+\alpha }}}+{\displaystyle \frac{1}{120}}\zeta (-5+\alpha )-{\displaystyle \frac{2}{15}}\zeta (-4+\alpha )+{\displaystyle \frac{19}{24}}\zeta (-3+\alpha )\hfill \\ \hfill & -{\displaystyle \frac{13}{6}}\zeta (-2+\alpha )+{\displaystyle \frac{27}{10}}\zeta (-1+\alpha )-{\displaystyle \frac{6}{5}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_6& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}[{\displaystyle \frac{1}{6^{1+\alpha }}}-{\displaystyle \frac{1}{720}}\zeta (-5+\alpha )+{\displaystyle \frac{1}{48}}\zeta (-4+\alpha )-{\displaystyle \frac{17}{144}}\zeta (-3+\alpha )\hfill \\ \hfill & +{\displaystyle \frac{5}{16}}\zeta (-2+\alpha )-{\displaystyle \frac{137}{360}}\zeta (-1+\alpha )+{\displaystyle \frac{1}{6}}\zeta \left(\alpha \right)],\hfill \\ \hfill b_j& ={\displaystyle \frac{1}{\Gamma (-\alpha )}}{\displaystyle \frac{1}{j^{1+\alpha }}}(j\ge 7).\hfill \end{array}$$

(41)

Example 1. 

To confirm the theoretical convergence order of the above high-order discrete operators, we test a compatible configuration without temporal singularity. The initial condition is set to zero, $u(x,0)=0$, and the source term is taken as $f(x,t)=t^{k}sin\left(\pi x\right)$. Under this choice, the first $k-1$ temporal derivatives of f vanish identically at $t=0$. Consequently, the exact solution possesses sufficient regularity at the origin and the scheme attains its formal convergence order without any starting correction.

All computations are performed with the Multiprecision Computing Toolbox (MP) at 100 decimal digits of precision, in order to circumvent the truncation limit imposed by standard IEEE 754 double-precision arithmetic. The reference solution is computed with the fully corrected scheme on a uniform mesh with ${\tau }_{\mathrm{ref}}=2^{-14}$ for all k.

For clarity of presentation, Table 1 reports the results on the test grids $\tau \in \{2^{-3},2^{-4},\dots ,2^{-9}\}$ together with the corresponding convergence rates. On the finest grids, the errors approach the double-precision machine floor ($\sim {10}^{-15}$); points falling below ${10}^{-12}$ are excluded from the rate computation to avoid contamination by round-off. Notably, for the $k=7$ scheme, the extreme decay rate limits the valid data to only a few points before hitting the noise floor, yet the estimated rates remain theoretically consistent.

The convergence rates in Table 1 confirm that, once the singular barrier at the initial time is removed, all higher-order Zeta discrete operators attain the formal convergence order $\mathcal{O}\left({\tau }^{k-\alpha }\right)$ under sufficiently fine grids and high-precision arithmetic. This confirms that the discrete operators themselves attain the formal convergence order, independently of any starting correction. In future work, we will investigate generalized multi-point starting correction strategies for these higher-order ($k=5,6,7$) discrete operators under incompatible source terms and nonsmooth initial data.

Remark 1. 

For $m=0$ ($Re\left(s\right)=\alpha +1>1$), $\zeta (\alpha +1)$ is computed directly from the Dirichlet series, which incurs a truncation error bounded by $N^{-\alpha }/\alpha$ for N terms. For $m=2,3$, the functional equation (59) reflects the arguments to $Re(1-s)>1$, enabling evaluation via the absolutely convergent reflected series. For $m=1$, the argument maps to $1-\alpha$, which remains in the critical strip $(0,1)$; here, evaluation typically relies on the Euler–Maclaurin summation or the alternating Dirichlet eta function. All four constants are evaluated exactly once at $\mathcal{O}\left(1\right)$ cost independent of the time-step count, with robust implementations readily provided by standard numerical libraries.

3. High-Order CQ Schemes

The definition of the Caputo derivative can be rewritten to establish the relationship:

$${{}_{0}D}_t^{\alpha }(u(x,t)-u_0\left(x\right))={{}_0{}^{C}D}_t^{\alpha }u(x,t).$$

(42)

By introducing the homogenized variable $V(x,t)=u(x,t)-u_0\left(x\right)$, we obtain the equivalent system with the Riemann-Liouville derivative:

$$\begin{array}{cccc}\hfill {{}_{0}D}_t^{\alpha }V(x,t)+AV(x,t)& =-A{u}_0\left(x\right)+f(x,t)\hfill & & in\Omega \times (0,T],\hfill \end{array}$$

(43a)

$$\begin{array}{cccc}\hfill V(x,t)& =0\hfill & & on\partial \Omega \times (0,T],\hfill \end{array}$$

(43b)

$$\begin{array}{cccc}\hfill V(x,0)& =0\hfill & & in\Omega ,\hfill \end{array}$$

(43c)

with the left-sided Riemann-Liouville fractional derivative of order $\alpha \in (0,1)$ defined by

$${{}_{0}D}_t^{\alpha }V(x,t):={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\displaystyle \frac{\partial }{\partial t}}{\int }_0^{t}V(x,s){(t-s)}^{-\alpha }ds.$$

(44)

To analyze this system, we utilize the continuous Laplace transform and its associated properties. Since $V(x,0)=0$, the Caputo and Riemann-Liouville derivatives coincide. Assuming that the exact solution $u\left(t\right)$ is analytically extendable to the sector $|argz|<\pi /2$, the Laplace transform of the fractional derivative is well-defined as

$$\widehat{{}_0{}^{C}D_t^{\alpha }u}\left(z\right)=z^{\alpha }\widehat{u}\left(z\right)-z^{\alpha -1}u\left(0\right).$$

(45)

where the frequency variable z belongs to the sector ${\Sigma }_{{\theta }_0}$ with $\pi /2<{\theta }_0<\pi$. We define the spatial operator A as a self-adjoint positive definite second-order elliptic operator with dense domain $\mathcal{D}\left(A\right)=H_0^1\left(\Omega \right)\cap H^2\left(\Omega \right)$ in $L_2\left(\Omega \right)$. As a strongly elliptic operator, A is sectorial: there exists an angle ${\theta }_0\in (\pi /2,\pi )$ such that the resolvent estimate

$${\parallel (zI+A)}^{-1}\parallel \le C{\left|z\right|}^{-1}$$

(46)

holds uniformly for all z in the sector ${\Sigma }_{{\theta }_0}=\{z\in \mathbb{C}\setminus \left\{0\right\}:|argz|<{\theta }_0\}$. Consequently, A generates a bounded analytic semigroup, which implies that the solution $u\left(t\right)$ can be analytically extended to a sector containing the positive real axis. This justifies the use of the Laplace transform and contour integration in the complex plane.

We now fix an angle $\theta \in (\pi /2,{\theta }_0)$ sufficiently close to $\pi /2$. For $0<\alpha <1$ and any $z\in {\Sigma }_{\theta }=\{z:|argz|<\theta \}$, we have $arg\left(z^{\alpha }\right)=\alpha \theta <\theta <{\theta }_0$; thus $z^{\alpha }\in {\Sigma }_{{\theta }_0}$. Replacing z by $z^{\alpha }$ in the resolvent estimate, we obtain

$$\parallel {(z^{\alpha }I+A)}^{-1}{\parallel \le C\left|z\right|}^{-\alpha },$$

(47)

valid for all $z\in {\Sigma }_{\theta }$, where C depends only on $\theta$ and $\alpha$. From the operator identity $A{(z^{\alpha }I+A)}^{-1}=I-z^{\alpha }{(z^{\alpha }I+A)}^{-1}$, which follows directly from the commutativity of A with its resolvent, we obtain a key relation for subsequent contour integral estimates.

Applying the Laplace transform to (Section 3) yields $z^{\alpha }\widehat{V}\left(z\right)+A\widehat{V}\left(z\right)=-z^{-1}A{u}_0+\widehat{f}\left(z\right)$. Let ${\Gamma }_{\theta ,\kappa }$ be the sectorial contour defined by

$${\Gamma }_{\theta ,\kappa }=\{z\in \mathbb{C}:|z|=\kappa ,|argz|\le \theta \}\cup \{z\in \mathbb{C}:z=\rho e^{\pm i\theta },\rho \ge \kappa \},$$

(48)

oriented with an increasing imaginary part, where $\kappa >0$. Then, by taking the inverse Laplace transform along this contour, the function $V\left(t\right)$ can be represented as

$$V\left(t\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\theta ,\kappa }}{e}^{zt}\left(-K\left(z\right)A{u}_0+zK\left(z\right)\widehat{f}\left(z\right)\right)dz,$$

(49)

with the kernel function defined by

$$K\left(z\right):=z^{-1}{(z^{\alpha }I+A)}^{-1}.$$

(50)

To approximate the fractional derivative, we employ a high-order convolution quadrature approximation characterized by the generating function. For the sake of simplicity in the error analysis, we now focus on the homogeneous case ($f\equiv 0$). Following the continuous formulation introduced earlier, we directly seek the fully discrete solution $V^n$ at time level $t_n$ with the initial condition $V^0=0$. The fully discrete evolution scheme satisfies

$${\displaystyle \sum _{j=0}^n}{b}_j{V}^{n-j}+{\tau }^{\alpha }A{V}^n=-{\tau }^{\alpha }A{v}_0,$$

(51)

where $v_0$ denotes the discrete initial data. Multiplying both sides of the equation by ${\xi }^n$ and summing from $n=1$ to ∞ yields

$${\displaystyle \sum _{n=1}^{\infty }}\left({\displaystyle \sum _{j=0}^n}{b}_j{V}^{n-j}\right){\xi }^n+{\tau }^{\alpha }A\tilde{V}\left(\xi \right)=-{\displaystyle \frac{\xi }{1-\xi }}{\tau }^{\alpha }A{v}_0,$$

(52)

where $\tilde{V}\left(\xi \right)={\sum }_{n=1}^{\infty }{V}^n{\xi }^n$ and $\tilde{b}\left(\xi \right)={\sum }_{j=0}^{\infty }{b}_j{\xi }^j$ is the generating function. Using the fact $V^0=0$ and the discrete convolution property of generating functions, the first term on the left-hand side can be directly expressed as $\tilde{b}\left(\xi \right)\tilde{V}\left(\xi \right)$. This leads to the frequency-domain algebraic equation:

$$\tilde{b}\left(\xi \right)\tilde{V}\left(\xi \right)+{\tau }^{\alpha }A\tilde{V}\left(\xi \right)=-{\displaystyle \frac{\xi }{1-\xi }}{\tau }^{\alpha }A{v}_0.$$

(53)

Extracting the common factor, the generating function of the fully discrete solution can be represented in the form of a resolvent operator:

$$\tilde{V}\left(\xi \right)=-{\displaystyle \frac{\xi }{1-\xi }}{\left({\tau }^{-\alpha }\tilde{b}\left(\xi \right)I+A\right)}^{-1}A{v}_0.$$

(54)

It is easy to verify that the function $\tilde{V}\left(\xi \right)$ is analytic at $\xi =0$. Hence, Cauchy’s integral theorem implies that for sufficiently small $\rho$, there holds

$$V^n={\displaystyle \frac{1}{2\pi i}}{\int }_{\left|\xi \right|=\rho }{\displaystyle \frac{1}{{\xi }^{n+1}}}\tilde{V}\left(\xi \right)d\xi =-{\displaystyle \frac{1}{2\pi i}}{\int }_{\left|\xi \right|=\rho }{\displaystyle \frac{1}{1-\xi }}{\displaystyle \frac{1}{{\xi }^n}}{\left({\tau }^{-\alpha }\tilde{b}\left(\xi \right)I+A\right)}^{-1}A{v}_{0}d\xi .$$

(55)

Upon changing the variable $\xi =e^{-z\tau }$, we obtain

$$V^n=-{\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }^{\tau }}{e}^{z{t}_n}{\displaystyle \frac{\tau }{e^{z\tau }-1}}{\left({\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)I+A\right)}^{-1}A{v}_{0}dz,$$

(56)

where the contour ${\Gamma }^{\tau }:=\{z=-ln\left(\rho \right)/\tau +iy:|y|\le \pi /\tau \}$ corresponds to the counterclockwise orientation in the $\xi$-plane. By continuously deforming the contour to ${\Gamma }_{\tau }:=\{z\in {\Gamma }_{\theta ,\kappa }:|\Im \left(z\right)|\le \pi /\tau \}$ and using the periodicity of the exponential function (which ensures that the integrals along the horizontal segments $\mathrm{Im}\left(z\right)=\pm \pi /\tau$ perfectly cancel each other), we obtain

$$V^n=-{\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}{\displaystyle \frac{\tau }{e^{z\tau }-1}}{\left({\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)I+A\right)}^{-1}A{v}_{0}dz.$$

(57)

This integral representation forms the foundation for the subsequent error analysis and kernel function estimation.

Lemma 3.([10,11]) The term $\zeta \left(s\right)$ denotes the Riemann zeta function. For $Re\left(s\right)>1$, it is defined by the absolutely convergent Dirichlet series

$$\zeta \left(s\right)={\displaystyle \sum _{n=1}^{\infty }}{\displaystyle \frac{1}{n^s}}.$$

(58)

It extends analytically to a meromorphic function on $\mathbb{C}$ whose only singularity is a simple pole at $s=1$. For $s\notin \mathbb{Z}$, the analytic continuation satisfies the functional equation

$$\zeta \left(s\right)=2^s{\pi }^{s-1}sin\left({\displaystyle \frac{\pi s}{2}}\right)\Gamma (1-s)\zeta (1-s).$$

(59)

Consequently, for any fractional parameter $\alpha \in (0,1)$ and integer $m\in \{0,1,2,3\}$, the evaluations $\zeta (\alpha +1-m)$ are finite, well-defined real constants. These constants uniquely determine the high-order correction coefficients $c_j$ for the generating function $\tilde{b}\left(\xi \right)$, as explicitly derived earlier.

Lemma 4. 

Let $\alpha \in (0,1)$ and define the discrete fractional derivative operator associated with the fully discrete $4-\alpha$ order Zeta scheme as $z_{\tau }^{\left(\alpha \right)}={\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)$. Assume that there exist a positive constant $c_0>0$ and an angle ${\theta }_0\in (\pi /2,\pi )$ independent of τ, such that for any $z\ne 0$ with $|argz|\le {\theta }_0$, the discrete operator satisfies $|z_{\tau }^{\left(\alpha \right)}|\ge c_0{\left|z\right|}^{\alpha }$ and $|arg\left(z_{\tau }^{\left(\alpha \right)}\right)|\le {\theta }_0$. Furthermore, as verified numerically in Figure 1, we assume $z_{\tau }^{\left(\alpha \right)}$ has no zeros. Under these hypotheses, for all $z\in {\Gamma }_{\tau }$, the spectral equivalence holds:

$${c\left|z\right|}^{\alpha }\le |z_{\tau }^{\left(\alpha \right)}{|\le C|z|}^{\alpha }.$$

(60)

Proof. 

By the hypothesis of the lemma, we directly obtain the lower bound ${c\left|z\right|}^{\alpha }\le \left|z_{\tau }^{\left(\alpha \right)}\right|$ for all $z\in {\Gamma }_{\tau }$. To establish the upper bound, we first evaluate the limit as $z\tau \to 0$. Multiplying the numerator and denominator by ${\tau }^{\alpha }$, we have:

$$\underset{z\tau \to 0}{lim}{\displaystyle \frac{|z_{\tau }^{\left(\alpha \right)}|}{{\left|z\right|}^{\alpha }}}=\underset{z\tau \to 0}{lim}{\displaystyle \frac{{|\left(z\tau \right)}^{\alpha }+c_4{\left(z\tau \right)}^4+c_5{\left(z\tau \right)}^5+\dots |}{{\left|z\tau \right|}^{\alpha }}}=1.$$

(61)

Therefore, there exists a constant $r_0>0$ such that ${\displaystyle \frac{|z_{\tau }^{\left(\alpha \right)}|}{{\left|z\right|}^{\alpha }}}\le C$ holds for $0<\left|z\tau \right|<r_0$ and $\forall z\in {\Gamma }_{\tau }$.

On the other hand, when $\left|z\tau \right|\ge r_0$, we have $\theta \in ({\displaystyle \frac{\pi }{2}},{\theta }_0)$ with ${\theta }_0$ being sufficiently close to ${\displaystyle \frac{\pi }{2}}$. In this case, the modulus satisfies:

$$|\Im \left(z\right)|\tau \le \left|z\right|\tau \le {\displaystyle \frac{\pi }{sin\theta }}\le \sqrt{2}\pi \le M_0.$$

(62)

Combining this bounded domain with the empirical boundedness of the generating function verified in Figure 1, we obtain the upper bound:

$$|z_{\tau }^{\left(\alpha \right)}{|\le C|z|}^{\alpha },\mathrm{for}{r}_0\le \left|z\tau \right|\le {\displaystyle \frac{\pi }{sin\theta }},\forall z\in {\Gamma }_{\tau }.$$

(63)

This completes the proof.    □

Lemma 5. 

For the discrete operator $z_{\tau }^{\left(\alpha \right)}$ defined in Lemma 4, and for all z on the low-frequency contour portion ${\Gamma }_{\mathrm{low}}:=\{z\in {\Gamma }_{\tau }:\left|z\tau \right|\le \pi \}$, there exists a positive constant C independent of τ such that the truncation error satisfies

$$|z_{\tau }^{\left(\alpha \right)}-z^{\alpha }|\le C{\tau }^{4-\alpha }{\left|z\right|}^4.$$

(64)

Proof. 

From the frequency-domain discrete scheme, the truncation error can be expanded via the generating function. By the high-order construction of $\tilde{b}\left(\xi \right)$ matching the fractional derivative symbol, the low-order error terms are precisely cancelled. Thus, for $\left|z\tau \right|\le \pi$, we have:

$$|z_{\tau }^{\left(\alpha \right)}-z^{\alpha }|=\left|{\displaystyle \frac{\tilde{b}\left(e^{-z\tau }\right)}{{\tau }^{\alpha }}}-z^{\alpha }\right|=\left|c_4{z}^4{\tau }^{4-\alpha }+c_5{z}^5{\tau }^{5-\alpha }+\dots \right|=\mathcal{O}\left(z^4{\tau }^{4-\alpha }\right).$$

(65)

Since $\left|z\tau \right|$ is bounded by $\pi$, the remainder of the power series is uniformly bounded, leading directly to the bound $|z_{\tau }^{\left(\alpha \right)}-z^{\alpha }|\le C{\tau }^{4-\alpha }{\left|z\right|}^4$. This completes the proof.    □

Theorem 2. 

Let $\alpha \in (0,1)$ and let $u\in W^{4,1}(0,T)$, extended by zero to $t<0$. Here $W^{4,1}(0,T)$ is the Sobolev space of functions whose weak derivatives up to order 4 are integrable on $(0,T)$; the zero extension implies $u^{\left(k\right)}\left(0\right)=0$ for $k=0,1,2,3$. Assume the Laplace transform $\widehat{u}\left(z\right)$ is initially analytic for $Re\left(z\right)>0$ and admits an analytic continuation to a domain encompassing the integration contour ${\Gamma }_{\tau }$. Assume further that the rapid decay of $\widehat{u}\left(z\right)$ ensures ${\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^4|\widehat{u}\left(z\right)\left|\right|dz|<\infty$. The truncation error of the $4-\alpha$ order discrete operator satisfies

$$|D_{\tau }^{\alpha }u\left(t_n\right){-}_0^C{D}_t^{\alpha }u\left(t_n\right)|\le M{\tau }^{4-\alpha }{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^4|\widehat{u}\left(z\right)\left|\right|dz|,$$

(66)

where the constant $M>0$ is independent of τ and n. The discrete weights $\left\{b_j\right\}$ are the corresponding coefficients of the generating function $\tilde{b}\left(\xi \right)$, whose explicit formulations incorporated with the high-order corrections have been provided. Here, the $4-\alpha$ order discrete fractional approximation operator $D_{\tau }^{\alpha }$ at $t_n$ is defined as:

$$D_{\tau }^{\alpha }u\left(t_n\right)={\tau }^{-\alpha }{\displaystyle \sum _{j=0}^n}{b}_{j}u\left(t_{n-j}\right).$$

(67)

Proof. 

By the causality condition, $u\left(t\right)$ and its derivatives up to order 3 vanish at $t=0$, ensuring $\mathcal{L}\left[{}_0{}^C\mathrm{D}_t^{\alpha }u\right]\left(z\right)=z^{\alpha }\widehat{u}\left(z\right)$. The exact Caputo derivative and the discrete approximation operator evaluated at $t_n$ are expressed as:

$${}_0{}^C\mathrm{D}_t^{\alpha }u\left(t_n\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}{z}^{\alpha }\widehat{u}\left(z\right)dz,D_{\tau }^{\alpha }u\left(t_n\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}{\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)\widehat{u}\left(z\right)dz.$$

(68)

Subtracting the exact derivative from the discrete operator yields the absolute truncation error:

$$|D_{\tau }^{\alpha }u\left(t_n\right){-}_0^C{D}_t^{\alpha }u\left(t_n\right)|\le {\displaystyle \frac{1}{2\pi }}{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}\left|\right|E(z,\tau )\left|\right|\widehat{u}\left(z\right)\left|\right|dz|,$$

(69)

where $E(z,\tau )={\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)-z^{\alpha }$. We partition the contour into ${\Gamma }_{\mathrm{low}}=\{z\in {\Gamma }_{\tau }:\left|z\tau \right|\le \pi \}$ and ${\Gamma }_{\mathrm{high}}=\{z\in {\Gamma }_{\tau }:\pi /sin\theta >\left|z\tau \right|>\pi \}$.

On ${\Gamma }_{\mathrm{low}}$, Lemma 5 directly yields

$${\int }_{{\Gamma }_{\mathrm{low}}}|e^{z{t}_n}\left|\right|E(z,\tau )\left|\right|\widehat{u}\left(z\right)\left|\right|dz|\le C_1{\tau }^{4-\alpha }{\int }_{{\Gamma }_{\mathrm{low}}}|e^{z{t}_n}{\left|\right|z|}^4|\widehat{u}\left(z\right)\left|\right|dz|.$$

(70)

On ${\Gamma }_{\mathrm{high}}$, the condition $\left|z\tau \right|>\pi$ gives ${\tau }^{-1}<\left|z\right|/\pi$. Furthermore, for all $z\in {\Gamma }_{\mathrm{high}}$, the mapped variable $\xi =e^{-z\tau }$ forms a compact set strictly bounded away from the singularity at $\xi =1$. Recall from its explicit formulation (as defined in Lemma 3.3) that the generating function $\tilde{b}\left(\xi \right)$ consists of polynomials and the polylogarithm ${Li}_{1+\alpha }\left(\xi \right)$, which admits an analytic continuation to the split complex plane $\mathbb{C}\setminus [1,\infty )$. Since $\xi \notin [1,\infty )$ for all $z\in {\Gamma }_{\mathrm{high}}$, $\tilde{b}\left(\xi \right)$ is analytic in this compact region. Therefore, the extreme value theorem guarantees its uniform boundedness, i.e., $|{\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)|\le C_2{\tau }^{-\alpha }$. By the triangle inequality, $\left|E(z,\tau )\right|\le |{\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right){|+|z|}^{\alpha }\le C_2{\tau }^{-\alpha }+{\left|z\right|}^{\alpha }$. We bound the two terms separately. Since ${\tau }^{-1}<\left|z\right|/\pi$, we have

$${\tau }^{-\alpha }\le {\pi }^{-4}{\tau }^{4-\alpha }{\left|z\right|}^4.$$

(71)

Raising ${\tau }^{-1}<\left|z\right|/\pi$ to the power $4-\alpha$ gives

$${\left|z\right|}^{\alpha }\le {\pi }^{\alpha -4}{\tau }^{4-\alpha }{\left|z\right|}^4.$$

(72)

Combining (71) and (72), the error symbol on ${\Gamma }_{\mathrm{high}}$ is bounded by $C_3{\tau }^{4-\alpha }{\left|z\right|}^4$, where $C_3=C_2{\pi }^{-4}+{\pi }^{\alpha -4}$. The rapid decay of $\widehat{u}\left(z\right)$ ensures integral convergence

$${\int }_{{\Gamma }_{\mathrm{high}}}|e^{z{t}_n}\left|\right|E(z,\tau )\left|\right|\widehat{u}\left(z\right)\left|\right|dz|\le C_3{\tau }^{4-\alpha }{\int }_{{\Gamma }_{\mathrm{high}}}|e^{z{t}_n}{\left|\right|z|}^4|\widehat{u}\left(z\right)\left|\right|dz|.$$

(73)

Summing (70) and (73), and setting $M=max(C_1,C_3)/\left(2\pi \right)$, concludes the proof.    □

Remark 2. 

The decay assumption on $\widehat{u}\left(z\right)$ stated earlier ensures the integral remains uniformly bounded as $\tau \to 0$. This guarantees the theoretical $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ asymptotic rate. For practical computations with any fixed $\tau >0$, the truncation $\left|z\tau \right|\le \pi /sin\theta$ naturally restricts ${\Gamma }_{\tau }$ to a compact set. The continuous integrand thus trivially yields a finite integral.

Remark 3. 

The preceding error estimate relies on the assumption that the initial data and its derivatives up to order $k-1$ vanish at $t=0$. If any $u^{\left(m\right)}\left(0\right)\ne 0$ for $m=0,1,\dots ,k-1$, the exact solution typically exhibits weak singularities at the initial time, which leads to a severe order reduction in the standard convolution quadrature. Instead of modifying the last few terms of the convolution sum, we are inspired by Jin et al. [12], who introduced starting correction terms in the time domain to construct high-order schemes. Motivated by this, we propose a novel approach: directly matching and constructing the starting correction terms in the frequency domain to neutralize the initial singularities. The specific construction and detailed rigorous analysis of this modified scheme will be presented in the subsequent section.

4. Corrected Schemes for Homogeneous and Inhomogeneous Problem

In this section, we derive optimal error estimates for the fully discrete scheme. The analysis is based on the representations of the semidiscrete solution $V\left(t_n\right)$ and the fully discrete solution $V^n$. Subtracting the two gives:

$$V\left(t_n\right)-V^n=I+II,$$

(74)

where the terms I and $II$ are given by

$$I={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\tau }}{e}^{z{t}_n}{K}_1\left(z\right)vdz$$

(75)

and

$$II={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}(K_1\left(z\right)-K_2\left(z\right))vdz.$$

(76)

Here the integral kernel functions $K_1\left(z\right)$ and $K_2\left(z\right)$ are defined by

$$K_1\left(z\right):=-z^{-1}{(z^{\alpha }I+A)}^{-1}A=:z^{-1}{B}_1\left(z\right)$$

(77)

$$K_2\left(z\right):=-{\displaystyle \frac{\tau }{e^{z\tau }-1}}{\left({\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)I+A\right)}^{-1}A=:{\mu }^{-1}\left(z\right)B_2\left(z\right),$$

(78)

where

$$\mu \left(z\right)={\tau }^{-1}(e^{z\tau }-1),$$

(79)

and

$$B_2\left(z\right)=-{\left(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A\right)}^{-1}A.$$

(80)

Since $e^{z{t}_n}$ is uniformly bounded and decays exponentially on ${\Gamma }_{\tau }$, it suffices to bound $\parallel K_1\left(z\right)-K_2\left(z\right)\parallel$. We first derive estimates for the auxiliary function and the high-order discrete operator.

4.1. Error Analysis of the Uncorrected Scheme

Lemma 6. 

Let $\mu \left(z\right)={\tau }^{-1}(e^{z\tau }-1)$. Then for all $z\in {\Gamma }_{\tau }$, there hold for some $c_1,c_2>0$

$$|\mu \left(z\right)-z|\le {c\left|z\right|}^2\tau \mathrm{and}{c}_1\left|z\right|\le \left|\mu \left(z\right)\right|\le c_2\left|z\right|.$$

(81)

Proof. 

We note that $\left|z\tau \right|\le \pi /sin\theta$ for $z\in {\Gamma }_{\tau }$. By means of Taylor expansion, the first assertion follows directly:

$$|\mu \left(z\right)-z|\le {\left|z\right|}^2\tau \left|{\displaystyle \sum _{j=0}^{\infty }}{\displaystyle \frac{{\left(z\tau \right)}^j}{(j+2)!}}\right|\le c{\left|z\right|}^2\tau \forall z\in {\Gamma }_{\tau }.$$

(82)

Next we consider the second claim. The upper bound of $\mu \left(z\right)$ is trivial, so it suffices to verify the lower bound. Split the contour ${\Gamma }_{\tau }$ into three disjoint parts: ${\Gamma }_{\tau }={\Gamma }_{\tau }^{+}\cup {\Gamma }_{\tau }^c\cup {\Gamma }_{\tau }^{-}$. For $z\in {\Gamma }_{\tau }^{+}$, $z\tau =\rho e^{i\theta }$ with $\rho \in (0,\pi /sin\theta )$. By $cos(\rho sin\theta )\le 1$, we obtain

$$\left|{\displaystyle \frac{e^{z\tau }-1}{z\tau }}\right|={\displaystyle \frac{|e^{\rho cos\theta }cos(\rho sin\theta )-1+i{e}^{\rho cos\theta }sin(\rho sin\theta )|}{\rho }}\ge {\displaystyle \frac{1-e^{\rho cos\theta }}{\rho }}\ge c_1>0$$

(83)

due to the positivity and monotonicity of $(1-e^{\rho cos\theta })/\rho$ as a function of $\rho$ over the interval $(0,\pi /sin\theta ]$. The case of $z\in {\Gamma }_{\tau }^{-}$ follows analogously. Last, we consider $z\in {\Gamma }_{\tau }^c$, the circular arc. In this case, by means of Taylor expansion, we have $\mu \left(z\right)=z(1+\mathcal{O}(z\tau \left)\right)$. From this and the fact that $\left|z\tau \right|<1$, it follows directly that $\left|\mu \left(z\right)\right|\ge c_3\left|z\right|$. This completes the proof of the lemma. □

To estimate the spatial resolvent term $B_2\left(z\right)$, we define the high-order discrete kernel $z_{\tau }^{\left(\alpha \right)}\left(z\right)={\tau }^{-\alpha }\tilde{b}\left(e^{-z\tau }\right)$.

   □

Lemma 7. 

Let $B_1\left(z\right)=-{(z^{\alpha }I+A)}^{-1}A$ and $B_2\left(z\right)=-{(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}A$. Then for any $z\in {\Gamma }_{\tau }$, there holds:

$$\left|\right|B_1\left(z\right)-B_2\left(z\right)\left|\right|\le {C\left|z\right|}^{4-\alpha }{\tau }^{4-\alpha }.$$

(84)

Proof. 

Using the operator identity ${(z^{\alpha }I+A)}^{-1}A=I-z^{\alpha }{(z^{\alpha }I+A)}^{-1}$, we can rewrite the operators as:

$$B_1\left(z\right)=-I+z^{\alpha }{(z^{\alpha }I+A)}^{-1},B_2\left(z\right)=-I+z_{\tau }^{\left(\alpha \right)}\left(z\right){(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}.$$

(85)

By subtracting them and inserting cross terms, we obtain

$$\begin{array}{cc}\hfill \left|\right|B_1\left(z\right)-B_2\left(z\right)\left|\right|& \le |z^{\alpha }-z_{\tau }^{\left(\alpha \right)}\left(z\right)|·||{(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}\left|\right|\hfill \\ \hfill & +{\left|z\right|}^{\alpha }·\left|\right|{(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}(z_{\tau }^{\left(\alpha \right)}\left(z\right)-z^{\alpha }){(z^{\alpha }I+A)}^{-1}\left|\right|.\hfill \end{array}$$

(86)

Based on the resolvent estimate of the sectorial operator $\left|\right|{(z^{\alpha }I+A)}^{-1}\left|\right|\le {c\left|z\right|}^{-\alpha }$ and the bound provided that $\left|\right|{(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}\left|\right|\le c|z_{\tau }^{\left(\alpha \right)}{\left(z\right)|}^{-1}\le C{\left|z\right|}^{-\alpha }$, we deduce:

$$\begin{array}{cc}\hfill \parallel B_1\left(z\right)-B_2\left(z\right)\parallel & \le {C\left|z\right|}^{-\alpha }\left|z_{\tau }^{\left(\alpha \right)}\left(z\right)-z^{\alpha }\right|\hfill \\ \hfill & +{C\left|z\right|}^{\alpha }{\left|z\right|}^{-\alpha }\left|z_{\tau }^{\left(\alpha \right)}\left(z\right)-z^{\alpha }\right|{\left|z\right|}^{-\alpha }\hfill \\ \hfill & \le {C\left|z\right|}^{-\alpha }\left|z_{\tau }^{\left(\alpha \right)}\left(z\right)-z^{\alpha }\right|.\hfill \end{array}$$

(87)

we finally arrive at $\left|\right|B_1\left(z\right)-B_2\left(z\right)\left|\right|\le {C\left|z\right|}^{-\alpha }·{\left|z\right|}^K{\tau }^{4-\alpha }=C{\left|z\right|}^{4-\alpha }{\tau }^{4-\alpha }$.

   □

Lemma 8. 

Let $K_1\left(z\right)=z^{-1}{(z^{\alpha }I+A)}^{-1}A$ and $K_2\left(z\right)={\mu }^{-1}\left(z\right){(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}A$. Then for any $z\in {\Gamma }_{\tau }$, there holds

$$\left|\right|K_1\left(z\right)-K_2\left(z\right)\left|\right|\le {C\left|z\right|}^{-\alpha }\tau .$$

(88)

Proof. 

By subtracting the operators and inserting a cross term, we decompose the norm as:

$$\left|\right|K_1\left(z\right)-K_2\left(z\right)\left|\right|\le {\displaystyle \frac{|z-\mu (z\left)\right|}{\left|z\right|\left|\mu \right(z\left)\right|}}\left|\right|B_1\left(z\right)\left|\right|+{\left|\mu \left(z\right)\right|}^{-1}\left|\right|B_1\left(z\right)-B_2\left(z\right)\left|\right|.$$

(89)

From Lemma 3.1, $|z-\mu \left(z\right)|\le {c\left|z\right|}^2\tau$ and $\left|\mu \left(z\right)\right|\ge c_1\left|z\right|$. Together with the sectorial resolvent estimate $\parallel B_1{\left(z\right)\parallel \le C\left|z\right|}^{-\alpha }$, the first term is bounded by ${C\left|z\right|}^{-\alpha }\tau$. For the second term, we apply the high-order bound $\parallel B_1\left(z\right)-B_2{\left(z\right)\parallel \le C\left|z\right|}^{3-\alpha }{\tau }^{4-\alpha }$ from Lemma 4.2. Note $\left|z\tau \right|\le {\displaystyle \frac{\pi }{sin\theta }}\le C_0$ uniformly on ${\Gamma }_{\tau }$, with $\alpha \in (0,1)$. Thereforce:

$$\left|\right|K_1\left(z\right)-K_2\left(z\right)\left|\right|\le {C\left|z\right|}^{-\alpha }\tau +{C\left|z\right|}^{-\alpha }{\left(\right|z\left|\tau \right)}^{4-\alpha }\le C{\left|z\right|}^{-\alpha }\tau .$$

(90)

This completes the proof.

   □

Theorem 3. 

Let $V\left(t_n\right)$ and $V^n$ be the semidiscrete and fully discrete homogeneous solutions, respectively. Then the total temporal discretization error satisfies $E=V\left(t_n\right)-V^n\le \mathcal{O}\left(t_n^{\alpha -1}\tau \right)$.

Proof. 

The total error is decomposed into $I+II$. By the preceding lemmas, we estimate $II$ directly:

$$\left|\right|II\left|\right|\le {\displaystyle \frac{1}{2\pi }}{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}|·||K_1\left(z\right)-K_2\left(z\right)\left|\right|·\left|\right|v\left|\right|\left|dz\right|\le C\tau \left|\right|v\left|\right|{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{-\alpha }\left|dz\right|.$$

(91)

Using the change of variables $s=z{t}_n$ ($dz=t_n^{-1}ds$), we have:

$${\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{-\alpha }\left|dz\right|=t_n^{\alpha -1}{\int }_{{\Gamma }_{\tau }}|e^s{\left|\right|s|}^{-\alpha }\left|ds\right|.$$

(92)

Since the integral ${\int }_{{\Gamma }_{\tau }}|e^s{\left|\right|s|}^{-\alpha }\left|ds\right|$ converges to a finite constant independent of $t_n$, we obtain $\left|\right|II\left|\right|\le C{t}_n^{\alpha -1}\tau$.

Similarly, the error for term I is bounded by:

$$\left|\right|I\left|\right|\le C\tau \left|\right|v\left|\right|{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{-\alpha }\left|dz\right|.$$

(93)

Applying the identical substitution $s=z{t}_n$, we compute:

$$\left|\right|I\left|\right|\le C\tau \left|\right|v\left|\right|t_n^{\alpha -1}{\int }_{{\Gamma }_{\tau }}|e^s{\left|\right|s|}^{-\alpha }\left|ds\right|\le C{t}_n^{\alpha -1}\tau .$$

(94)

By the triangle inequality $\left|\right|E\left|\right|\le \left|\right|I\left|\right|+\left|\right|II\left|\right|$, the proof is complete.    □

Remark 4. 

The proof shows that although the CQ operator $z_{\tau }^{\left(\alpha \right)}-z^{\alpha }$ in the high-order scheme attains $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ accuracy, the overall truncation error is dominated by the $\mathcal{O}\left(\tau \right)$ error between z and $\mu \left(z\right)$ in the first term. This forms the algebraic bottleneck for further accuracy improvement.

4.2. Corrected Schemes for Nonsmooth Problems

In the generating function framework, the discrete time-stepping operator is encoded in the function

$$\mu \left(z\right)={\displaystyle \frac{e^{z\tau }-1}{\tau }},$$

(95)

which appears in the contour integral representation of the numerical solution. Let $K_1\left(z\right)=z^{-1}{(z^{\alpha }I+A)}^{-1}A$ be the continuous solution kernel and $K_2\left(z\right)={\mu }^{-1}\left(z\right){(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}A$ the kernel induced by the standard L1 discretisation. As established in Lemma 4.3, the discrepancy between these two kernels satisfies

$$\parallel K_1\left(z\right)-K_2{\left(z\right)\parallel \le C\left|z\right|}^{-\alpha }\tau ,\forall z\in {\Gamma }_{\tau }.$$

(96)

The root cause of this saturation is the first-order accuracy of $\mu \left(z\right)$ as an approximation to z, namely $|\mu \left(z\right)-z|=\mathcal{O}\left(\right|z{|}^2\tau )$. To break this barrier we replace $\mu \left(z\right)$ by a corrected prefactor $\overline{\mu }\left(z\right)$ that approximates z to higher order, thereby realigning the accuracy of the two kernel layers in the contour integral representation.

The construction starts from the elementary identity

$$z\tau =ln\left(1+(e^{z\tau }-1)\right)=ln(1+\tau \mu \left(z\right)).$$

(97)

Expanding the logarithm with the Mercator series $ln(1+x)={\sum }_{j=1}^{\infty }{\displaystyle \frac{{(-1)}^{j-1}}{j}}{x}^j$ and substituting $x=\tau \mu \left(z\right)$ yields the formal power series representation

$$z={\displaystyle \frac{1}{\tau }}{\displaystyle \sum _{j=1}^{\infty }}{\displaystyle \frac{{(-1)}^{j-1}}{j}}{\left(\tau \mu \left(z\right)\right)}^j={\displaystyle \sum _{j=1}^{\infty }}{\displaystyle \frac{{(-1)}^{j-1}{\tau }^{j-1}}{j}}{\mu }^j\left(z\right).$$

(98)

Truncating the infinite series at the K-th term defines the order-K correction prefactor

$$\overline{\mu }\left(z\right)={\displaystyle \sum _{j=1}^K}{\displaystyle \frac{{(-1)}^{j-1}{\tau }^{j-1}}{j}}{\mu }^j\left(z\right).$$

(99)

By construction the residual satisfies

$$\overline{\mu }\left(z\right)-z=\mathcal{O}\left({\tau }^K\mu {\left(z\right)}^{K+1}\right)=\mathcal{O}\left(z^{K+1}{\tau }^K\right),$$

(100)

because $\mu \left(z\right)=z+\mathcal{O}\left(z^2\tau \right)$. The integer K controls the formal accuracy: larger K eliminates more low-order error terms in the kernel expansion.

Setting $K=4$ in (99) gives the explicit expression

$$\overline{\mu }\left(z\right)=\mu \left(z\right)-{\displaystyle \frac{\tau }{2}}{\mu }^2\left(z\right)+{\displaystyle \frac{{\tau }^2}{3}}{\mu }^3\left(z\right)-{\displaystyle \frac{{\tau }^3}{4}}{\mu }^4\left(z\right).$$

(101)

For the subsequent analysis it is convenient to introduce the scaled variable

$$w=z\tau ,$$

(102)

so that $\tau \mu \left(z\right)=e^w-1$ and consequently

$$\overline{\mu }\left(z\right)={\displaystyle \frac{1}{\tau }}{\Phi }_K\left(w\right),{\Phi }_K\left(w\right)={\displaystyle \sum _{j=1}^K}{\displaystyle \frac{{(-1)}^{j-1}}{j}}{(e^w-1)}^j.$$

(103)

The function ${\Phi }_K\left(w\right)$ is precisely the K-th order Taylor polynomial of $ln(1+(e^w-1))=w$ centred at $w=0$. For $K=4$, factoring out the common factor $(e^w-1)$ yields

$${\Phi }_4\left(w\right)=(e^w-1)P_4\left(e^w\right),$$

(104)

where $P_4\left(Y\right)$ is the cubic polynomial

$$\begin{array}{cc}\hfill P_4\left(Y\right)& =1-{\displaystyle \frac{Y-1}{2}}+{\displaystyle \frac{{(Y-1)}^2}{3}}-{\displaystyle \frac{{(Y-1)}^3}{4}}\hfill \\ \hfill & ={\displaystyle \frac{1}{12}}\left(25-23Y+13Y^2-3Y^3\right).\hfill \end{array}$$

(105)

On the rays of the integration contour ${\Gamma }_{\tau }^{\pm }$ defined in Lemma 10 below, we have $Re\left(z\tau \right)<0$, hence $|e^{z\tau }|<1$. For the error analysis it is essential that $\overline{\mu }\left(z\right)$ does not vanish on the contour. By the factorisation (104) this is guaranteed if $P_4\left(Y\right)$ has no root inside the closed unit disk $\left|Y\right|\le 1$. The following classical result [15,16] provides a convenient sufficient condition.

Lemma 9. 

[15] [16] Let $Q_0\left(Y\right)={\sum }_{j=0}^n{a}_j^{\left(0\right)}{Y}^j$ be a polynomial of degree n with complex coefficients. For $k=0,1,\dots ,n-1$, define the sequence

$$Q_{k+1}\left(Y\right)=\overline{a_0^{\left(k\right)}}{Q}_k\left(Y\right)-a_{n-k}^{\left(k\right)}{Q}_k^{*}\left(Y\right),$$

(106)

where $Q_k^{*}\left(Y\right)={\sum }_{j=0}^{n-k}\overline{a_{n-k-j}^{\left(k\right)}}{Y}^j$ is the reciprocal polynomial of $Q_k$, and $a_0^{\left(k\right)}$, $a_{n-k}^{\left(k\right)}$ denote respectively the constant term and the leading coefficient of $Q_k$. If the strict inequalities

$$|a_0^{\left(k\right)}|>|a_{n-k}^{\left(k\right)}|,k=0,1,\dots ,n-1,$$

(107)

hold for every step, then all roots of $Q_0\left(Y\right)$ lie strictly outside the unit disk $\{Y\in \mathbb{C}:|Y|\le 1\}$.

We apply Lemma 9 to $Q\left(Y\right)=12P_4\left(Y\right)=-3Y^3+13Y^2-23Y+25$ (clearing the denominator does not affect the root locations). The coefficients are

$$a_0^{\left(0\right)}=25,a_1^{\left(0\right)}=-23,a_2^{\left(0\right)}=13,a_3^{\left(0\right)}=-3.$$

(108)

For $k=0$, $|a_0^{\left(0\right)}|=25>3=|a_3^{\left(0\right)}|$. The reciprocal polynomial is $Q_0^{*}\left(Y\right)=-3+13Y-23Y^2+25Y^3$, hence

$$Q_1\left(Y\right)=25Q_0\left(Y\right)-(-3)Q_0^{*}\left(Y\right)=256Y^2-536Y+616.$$

(109)

For $k=1$, $|a_0^{\left(1\right)}|=616>256=|a_2^{\left(1\right)}|$. With $Q_1^{*}\left(Y\right)=256-536Y+616Y^2$,

$$Q_2\left(Y\right)=616Q_1\left(Y\right)-256Q_1^{*}\left(Y\right)=-192960Y+313920.$$

(110)

For $k=2$, $|a_0^{\left(2\right)}|=313920>192960=|a_1^{\left(2\right)}|$. All three inequalities are satisfied. Lemma 9 therefore guarantees that every root of $Q\left(Y\right)$, and consequently of $P_4\left(Y\right)$, lies strictly outside $\left|Y\right|\le 1$. Thus for all z on the rays ${\Gamma }_{\tau }^{\pm }$ where $|e^{z\tau }|<1$, we have $P_4\left(e^{z\tau }\right)\ne 0$. Since $e^{z\tau }\ne 1$ on these rays, (104) implies ${\Phi }_4\left(z\tau \right)\ne 0$, which is equivalent to $\overline{\mu }\left(z\right)\ne 0$ by (103).

Lemma 10. 

Let Γ be the Hankel-type contour consisting of the two rays

$${\Gamma }^{\pm }=\{r{e}^{\pm i\theta }:r_0\le r\le \pi /sin\theta \},\theta \in (\pi /2,\pi ),$$

(111)

and a connecting circular arc ${\Gamma }^c$ of radius $r_0>0$. Define ${\Gamma }_{\tau }=\{z\in \mathbb{C}:z\tau \in \Gamma \}$. For $K=4$ and every $z\in {\Gamma }_{\tau }$, there exist constants $c,C>0$, independent of τ, such that

$$c\left|z\right|\le |\overline{\mu }\left(z\right)|\le C|z|,$$

(112)

and the approximation error satisfies

$$|\overline{\mu }{\left(z\right)-z|\le C|z|}^5{\tau }^4.$$

(113)

Proof. 

Set $w=z\tau \in \Gamma$. From (103), $\overline{\mu }\left(z\right)={\tau }^{-1}{\Phi }_4\left(w\right)$. Truncating the Mercator series at $j=4$ gives

$$w=ln\left(1+(e^w-1)\right)={\Phi }_4\left(w\right)+{\displaystyle \frac{1}{5}}{(e^w-1)}^5+\mathcal{O}\left({(e^w-1)}^6\right).$$

(114)

Since $e^w-1=w+\mathcal{O}\left(w^2\right)$, we have ${(e^w-1)}^5=w^5+\mathcal{O}\left(w^6\right)$, and therefore

$${\Phi }_4\left(w\right)=w-{\displaystyle \frac{1}{5}}{w}^5+\mathcal{O}\left(w^6\right).$$

(115)

Hence ${\Phi }_4\left(w\right)-w=\mathcal{O}\left(w^5\right)$, which directly implies $|\overline{\mu }\left(z\right)-z|={\tau }^{-1}|{\Phi }_4{\left(w\right)-w|\le C|z|}^5{\tau }^4$, proving (113).

We now turn to (112). On the rays ${\Gamma }_{\tau }^{\pm }$, $Re\left(w\right)<0$ and $|e^w|<1$. As shown above, $P_4\left(e^w\right)\ne 0$ and $e^w\ne 1$, so ${\Phi }_4\left(w\right)\ne 0$. Define $F\left(w\right)={\Phi }_4\left(w\right)/w$ for $w\ne 0$, and $F\left(0\right)={\Phi }_4^{\prime }\left(0\right)=1$. The function F is analytic on the compact set ${\Gamma }^{\pm }$ and never vanishes there. By the extreme value theorem, there exist $c_{\theta },C_{\theta }>0$ with

$$c_{\theta }\le \left|F\left(w\right)\right|\le C_{\theta }\forall w\in {\Gamma }^{\pm }.$$

(116)

Multiplying by $\left|w\right|/\tau =\left|z\right|$ yields (112) on ${\Gamma }_{\tau }^{\pm }$.

On the circular arc ${\Gamma }_{\tau }^c$, $\left|w\right|=r_0$. Because ${lim}_{w\to 0}{\Phi }_4\left(w\right)/w=1$, we may select $r_0\in (0,1)$ sufficiently small so that $|{\Phi }_4\left(w\right)/w-1|\le 1/2$ for all $\left|w\right|\le r_0$. The triangle inequality then gives $\frac{1}{2}\left|w\right|\le |{\Phi }_4\left(w\right)|\le \frac{3}{2}\left|w\right|$, which translates to $c\left|z\right|\le |\overline{\mu }\left(z\right)|\le C|z|$ on ${\Gamma }_{\tau }^c$. Combining the estimates on the three pieces of ${\Gamma }_{\tau }$ completes the proof.    □

The same construction principle applies to any fixed order $K\ge 1$. The polynomial $P_K\left(Y\right)$, defined implicitly by ${\Phi }_K\left(w\right)=(e^w-1)P_K\left(e^w\right)$, can be verified to satisfy the Schur–Cohn criterion for each fixed K by a finite number of elementary operations, guaranteeing that ${\Phi }_K\left(w\right)\ne 0$ on the contour. The resulting prefactor then achieves $|\overline{\mu }\left(z\right)-z|=\mathcal{O}\left(z^{K+1}{\tau }^K\right)$, providing a systematic route to arbitrarily high-order starting correction schemes.

Lemma 11. 

Let ${\tilde{K}}_2\left(z\right)={\overline{\mu }}^{-1}\left(z\right){(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}A$. Assuming the fractional order is $0<\alpha <1$ and the time step satisfies $0<\tau \le T$ for a fixed final time T, for any $z\in {\Gamma }_{\tau }$, there exists a constant $C>0$ such that

$$\left|\right|K_1\left(z\right)-{\tilde{K}}_2\left(z\right)\left|\right|\le C{\tau }^{4-\alpha }{\left(\right|z|}^{2-\alpha }+{\left|z\right|}^{3-\alpha }).$$

(117)

Proof. 

By subtracting the operators and inserting a cross term, the norm decomposes as:

$$\left|\right|K_1\left(z\right)-{\tilde{K}}_2\left(z\right)\left|\right|\le {\displaystyle \frac{|z-\overline{\mu }\left(z\right)|}{\left|z\right||\overline{\mu }\left(z\right)|}}\left|\right|B_1\left(z\right)\left|\right|+|\overline{\mu }{\left(z\right)|}^{-1}\left|\right|B_1\left(z\right)-{\tilde{B}}_2\left(z\right)\left|\right|.$$

(118)

From the mapping bounds : $|\overline{\mu }{\left(z\right)-z|\le C|z|}^5{\tau }^4$ and $|\overline{\mu }\left(z\right)|\ge c|z|$. Given the spatial resolvent norm $\left|\right|B_1\left(z\right)\left|\right|\le {C\left|z\right|}^{-\alpha }$, the first term yields:

$${\displaystyle \frac{|z-\overline{\mu }\left(z\right)|}{\left|z\right||\overline{\mu }\left(z\right)|}}\left|\right|B_1\left(z\right)\left|\right|\le {\displaystyle \frac{{C\left|z\right|}^5{\tau }^4}{{c\left|z\right|}^2}}·{C\left|z\right|}^{-\alpha }=C{\tau }^4{\left|z\right|}^{3-\alpha }.$$

(119)

From the spatial approximation bound, $\left|\right|B_1\left(z\right)-{\tilde{B}}_2\left(z\right)\left|\right|\le C{\tau }^{4-\alpha }{\left|z\right|}^{3-\alpha }$. For the second term, we have

$$|\overline{\mu }{\left(z\right)|}^{-1}\left|\right|B_1\left(z\right)-{\tilde{B}}_2\left(z\right)\left|\right|\le c^{-1}{\left|z\right|}^{-1}·C{\tau }^{4-\alpha }{\left|z\right|}^{4-\alpha }=C{\tau }^{4-\alpha }{\left|z\right|}^{3-\alpha }.$$

(120)

Summing these bounds completes the proof.

   □

Theorem 4. 

Let ${\tilde{V}}_h^n$ be the fully discrete homogeneous solution of the corrected scheme.For nonsmooth initial data $v\in L^2\left(\Omega \right)$, the $L^2$-norm temporal discretization error $\tilde{E}=V\left(t_n\right)-{\tilde{V}}^n$ satisfies:

$$\left|\right|\tilde{E}\left|\right|\le C\left(t_n^{-4}{\tau }^4+t_n^{\alpha -4}{\tau }^{4-\alpha }\right)\left|\right|v\left|\right|.$$

(121)

Proof. 

The total error in the $L^2$-norm is decomposed into the outer and inner contour integral components $I+II$.

For the inner contour integral $II$, substituting the kernel difference bound from Lemma 4.2 yields

$$\left|\right|II\left|\right|\le {\displaystyle \frac{1}{2\pi }}{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}|·||K_1\left(z\right)-{\tilde{K}}_2\left(z\right)\left|\right|·\left|\right|v\left|\right|\left|dz\right|\le C{\tau }^{4-\alpha }\left|\right|v\left|\right|{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{3-\alpha }\left|dz\right|.$$

(122)

Utilizing the properties of the Gamma function, the above integral can be evaluated as

$${\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{3-\alpha }\left|dz\right|\le C{\int }_0^{\infty }{e}^{-r{t}_n|cos\theta |}{r}^{3-\alpha }dr=C{t}_n^{\alpha -3}.$$

(123)

This derives $\left|\right|II\left|\right|\le C{t}_n^{\alpha -4}{\tau }^{4-\alpha }\left|\right|v\left|\right|$.

For the outer contour integral I, on ${\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\tau }$ we have $\left|z\right|=r\ge {\displaystyle \frac{\pi }{\tau sin\theta }}$. Extracting the Kth-order temporal scaling term gives $r^{-4}\le C{\tau }^4$. Combining with the $L^2$ resolvent estimate $\parallel K_1{\left(z\right)\parallel \le C\left|z\right|}^{-1}$ and matching orders yields

$$\left|\right|I\left|\right|\le C\left|\right|v\left|\right|{\int }_{{\displaystyle \frac{\pi }{\tau sin\theta }}}^{\infty }{e}^{-r{t}_n|cos\theta |}{r}^{-1}dr=C\left|\right|v\left|\right|{\int }_{{\displaystyle \frac{\pi }{\tau sin\theta }}}^{\infty }{e}^{-r{t}_n|cos\theta |}{r}^3·r^{-4}dr.$$

(124)

The lower integration limit to 0 gives a strict upper bound:

$$\left|\right|I\left|\right|\le C{\tau }^4\left|\right|v\left|\right|{\int }_0^{\infty }{e}^{-r{t}_n|cos\theta |}{r}^{3}dr=C{t}_n^{-4}{\tau }^4\left|\right|v\left|\right|.$$

(125)

By the triangle inequality $\left|\right|\tilde{E}\left|\right|\le \left|\right|I\left|\right|+\left|\right|II\left|\right|$, and recognizing that for the given fractional parameter $\alpha \in (0,1)$ the term ${\tau }^{4-\alpha }$ dominates the overall truncation error, the theorem is proved. □

   □

4.3. Corrected Scheme for Smooth Cases

Having discussed the $L^2$ space estimates for non-smooth initial data, we now address the problem under smooth initial conditions. Here, we define $u\left(0\right)=v\in D\left(A\right)=H^2\left(\Omega \right)\cap H_0^1\left(\Omega \right)$.

Lemma 12. 

Let the continuous and discrete smooth kernel functions be defined as $K_1^s\left(z\right)=-z^{-1}{(z^{\alpha }I+A)}^{-1}$ and $K_2^s\left(z\right)=-{\overline{\mu }}^{-1}\left(z\right){(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}$, respectively. For any $z\in {\Gamma }_{\tau }$, there exists a constant $C>0$ such that:

$$\left|\right|K_1^s\left(z\right)-K_2^s\left(z\right)\left|\right|\le C{\tau }^{4-\alpha }{\left|z\right|}^{3-2\alpha }.$$

(126)

Proof. 

Utilizing operator difference, the norm is bounded by:

$$\left|\right|K_1^s\left(z\right)-K_2^s\left(z\right)\left|\right|\le |z^{-1}-{\overline{\mu }}^{-1}\left(z\right)|·||B_1^s\left(z\right)\left|\right|+|{\overline{\mu }}^{-1}\left(z\right)|·||B_1^s\left(z\right)-B_2^s\left(z\right)\left|\right|,$$

(127)

where $B_1^s\left(z\right)={(z^{\alpha }I+A)}^{-1}$ and $B_2^s\left(z\right)={(z_{\tau }^{\left(\alpha \right)}\left(z\right)I+A)}^{-1}$. For the first term, from previous bounds, it is easy to obtain

$$|z^{-1}-{\overline{\mu }}^{-1}\left(z\right)|·||B_1^s\left(z\right)\left|\right|\le {C\left|z\right|}^{3-\alpha }{\tau }^5.$$

(128)

For the second term, applying the resolvent identity along with the operator estimates yields:

$$\begin{array}{cc}\hfill |{\overline{\mu }}^{-1}\left(z\right)|·||B_1^s\left(z\right)-B_2^s\left(z\right)\left|\right|& \le {C\left|z\right|}^{-1}·\left|\right|B_1^s\left(z\right)\left|\right|·\left|\right|B_2^s\left(z\right)\left|\right|·|z_{\tau }^{\left(\alpha \right)}\left(z\right)-z^{\alpha }|\hfill \\ \hfill & \le {C\left|z\right|}^{-1}·{\left|z\right|}^{-\alpha }·{\left|z\right|}^{-\alpha }·{\left|z\right|}^K{\tau }^{4-\alpha }\hfill \\ \hfill & =C{\tau }^{4-\alpha }{\left|z\right|}^{3-2\alpha }.\hfill \end{array}$$

Combining both terms strictly derives $\left|\right|K_1^s\left(z\right)-K_2^s\left(z\right)\left|\right|\le C{\tau }^{4-\alpha }{\left|z\right|}^{3-2\alpha }$.    □

Theorem 5. 

Let $\tilde{V}\left(t_n\right)$ be the numerical solution under the corrected scheme. For smooth initial data $v\in D\left(A\right)$ with $U^0=v=R_{h}v$, the total error ${\tilde{E}}^s=V\left(t_n\right)-{\tilde{V}}^n$ measured in the $L^2$-norm satisfies:

$$\left|\right|{\tilde{E}}^s\left|\right|\le C\left(t_n^{2\alpha -4}{\tau }^4+t_n^{2\alpha -4}{\tau }^{4-\alpha }\right)\left|\right|Av\left|\right|.$$

(129)

Proof. 

The error is decomposed into the outer and inner contour integrals $I^s+I{I}^s$:

$${\tilde{E}}^s={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\tau }}{e}^{z{t}_n}{K}_1^s\left(z\right)Avdz+{\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}(K_1^s\left(z\right)-K_2^s\left(z\right))Avdz.$$

(130)

Analogous to the proof strategy for $\left|\right|\tilde{E}\left|\right|$ in Theorem 4.5, the inner integral $I{I}^s$ can be bounded as

$$\begin{array}{cc}\hfill \left|\right|I{I}^s\left|\right|& \le {\displaystyle \frac{1}{2\pi }}{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}|·||K_1^s\left(z\right)-K_2^s\left(z\right)\left|\right|·\left|\right|Av\left|\right|\left|dz\right|\hfill \\ \hfill & \le C{\tau }^{4-\alpha }\left|\right|Av\left|\right|{\int }_{{\Gamma }_{\tau }}|e^{z{t}_n}{\left|\right|z|}^{4-2\alpha }\left|dz\right|\hfill \\ \hfill & \le C{t}_n^{2\alpha -4}{\tau }^{4-\alpha }\left|\right|Av\left|\right|.\hfill \end{array}$$

Similarly, for the outer integral $I^s$, extracting the 4th-order scaling term and using smooth-state estimates gives

$$\parallel I^s\parallel \le C{t}_n^{2\alpha -4}{\tau }^4\parallel Av\parallel .$$

(131)

Combining bounds for $I^s$ and $I{I}^s$ via the triangle inequality, the scheme attains the optimal temporal convergence order $\mathcal{O}\left({\tau }^{4-\alpha }\right)$.    □

Remark 5. 

Theorems 4 and 5 reveal the core difference in error behavior of the corrected scheme between nonsmooth initial data $v\in L^2\left(\Omega \right)$ and smooth initial data $v\in D\left(A\right)$. With higher spatial regularity, the singularity of the temporal discretization error as $t_n\to 0$ weakens from $t_n^{\alpha -4}$ to $t_n^{2\alpha -4}$. The singularity persists even for smooth initial data, reflecting the limited smoothing effect of the fractional derivative operator.

The high-order discrete symbol $\overline{\mu }\left(z\right)$ is constructed to approximate the continuous fractional operator $z^{-1}$ with $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ accuracy. For the homogeneous problem, the equivalent constant source term has a low-frequency truncation deficit in the discrete Laplace domain. To correct this, we perform high-order asymptotic expansion of the discrete generating function as $\omega =z\tau \to 0$.

Table 2. Starting correction coefficients ${\overline{c_j}}^{\left(K\right)}$ for the K-th order corrected $\overline{\mu }$ schemes.

K	${\overline{{\mathit{c}}_1}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_2}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_3}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_4}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_5}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_6}}^{\left(\mathit{K}\right)}$
$K=4$	${\displaystyle \frac{31}{24}}$	$-{\displaystyle \frac{7}{6}}$	${\displaystyle \frac{3}{8}}$
$K=5$	${\displaystyle \frac{1181}{720}}$	$-{\displaystyle \frac{177}{80}}$	${\displaystyle \frac{341}{240}}$	$-{\displaystyle \frac{251}{720}}$
$K=6$	${\displaystyle \frac{2837}{1440}}$	$-{\displaystyle \frac{2543}{720}}$	${\displaystyle \frac{17}{5}}$	$-{\displaystyle \frac{1201}{720}}$	${\displaystyle \frac{95}{288}}$
$K=7$	${\displaystyle \frac{138241}{60480}}$	$-{\displaystyle \frac{309047}{60480}}$	${\displaystyle \frac{198251}{30240}}$	$-{\displaystyle \frac{145877}{30240}}$	${\displaystyle \frac{23077}{12096}}$	$-{\displaystyle \frac{19087}{60480}}$

Table 2. Starting correction coefficients ${\overline{c_j}}^{\left(K\right)}$ for the K-th order corrected $\overline{\mu }$ schemes.

K	${\overline{{\mathit{c}}_1}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_2}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_3}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_4}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_5}}^{\left(\mathit{K}\right)}$	${\overline{{\mathit{c}}_6}}^{\left(\mathit{K}\right)}$
$K=4$	${\displaystyle \frac{31}{24}}$	$-{\displaystyle \frac{7}{6}}$	${\displaystyle \frac{3}{8}}$
$K=5$	${\displaystyle \frac{1181}{720}}$	$-{\displaystyle \frac{177}{80}}$	${\displaystyle \frac{341}{240}}$	$-{\displaystyle \frac{251}{720}}$
$K=6$	${\displaystyle \frac{2837}{1440}}$	$-{\displaystyle \frac{2543}{720}}$	${\displaystyle \frac{17}{5}}$	$-{\displaystyle \frac{1201}{720}}$	${\displaystyle \frac{95}{288}}$
$K=7$	${\displaystyle \frac{138241}{60480}}$	$-{\displaystyle \frac{309047}{60480}}$	${\displaystyle \frac{198251}{30240}}$	$-{\displaystyle \frac{145877}{30240}}$	${\displaystyle \frac{23077}{12096}}$	$-{\displaystyle \frac{19087}{60480}}$

Table 3. Starting correction coefficients $d_{\ell ,n}^{\left(K\right)}$ for high-order Zeta schemes.

K	ℓ	${\mathit{d}}_{\mathit{\ell },1}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },2}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },3}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },4}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },5}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },6}^{\left(\mathit{K}\right)}$
$K=4$	$\ell =1$	${\displaystyle \frac{31}{24}}$	$-{\displaystyle \frac{7}{3}}$	${\displaystyle \frac{9}{8}}$
	$\ell =2$	${\displaystyle \frac{31}{48}}$	$-{\displaystyle \frac{7}{3}}$	${\displaystyle \frac{27}{16}}$
	$\ell =3$	${\displaystyle \frac{31}{144}}$	$-{\displaystyle \frac{14}{9}}$	${\displaystyle \frac{27}{16}}$
$K=5$	$\ell =1$	${\displaystyle \frac{1181}{720}}$	$-{\displaystyle \frac{177}{40}}$	${\displaystyle \frac{341}{80}}$	$-{\displaystyle \frac{251}{180}}$
	$\ell =2$	${\displaystyle \frac{1181}{1440}}$	$-{\displaystyle \frac{177}{40}}$	${\displaystyle \frac{1023}{160}}$	$-{\displaystyle \frac{251}{90}}$
	$\ell =3$	${\displaystyle \frac{1181}{4320}}$	$-{\displaystyle \frac{59}{20}}$	${\displaystyle \frac{1023}{160}}$	$-{\displaystyle \frac{502}{135}}$
	$\ell =4$	${\displaystyle \frac{1181}{17280}}$	$-{\displaystyle \frac{59}{40}}$	${\displaystyle \frac{3069}{640}}$	$-{\displaystyle \frac{502}{135}}$
$K=6$	$\ell =1$	${\displaystyle \frac{2837}{1440}}$	$-{\displaystyle \frac{2543}{360}}$	${\displaystyle \frac{51}{5}}$	$-{\displaystyle \frac{1201}{180}}$	${\displaystyle \frac{475}{288}}$
	$\ell =2$	${\displaystyle \frac{2837}{2880}}$	$-{\displaystyle \frac{2543}{360}}$	${\displaystyle \frac{153}{10}}$	$-{\displaystyle \frac{1201}{90}}$	${\displaystyle \frac{2375}{576}}$
	$\ell =3$	${\displaystyle \frac{2837}{8640}}$	$-{\displaystyle \frac{2543}{540}}$	${\displaystyle \frac{153}{10}}$	$-{\displaystyle \frac{2402}{135}}$	${\displaystyle \frac{11875}{1728}}$
	$\ell =4$	${\displaystyle \frac{2837}{34560}}$	$-{\displaystyle \frac{2543}{1080}}$	${\displaystyle \frac{459}{40}}$	$-{\displaystyle \frac{2402}{135}}$	${\displaystyle \frac{59375}{6912}}$
	$\ell =5$	${\displaystyle \frac{2837}{172800}}$	$-{\displaystyle \frac{2543}{2700}}$	${\displaystyle \frac{1377}{200}}$	$-{\displaystyle \frac{9608}{675}}$	${\displaystyle \frac{59375}{6912}}$
$K=7$	$\ell =1$	${\displaystyle \frac{138241}{60480}}$	$-{\displaystyle \frac{309047}{30240}}$	${\displaystyle \frac{198251}{10080}}$	$-{\displaystyle \frac{145877}{7560}}$	${\displaystyle \frac{115385}{12096}}$	$-{\displaystyle \frac{19087}{10080}}$
	$\ell =2$	${\displaystyle \frac{138241}{120960}}$	$-{\displaystyle \frac{309047}{30240}}$	${\displaystyle \frac{198251}{6720}}$	$-{\displaystyle \frac{145877}{3780}}$	${\displaystyle \frac{576925}{24192}}$	$-{\displaystyle \frac{19087}{3360}}$
	$\ell =3$	${\displaystyle \frac{138241}{362880}}$	$-{\displaystyle \frac{309047}{45360}}$	${\displaystyle \frac{198251}{6720}}$	$-{\displaystyle \frac{145877}{2835}}$	${\displaystyle \frac{2884625}{72576}}$	$-{\displaystyle \frac{19087}{1680}}$
	$\ell =4$	${\displaystyle \frac{138241}{1451520}}$	$-{\displaystyle \frac{309047}{90720}}$	${\displaystyle \frac{198251}{8960}}$	$-{\displaystyle \frac{145877}{2835}}$	${\displaystyle \frac{14423125}{290304}}$	$-{\displaystyle \frac{19087}{1120}}$
	$\ell =5$	${\displaystyle \frac{138241}{7257600}}$	$-{\displaystyle \frac{309047}{226800}}$	${\displaystyle \frac{594753}{44800}}$	$-{\displaystyle \frac{583508}{14175}}$	${\displaystyle \frac{14423125}{290304}}$	$-{\displaystyle \frac{57261}{2800}}$
	$\ell =6$	${\displaystyle \frac{138241}{43545600}}$	$-{\displaystyle \frac{309047}{680400}}$	${\displaystyle \frac{594753}{89600}}$	$-{\displaystyle \frac{1167016}{42525}}$	${\displaystyle \frac{72115625}{1741824}}$	$-{\displaystyle \frac{57261}{2800}}$

Table 3. Starting correction coefficients $d_{\ell ,n}^{\left(K\right)}$ for high-order Zeta schemes.

K	ℓ	${\mathit{d}}_{\mathit{\ell },1}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },2}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },3}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },4}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },5}^{\left(\mathit{K}\right)}$	${\mathit{d}}_{\mathit{\ell },6}^{\left(\mathit{K}\right)}$
$K=4$	$\ell =1$	${\displaystyle \frac{31}{24}}$	$-{\displaystyle \frac{7}{3}}$	${\displaystyle \frac{9}{8}}$
	$\ell =2$	${\displaystyle \frac{31}{48}}$	$-{\displaystyle \frac{7}{3}}$	${\displaystyle \frac{27}{16}}$
	$\ell =3$	${\displaystyle \frac{31}{144}}$	$-{\displaystyle \frac{14}{9}}$	${\displaystyle \frac{27}{16}}$
$K=5$	$\ell =1$	${\displaystyle \frac{1181}{720}}$	$-{\displaystyle \frac{177}{40}}$	${\displaystyle \frac{341}{80}}$	$-{\displaystyle \frac{251}{180}}$
	$\ell =2$	${\displaystyle \frac{1181}{1440}}$	$-{\displaystyle \frac{177}{40}}$	${\displaystyle \frac{1023}{160}}$	$-{\displaystyle \frac{251}{90}}$
	$\ell =3$	${\displaystyle \frac{1181}{4320}}$	$-{\displaystyle \frac{59}{20}}$	${\displaystyle \frac{1023}{160}}$	$-{\displaystyle \frac{502}{135}}$
	$\ell =4$	${\displaystyle \frac{1181}{17280}}$	$-{\displaystyle \frac{59}{40}}$	${\displaystyle \frac{3069}{640}}$	$-{\displaystyle \frac{502}{135}}$
$K=6$	$\ell =1$	${\displaystyle \frac{2837}{1440}}$	$-{\displaystyle \frac{2543}{360}}$	${\displaystyle \frac{51}{5}}$	$-{\displaystyle \frac{1201}{180}}$	${\displaystyle \frac{475}{288}}$
	$\ell =2$	${\displaystyle \frac{2837}{2880}}$	$-{\displaystyle \frac{2543}{360}}$	${\displaystyle \frac{153}{10}}$	$-{\displaystyle \frac{1201}{90}}$	${\displaystyle \frac{2375}{576}}$
	$\ell =3$	${\displaystyle \frac{2837}{8640}}$	$-{\displaystyle \frac{2543}{540}}$	${\displaystyle \frac{153}{10}}$	$-{\displaystyle \frac{2402}{135}}$	${\displaystyle \frac{11875}{1728}}$
	$\ell =4$	${\displaystyle \frac{2837}{34560}}$	$-{\displaystyle \frac{2543}{1080}}$	${\displaystyle \frac{459}{40}}$	$-{\displaystyle \frac{2402}{135}}$	${\displaystyle \frac{59375}{6912}}$
	$\ell =5$	${\displaystyle \frac{2837}{172800}}$	$-{\displaystyle \frac{2543}{2700}}$	${\displaystyle \frac{1377}{200}}$	$-{\displaystyle \frac{9608}{675}}$	${\displaystyle \frac{59375}{6912}}$
$K=7$	$\ell =1$	${\displaystyle \frac{138241}{60480}}$	$-{\displaystyle \frac{309047}{30240}}$	${\displaystyle \frac{198251}{10080}}$	$-{\displaystyle \frac{145877}{7560}}$	${\displaystyle \frac{115385}{12096}}$	$-{\displaystyle \frac{19087}{10080}}$
	$\ell =2$	${\displaystyle \frac{138241}{120960}}$	$-{\displaystyle \frac{309047}{30240}}$	${\displaystyle \frac{198251}{6720}}$	$-{\displaystyle \frac{145877}{3780}}$	${\displaystyle \frac{576925}{24192}}$	$-{\displaystyle \frac{19087}{3360}}$
	$\ell =3$	${\displaystyle \frac{138241}{362880}}$	$-{\displaystyle \frac{309047}{45360}}$	${\displaystyle \frac{198251}{6720}}$	$-{\displaystyle \frac{145877}{2835}}$	${\displaystyle \frac{2884625}{72576}}$	$-{\displaystyle \frac{19087}{1680}}$
	$\ell =4$	${\displaystyle \frac{138241}{1451520}}$	$-{\displaystyle \frac{309047}{90720}}$	${\displaystyle \frac{198251}{8960}}$	$-{\displaystyle \frac{145877}{2835}}$	${\displaystyle \frac{14423125}{290304}}$	$-{\displaystyle \frac{19087}{1120}}$
	$\ell =5$	${\displaystyle \frac{138241}{7257600}}$	$-{\displaystyle \frac{309047}{226800}}$	${\displaystyle \frac{594753}{44800}}$	$-{\displaystyle \frac{583508}{14175}}$	${\displaystyle \frac{14423125}{290304}}$	$-{\displaystyle \frac{57261}{2800}}$
	$\ell =6$	${\displaystyle \frac{138241}{43545600}}$	$-{\displaystyle \frac{309047}{680400}}$	${\displaystyle \frac{594753}{89600}}$	$-{\displaystyle \frac{1167016}{42525}}$	${\displaystyle \frac{72115625}{1741824}}$	$-{\displaystyle \frac{57261}{2800}}$

5. Error Analysis for the Inhomogeneous Problem

We now extend the corrected scheme to the inhomogeneous equation,where $g\left(t\right)=f\left(t\right)-Av$ and $v=u\left(0\right)$. The homogeneous problem treated in the previous sections corresponds to the special case $g\left(t\right)\equiv -Av$, i.e., a constant source.

The correction coefficients ${\overline{c}}_n^{\left(4\right)}(n=1,2,3)$ were determined in Section 4.1 by requiring that the homogeneous scheme be exact for the leading singular modes of the continuous solution. For the inhomogeneous equation the source $g\left(t\right)$ is no longer constant; its variation in time introduces additional low-order errors that must also be compensated. The key observation is that, by linearity, the correction for a variable source can be obtained by Taylor-expanding $g\left(t\right)$ about $t=0$ and correcting each monomial term individually.

Assume that g is sufficiently smooth and write

$$g\left(t\right)=g\left(0\right)+t{g}^{\prime }\left(0\right)+{\displaystyle \frac{t^2}{2}}{g}^{\prime \prime }\left(0\right)+{\displaystyle \frac{t^3}{6}}{g}^{\prime \prime \prime }\left(0\right)+R\left(t\right),$$

(132)

where the remainder $R\left(t\right)=\mathcal{O}\left(t^4\right)$ does not affect the leading-order error and requires no correction. Evaluating (132) at the discrete time levels $t_n=n\tau$ gives

$$g\left(t_n\right)=g\left(0\right)+n\tau g^{\prime }\left(0\right)+{\displaystyle \frac{n^2{\tau }^2}{2}}{g}^{\prime \prime }\left(0\right)+{\displaystyle \frac{n^3{\tau }^3}{6}}{g}^{\prime \prime \prime }\left(0\right)+R\left(t_n\right).$$

(133)

The constant term $g\left(0\right)$ is precisely of the same form as the homogeneous source $-Av$; its correction is therefore identical to the homogeneous case. The terms involving $g^{\prime }\left(0\right),g^{\prime \prime }\left(0\right),g^{\prime \prime \prime }\left(0\right)$ are new and require separate treatment.

Consider the contribution of the ℓ-th derivative term ${\displaystyle \frac{t^{\ell }}{\ell !}}{g}^{(\ell )}\left(0\right)$ to the right-hand side. At time $t_n$ this term evaluates to ${\displaystyle \frac{n^{\ell }{\tau }^{\ell }}{\ell !}}{g}^{(\ell )}\left(0\right)$. In the continuous equation, the response of the system to a source proportional to $t^{\ell }$ is $\ell !$ times weaker than its response to a constant source of the same pointwise magnitude. Consequently, the correction needed for this term inherits the same relative factor ${\overline{c}}_n^{\left(4\right)}$ as the constant-source correction, weighted by the Taylor coefficient ${\displaystyle \frac{n^{\ell }}{\ell !}}$. Formally, matching the discrete correction term ${\tau }^{\ell }{d}_{\ell ,n}^{\left(4\right)}{g}^{(\ell )}\left(0\right)$ with the Taylor contribution ${\displaystyle \frac{{\left(n\tau \right)}^{\ell }}{\ell !}}{\overline{c}}_n^{\left(4\right)}{g}^{(\ell )}\left(0\right)$ yields the relation

$$d_{\ell ,n}^{\left(4\right)}={\displaystyle \frac{n^{\ell }}{\ell !}}{\overline{c}}_n^{\left(4\right)},\ell =1,2,3,n=1,2,3.$$

(134)

For $n\ge 4$, ${\overline{c}}_n^{\left(4\right)}=0$ and consequently $d_{\ell ,n}^{\left(4\right)}=0$ for all ℓ; no correction is applied beyond the third time step.

Theorem 6. 

Let $V^n$ be the numerical solution of the corrected fully discrete scheme and $V\left(t_n\right)$ the exact continuous solution. Suppose $f\in C^4(0,T;L^2\left(\Omega \right))$ (which implies $g_h\in C^3([0,T];L^2\left(\Omega \right))$) and ${\int }_0^{t_n}{(t_n-s)}^{2\alpha -1}\parallel g^{\left(4\right)}\left(s\right)\parallel ds<\infty$ .

There exists a constant $C>0$, independent of τ and h, subject to the initial condition $u\left(0\right)=v\in D\left(A\right)$ (such that $Av\in L^2\left(\Omega \right)$), such that for all $t_n>0$

$$\begin{array}{cc}\hfill & \parallel V^n-V\left(t_n\right)\parallel \hfill \\ & \le C{\tau }^{4-\alpha }\left(t_n^{2\alpha -4}\parallel -Av+g\left(0\right)\parallel +{\displaystyle \sum _{\ell =1}^3}{t}_n^{2\alpha -5+\ell }\parallel g^{(\ell )}\left(0\right)\parallel +{\int }_0^{t_n}{(t_n-s)}^{2\alpha -1}\parallel g^{\left(4\right)}\left(s\right)\parallel ds\right).\hfill \end{array}$$

(135)

Proof. 

The Taylor expansion of the total source term $g\left(t\right)=-Av+f\left(t\right)$ at $t=0$ gives

$$g\left(t\right)={\displaystyle \sum _{\ell =0}^2}{\displaystyle \frac{t^{\ell }}{\ell !}}{g}^{(\ell )}\left(0\right)+R_3\left(t\right),$$

(136)

where $R_3\left(t\right)={\displaystyle \frac{t^3}{3!}}{g}^{\left(3\right)}\left(0\right)+{\int }_0^t{\displaystyle \frac{{(t-s)}^3}{3!}}{g}^{\left(4\right)}\left(s\right)ds$.

Using the Laplace transform identity $\mathcal{L}\{t^{\ell }/\ell !\}=z^{-\ell -1}$, we obtain

$${\widehat{g}}_h\left(z\right)={\displaystyle \sum _{\ell =0}^2}{z}^{-\ell -1}{g}^{(\ell )}\left(0\right)+{\widehat{R}}_3\left(z\right),$$

(137)

where ${\widehat{R}}_3\left(z\right)=z^{-4}{g}^{\left(3\right)}\left(0\right)+\mathcal{L}\left\{{\int }_0^t{\displaystyle \frac{{(t-s)}^3}{3!}}{g}^{\left(4\right)}\left(s\right)ds\right\}$.

The exact continuous solution is given by the inverse Laplace transform

$$V\left(t_n\right)={\displaystyle \frac{1}{2\pi i}}{\int }_{\Gamma }{e}^{z{t}_n}{(z^{\alpha }I+A)}^{-1}\left[{\displaystyle \sum _{\ell =0}^2}{z}^{-\ell -1}{g}^{(\ell )}\left(0\right)+{\widehat{R}}_4\left(z\right)\right]dz.$$

(138)

Applying the discrete Laplace transform ${\widehat{V}}_{\tau }\left(z\right)={\sum }_{n=1}^{\infty }{V}^n{\zeta }^n$ (where $\xi =e^{-z\tau }$) to the fully discrete scheme yields $(z_{\tau }^{\left(\alpha \right)}I+A){\widehat{V}}_{\tau }\left(z\right)={\widehat{g}}_{\tau }\left(z\right)$. To incorporate the starting corrections, we introduce the generating functions ${\mathcal{M}}_0\left(\zeta \right)={\displaystyle \frac{\zeta }{1-\zeta }}+{\sum }_{n=1}^3{\overline{c_n}}^{\left(4\right)}{\zeta }^n$ and ${\mathcal{M}}_{\ell }\left(\zeta \right)={\sum }_{n=1}^{\infty }{\displaystyle \frac{n^{\ell }}{\ell !}}{\zeta }^n+{\sum }_{n=1}^3{d}_{\ell ,n}^{\left(4\right)}{\zeta }^n$, which satisfy the asymptotic relation ${\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(e^{-z\tau }\right)=z^{-\ell -1}+\mathcal{O}({\tau }^4{\left|z\right|}^{3-\ell })$. Denoting the discrete transform of the remainder as ${\widehat{R}}_{3,\tau }\left(z\right)={\sum }_{n=1}^{\infty }{R}_3\left(t_n\right){\zeta }^n$, the expanded source term in the frequency domain is

$${\widehat{g}}_{\tau }\left(z\right)={\displaystyle \sum _{\ell =0}^2}{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)g^{(\ell )}\left(0\right)+\tau {\widehat{R}}_{3,\tau }\left(z\right).$$

(139)

The frequency-domain solution is

$${\widehat{V}}_{\tau }\left(z\right)={(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}\left[{\displaystyle \sum _{\ell =0}^2}{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)g^{(\ell )}\left(0\right)+\tau {\widehat{R}}_{3,\tau }\left(z\right)\right].$$

(140)

Applying the inverse discrete Laplace transform over the contour ${\Gamma }_{\tau }$, the numerical solution is recovered as

$$V^n={\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}{(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}\left[{\displaystyle \sum _{\ell =0}^2}{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)g^{(\ell )}\left(0\right)+\tau {\widehat{R}}_{3,\tau }\left(z\right)\right]dz.$$

(141)

Thus,the error between the discrete solution and the exact solution is

$$|V^n-V\left(t_n\right)|=\left[{\displaystyle \sum _{\ell =0}^2}{I}_{\ell }+I_R\right]$$

For each $\ell \in (0,1\dots 3)$, $I_{\ell }$ is split into the tail error outside ${\Gamma }_{\tau }$ and the approximation error inside ${\Gamma }_{\tau }$

$$\begin{array}{cc}\hfill I_{\ell }=& {\displaystyle \frac{1}{2\pi i}}{\int }_{\Gamma \setminus {\Gamma }_{\tau }}{e}^{z{t}_n}{(z^{\alpha }I+A)}^{-1}{z}^{-\ell -1}{g}^{(\ell )}\left(0\right)dz\hfill \\ \hfill & +{\displaystyle \frac{1}{2\pi i}}{\int }_{{\Gamma }_{\tau }}{e}^{z{t}_n}\left[{(z^{\alpha }I+A)}^{-1}{z}^{-\ell -1}-{(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)\right]g^{(\ell )}\left(0\right)dz.\hfill \end{array}$$

(142)

The tail error is bounded by $C{\tau }^{4-\alpha }{t}_n^{2\alpha -4+\ell }\parallel g^{(\ell )}\left(0\right)\parallel$. Inside ${\Gamma }_{\tau }$, using the operator difference bound from lemma 4.7 $\parallel {(z^{\alpha }I+A)}^{-1}-{(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}\parallel \le C{\tau }^{4-\alpha }{\left|z\right|}^{4-2\alpha }$ yields

$$\begin{array}{cc}\hfill & ∥{(z^{\alpha }I+A)}^{-1}{z}^{-\ell -1}-{(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)∥\hfill \\ \hfill & \le ∥{(z^{\alpha }I+A)}^{-1}-{(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}∥{\left|z\right|}^{-\ell -1}+\parallel {(z_{\tau }^{\left(\alpha \right)}I+A)}^{-1}\parallel ∥z^{-\ell -1}-{\tau }^{\ell +1}{\mathcal{M}}_{\ell }\left(\zeta \right)∥\hfill \\ \hfill & \le C{\tau }^{4-\alpha }{\left|z\right|}^{3-\ell -2\alpha }.\hfill \end{array}$$

(143)

By analogy with the proof of $\left|\right|I\left|\right|$ in Theorem 4.2, integrating this bound over the ray portions of ${\Gamma }_{\tau }$ and applying the change of variables $\rho =r{t}_n$ ($dr=d\rho /t_n$) yields:

$$C{\tau }^{4-\alpha }{\int }_{1/t_n}^{{\displaystyle \frac{\pi }{\tau sin\theta }}}{e}^{-cr{t}_n}{r}^{3-\ell -2\alpha }dr=C{\tau }^{4-\alpha }{t}_n^{2\alpha -4+\ell }{\int }_1^{\infty }{e}^{-c\rho }{\rho }^{3-\ell -2\alpha }d\rho \le C{\tau }^{4-\alpha }{t}_n^{2\alpha -4+\ell }.$$

(144)

For the remainder $I_R$, we decompose $R_3\left(t\right)$ into a polynomial part and a convolution part: $R_3\left(t\right)=R_3^1\left(t\right)+R_3^2\left(t\right)$, where $R_3^1\left(t\right)={\displaystyle \frac{t^3}{3!}}{g}^{\left(3\right)}\left(0\right)$ and $R_3^2\left(t\right)={\int }_0^t{\displaystyle \frac{{(t-s)}^3}{3!}}{g}^{\left(4\right)}\left(s\right)ds$. Correspondingly, the error is split as $I_R=I_3+I_c$.

For the polynomial term $R_3^1\left(t\right)$, integration over the $z^{-4}$ pole ($\ell =3$) yields[3]:

$$\parallel I_3\parallel \le C{\tau }^{4-\alpha }{t}_n^{2\alpha -1}\parallel g^{\left(3\right)}\left(0\right)\parallel .$$

(145)

For the convolution part $I_c$

$$\parallel I_c\parallel \le C{\tau }^{4-\alpha }{\int }_0^{t_n}{(t_n-s)}^{2\alpha -1}\parallel g^{\left(4\right)}\left(s\right)\parallel ds.$$

(146)

Combining the bounds for $I_0$ to $I_3$ and $I_c$ completes the proof. Similarly, by making relevant smoothness assumptions on the nonlinear term f.    □

6. Numerical Simulations

The numerical implementation of the high-order Zeta scheme decouples the precomputation of coefficients from the time evolution to ensure reproducibility and computational efficiency.

The Caputo integration kernel ${(t-s)}^{-\alpha }$ induces a weak singularity at $t=0$. For solutions lacking sufficient vanishing initial derivatives, this restricts temporal regularity and severely reduces the convergence order of high-order schemes, necessitating starting corrections.

This section empirically validates the temporal convergence. Spatial discretization employs linear finite elements with a consistent mass matrix ($N_x=64$). To isolate temporal errors, numerical solutions are compared against a fully corrected reference solution on an ultra-fine mesh (${\tau }_{\mathrm{ref}}=2^{-14}$). We report $L_2$-norm errors at $T=1.0$ and local convergence orders ${log}_2(e_{2\tau }/e_{\tau })$. Rate calculations strictly omit errors below $5\times {10}^{-14}$ (to suppress machine precision contamination) and the coarsest grid pair $2^{-3}\to 2^{-4}$ (to avoid pre-asymptotic order inflation). The average convergence order over valid pairs is reported for each $\alpha$.

The experiments are structured progressively: verifying baseline homogeneous cases, introducing source incompatibility, testing polynomial compatibility, examining source regularity limits, and establishing the empirical regularity threshold.

6.1. Numerical Implementation and Time Evolution

The numerical implementation of the high-order Zeta scheme decouples the precomputation of coefficients from the time evolution to ensure reproducibility and computational efficiency.

The spatial domain is discretized via linear finite elements on a uniform mesh of size h. To strictly preserve the high-order accuracy, we utilize the exact consistent mass matrix rather than mass lumping. This yields the standard tridiagonal mass matrix $\mathbf{M}$ and stiffness matrix $\mathbf{S}$:

$$\mathbf{M}={\displaystyle \frac{h}{6}}\left[\begin{array}{cccc}4& 1& & \\ 1& \ddots & \ddots & \\ & \ddots & \ddots & 1\\ & & 1& 4\end{array}\right],\mathbf{S}={\displaystyle \frac{1}{h}}\left[\begin{array}{cccc}2& -1& & \\ -1& \ddots & \ddots & \\ & \ddots & \ddots & -1\\ & & -1& 2\end{array}\right].$$

(147)

Both matrices are assembled in sparse formats, and the initial condition $u_0$ is projected onto this discrete space.

The convolution weights $b_j$ and the starting correction coefficients ${\overline{c}}_n$ are computed once prior to the time loop, following their exact analytical definitions established in Section 2 and Section 3. For our baseline 4th-order experiments ($K=4$), the starting corrections are applied strictly to the first 3 time steps and take the exact values ${\overline{c}}_1^{\left(4\right)}={\displaystyle \frac{31}{24}}$, ${\overline{c}}_2^{\left(4\right)}=-{\displaystyle \frac{7}{6}}$, and ${\overline{c}}_3^{\left(4\right)}={\displaystyle \frac{3}{8}}$.

For the inhomogeneous equation, the correction also involves the source term $g_h\left(t\right)=\mathbf{M}f\left(t\right)-\mathbf{S}{u}_0$ and its time derivatives at $t=0$. Matching the Taylor expansion of $g_h\left(t\right)$ with the homogeneous correction factors yields the auxiliary derivative-correction coefficients:

$$d_{\ell ,n}^{\left(4\right)}={\displaystyle \frac{n^{\ell }}{\ell !}}{\overline{c}}_n^{\left(4\right)},\ell =1,\dots ,K-1.$$

(148)

For high-order schemes ($K\ge 5$), calculating these coefficients using standard IEEE 754 double-precision arithmetic leads to a severe loss of significant digits due to the ill-conditioned Vandermonde matrices. To address this, the coefficient generation is executed using the Advanpix Multiprecision Computing Toolbox (MP) for MATLAB at 100 decimal digits. After accurately generating $b_j$ and ${\overline{c}}_n$, they are cast back to double-precision arrays. This strategy confines the rounding errors entirely to the initialization phase without inflating the computational cost of the time-stepping loop.

Since the implicit discrete operator is time-independent, the coefficient matrix $\mathbf{A}=\mathbf{S}+{\tau }^{-\alpha }{b}_0\mathbf{M}$ is assembled and LU-factorized only once. For each time step $n=1,\dots ,N_t$, the fully discrete scheme reads:

$$\mathbf{A}{V}^n=\mathbf{M}f\left(t_n\right)-\mathbf{S}{u}_0-{\tau }^{-\alpha }\mathbf{M}{\displaystyle \sum _{j=1}^{n-1}}{b}_{n-j}{V}^j+C^{\left(n\right)},$$

(149)

where $V^n$ represents the homogenized discrete solution, and the correction vector $C^{\left(n\right)}$ is active only for $n\le K-1$:

$$C^{\left(n\right)}={\overline{c}}_n^{\left(K\right)}{g}_h\left(0\right)+{\displaystyle \sum _{\ell =1}^{K-1}}{\tau }^{\ell }{d}_{\ell ,n}^{\left(K\right)}{g}_h^{(\ell )}\left(0\right),1\le n\le K-1,$$

(150)

with $g_h\left(0\right)=\mathbf{M}f\left(0\right)-\mathbf{S}{u}_0$. For $n\ge K$, the right-hand side reduces to $\mathbf{M}f\left(t_n\right)-\mathbf{S}{u}_0$. The solution $V^n$ is obtained via rapid forward and backward substitution using the precomputed LU factors.

Errors are evaluated against the reference solution $V_{\mathrm{ref}}$ measured in the mass-weighted $L^2$-norm: ${\parallel e\parallel }_{L^2}=\sqrt{e^T\mathbf{M}e}$.

6.2. Numerical Examples

Numerical experiments test the proposed scheme across six scenarios. The spatial domain is $\Omega =(0,1)$ with final time $T=1.0$. We apply piecewise linear finite elements on a uniform mesh with $h=1/64$. The nonsmooth initial data $u_0\left(x\right)={\chi }_{(0,0.5)}$ is discretized via a lumped-mass $L^2$-projection. Average rates are computed via log-linear regression, excluding the coarsest step $\tau =2^{-3}$ to avoid pre-asymptotic pollution. The reference solution is generated by the Zeta scheme with $K=4$ at $\tau =2^{-14}$. Furthermore, a systematic comparison against the corrected Gao and Cao schemes of [2] is included to evaluate pre-asymptotic stability.

Example 2. 

Consider the homogeneous equation:

$${}_0{}^{C}D_t^{\alpha }u(x,t)-\Delta u(x,t)=0,0<x<1,t>0,$$

$$u(0,t)=u(1,t)=0,$$

$$u(x,0)=u_0\left(x\right),$$

where Case (I) takes $u_0\left(x\right)=sin\left(\pi x\right)$, and Case (II) takes $u_0\left(x\right)={\chi }_{(0,0.5)}$.

Table 4 lists the errors and convergence orders. Uncorrected steps (Scheme A) limit the rate to $\mathcal{O}\left(\tau \right)$ due to the singularity at $t=0$. Scheme B recovers the theoretical $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ rate. The correction remains effective under nonsmooth initial data.

Example 3. 

Consider inhomogeneous problems with smooth sources:

$${}_0{}^{C}D_t^{\alpha }u(x,t)-\Delta u(x,t)=f(x,t),0<x<1,t>0,$$

$$u(0,t)=u(1,t)=0,$$

$$u(x,0)=sin\left(\pi x\right),$$

where Case (I) sets $f(x,t)=(1+{\pi }^2)sin\left(\pi x\right)$ (zero-order mismatch), and Case (II) sets $f(x,t)=e^{t}sin\left(\pi x\right)$.

Table 5 shows the results. For Case (I), $g\left(t\right)=f\left(t\right)-\Delta u_0=sin\left(\pi x\right)$ is constant. Only ${\overline{c}}_n$ activates, eliminating the constant offset. For Case (II), Scheme B corrects the origin singularity and restores the $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ rate.

Example 4. 

We test the effect of compatibility using $f(x,t)=t^{4}sin\left(\pi x\right)$:

$${}_0{}^{C}D_t^{\alpha }u(x,t)-\Delta u(x,t)=t^{4}sin\left(\pi x\right),0<x<1,t>0,$$

$$u(0,t)=u(1,t)=0,$$

$$u(x,0)=u_0\left(x\right).$$

Case (I) uses $u_0\left(x\right)=0$ (compatible), and Case (II) uses $u_0\left(x\right)=sin\left(\pi x\right)$ (incompatible).

As shown in Table 6, Case (I) yields $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ naturally because $f(x,t)$ and its derivatives vanish at $t=0$. In Case (II), $u_0\ne 0$ breaks compatibility, reducing Scheme A to $\mathcal{O}\left(\tau \right)$. Scheme B recovers full accuracy.

Example 5. 

We consider a low-regularity source to observe order limits:

$${}_0{}^{C}D_t^{\alpha }u(x,t)-\Delta u(x,t)=t^{1+\alpha }sin\left(\pi x\right),0<x<1,t>0,$$

$$u(0,t)=u(1,t)=0,$$

$$u(x,0)=0.$$

Table 7 shows that low source regularity limits convergence to $\mathcal{O}\left({\tau }^{2+\alpha }\right)$. Correction terms vanish automatically, confirming that starting steps cannot bypass the regularity barrier of the underlying PDE solution.

Example 6. 

To track how regularity influences uncorrected steps, we test manufactured solutions $u(x,t)=t^{\sigma }sin\left(\pi x\right)$ with varying $\sigma >0$. With $u(x,0)=0$, the source is:

$$f(x,t)=\left[{\displaystyle \frac{\Gamma (\sigma +1)}{\Gamma (\sigma +1-\alpha )}}{t}^{\sigma -\alpha }+{\pi }^2{t}^{\sigma }\right]sin\left(\pi x\right).$$

This test uses Scheme A exclusively.

Table 8 records the rates. Small values of $\sigma$ restrict convergence. Rates grow with $\sigma$ and hit the formal $4-\alpha$ limit once $\sigma \ge 4-\alpha$.

Remark 6. 

To isolate the temporal truncation error, an operator test evaluates the Caputo derivative of the supersmooth function $u\left(t\right)=t^6$ at $T=1.0$, with temporal grid size N up to 3000. Figure 2 displays the absolute errors across different α for three schemes: the $\mathcal{O}\left({\tau }^{2-\alpha }\right)$ Liao L1 scheme [20], the $\mathcal{O}\left({\tau }^{3-\alpha }\right)$ Alikhanov L2-1_σ scheme [21], and the proposed $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ Zeta scheme. Under ideal smoothness, the Zeta scheme achieves the designed $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ accuracy and shows notably smaller errors than the lower-order methods for all tested fractional orders.

Example 7. 

To evaluate absolute accuracy and numerical stability, we compare the proposed Zeta scheme against the Gao $\mathcal{O}\left({\tau }^{3-\alpha }\right)$ and Cao $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ schemes of [2] using the homogeneous fractional diffusion equation with nonsmooth initial data:

$${}_0{}^{C}D_t^{\alpha }u(x,t)-\Delta u(x,t)=0,0<x<1,t>0,$$

$$u(0,t)=u(1,t)=0,$$

$$u(x,0)={\chi }_{[0,0.5]}.$$

The spatial domain is discretized via a lumped-mass $L^2$-projection onto linear finite elements with a uniform mesh size $h=1/64$. The final time is $T=1$ and the test grid is $\tau \in \{2^{-3},2^{-4},\dots ,2^{-7}\}$. To avoid the systematic bias inherent in self-referential measurements, the reference solution is strictly computed using the Zeta scheme ($K=6$) at $\tau =2^{-14}$, establishing an independent baseline accurate to approximately $\mathcal{O}\left({10}^{-12}\right)$.

Table 9 reports the absolute $L^2$ errors and convergence rates. Evaluated against the independent reference, all schemes exhibit high-order convergence. The Gao and Cao schemes maintain their expected $\mathcal{O}\left({\tau }^{3-\alpha }\right)$ and $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ trajectories, respectively.

The Zeta scheme ($K=4$) achieves comparable accuracy to the Cao scheme, with absolute errors reaching $\mathcal{O}\left({10}^{-11}\right)$. The structural advantage of the Zeta scheme is that it naturally assigns exact spectral coefficients via the Euler–Maclaurin expansion, bypassing the need for customized starting corrections. This intrinsically preserves spectral consistency and provides a unified, stable framework for nonsmooth data.

6.3. Summary

In summary, the uncorrected scheme degrades to first-order accuracy whenever the exact solution exhibits a weak singularity at $t=0$, whether due to a non-zero initial condition or an incompatibility between the initial data and the source. The implementation of starting corrections successfully restores the $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ convergence rate. This restoration remains effective regardless of the spatial regularity of the initial data.

The two correction mechanisms, ${\overline{c}}_n^4$ and $d_{\ell ,n}^4$, function independently. As confirmed by the zero-order mismatch experiment, ${\overline{c}}_n$ alone is sufficient to recover the full homogeneous rate when only the constant offset is non-zero. However, these corrections cannot bypass the limits imposed by the source term’s own regularity. If the source lacks sufficient smoothness, the convergence rate becomes bounded by the continuous solution’s regularity; for instance, a source $f\sim t^{1+\alpha }$ restricts the rate to approximately $\mathcal{O}\left({\tau }^{2+\alpha }\right)$. The regularity threshold experiment precisely maps this boundary, showing that without correction, the scheme only attains its formal order if the solution satisfies $u\sim t^{\sigma }$ with $\sigma \ge 4-\alpha$ near $t=0$.

Finally, the numerical stability of the Zeta scheme stems directly from its generating function structure. Instead of relying on algebraic coefficient cancellation, the method utilizes the Euler–Maclaurin expansion to assign exact spectral coefficients. This specific design preserves spectral consistency and guarantees high-order convergence for nonsmooth data, entirely avoiding the need for ad-hoc initialization adjustments.

7. Conclusion

We constructed a high-order fully discrete Zeta scheme for time-fractional partial differential equations and established an optimal $\mathcal{O}\left({\tau }^{4-\alpha }\right)$ convergence rate for the fourth-order discretization. Numerical tests confirm this theoretical bound. The results demonstrate that the starting corrections successfully restore the optimal order by resolving the initial weak singularities. In nonhomogeneous scenarios, the Zeta scheme yields smaller errors than the reference method, an advantage that becomes particularly evident for small $\alpha$. The generating function framework provides a clear path for extensions to nonlinear problems and variable-order models, as the discrete weights depend smoothly on $\alpha$.