Numerical Approximation of Fractional-Order Volterra Integrodifferential Equation

Institute of Computing Science and Technology, Guangzhou University, Guangzhou 510006, China Department of Mathematics, Islamia College Peshawar, Khyber Pakhtoon Khwa, Pakistan Department of Mathematics, Division of Science and Technology, University of Education, Lahore, Pakistan Department of Mathematics, Huzhou University, Huzhou 313000, China Hunan Provincial Key Laboratory of Mathematical Modeling and Analysis in Engineering, Changsha University of Science & Technology, Changsha 410114, China


Introduction
Fractional calculus has a large number of applications in engineering and mathematical sciences [1][2][3][4][5]. Problems from engineering and other sciences which involve derivatives or integrals of noninteger orders are large in number and still growing. That is why the research community has showed great interest in the area of fractional calculus. One can find various applications of fractional calculus in numerous phenomena such as frequency-dependent damping behavior of viscoelastic materials [6], control theory [7], economics [8], mass and heat diffusion processes, electromagnetic theory, biological species [9], robotics, and many other engineering problems [10]. Accurate models of systems under consideration can be obtained using fractional differential and integrodifferential equations [11]. Many remarkable works on fractional calculus are available in literature for the approxi-mation of the solution fractional differential or integrodifferential equations, for example, the sinc-collocation method [12], Legendre collocation method [13], Laguerre polynomials [14], Adomian decomposition method [15], Variational iteration method [16][17][18], and their references.
Numerous essential phenomena in nature are resolved using the solutions of fractional integrodifferential equations [11,19]. For example, one can find applications of fractional integrodifferential equations (FIDEs) in electromagnetics [20], heat conduction [21], etc. Due to the wide range of applications, FIDEs have attracted researchers. Therefore, a large number of analytic and numerical methods have been developed for finding the solutions of FIDEs [22][23][24][25][26]. For example, the authors [27] utilized the fractional differential transform method for the solution of fractional integrodifferential equations. A Legendre wavelet method [28] is proposed for the solution of FIDEs. In [29], the analytic solution of FIDEs is obtained using the homotopy analysis method. The author in [30] obtained the analytic solution of FIDEs using the fractional residual power series method.
In many cases, the analytic solutions of fractional differential or integrodifferential equations are hard to obtain, so numerical methods must be used. In this connection, the researchers have developed numerous numerical methods. For instance, the authors in [19,31] have solved FIDEs using the reproducing kernel Hilbert space method. A hybrid collocation method [32] is proposed to solve the FIDEs. The authors in [33] solved the first-order Volterra equation by the collocation method. The author in [34] utilized the shifted Chebyshev polynomial and least squares method for the approximation of the solution of FIDEs. Mahdy and Shwayyea [35] obtained the approximate solution of FIDEs using the shifted Laguerre polynomial pseudospectral method and least squares method. The authors in [36] approximated the solution of Volterra-type FIDEs via Bernoulli wavelet approximation. In [14], Laguerre polynomials are utilized to approximate FIDEs of Volterra type. More work on the solution of FIDEs can be found in [31,33,37,38] and their references.
There are various definitions of fractional derivatives such as Caputo-Liouville's and Riemann-Liouville's [11,39]. These derivatives have certain disadvantages due to nonlocal and nonsingular kernel functions involved in these derivatives. In order to avoid these difficulties, it is better to use the ABC derivative. The ABC derivative contains a nonsingular kernel and therefore can model a phenomenon which cannot be handled by fractional-order derivatives having singular kernels [40]. Some recent articles on the application of the ABC derivative can be found in References [41][42][43][44][45]. In the present work, we consider a FIDEs of Volterra type with the ABC derivative of the form

Preliminaries
Definition 1. The ABC fractional derivative of order 0 < β < 1, of ζ ∈ S 1 ða, bÞ with base point a at t ∈ ða, bÞ, is defined as [46] ABC 0 where S 1 is a Sobolev-space of first-order fitted with L 2 -norm defined over the Ω ⊂ R defined as follows: where MðβÞ is defined as follows: also, E β ðtÞ is a one-parameter Mittage-Leffler (ML) function defined as Definition 2. The ABC fractional integral of order 0 < β < 1 of ζ ∈ S 1 ða, bÞ with base point a at t ∈ ða, bÞ is defined as [46] Definition 3. The Laplace transform of a piecewise continues function ζðtÞ is defined as Definition 4. If 0 < β ≤ 1, then the Laplace transform of the ABC derivative is defined by

