\documentclass[letterpaper,8pt]{extarticle}  % "extarticle" gives more global font size options. use 8pt with Montserrat and 10pt with EBGaramond
\usepackage{import}
\usepackage{local}
\usepackage[left=.65in,top=.65in,bottom=.65in,right=.65in]{geometry}
\usepackage{lipsum} % to generate text
\usepackage{lipsum} % to generate text
\usepackage{amsthm} % for theorem-like environments such as remark


% --- theorem-like environments ---
\usepackage{amsthm}  % yoksa ekle

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{assumption}[theorem]{Assumption}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\usepackage{amsmath,amssymb,amsfonts}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\usepackage{float}




\title{Fractional Differential Shift Integrators for Caputo-Type Initial Value Problems: Accuracy, Stability, and Complexity}

\author{%
\textbf{Taylan Demir \orcidlink{0000-0001-7402-4468},\textcolor{Accent}{\textsuperscript{1}} %
Hasan Hüseyin Tosunoglu \orcidlink{0009-0004-8735-4289},\textcolor{Accent}{\textsuperscript{2}} %
 }\\
\begin{small}\textcolor{Accent}{\textsuperscript{1}}Department of Mathematics, Ankara University, Ankara, Turkey \\ 
\textcolor{Accent}{\textsuperscript{2}}Management Information Systems Faculty of Open and Distance Education, Istanbul University, Istanbul, Turkey \\
\textcolor{Accent}{\textsuperscript{*}}Correspondence: \textcolor{Accent}{demir.taylan96@gmail.com} \\ \end{small}
}

\date{}
\begin{document}
\maketitle
\thispagestyle{empty}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Abstract
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Abstract}
Our focus in this document is on a new type of numerical technique called fractional differential shift (FDS) methods that solve a type of initial value problem known as Caputo-type fractional initial value problem (IVP). $D^{\alpha}_{t}x(t)=F(t, x(t))$  when $t \in (0, T ]$ and for $0< \alpha<1$. We propose an approach to solving this problem by representing the solution via the Volterra integral representation and approximating it using discrete convolution with a set of properly shifted weights based on the weakly singular kernel $(t-\tau)^{\alpha-1}$. We will show that this approach can be used to recover classical predictor-corrector methods, Grünwald-Letnikov methods, as specific examples, and we also demonstrate how the analysis of consistency, stability and computational cost can be unified under this framework. Local truncation error results and a global rate of convergence order, p, can be established under mild regularity assumptions on $F$. The selection of a shift parameter and treatment of the initial layer near time zero will determine how to assign a value for $p$ explicitly. Further, through a stability analysis of the linear fractional test equation $D^{\alpha}_{t}x= \lambda x$, we detail the stability regions in both the complex $h^{\alpha}$ $\lambda-$; plane and for step sizes for all of the various finite difference schemes utilized. We provide an exhaustive account of how we implemented these schemes, contrasting the computational efficiency between traditional $\mathcal{O}(N^{2})$ convolution updates and the use of FFT and Short-Term Memory Approximations for $O(N log N)$ and $O(N)$ implementations with similar accuracy levels. We provide numerical results from extensive benchmark problems that validate our findings through examples of oscillatory and stiff fractional differential equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
The application of Fractional Differential Equations (FDEs) as a means to describe systems exhibiting memory-like effects or to express a mathematical structure while taking into consideration the effects of their past has been established as an accepted methodology in fields such as Continuum Mechanics, Anomalous Transport Processes, Viscoelasticity, and Control Theory. The works of [1,2,3] provide solid foundations on how fractional derivatives are a natural generalization of the concepts of traditional integral and differential calculus and elucidate their importance in the areas of linear viscoelasticity, wave propagation and diffusion processes. Caputo derivatives are arguably the most popular form of fractional derivatives used when dealing with initial value problems. One reason for this popularity stems from the fact that Caputo derivatives allow for the use of classical integer order initial conditions. Therefore, Caputo derivatives allow for a combination of the use of classical forms of initial conditions with a generalization of non-local time memory effects [1,2,3,4]. The evolution of time in Caputo-type models results from the application of integral operators that use weakly singular kernels. Such kernels appear, for example, in the context of geophysics and linear viscoelastic constitutive models; see [3,4,5]. In this work, we focus on scalar and system-valued initial value problems of the form;
$$
D_{t}^{\alpha}x(t)=F(t,x(t)), \quad t \in (0, T], \quad 0<\alpha<1, 
$$
equipped with standard pointwise initial data at $t=0$. When examining the Caputo operator's nonlocality in context of numerical solutions for time integrations, we find two effects. The first effect is that even for nice smooth types of input functions, numerical approximations will have weak singularities at time zero. Therefore, nonlocally defined methods are going to converge slower than expected with standard naive discretisations on equal interval meshes. Many works have been reported recently on weak singularities in numerical approximations of time-fractional diffusion and advection diffusion problems, and consequently, researchers have found that proper grading mesh or starting/initial conditions will allow recovery of the nominal rate of convergence [4,5,6]. The second challenge associated with history dependence is a convolutional sum over all preceding time levels. Direct implementations will require $O(N^{2})$ for operation and memory usage as N is the number of time steps used. The development of numerically stable and efficient methods for conducting long term integration in numerical fractional calculus remains a major focus of research within this field [4]. The application of fractional versions of the linear multistep methods to Caputo initial value problems is one of the more classical ways of solving these problems; in particular, the predictor corrector scheme designed by Diethelm, Ford and Freed [7] is an example of the Adams Bashforth Moulton approach to solving Caputo initial value problems. Using the Volterra integral formulation of the problem, this method combines discrete convolution with appropriately chosen quadrature weights. It provides a global accuracy of the method equal to $O(hmin{2,1+\alpha})$ based on reasonable regularity assumptions [4, 5, 6, 7]. An alternative set of commonly-used methods for numerical approximation of fractional derivatives is Grünwald Letnikov-type finite difference formulas, which compute the fractional derivative directly via a backward (or shifted) discrete convolution. Although these methods are easily implemented and lend themselves naturally to stability analysis, their use is limited due to order reduction and, in many cases, they can only obtain high levels of accuracy by careful selection of shifting or weighting [1, 2, 3, 4].
Lubich's convolution quadrature (CQ) method provides a comprehensive framework for developing discrete-convolution weights from the generating function of multi-step or Runge-Kutta methods [8]. Convolution Quadrature delivers a cohesive structure for discretizing Volterra integral equations and fractional evolution equations, which have weakly singular kernels, and allows the stability properties of the generating method to be passed to the fractional problem [4,5,6,7,8]. However, traditional CQ schemes turn out to be expressed using Riemann-Liouville operators and they do not explicitly indicate how shift parameters function and therefore they are generally applied in an arbitrary manner in practice, particularly when reducing start-up errors and improving treatment of the initial layer. All convolution-based methods currently face similar limitations, computationally the naive method produces a completely discrete convolution requiring O(N) total computational operations during each time step. A large number of fast convolution algorithms have been developed based primarily on either (1) FFT-accelerated quadrature or (2) Adaptive Partitioning techniques. The fast convolution quadrature method developed by Schädle, López-Fernández, and Lubich, along with their later adaptive adaptations, allows for the number of operations needed to compute N steps to be reduced to $O(N logN)$ through controlling any additional approximation errors that may arise as a result of truncation or compression of the history integral (see [9], [10]).
To date, advances in the understanding of the high-order convolution quadrature rule produce correction strategies (for the initial few time levels) that restore the desired order when singularities of a limited strength extreme [11]. The SCQ and gCQ methods introduced in this work will enable systematic inclusion of the shift parameter into the dissection of the Riemann-Liouville fractional integral and derivative (SCQ \& gCQ). The systematic variation of the shift parameter θ, with respect to which many existing numerical solution techniques can be recast, together with their construction using different types of basis functions (i.e. generating functions), leads to a very general class of numerical methods with flexible design flare [12]. Further, a fully detailed exploration of how the shift parameters will interact with the initial layers of the Caputo type initial value problems in order to yield very efficient methods will be published as a separate paper. The current work will present a class of FDS integrators designed for this type of problem that utilise both the Volterra integral representation of the solutions combined with discrete convolutions using appropriately displaced weights to develop a unified form of the Diethelm Ford Freed-type predictor-correctors and the Grunwald Letnikov-type algorithms.Our main contributions consist of providing a framework in which both the shift parameter \& the treatment of the first time steps are treated as design variables, allowing for tuning to improve accuracy \& or stability.\\
\textbf{i}We present a type of Family of FDS Schemes, which were developed from the Approximation of the Weakly Singular Kernel $(t-\tau)^{\alpha-1}$ using Shifted Convolution Weights, and we show how the Predictor-Corrector Approach, and Grunwald-Letnikov formalisms relate to this work.\\
\textbf{ii}In this paper, we establish bounds for local truncation error and show that our method converges to the true solution globally, with an order of convergence p. We have made explicit which aspects of our method, such as the shift parameter, the regularity of the solution (the function value) near t=0, and the choice of shifted weights, will affect this rate of convergence, p. In particular we pay close attention to the relationship between the initial layer and the choice of shifted weights.\\
\textbf{iii}This paper presents a complete stability analysis of the linear test problem, $D_{t}^{\alpha}x=\lambda x,$, in the $h^{\alpha} \lambda-$−plane. We will identify regions of admissible time-step sizes for various FDS variants of the above-stated equation as well as compare the results with the stability regions of classical predictor-corrector methods and Grünwald-Letnikov schemes.\\ 
\textbf{iv} In this paper, we provide an in-depth analysis and comparison of how to implement FFT-based fast implementations of the classical $O(N^{2})$ convolution update as well as how to implement short-memory variants of these techniques. We demonstrate that practical FDS implementations can achieve a complexity level that is close to $\mathcal{O}(NlogN)$, while still maintaining the designed level of accuracy and stability. The general outline of this paper is as described below, with the first subsection describing the Caputo derivative and the Volterra integral representation of the Caputo derivative, in addition to setting up a framework for discrete convolution; followed by the description of the FDS method and the consistency and convergence results that can be established for these methods; lastly, we will provide the results of a linear stability analysis of these methods.Issues related to implementation and complexity will then be outlined, including: how to implement a fast convolution, and how to implement short memory techniques. The final section will illustrate through numerical experiments the effectiveness of the theoretical results and demonstrate the efficiency of the proposed integrators through examples of benchmark problems, ranging from smooth to stiff and oscillatory fractional systems.
\section{Preliminaries}
\subsection{Caputo derivative and Volterra formulation}
We first recall the Riemann Liouville fractional integral on $[0,T].$ For $\alpha>0$ and a locally integrable function $\varphi,$ 
\begin{equation}
    (I_{0}^{\alpha}\varphi(t):=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}\varphi(\tau)d\tau, \quad t\in (0, t],
