Exponentially Convergent Numerical Method for Abstract Cauchy Problem with Fractional Derivative of Caputo Type

TL;DR

This work addresses the numerical solution of the abstract Cauchy problem $\partial_t^\alpha u + A u = f$ with Caputo derivative for $\alpha\in(0,2)$ and sectorial $A$, by developing an exponentially convergent, contour-integral-based method for the mild solution. A corrected propagator representation together with a time-independent hyperbolic contour and sinc-quadrature enables uniform exponential convergence across all $\alpha\in(0,2)$ and finite times, while a multi-level parallel algorithm evaluates independent resolvent problems to obtain both homogeneous and inhomogeneous parts with rigorous error bounds. The approach handles unbounded operator coefficients and data known only on $[0,T]$, and remains robust to spectral placement of $A$; with detailed a priori error estimates, it demonstrates significant accuracy gains over time-stepping methods when parallel resources are available. Practical limitations include analyticity requirements on the right-hand side and moderate time horizons, but numerical experiments confirm the theoretical rates and illustrate the method’s effectiveness across a broad range of spectral parameters.

Abstract

We present an exponentially convergent numerical method to approximate the solution of the Cauchy problem for the inhomogeneous fractional differential equation with an unbounded operator coefficient and Caputo fractional derivative in time. The numerical method is based on the newly obtained solution formula that consolidates the mild solution representations of sub-parabolic, parabolic and sub-hyperbolic equations with sectorial operator coefficient $A$ and non-zero initial data. The involved integral operators are approximated using the sinc-quadrature formulas that are tailored to the spectral parameters of $A$, fractional order $α$ and the smoothness of the first initial condition, as well as to the properties of the equation's right-hand side $f(t)$. The resulting method possesses exponential convergence for positive sectorial $A$, any finite $t$, including $t = 0$ and the whole range $α\in (0,2)$. It is suitable for a practically important case, when no knowledge of $f(t)$ is available outside the considered interval $t \in [0, T]$. The algorithm of the method is capable of multi-level parallelism. We provide numerical examples that confirm the theoretical error estimates.

Figures (10)

• Figure 1: Schematic plot of the complex neighborhood $D \equiv D_d$ of $\mathbb R$ where the parametrized integrand $\mathcal{F}_\alpha(t,\xi)$ remains analytic and exponentially decaying for any $t \in [0, T]$ (a) along with the image of $D_d$ under the mapping $v \to z(v)$ defined by $\Gamma_I$ (b) and the region $z^\alpha(v)$, $v \in D_d$ (c). The "forbidden" regions of complex plane are indicated by "beige" color). ($\alpha =1.3$, $\rho_s = \pi$, $\varphi_s = {\pi}/{6}$).
• Figure 1: Exact solution $u(t,0.5)$ of problem \ref{['eq:FCP_DE']}, \ref{['eq:FCP_BC']} with $f(t) = 0$ and $A$, $u_0$, $u_1$ defined by \ref{['eq:FCP_ex1_hom_A']}, \ref{['eq:FCP_ex1_hom_IV']} ($L =1$, $k_0 = 1$, $k_1 = 4$, $a = 1$): (a) the case $\alpha = 0.1, 0.3, 0.5, 0.7, 1$; (b) the case $\alpha = 1, 1.2, 1.5, 1.7, 1.9$.
• Figure 2: Error $\mathcal{E}_{\mathrm{h}}(t,0.5)$ of the approximate solution $\widetilde{u}_{\mathrm{h}}^N$ to problem \ref{['eq:FCP_DE']}, \ref{['eq:FCP_BC']} with $f(t) = 0$, $A$, $u_0$, $u_1$ being defined by \ref{['eq:FCP_ex1_hom_A']}, \ref{['eq:FCP_ex1_hom_IV']} and $L =1$, $k_0 = 1$, $k_1 = 4$, $a = 1$. Graphs from the top row of subplots are for $\alpha = 0.1, 0.3, 0.5, 0.7, 1$ and (a)$N = 32$; (b)$N = 64$(c); $N = 128$. Graphs from the bottom row of plots correspond to $\alpha = 1, 1.2, 1.5, 1.7, 1.9$ and (d)$N = 32$; (e)$N = 64$; (f)$N = 128$.
• Figure 3: Sup-norm error $\mathrm{err}_{\mathrm{h}}$ of the approximate solution $\widetilde{u}_{\mathrm{h}}^N$ to problem \ref{['eq:FCP_DE']}, \ref{['eq:FCP_BC']} with $f(t) = 0$, $A$, $u_0$, $u_1$ being defined by \ref{['eq:FCP_ex1_hom_A']}, \ref{['eq:FCP_ex1_hom_IV']} and $L =1$, $k_0 = 1$, $k_1 = 4$. Graphs from sublots correspond to the different values of diffusivity constant: (a)$a = 1 \times 10^{-5}$; (b)$a = 0.1$; (c)$a = 1$; (d)$a = 10$;
• Figure 4: Exact solution $u(t,0.5)$ of problem \ref{['eq:FCP_DE']}, \ref{['eq:FCP_BC']} with $f(t) = \sin{\pi x} + t \sin{4\pi x}$, $u_0 = u_1 = 0$ and $A$, defined by \ref{['eq:FCP_ex1_hom_A']} with $a = 1$, $N_I = 256$: (a) the case $\alpha = 0.1, 0.3, 0.5, 0.7, 1$; (b) the case $\alpha = 1, 1.2, 1.5, 1.7, 1.9$.

Theorems & Definitions (23)

• Theorem 1.1: SytnykWohlmuth2023
• Proposition 3.1
• Proof 1
• Lemma 3.2
• Proof 2
• Remark 3.3
• Lemma 3.4
• Proof 3
• Lemma 3.5
• Theorem 3.6