Proposed Method
We provide a description of our proposed schemes. Our numerical scheme has three main steps. In our first step, we apply the Laplace transform to the given model, with which the problem is reduced to an algebraic equation. In the second step, we solve the reduced equation in Laplace space for the unknown. Finally, the solution of the original problem is obtained using numerical inversion of Laplace transform. In this work, the numerical approximation of Laplace transform is obtained using two schemes. One is the contour integration method (CIM), and second is the Stehfest method (SM). We discuss these methods in the next sections. First, we apply Laplace transform to Eq. (1); we get which can be rearranged as whereĜ Journal of Function Spaces In order to obtain the solution of problem (1), we apply Laplace inverse on (12), and the Laplace inverse is then approximated using the contour integration method and the Stehfest method, which are discussed as follows.

Contour Integration Method (CIM)
In CIM, we represent the solution ζðtÞ of (1) as Bromwich integral: where ϑ ∈ R and Γ is an initially appropriately chosen line Γ 0 in a complex plane parallel to the y axis, with Ims → ±∞. Then, in (14), ζðtÞ is the inverse transform of ζ b ðsÞ, satisfying the condition that all the singularities of ζ b ðsÞ lie to the left of Γ 0 . However, for our purposes, we assume that b ζðsÞ may be continued analytically in an appropriate way; we shall want to take for Γ 0 a deformed contour Γ in Σ ϑ ϕ = f0g ∪ fs ≠ 0 : | args|<ϕg, which has asymptotic behavior in the left complex plane, with Res → −∞ when Ims → ±∞, which force e st to decay towards both ends of Γ. In our work, we choose Γ as where Letting s = ς + ιζ, we see that (15) represents the left branch of the hyperbola given by where ζ = ±ðς − ϑ − ρÞ cot η are the asymptotes for (17) and and grows into the left half plane. From (15) and (14), we get The approximation of Eq. (18) using the trapezoidal rule with uniform step k is given as 3.1. Error Analysis. While solving the fractional-order problem defined in Eqs. (1)-(3), the first step is the application of Laplace transform to the given problem. The Laplace transform reduces the problem to a time-independent problem in Laplace space, and in this process, no error occurs. Next, the problem is solved in complex space for the unknown, and this process also incurs no error. In the final step, the solution of the given problem is retrieved using the Laplace inverse. We represent our solution as integral (18). We approximate this integral by the trapezoidal rule. During this process, integral (18) converges at different time rates which depends on Γ. In the process of numerical approximation of the Bromwich integral (18), the orders of convergence greatly depend on the step k of the trapezoidal rule and the temporal domain ½t 0 , T. For useful results and strong convergence, a suitable temporal domain is chosen. The next theorem gives the order of quadrature error.

Stehfest Method (SM)
In the Stehfest method, the approximate value of ζðtÞ is given as where the weights w i are given by Solving (12) for the corresponding Laplace parameters s = ðln 2/tÞi, i = 1, 2, 3, ⋯, M: The solution of the original problem can be obtained using (21).

Error Analysis.
The authors [48,49] performed a large number of numerical experiments to study the effect of the parameter on numerical accuracy, and they summarized there experimental work as the following: "If j significant digits are desired, then let M be the positive integer d1:1je. Given M, set the system precision at d2:2Me. Given M and the system precision, calculate the weights w i , 1 ≤ i ≤ 2M, using (22). Then, given the transform b ζðsÞ and the argument 3 Journal of Function Spaces              Journal of Function Spaces t, calculate ζðtÞ in (21)." According to these conclusions, the error is 10 −j+1 ≤ ð b ζ − ζÞ/ζ ≤ 10 −j , where M = d1:1je [50].

Results and Discussion
In order to check the accuracy and efficiency of the two numerical schemes. We consider linear and nonlinear and system of fractional-order VIDEs. In our numerical computations, the accuracy is achieved using the optimal parameters. The optimal values of the parameters utilized for CIM in our work are θ = 0:1, τ = t 0 /T, η = 0:15410, r = 2πr, r = 0:13870, ½t 0 , T = ½0:5,5, ϑ = 2:0. The MATLAB command γ = −N : k : N is used to generate the quadrature nodes for CIM. The parameter used in our numerical computation for the Stehfest method is M.

