Error analysis of arbitrarily high-order stepping schemes for fractional integro-differential equations with weakly singular kernels

1 Introduction

Fractional integro-differential equations (FIDEs) play a pivotal role in modeling numerous complex phenomena across various scientific disciplines, showcasing their significance and broad applicability. These equations offer a powerful framework for capturing memory effects and long-range interactions, which are prevalent in diverse fields such as physics, biology, finance, and engineering. In physics, FIDEs are utilized to describe anomalous diffusion processes, viscoelastic materials, and fractional quantum mechanics. In biology, they are employed to model population dynamics, epidemic spread, and biological transport phenomena. Moreover, in finance, FIDEs are instrumental in analyzing financial markets with memory effects and non-local interactions, contributing to risk assessment and option pricing. Additionally, FIDEs find applications in control theory, signal processing, image processing, and many other areas where systems exhibit complex behaviors beyond the reach of classical differential equations. The ability of FIDEs to capture intricate dynamics and non-local interactions makes them indispensable in understanding and predicting real-world phenomena, underscoring their significance and wide-ranging applicability. For additional applications, interested readers can refer to previous studies [1–7].

Recent developments have seen the emergence of various numerical and analytical strategies for solving FIDEs. These strategies encompass a range of numerical approaches, including collocation methods [8–10], reproducing kernel methods [11,12], wavelet methods [13–15], spectral methods [16,17], finite difference schemes [18,19], differential transform methods [20], and convolution quadrature (CQ) techniques [21]. Among these, the Runge–Kutta convolution quadrature (RKCQ) techniques stand out as widely employed [22–25]. Additionally, a variety of analytical schemes have been proposed, as seen in previous studies [26–28].

Limited numerical methods are specifically developed for the most general form of FIDEs, and existing methods may not be applicable for certain types of FIDEs. Moreover, commonly used methods resort to standard numerical integration techniques, resulting in low-order methods. This motivates our endeavor to develop an efficient high-order numerical scheme for solving nonlinear FIDEs with weakly singular kernels of the form

where β > 0 , y 0 , y ( t ) ∈ R , and the functions p : [ 0 , T ] × R → R , and K : [ 0 , T ] × [ 0 , T ] × R → R are sufficiently smooth to ensure the existence of a unique solution. Here, D t α represents the Caputo differential operator of order α ( 0 < α < 1 ), defined as

The preference for Caputo’s derivative over alternative fractional derivatives, such as Riemann–Liouville, is motivated by its ability to offer a more accessible and transparent physical interpretation to initial values [29].

The nonlocal character of the fractional integral and the singularity of its kernel pose substantial challenges in the numerical treatment of FIDEs, surpassing the complexities associated with standard differential equations. Direct approaches to discretizing FIDEs involving a fractional operator with respect to time require advancing the solution to t + γ , where γ > 0 , by approximating fractional α -integrals of the form I α f ( t + γ ) , where

with k ( t ) = t α − 1 ∕ Γ ( α ) and α ∈ ( 0 , 1 ) . A standard discretization of (2) is obtained by CQ based on a Runge–Kutta scheme

Here, the convolution weights ω n can be expressed as

where t n = n h with a time step h > 0 , and e n ( ⋅ ) is a function dependent on the specific Runge–Kutta scheme [30].

This requires that the entire solution history is stored and used throughout the computation. This can be expensive in terms of both computational and memory costs. Precisely, to compute up to time T = N h , using Formula (3) requires O ( N ) memory and O ( N 2 ) arithmetic operations. To address these challenges, various strategies have been proposed, aiming to alleviate the associated computational and memory costs in the numerical treatment of FIDEs. For a qualitative discussion comparing various methods and their impact on memory and computational efficiency, we refer the reader to previous studies [31–34].

Another approach is to seek a suitable sum of exponentials to approximate the kernel function (2), i.e.,

where w i and s i , i = 1 , … , N ε , are in general complex numbers, and ε > 0 is a given precision. The key is to determine w i and s i such that the desired accuracy up to O ( ε ) can be achieved. Over the past decades, the sum of exponentials methods have been proven effective in many applications of scientific computing.

In our approach, the fractional integral I α f ( t + γ ) is split into a local term

and a history term of the form k γ ∗ f ( t ) with k γ ( t ) = k ( t + γ ) . Inspired by Yan et al. [35], we develop a sum of exponential scheme that gives the approximation ∣ k ( t ) − S ( t ) ∣ < ε , for t ∈ [ δ , T ] , with