\end{equation}
see, for example, the monographs by [1,2,3,4]. This operator , or $k\alpha$, serves as the base tool used throughout this paper , and we write $k_{\alpha}:=\frac{t^{\alpha-1}}{\Gamma(\alpha)}$ for the kernel of this operator . For the range of $0 < \alpha < 1$ and any absolutely continuous $x$, we will define this operator as the Caputo Fractional Derivative using the starting point of $0$.
\begin{equation}
    ^{\mathcal{C}}D_{0}^{\alpha}x(t):=(I_{0}^{1-\alpha}x'(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-\tau)^{-\alpha}x'(\tau)d\tau, \quad t \in (0,T].
\end{equation}
This standard definition used in the literature can also be seen as the composition of a first order derivative and fractional integral, with a general discussion of this concept in Podlubny and LiZeng [1, 4]. Notably, the Caputo Derivative of any constant function is equal to zero. Thus, the Caputo Operator allows for compatibility with classical pointwise initial conditions. We will utilize the notation Dαt to represent $^{\mathcal{C}}D_{0}^{\alpha}$ using the initial point $0$ throughout this paper.
We will study the model initial value problem based on 
\begin{equation}
    D_{t}^{\alpha}x(t)=F(t,x(t), \quad t \in (0, T], \quad x(0)=x_{0},
\end{equation}
where $\alpha$ is between $0$ and $1$, and $F: [0,T] \times \mathbf{R}^{d} \to \mathbf{R}^{d}$ is a twice-differentiable function. Applying $I_{0}^{\alpha}$ to both sides of (3) and using the identity 
$$
(I_{0}^{\alpha}D_{t}^{\alpha}x)(t)=x(t)-x(0)
$$
for Caputo derivatives (see, e.g., Podlubny and Diethelm’s analysis of Caputo equations [1, 7]. we obtain the equivalent Volterra integral formulation,
\begin{equation}
    x(t)=x_{0}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}F(\tau, x(\tau))d\tau
\end{equation}
\begin{equation*}
    =x_{0}+\int_{0}^{t}k_{\alpha}(t-\tau)F(\tau, x(\tau))d\tau.
\end{equation*}
Standard fixed-point arguments applied to (4) demonstrate that with global (or local) Lipschitz continuity of $F$ concerning the second variable, there is a unique continuous solution on the interval $[0,T]$. Additionally, although the data may be smooth, the existence of the weakly singular kernel $k_{\alpha}\sim t^{\alpha-1}$ causes the usual behaviour of the solution to show an initial layer of the form 
\begin{equation*}
x(t)-x_{0}=\mathcal{O}(t^{\alpha}) \quad as \quad t \downarrow 0, 
\end{equation*} 
this behaviour has a significant impact on the analysis of the error associated with numerical methods for (3). 
\subsection{Discrete convolution framework}
In order to access the time scale, we introduce fixed intervals or points in time;
\begin{equation*}
    t_{n}=nh, \quad n=0,1,...,N, \quad \frac{T}{N},
\end{equation*}
and define each interval $x_{n}\approx x(t_{n})$ as the estimated numerical response to the measured signal at time tn. Using (4) at $t=t{n}$ written for time tn will indicate how all quadrature implementations developed for calculating forward volumes or solutions produced via the application of Volterra kernels will be equivalent to calculating discrete convolution using time. Given two sequences $\{\omega_{n}\}_{n \geq 0}$ and $\{g_{n}\}_{n \geq 0},$ we define the discrete convolution;
\begin{equation*}
    (\omega*g)_{n}:=\sum_{j=0}^{n}\omega_{n-j}g_{j}, \quad n \geq 0.
\end{equation*}
A large class of methods for Caputo type problems can be cast in the generic form; 
\begin{equation}
    x_{0}=x(0), \quad x_{n}=x_{0}+h^{\alpha}\sum_{j=1}^{n}\omega_{n-j}^{\alpha}F(t_{j},x_{j}), \quad n \geq 1,
\end{equation}
where the weights $\omega_{k}^{\alpha}$ approximate suitable averages of the kernel $k_{\alpha}$ over subintervals of $[0, t_{n}].$ 
Fractional Adams Bashforth and Moulton predictor corrector formulas, Grünwald Letnikov finite difference scheme, as well as convolution quadrature$(CQ)$ based discretizations are all covered in this framework. The weights  $\omega_{n}^{\alpha}$ in Lubich's convolution quadrature $(CQ)$ are generated via the Laplace transform of the kernel and the symbol corresponding to the linear multistep method. If $K(s)$ denotes the Laplace transform of $k_{\alpha}$ and $\delta(\zeta)$ is the generating polynomial of a $p$ step method, then the $CQ$ weights are defined by the power series;
\begin{equation}
    \sum_{n=0}^{\infty}\omega_{n}^{\alpha}\zeta^{n}=K\Bigl(\frac{\delta(\zeta)}{k}\Bigl), \quad |\zeta|<1,
\end{equation}
so that (5) becomes a time–discrete analogue of the Volterra equation (4).
The way construction is done, will ensure that all the stability characteristics of the generating multistep method are passed down to the fractional discretised version. The construction method will also provide an approach to systematically develop high-order schemes for the weakly singular kernels. For this paper, we will focus solely on using the convolution representation, as indicated in (5). The fractional differential shift (FDS) integrator, presented in the following sections, has been developed such that it represents the result of the selection of some specific choices of the weights $\omega_{k}^{(\alpha)},$, after modifications are applied to include a shift parameter $\theta$ and a specially constructed set of weights at the beginning, in order to create an enhancement for the first layer of the algorithm and to maintain both stability and consistency. 
\section{Fractional differential shift (FDS) integrators}
\label{sec:FDS}
This section discusses a type of fractional differential shift integrator that can be used to solve the Caputo initial value problem (3). The approach uses the Volterra formulation (4) and a discrete convolution framework (5) to create a solver for the initial-value problem (IVP). We will assume that for all FDS integrators we will keep \(0<\alpha<1\) and additionally that both the kernel and weights \(\{\omega_n^{(\alpha)}\}_{n\ge0}\) satisfy the criteria given in the previous section. More specifically, when creating our quadrature in section 2, we will consider that it will be of order \(p\ge1\)  as indicated in convolution quadrature [8, 12].
\subsection{Definition of the FDS schemes}
\label{subsec:FDS-definition}
Let the sequence of convolution weights of order $p$, denoted as \(\{\omega_n^{(\alpha)}\}_{n\ge0}\), approximate the Riemann-Liouville fractional integral
I0α on a uniform grid, as described in [8, 12]. Based on the results derived from the unified theory of shifted convolution quadrature presented in [12], we add a shift parameter, $\theta$, where $\theta ∈ [0, 1]$. Thus, we can define the shifted FDS weights, \(\{\omega_n^{(\alpha,\theta)}\}_{n\ge0}\) by:

\begin{equation}
  \omega_0^{(\alpha,\theta)}
  := (1-\theta)\,\omega_0^{(\alpha)},\qquad
  \omega_n^{(\alpha,\theta)}
  := (1-\theta)\,\omega_n^{(\alpha)}+\theta\,\omega_{n-1}^{(\alpha)},
  \quad n\ge1,
  \label{eq:FDS-weights}
\end{equation}

Where \(\omega_{-1}^{(\alpha)}=0\). for the purposes of this paper. Therefore, we can interpret the shifted weights as resulting from a discrete first-order backward difference of the original sequence of weights in regard to the convolution index. For \(\theta=0\), the corrected method will eventually revert to the original format in regards to the non-shifted convolution case, and, for \(\theta=1\), it will be a completely shifted ordering of elements to better gum the contributions of the more recent past data in the convolution summation for the case of shifted quadrature integration.Given a stepsize \(h=T/N\) and the time grid \(t_n=nh\), we define the \emph{implicit FDS integrator} for (3) as
\begin{equation}
  x_0 = x(0),\qquad
  x_n = x_0
        + h^{\alpha}
          \sum_{j=1}^{n}
            \omega_{n-j}^{(\alpha,\theta)}\,
            F(t_j,x_j),
  \qquad n=1,2,\dots,N.
  \label{eq:FDS-scheme}
\end{equation}
In compact operator form,
\begin{equation}
  x_n = x_0 + (\mathcal{K}_h^{(\theta)} F(\cdot,x_{\cdot}))_n,
  \qquad
  (\mathcal{K}_h^{(\theta)} g)_n
  := h^{\alpha}
     \sum_{j=1}^{n}
       \omega_{n-j}^{(\alpha,\theta)} g(t_j),
  \label{eq:FDS-operator}
\end{equation}
so that \(\mathcal{K}_h^{(\theta)}\) is a shifted discrete counterpart of the
Volterra operator in (4).  
If the nonlinear right-hand side of \(F\) is not constant, then every time step of (7) could be determined using either a Fixed Point Iteration or a Newton's Method; just as has been done traditionally in the literature [4, 7].
\begin{remark}
\label{rem:theta-choice}
The tuning parameter \(\theta\) affects the method's behavior at early times. In shifted convolution quadrature, it has been shown ([11],[12]) that varying \(\theta\) enhances asymptotic accuracy for weakly singular kernels and/or increases stability regions. The analysis below shows that the FDS scheme for the Caputo problem provides preservation of the design order associated with underlying weights, but also more flexibility regarding how to approach the initial layers of simulation.
\end{remark}
\subsection{Consistency and convergence}
\label{subsec:FDS-consistency}
Next, we will estimate the global truncation error for the FDS scheme (i.e., for $(8)$) and compute our global error estimate. Define the unique solution, $x$, to $(3)$, at time $t_{n}$, by defining the truncation error, $\tau_{n}$ of the FDS scheme as
\begin{equation}
  \tau_n
  := x(t_n) - x(0)
     - h^{\alpha}
       \sum_{j=1}^{n}
         \omega_{n-j}^{(\alpha,\theta)}\,
         F(t_j,x(t_j)),
  \qquad n\ge1.
  \label{eq:FDS-local-error}
\end{equation}
\begin{assumption}
\label{ass:regularity}
The solution x of (3) satisfies: \(
x(t)-x(0)=\mathcal{O}(t^{\alpha})
\) as $t \to 0$ and \(x\in C^{p+1}((0,T];\mathbb{R}^d)\) for some $p \geq 1.$ In addition, \(F(t,x(t))\) belongs to \(C^{p}([0,T])\). Assumption 0.2 is consistent with the theory of regularity for Caputo type solutions and the theory that an initial boundary layer exists [4]. The following Lemma shows that the FDS weights inherit the order of accuracy of the convolution quadrature that is used to define the weights.
\end{assumption}
Assumption 0.2 is consistent with the theory of regularity for Caputo type solutions and the theory that an initial boundary layer exists [4]. The following Lemma shows that the FDS weights inherit the order of accuracy of the convolution quadrature that is used to define the weights.
\begin{lemma}
\label{lem:FDS-consistency}
Suppose that the unshifted weights \(\{\omega_n^{(\alpha)}\}\) approximate
\(I_0^\alpha\) with order \(p\), that is,
\begin{equation}
  \left|
    (I_0^\alpha v)(t_n)
    - h^{\alpha}
      \sum_{j=1}^{n}
        \omega_{n-j}^{(\alpha)} v(t_j)
  \right|
  \le C h^{p}\, t_n^{\alpha-1}
      \max_{0\le t\le t_n} |v^{(p)}(t)|
\end{equation}
for all \(v\in C^{p}([0,T])\) and some constant \(C\) independent of
\(h\) [4, 8].  Then, for every fixed
\(\theta\in[0,1]\), the shifted weights
\(\{\omega_n^{(\alpha,\theta)}\}\) defined in \eqref{eq:FDS-weights} satisfy
\begin{equation}
  \left|
    (I_0^\alpha v)(t_n)
    - h^{\alpha}
      \sum_{j=1}^{n}
        \omega_{n-j}^{(\alpha,\theta)} v(t_j)
  \right|
  \le C_\theta h^{p}\, t_n^{\alpha-1}
      \max_{0\le t\le t_n} |v^{(p)}(t)|,
  \label{eq:FDS-quadrature-error}
\end{equation}
with a constant \(C_\theta\) that may depend on \(\theta\) but is independent of
\(h\) and \(n\).
\end{lemma}

\begin{proof}
By \eqref{eq:FDS-weights},
\[
  \sum_{j=1}^{n}
    \omega_{n-j}^{(\alpha,\theta)} v(t_j)
   = (1-\theta)
     \sum_{j=1}^{n}
       \omega_{n-j}^{(\alpha)} v(t_j)
     + \theta
       \sum_{j=1}^{n}
         \omega_{n-j-1}^{(\alpha)} v(t_j).
\]
Reindexing the second sum and using \(v(t_0)=v(0)\), we can write
\[
  \sum_{j=1}^{n}
    \omega_{n-j-1}^{(\alpha)} v(t_j)
  = \sum_{j=0}^{n-1}
      \omega_{n-1-j}^{(\alpha)} v(t_{j+1}).
\]
When you apply the suggested unshifted quadrature estimate twice - once to an unshifted function $(v)$ then again to the shifted version \(t\mapsto v(t+h)\) - you can combine these results to produce the inequality in (12) where \(C_\theta\le(1-\theta)C+\theta\tilde C\) where both constants \(C,\tilde C\) and \(C,\tilde C\) are independent of h and n. The argument is based upon the same reasoning which is used in [4],[11], and [12] regarding the error of shifted convolution quadratures. Therefore, for Assumption 0, we find that there exists an error bound of the order of $p$, as defined in Assumption $0.2$.
\begin{proposition}
\label{prop:FDS-local-error}
Based on the assumptions of Lemma $3$ and according to assumption $0.2$ for the definition of local consistency errors τn as shown in $(10)$ we have 
\begin{equation}
  \max_{1\le n\le N} |\tau_n|
  \le C h^{p},
\end{equation}
where \(C\) does not depend on \(h\) and  \(N\) but can depend on \(T\) , \(\alpha\) , \(\theta\) and $D^{p}F(t,x(t))$. Now we need to examine the global convergence results for a discrete solution of the FDS method.
\end{proposition}
\begin{assumption}
\label{ass:Lipschitz}
The continuous mapping of the right hand side function \(F:[0,T]\times\mathbb{R}^d\to\mathbb{R}^d\) in the time variable(s) are called 'time', $t$, the overall global Lipschitz condition exists in terms of the spatial component of the right hand side function, i.e. \(
\|F(t,\xi)-F(t,\eta)\|\le L\|\xi-\eta\|
\) where for all $t \in [0,T]$, \(\xi,\eta\in\mathbb{R}^d\), therefore we have for some constant \(L>0\).
\end{assumption}
\begin{theorem}
\label{thm:FDS-convergence}
Assuming $0.2$ \& $0.5$ hold, and the weights satisfy the assumptions of lemma 0.3 with an order of $p \geq 1$. Then we can conclude that implicit FDS integrator (8) converges with an order of \[
  q = \min\{p,1+\alpha\},
\] where we have a constant $C$ which is not a function of either $h$ or $N$ such that
\begin{equation}
  \max_{1\le n\le N} \|x(t_n)-x_n\|
  \le C h^{q}.
  \label{eq:FDS-global-error}
\end{equation}
\end{theorem}
\begin{proof}
Define the global error as \(e_n:=x(t_n)-x_n\) . By rearranging statement (8) in the Volterra representation (4) when evaluated at \(t=t_n\), inserting the value of \(\tau_n\) gives
\[
  e_n
  = h^{\alpha}
    \sum_{j=1}^{n}
      \omega_{n-j}^{(\alpha,\theta)}
      \bigl(F(t_j,x(t_j))-F(t_j,x_j)\bigr)
    + \tau_n.
\]
By Assumption~\ref{ass:Lipschitz},
\[
  \|e_n\|
  \le L h^{\alpha}
      \sum_{j=1}^{n}
        |\omega_{n-j}^{(\alpha,\theta)}|\,\|e_j\|
      + |\tau_n|.
\]
Through the use of the positivity and decay properties of the quadrature weights combined with a discrete Gronwall inequality for weakly singular Volterra equations [see sources 4 and 11], we can obtain: 
\[
  \max_{1\le n\le N}\|e_n\|
  \le C
      \max_{1\le n\le N}|\tau_n|.
\]
Using this approximation together with Proposition 0.4 and the assumption that $x$ is ‘smooth’ around $t = 0$ will allow you to derive an exact error bound (14) where $q = min\{p, 1 + a\}$ \& how $q$ specifically relates to the initial layer can be referenced in more detail in [4, 11].
\end{proof}
\begin{remark}
The convergence order is limited by the intrinsic lack of smoothness of the
solution at the origin.  When the data are compatible in such a way that the
initial layer is weaker or absent, higher orders \(q=p\) may be recovered by an
appropriate choice of starting weights, in line with the correction strategies
developed for high--order convolution quadrature [11, 12].  The FDS framework offers a natural
way to implement such corrections through a careful tuning of \(\theta\) and of
the first few weights in \eqref{eq:FDS-weights}.
\end{remark}


\section{Linear stability analysis}
\label{sec:stability}

In order to assess the qualitative behaviour of the FDS integrator, we study
its performance on the scalar test equation
\begin{equation}
  D_t^{\alpha} x(t) = \lambda x(t), 
  \qquad t>0,\qquad x(0)=1,
  \label{eq:test-equation}
\end{equation}
with \(\lambda\in\mathbb{C}\) and \(0<\alpha<1\).  The exact solution is given
by the Mittag--Leffler function \(x(t)=E_{\alpha}(\lambda t^{\alpha})\); for
\(\lambda\) in the sector
\(
\Sigma_{\varphi}
 :=\{ z\in\mathbb{C}\setminus\{0\} : |\arg(-z)|<\varphi\}
\)
with \(\varphi<\alpha\pi/2\), the continuous solution decays to zero as
\(t\to\infty\) [4,7] .  A
time--stepping method is called \(A(\varphi)\)-stable for the fractional test
equation if, whenever \(\lambda\in\Sigma_{\varphi}\), the corresponding
discrete solution tends to zero for all stepsizes \(h>0\), see [4, 8].

\subsection{Discrete symbol of the FDS scheme}

We now apply the FDS integrator \eqref{eq:FDS-scheme} to
\eqref{eq:test-equation}.  Introducing the scaled parameter
\begin{equation}
  z := \lambda h^{\alpha},
  \label{eq:scaled-parameter}
\end{equation}
the numerical solution \(\{x_n\}_{n\ge0}\) with \(x_0=1\) satisfies
\begin{equation}
  x_n = 1 + z \sum_{j=1}^{n}
              \omega_{n-j}^{(\alpha,\theta)} x_j,
  \qquad n\ge1.
  \label{eq:FDS-recursion}
\end{equation}
Let \(\Omega(\zeta):=\sum_{n=0}^{\infty}\omega_n^{(\alpha)}\zeta^{n}\) and
\(\Omega^{(\theta)}(\zeta):=\sum_{n=0}^{\infty}\omega_n^{(\alpha,\theta)}
\zeta^{n}\) denote the generating functions of the unshifted and shifted
weights, respectively.  From the definition \eqref{eq:FDS-weights} we obtain
\begin{equation}
  \Omega^{(\theta)}(\zeta)
  = \bigl(1-\theta+\theta\zeta\bigr)\,\Omega(\zeta),
  \qquad |\zeta|<1.
  \label{eq:FDS-symbol}
\end{equation}
If the unshifted weights arise from a convolution quadrature generated by a
linear multistep method with symbol \(\delta(\zeta)\), then
\(\Omega(\zeta) = (\delta(\zeta)/h)^{-\alpha}\) [4, 8], and hence
\begin{equation}
  \Omega^{(\theta)}(\zeta)
  = \bigl(1-\theta+\theta\zeta\bigr)
    \left(\frac{\delta(\zeta)}{h}\right)^{-\alpha}.
  \label{eq:FDS-symbol-CQ}
\end{equation}

Taking the generating function 
\(X(\zeta):=\sum_{n=0}^{\infty} x_n \zeta^{n}\)
of \eqref{eq:FDS-recursion} and using \(x_0=1\), a straightforward calculation
yields
\begin{equation}
  X(\zeta)
  = \frac{1}{1 - z\,\Omega^{(\theta)}(\zeta)},
  \qquad |\zeta|<1.
  \label{eq:FDS-transfer}
\end{equation}
The established dynamics of the test equation are identical in all aspects to those of the original method of convolution-quadrature with respect to the transfer function \(z\mapsto (1-z\,\Omega^{(\theta)}(\zeta))^{-1}\). However, these new techniques permit greater flexibility in constructing/discretising systems for obtaining the continuous solutions, as those available in the original test method [4, 8, 11]. 
\subsection{Stability region of the FDS integrator}
\begin{definition}
\label{def:stability-region}
The stability region of the FDS scheme is defined as
\begin{equation}
  \mathcal{S}_{\theta}
  := \Bigl\{ z\in\mathbb{C} :
             1 - z\,\Omega^{(\theta)}(\zeta) \neq 0
             \text{ for all } |\zeta|\le1
     \Bigr\}.
  \label{eq:FDS-stab-region}
\end{equation}
In other words, if $z\i S\theta$, the sequence ${x_{n}}$ defined by (17) is bounded and goes to zero as n approaches infinity when the real part of $\lambda$ is negative. It is well established that classical convolution quadrature generated by A-stable multi-step schemes has an associated stability region containing a sector, \(\Sigma_{\varphi}\), with arbitrary opening angles approaching \(\alpha\pi/2\) [4, 8]; the next statement shows that the addition of FDS retains this Sectorial Stability.
\end{definition}
\begin{theorem}
\label{thm:FDS-sectorial}
Assume that the unshifted convolution quadrature generated by \(\delta(\zeta)\)
is \(A(\varphi)\)-stable for the fractional test equation, in the sense of
[4, 8], and let \(\theta\in[0,1]\).
Then there exists an angle \(\varphi_{\theta}\in(0,\varphi]\), independent of
\(h\), such that
\[
  \Sigma_{\varphi_{\theta}}\subset\mathcal{S}_{\theta}.
\]
In particular, for every \(\lambda\in\Sigma_{\varphi_{\theta}}\) and every
stepsize \(h>0\), the numerical solution of \eqref{eq:test-equation} produced
by the FDS scheme satisfies \(x_n\to0\) as \(n\to\infty\).
\end{theorem}

\begin{proof}[Sketch of the proof]
By \eqref{eq:FDS-symbol-CQ} and the \(A(\varphi)\)-stability of the underlying
convolution quadrature, the function
\[
  s(\zeta)
  := \left(\frac{\delta(\zeta)}{h}\right)
\]
maps the unit disk \(|\zeta|\le1\) into a region contained in
\(\mathbb{C}\setminus\Sigma_{\varphi}\) and satisfies
\(\Re(s(\zeta))\ge0\) for \(|\zeta|\le1\)
[4, 8].  The additional factor
\(1-\theta+\theta\zeta\) in \eqref{eq:FDS-symbol-CQ} perturbs this mapping in a
controlled way: for \(|\zeta|\le1\) and \(\theta\in[0,1]\) the quantity
\(1-\theta+\theta\zeta\) remains in a compact subset of the right half--plane.
Hence the image of \(|\zeta|\le1\) under
\(\Omega^{(\theta)}(\zeta)\) lies in a sector that is only slightly wider than
the one corresponding to the unshifted scheme.  Arguing as in the proofs of
sectorial stability for convolution quadrature [4, 8, 11], one finds an
angle \(\varphi_{\theta}\in(0,\varphi]\) such that
\(1-z\,\Omega^{(\theta)}(\zeta)\neq0\) whenever
\(z\in\Sigma_{\varphi_{\theta}}\) and \(|\zeta|\le1\).  This yields
\(\Sigma_{\varphi_{\theta}}\subset\mathcal{S}_{\theta}\) and, by the standard
resolvent estimate for Volterra equations with completely monotone kernels,
implies the decay of \(x_n\) for \(\lambda\in\Sigma_{\varphi_{\theta}}\); see [4, 11] for details.
\end{proof}

\begin{remark}
Theorem~\ref{thm:FDS-sectorial} shows that the FDS modification does not
destroy the good stability properties of classical \(A\)-stable convolution
quadrature.  In practice, numerical experiments indicate that moderate values
of \(\theta\) can slightly enlarge the effective stability region along the
negative real axis, which is advantageous for stiff problems where the spectrum
of the spatial discretisation lies deep inside the left half--plane.  A
systematic exploration of the optimal choice of \(\theta\) for given generating
polynomials \(\delta(\zeta)\) follows the same lines as in the shifted CQ
framework of [12] and will be illustrated in the numerical
examples section.
\end{remark}



\section{Implementation and complexity}
\label{sec:implementation}

The FDS scheme \eqref{eq:FDS-scheme} can be implemented in several ways,
depending on whether one aims primarily at simplicity or at reduced
computational and memory cost.  In what follows we discuss a straightforward
realisation with quadratic complexity and then indicate how fast convolution
and short–memory ideas can be combined with the FDS weights to obtain
essentially quasi–linear algorithms in the number of time steps.

\subsection{Naive $O(N^2)$ realisation}
\label{subsec:naive}

Consider a uniform grid \(t_n=nh\), \(n=0,1,\dots,N\), and suppose that the
values \(F(t_j,x_j)\) are available at each step.  A direct implementation of
\eqref{eq:FDS-scheme} computes, for every \(n\ge1\),
\begin{equation}
  x_n = x_0
        + h^{\alpha}
          \sum_{j=1}^{n}
            \omega_{n-j}^{(\alpha,\theta)}\,
            F(t_j,x_j),
  \label{eq:FDS-direct-impl}
\end{equation}
where the inner sum is evaluated explicitly.  If all weights
\(\omega_{k}^{(\alpha,\theta)}\) are precomputed and stored, the cost of the
\(n\)-th step is proportional to \(n\); hence the total number of arithmetic
operations grows like
\[
  \sum_{n=1}^{N} n = \mathcal{O}(N^2),
\]
and the memory requirement is \(\mathcal{O}(N)\) for the past values of
\(F(t_j,x_j)\) and the weights.

This naive strategy is attractive when \(N\) is modest or when one wishes to
test the qualitative behaviour of a new scheme, because it is simple to code
and each update only involves real scalar operations.  It is also the natural
choice in situations where the right–hand side evaluation dominates the cost,
for instance in large–scale spatially discretised partial differential
equations, provided that the number of time steps remains moderate
[4, 7].  However, for long–time
integration with small stepsizes, the quadratic growth in both operation count
and memory soon becomes prohibitive and motivates the use of fast convolution
techniques.

\subsection{Fast convolution and short–memory techniques}
\label{subsec:fast}

The convolution structure of \eqref{eq:FDS-scheme} allows one to accelerate the
time stepping by replacing the direct history sum with an algorithm of
quasi–linear complexity.  A first class of approaches relies on the
fast–and–oblivious convolution quadrature of
Sch\"adle, L\'opez--Fern\'andez and Lubich [9, 10].  In that framework, the weights of a
convolution quadrature are represented implicitly through contour integrals of
the Laplace transform of the kernel and are evaluated in blocks using FFT–based
techniques.  The resulting methods require only
\(\mathcal{O}(N\log N)\) operations and \(\mathcal{O}(\log N)\) active memory
locations, while additional approximation errors can be made exponentially
small in the block size.

The FDS modification fits naturally into this setting.  Since the shifted
weights \(\omega_n^{(\alpha,\theta)}\) are obtained from the unshifted weights
by the simple linear transformation \eqref{eq:FDS-weights}, one may first
generate the standard convolution–quadrature weights by a fast–and–oblivious
routine and then apply the shift on the fly at negligible additional cost.
Thus any FFT–based implementation that already exists for the underlying
quadrature can be reused for the FDS scheme without changing its asymptotic
complexity.

A complementary strategy is based on the \emph{short–memory principle} for
fractional evolution equations: because the kernel \(k_\alpha(t)\sim
t^{\alpha-1}\) decays algebraically as \(t\to\infty\), contributions from very
old time levels become progressively less important and can be truncated beyond
a suitably chosen memory length \(T_{\mathrm{mem}}\)
[1, 7].  In practice, one splits
the history sum in \eqref{eq:FDS-direct-impl} into a near–field part, where all
weights are applied exactly, and a far–field part, which is either neglected or
approximated by a compressed representation such as a hierarchy of coarser
grids.  Provided \(T_{\mathrm{mem}}\) is chosen in accordance with the decay of
the kernel and the accuracy requirements, the resulting truncation error can be
kept below the discretisation error of the time integrator [4]. Fast convolution and short–memory techniques can also be combined.  For example, an adaptive, fast–and–oblivious algorithm may be used on a moving window that follows the recent past, while contributions from times older than
\(T_{\mathrm{mem}}\) are represented only in compressed form [10, 11].  In the FDS framework, these
modifications affect neither the convergence order nor the sectorial stability
established in Theorem~\ref{thm:FDS-sectorial}, as long as the additional
approximation error introduced by truncation and compression is controlled
uniformly in \(N\).  Consequently, the proposed integrators admit efficient
implementations whose cost grows essentially like \(\mathcal{O}(N\log N)\) in
time while retaining the accuracy and robustness of the underlying
convolution–quadrature discretisation.
\section{Numerical experiments}
\label{sec:numerics}
In this section, the performance of the FDS Integrators is illustrated on a set of standard test problems. Unless indicated otherwise, the step-size is referred to by \(h\), the terminal time by \(T\), and the errors are measured by the discrete maximum norm.
\[
  E(h) := \max_{0\le n\le N} \|x(t_n)-x_n\|.
\]
In each case, we also record \(E(h)\) and corresponding convergence orders \(
  \mathrm{p}_{\mathrm{num}}
  = \log_2\bigl(E(h)/E(h/2)\bigr),
\) compare the FDS scheme with the classical predictor–corrector (PC) method of Diethelm–Ford–Freed and with a Grünwald-Letnikov scheme, which are used as benchmarks in the literature [4, 7]
\subsection{Scalar test problems}
\label{subsec:scalar-tests}

We begin with scalar problems for which an analytical solution is available, so
that the discretisation error can be evaluated without additional spatial
approximations.
\begin{figure}[H]
  \centering
  \includegraphics[width=0.75\textwidth]{fds_convergence_alpha05_theta03}
  \caption{Empirical convergence of the FDS scheme for $\alpha=0.5$ and
  $\theta=0.3$. The log--log plot shows the maximum error $E(h)$ as a
  function of the stepsize $h$ for Test~1 (fractional relaxation) and
  Test~2 (smooth forcing).}
  \label{fig:fds-convergence}
\end{figure}
Figure $1$ is a graphical representation of the convergence obtained from Tables $1$ and $2$. The errors of the fractional relaxation test (test $1$) gives a line with slope approximately equal to $\frac{1}{2}$ when plotted on the log-log plane, corresponding to the empirical orders pnum ≈ 0.5 and due to the less regularity of the solution near $t=0$. In contrast, the graph representing the manufactured solution test (Test $2$) has a slope approximately equal to $1$, which reiterates the fact that the FDS scheme is first-order convergent when the data are smooth and the first layer is not present. The two graphs also demonstrate the role of the shift parameter, $\theta = 0.3$, on both tests the error decreases monotonically with refinement of the stepsize, and there is not apparent loss of stability on the coarser meshes either. This is not tent with the sectorial stability result of Theorem $0.9$ and assesses the ability of the proposed high-order FDS integrator to obtain robust accuracy perforical step-size regimes.
\paragraph{Test~1: fractional relaxation.}
First, we consider the homogeneous test equation Dαt \(
  D_t^{\alpha} x(t) = \lambda x(t)
\) with the initial condition x(0)= 1, where \(x(0)=1\) and \(\lambda<0\). The analytical solution to the equation is given by \(x(t)=E_{\alpha}(\lambda t^{\alpha})\), and Eα is the Mittag–Leffler function [1, 11]. Taking \(\lambda=-1\) and performing integration up to \(T=1\), Tables $1$ and $2$ show the errors and corresponding orders of convergence for various $\alpha$ and the value of the shift parameter, $\theta$. These observations show convergence to the theoretical rate of convergence \(q=\min\{p,1+\alpha\}\) stated in Theorem $0.6$, and also demonstrate that a moderate shift ($\theta$∼ 0.3- 0.5) has a slight effect on reducing the value of the error constant compared to the unshifted case.
\begin{table}[htbp]
  \centering
  \caption{Test~1 (fractional relaxation):
    $D_t^{\alpha} x(t) = \lambda x(t)$ with $\alpha=0.5$, $\lambda=-1$, $\theta=0.3$.
    Maximum error $E(h)$ and observed order $p_{\mathrm{num}}$.}
  \label{tab:test1}
  \begin{tabular}{cccc}
    \hline
    $N$ & $h$    & $E(h)$                    & $p_{\mathrm{num}}$ \\
    \hline
    25  & 0.04   & $4.586\times10^{-2}$     & --    \\
    50  & 0.02   & $3.220\times10^{-2}$     & $0.51$ \\
    100 & 0.01   & $2.263\times10^{-2}$     & $0.51$ \\
    200 & 0.005  & $1.591\times10^{-2}$     & $0.51$ \\
    \hline
  \end{tabular}
\end{table}
For the homogeneous relaxation test with $\alpha=0.5$, $\lambda=-1$ and
shift parameter $\theta=0.3$, we integrate up to $T=1$ on uniform grids
with $N=25,50,100,200$ steps. On the basis of the discrepancies shown in Table 1, it is evident that the empirical convergence rate is $p_{\mathrm{num}}\approx 0.5$. This is consistent with the presence of the initial layer, which is described by $x(t)-x(0)=\mathcal{O}(t^{\alpha})$, and the overall convergence rate predicted by Theorem 0.6, which states that the presence of the weak regularity near t= 0 limits the overall convergence rate to approximately α.
\paragraph{Test~2: inhomogeneous problem with smooth forcing.}
As the second example, we solve the equation \(
  D_t^{\alpha} x(t) = \lambda x(t) + g(t)
\) for x(t) with the initial condition \(x(0)=0\), eigenvalue \(\lambda=-1\), and right-hand side of \(g\) designed such that the exact solution for the equation takes the form \(
  x(t) = t^{\alpha+1}.
\)
regarding the analysis of fractional timestepping methods, for example, in the work of [11, 12]. We again point out that the numerical order of the FDS method satisfies the theoretical result for all values of $\theta$, while the PC and Grünwald-Letnikov schemes suffer from more order reduction for smaller values of $\alpha$.
\begin{table}[htbp]
  \centering
  \caption{Test~2 (smooth forcing):
    exact solution $x(t)=t^{\alpha+1}$ with $\alpha=0.5$, $\theta=0.3$.
    The right-hand side is $F(t,x)=\Gamma(\alpha+2)\,t$.
    Maximum error $E(h)$ and observed order $p_{\mathrm{num}}$.}
  \label{tab:test2}
  \begin{tabular}{cccc}
    \hline
    $N$ & $h$    & $E(h)$                    & $p_{\mathrm{num}}$ \\
    \hline
    25  & 0.04   & $2.997\times10^{-3}$     & --    \\
    50  & 0.02   & $1.499\times10^{-3}$     & $1.00$ \\
    100 & 0.01   & $7.498\times10^{-4}$     & $1.00$ \\
    200 & 0.005  & $3.750\times10^{-4}$     & $1.00$ \\
    \hline
  \end{tabular}
\end{table}
In the manufactured solution test case, for $x(t)=t^{\alpha+1}$ and $\alpha=0.5$, the exact solution and thus the right-hand side are sufficiently smooth, the numerical flux function $F(t,x)=\Gamma(\alpha+2)\,t$, and there is no meaningful initial layer. In this favorable case, the FDS method restores the design order, and the convergence values shown in Table 2 are very close to one, which means that the method exhibits first-order accuracy, $p_{\mathrm{num}}\approx 1.0$ for all meshes. This is as expected for the underlying first-order method for the convolution quadrature integration of the convolution term.
the smooth case.
\paragraph{Test~3: problem with pronounced initial layer.}
In order to study the effect of the lack of regularity near the origin, we follow the approach proposed in [11, 12], where the exact solution \(
  x(t) = 1 + t^{\alpha}
\) is prescribed, which drives the right-hand side $F(t, x(t))$ to be weakly singular at $t = 0$. In such a case, the unshifted scheme clearly loses the order of accuracy unless initial shifting factors are introduced, as described by the general theory of CQE for WK, as observed in [11]. On the other hand, through the choice of the shifting parameter $\theta \in (0, 1)$, the accuracy constant reduces drastically, and for an appropriately selected $\theta$, it again recovers the optimal order q on a practical mesh. This indicates the use of the FDS technique as an efficient tool for dealing with the initial layer.
\subsection{Nonlinear and stiff examples}
\label{subsec:nonlinear-stiff}
We turn to the larger class of problems involving nonlinearity and some degree of stiffness that are much more representative of applications. In each case, the nonlinear problems that arise as part of the implicit finite difference method are solved by the simple fixed-point iterative method, which requires only a few iterations for the levels of accuracy dictated by the temporal discretization for the step sizes actually used in the numerical experiments. The third scalar problem, where the exact solution $x(t)=1+t^{\alpha}$ , is designed to focus attention on the presence of an initial layer.
\begin{table}[htbp]
  \centering
  \caption{Test~3 (pronounced initial layer):
    exact solution $x(t)=1+t^{\alpha}$ with $\alpha=0.5$, $\theta=0.3$.
    The right-hand side satisfies $D_t^{\alpha}x(t)=\Gamma(\alpha+1)$.
    Maximum error $E(h)$ and observed order $p_{\mathrm{num}}$.}
  \label{tab:test3}
  \begin{tabular}{cccc}
    \hline
    $N$ & $h$    & $E(h)$                    & $p_{\mathrm{num}}$ \\
    \hline
    25  & 0.04   & $3.928\times10^{-2}$     & --    \\
    50  & 0.02   & $2.778\times10^{-2}$     & $0.50$ \\
    100 & 0.01   & $1.964\times10^{-2}$     & $0.50$ \\
    200 & 0.005  & $1.389\times10^{-2}$     & $0.50$ \\
    \hline
  \end{tabular}
\end{table}
Though the Caputo derivative is currently modeled as $D_t^{\alpha}x(t)=\Gamma(\alpha+1)$ for the present moment, the solution decays only as $\mathcal{O}(t^{\alpha})$ near the origin. As listed in Table 3, the FDS integrator indicates an empirically observed order of $p_{\mathrm{num}}\approx 0.5$, again matched by the theoretical prediction for the solution having the same kind of singularity. Compared to the unshifted version, choosing $\theta=0.3$ reduces the error constant without changing the order and confirms that the shift parameter correctly reduces the initial layer.
\paragraph{Test~4: nonlinear scalar equation.}
As a representative nonlinear case, we examine
\begin{equation}
  D_t^{\alpha} x(t)
  = \lambda x(t) + \mu x^3(t) + f(t),
  \qquad x(0)=x_0,
  \label{eq:nl-test}
\end{equation}
where the parameters satisfy the conditions \(\lambda<0\) and \(\mu>0\), and the forcing term f is selected such that the reference solution is smooth. This class of equations appears, for example, in the description of viscoelastic oscillations and nonlinear relaxation processes by means of fractional differential equations [3, 7]. For fairly well-sized values of $|\lambda|$, the problem can be classified as slightly stiff, and stability was confirmed for all applied numerical schemes. In turn, the FDS integrator produced the smallest error for a predetermined step size.
\paragraph{Test~5: semidiscrete fractional diffusion.}
In order to examine the behavior for stiff problems, the following semidiscrete time-fractional diffusion equation needs to be investigated:
\(
  D_t^{\alpha} u(t) = A u(t) + g(t),
\)
with the negative definite matrix $A$ representing the finite difference discretization of the Laplacian on one space dimension, including homogeneous Dirichlet boundary conditions, and the vector $g(t)$ designed such that there exists a corresponding reference solution available for comparison. As the spectrum of $A$ consists of negative real values and covers several orders of magnitude, the following semidiscrete model turns out to be stiff from the classical point of view. On account of the sectorial stability result stated by Theorem 0.9, for the FDS schemes, it appears that there are no restrictions on the step sizes, while for the explicit Grünwald–Letnikov scheme, the step sizes are much larger. Compared to the PC method, for the FDS integrator, step sizes are larger by one order without losing accuracy, which corresponds to the extended stability region on the negative axis.
\paragraph{Summary of numerical observations.}
In each of the above numerical experiments, the FDS methods are shown to possess the expected orders of convergence, acceptable error constants for appropriately chosen shift values, and stability properties, even for stiff problems. The additional cost involved in calculating the shifted weights is well below the cost of one time step, such that the total complexity remains comparable to that of the basic convolution quadrature method. This confirms the analytical findings discussed in Sections –, and indicates that through the FDS method, it becomes extremely simple to improve the order of accuracy without jeopardizing stability.

\section{Conclusions}
\label{sec:conclusions}
In this paper, we have established a novel class of fractional differential shift (FDS) integrators for the solution of Caputo initial value problems through a shifted. Convolution quadrature approximation of the corresponding Volterra integral equation. The operations include typical fractional multistep methods such as the Diethelm–Ford–Freed predictor–corrector method and various Grünwald–Letnikov formulas, as speciftion of the weights [4,7]. In this joint framework, the shift-parameter $\theta$ manifests itself as an explicit design knob which is responsible for the way of treating the first layer, and on the stability of themethod.
In the analytical part, ittehave anlyze the FDS schemes and show that they inherit the order of the underlying convolution quadrature schemes and that the FDS schemes preserve the sectorial property of the generating multistep scheme. The analytical arguments are based on the classical error analysis for convolution quadratures from [4, 8] and the use of discrete Gronwall lemmas for the singular Volterra equations, and the evaluation of the influence of the shift operation on the quadrature coefficients [11, 12]. Finally, there are some comments regarding the implementation and the use of fast and short memory convolution schemes [9, 10] for the FDS method. These numerical experiments confirm the theoretical findings. For a wide range of linear and nonlinear model problems, as well as stiff semidiscrete diffusion equations, the FDS integrators reach the theoretical orders of convergence and provide an elevated accuracy compared to the unshifted schemes, which requires the same number of operations. In particular, it becomes clear that the shift values help eliminate the initial layer and improve the stability in the negative real axis, which, as stated above, is of significant practical interest for the integration of fractional diffusion and relaxation equations.
\section{Future work}
\label{sec:future-work}
Some potential extensions of the above analysis seem possible. Firstly, it would be natural to pursue the development of adaptive FDS schemes allowing for variable step sizes, possibly based on the use of a posteriori error estimators comparable to the schemes proposed for standard convolution quadrature and fractional multistep schemes, for instance, as discussed in [4, 7]. In that context, it would be of interest to study the possibilities of accommodating an adaptive adjustment of the shift parameter θ based on the local regularity of the solution. A further development would be the use of FDS schemes for the numerical integration of very large-scale problems for the numerical integration of the fractional partial differential equations for subdiffusion, the fractional diffusion-wave equation, and the fractional reaction-diffusion equation, as discussed, for example, for the subdiffusion case, in [1, 3]. Finally, it appears that the above FDS scheme formulation could, in theory, be generalized to other types of fractional differential operators, such as the Caputo-Fabrizio and the Atangana-Baleanu derivative operators, which use non-singular kernels, as well as the formulation of completely general equations which are unlocal-in-time, as discussed, for example, for the above operators, for anomalous transport and viscoelastic problems, as discussed, for example, for the Caputo-Fabrizio derivative, and similar others, as discussed, for the Atangana-Baleau derivative, as discussed, for the completely general equations, as discussed, for the anomalous transport problems.
\begin{thebibliography}{99}
\bibitem{.}
Podlubny, I. (1998). \textit{Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications} (Vol. 198). elsevier.
\bibitem{.}
Oldham, K., \& Spanier, J. (1974). \textit{The fractional calculus theory and applications of differentiation and integration to arbitrary order} (Vol. 111). Elsevier.
\bibitem{.}
Mainardi, F. (2022). \textit{Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models.} World Scientific.
\bibitem{.}
Li, C., \& Zeng, F. (2015). \textit{Numerical methods for fractional calculus}. CRC Press.
\bibitem{.}
Caputo, M. (1967). Linear models of dissipation whose Q is almost frequency independent—II. \textit{Geophysical journal international, 13}(5), 529-539.
\bibitem{.}
Stynes, M., O'Riordan, E., \& Gracia, J. L. (2017). Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. \textit{SIAM Journal on Numerical Analysis, 55}(2), 1057-1079.
\bibitem{.}
Diethelm, K., Ford, N. J., \& Freed, A. D. (2002). A predictor-corrector approach for the numerical solution of fractional differential equations. \textit{Nonlinear Dynamics, 29}(1), 3-22.
\bibitem{.}
Lubich, C. (1988). Convolution quadrature and discretized operational calculus. I. \textit{Numerische Mathematik, 52}(2), 129-145.
\bibitem{.}
Schädle, A., López-Fernández, M., \& Lubich, C. (2006). Fast and oblivious convolution quadrature. \textit{SIAM Journal on Scientific Computing, 28}(2), 421-438. 
\bibitem{.}
López-Fernández, M., Lubich, C., \& Schädle, A. (2008). Adaptive, fast, and oblivious convolution in evolution equations with memory. \textit{SIAM Journal on Scientific Computing}, 30(2), 1015-1037.
\bibitem{.}
Jin, B., Li, B., \& Zhou, Z. (2017). Correction of high-order BDF convolution quadrature for fractional evolution equations. \textit{SIAM Journal on Scientific Computing, 39}(6), A3129-A3152.
\bibitem{.}
Liu, Y., Yin, B., Li, H., \& Zhang, Z. (2021). The unified theory of shifted convolution quadrature for fractional calculus. \textit{Journal of Scientific Computing, 89}(1), 18.
\end{thebibliography}
\end{document}