Problem 1.
Here, we consider the fractional linear Volterra integrodifferential equation [31]: The exact solution of this problem is ζðtÞ = sin ðtÞ. In Table 1, the approximate solutions of ζðtÞ = sin ðtÞ for various values of β and t ∈ ½0, 1 using the two proposed methods are presented, and in Table 2, the absolute errors are shown. The plot of error function for various quadrature nodes N is        Figure 1(a), and the plot of the error function for t ∈ ½0, 1 is shown in Figure 1(b). The graphs of approximate solutions for β are presented in Figure 2(b), and Figure 2(a) displays the plot absolute error vs. error estimate for β = 1. The numerical results using the Stehfest method are obtained by using M = 18. All the plots are obtained using the CIM. It is evident that these methods can solve such type of equations efficiently.
The numerical results for various values of β and t ∈ ½0, 1 are presented in Table 3. The graphs of numerical solutions for various values of β are shown in Figure 3(a), whereas Figure 3(b) presents the graph of error estimate vs. absolute error for β = 1. For t ∈ ½0, 1 and β = 1, the error function is shown in Figure 4(a), and the error function for various N, and β = 1 is shown in Figure 4(b). A similar performance is observed as was observed for Problem 1.

Problem 3.
Here, we consider the fractional nonlinear Volterra integrodifferential equation [51]: The exact solution of this problem is ζðtÞ = e t . The approximate solutions for different values of β and t ∈ ½0, 1 are presented in Table 4. The Stehfest parameter used in this experiment is M = 22. The graphs of approximate solutions for various values of β are depicted in Figure 5(a), whereas Figure 3(b) presents the plot of error estimate vs. observed error for β = 1. In Figure 6(a), for β = 1 and t ∈ ½0, 1, the error function is shown, and in Figure 6(b) for different nodes N and for β = 1 is shown.

Problem 4.
Here, we consider the fractional system of nonlinear Volterra integrodifferential equations [51]: The exact solutions of this system are ζ 1 ðtÞ = t and ζ 2 ðtÞ = 1 + t. The approximate solution of ζ 1 ðtÞ for various values of β and t ∈ ½0, 1 and M = 22 is using the two numerical schemes N = 240, M = 22, corresponding to Problem 3 presented in Table 5, and the approximate solution of ζ 2 ðtÞ for various values of β and t ∈ ½0, 1 and M = 20 is presented in Table 6. In Table 7, the errors obtained using the two proposed schemes are shown. For β = 1 and various nodes N for ζ1ðtÞ, the error function is shown in Figure 7(a), and for ζ 2 ðtÞ, the error function is shown in Figure 7(b). The plots of numerical solutions of ζ 1 ðtÞ for different fractional-order β are displayed in Figure 8(a), and the plots of numerical solutions of ζ 2 ðtÞ for various values of β are displayed in Figure 8(b).
The absolute errors for β = 1, and quadrature nodes N using the CIM are presented in Table 9. The plots of numerical solutions of ζ 1 ðtÞ for different fractional-orders β are displayed in Figure 9(a), and the plots of numerical solutions of ζ 2 ðtÞ for various values β are displayed in Figure 9(b). In Figure 10(a), the plots of error function of ζ 1 ðtÞ for β = 1, and t ∈ ½0, 1 are displayed, and in Figure 10(b), the plots of error function of ζ 2 ðtÞ for β = 1, and t ∈ ½0, 1 are displayed. The plots of error function of ζ 1 ðtÞ for various nodes N and β = 1 and are displayed in Figure 11(a), and the plots of error function of ζ 2 ðtÞ for quadrature nodes N and β = 1 are displayed in Figure 11(b).

Conclusion
In this work, we have successfully applied the Laplace transform to fractional VIDEs. The numerical inversion of Laplace transform is performed using the CIM and SM methods. In our numerical experiments, we considered linear, nonlinear, and system of fractional VIDEs with the Table 9: The absolute error fractional-order β = 1 and different quadrature nodes N using the CIM numerical scheme, corresponding to Problem 5.      Figure 10: In (a), for β = 1 and t ∈ ½0, 1, the error function of ζ 1 ðtÞ corresponding to Problem 5 is shown. In (b), for β = 1 and t ∈ ½0, 1, the error function of ζ 2 ðtÞ corresponding to Problem 5 is shown.  Journal of Function Spaces ABC derivative. The obtained results led us to the conclusion that the two proposed methods are capable of solving such type of equations. However, it was observed that the CIM is more accurate and converges faster than SM. Hence, the two proposed numerical schemes are the best alternative methods for solving such type of equations.

Data Availability
All data required for this paper is included within this paper.

Conflicts of Interest
The authors do not have any competing interests.