The proposed scheme reduces the memory requirement to O ( log N ) and the computational cost to O ( N log N ) , with ε being the accuracy in the computation of the convolution weights. Thus, the approximation to the history term is given by the convolution S ∗ f ( t ) , where S is of the form (4). This approach has three main advantages. The first is that the convolution S ∗ f requires only local information of f to advance. The second advantage lies in the accurate approximation of k at a positive distance γ from its singularity using S with a relatively small number of terms N ε . Finally, only real arithmetic is required as w i and s i , i = 1 , … , N ε , are the real numbers in the proposed scheme.

To approximate the local part of the fractional integral, we employ an integral deferred correction (IDC) method [31,36–39]. This scheme involves two stages: a “propagation stage,” utilizing a low-order quadrature method to advance the solution over a certain time step, and a “correction stage,” modifying the solution for improved accuracy by solving a system of algebraic equations. The correction stage can be performed multiple times to achieve higher-order accuracy. IDC schemes can be more efficient than traditional higher-order methods, especially for problems with stiff or singular solutions [40].

The present fully discrete time-stepping scheme retains these properties while benefiting from the high-order convergence of the IDC method. The main result of this article, as articulated in Theorem 4, provides an estimate of the local truncation error of the proposed method. Specifically, for α ∈ ( 0 , 1 ) and P = ( γ 0 , … , γ q 1 ) ⊂ [ 0 , τ n ] being a set of q 1 + 1 distinct quadrature nodes, the numerical solution V k to FIDE (1), obtained by the proposed IDC scheme at the kth correction iteration, satisfies

where τ n = t n + 1 − t n and q 2 + α is the order of the method used to determine a provisional solution at the 0th correction stage. Theorem 4 states that the local order of the IDC method increases by at least α each correction iteration. The proposed scheme has the advantage of being easy to implement, as it does not require sophisticated memory management and the optimization of quadrature parameters is simple enough.

The rest of this article is organized as follows. In Section 2, we provide an overview of the proposed methods. Section 3 discusses the sum of exponentials method. High-order time-stepping methods to solve FIDEs and their error analysis are introduced and discussed in Sections 4 and 5. Finally, as an illustration of the efficiency of the method, some numerical results are presented in Section 6.

2 Overview

Throughout this article, we use the notations

Furthermore, for the functions p and K in Eq. (1), we assume that

1. p is bounded and continuous on I ( T ) × R ;

2. K is bounded and continuous on Ω T ;

3. p satisfies a Lipschitz condition with respect to the second argument;

4. K satisfies a Lipschitz condition with respect to the third argument.

It is well known that the initial value problem

is equivalent to the Volterra integral equation

in the sense that a continuous function is a solution of (5) if and only if it is a solution of (6).

Owing to (6), we conclude that the solution of (1) is equivalent to the solution of the integral equation,

Using Dirichlet’s formula

we obtain that

Using the variable substitution τ = s − ω t − ω , we can express (7) as

where

Observe that y ( t ) satisfies the integral equation

where

and

Assume that

Then, we obtain the estimate

where B ( ⋅ , ⋅ ) is the beta function. The last estimate and the fact that 1 − β > 0 imply that K ˜ ( t , ω , y ) is bounded and continuous in Ω T .

From the aforementioned discussion, one can easily obtain the following result.

3 Approximation of the history term

An efficient approach to approximating the kernel t α − 1 on the interval [ δ , T ] with a uniform absolute error ε has been proposed in previous studies [35,41], which is a sum of exponential scheme that we will call throughout this work the SOE approximation.

Based on Theorem 2, if δ = T N T and ε is fixed, then we can determine that N ε is of order O ( log N T ) if T ≫ 1 , or O ( log 2 N T ) if T ≈ 1 .

Using the SOE approximation (14) for the kernel ( t − s + γ ) α − 1 in H ( t , γ ) , we can express the history term (11) as

where

and

Thus, we require a technique to approximate (16). On the other hand, quadrature methods can be employed to approximate equation (15) effectively if f is known. We can use the data generated by the IDC at the nodes γ 0 , … , γ q 1 to approximate f at the quadrature nodes and treat it as a known function in order to apply a quadrature method.

4 Approximation of the local term

At t = t n , we consider (13)

in ( 0 , τ n ] , where H ( γ ) is given. Assuming V represents the numerical approximation to Y , with V k = V ( γ k ) for each k = 0 , … , j − 1 , we can compute V j by solving the following equation:

where H j = H ( γ j ) and ω k , j ∈ { γ 0 , … , γ j } are the weights that characterize the scheme.