1. Introduction
Fractional calculus has been playing a very important role in scientific computations. Scientists are able to describe and model many physical phenomena with fractional-order differential equations. As a result, fractional-order differential operators are widely used to solve systems by developing more accurate models [
1,
2,
3,
4]. The nonlocal property of the fractional-order operators makes them more efficient for modeling the various problems of physics, fluid dynamics and their related disciplines [
1,
5,
6,
7,
8,
9]. For example, consider a thin rigid plate of mass
and area
R immersed in a Newtonian fluid of infinite extent and connected by a massless spring of stiffness
K to a fixed point. A force
is applied to the plate. Assume that the spring has no effect on the fluid and that the area of the plate is large enough to produce the fluid adjacent to the plate, whereas stresses
can be defined by the following relation:
where
x is the distance of a point in the fluid from the spring to the submerged plate. By some assumptions discussed in [
10], the dynamics of the system are given by
where
. Equation (
2), with some assumptions considered in [
10], takes the following form of a Bagley–Torvik-type problem solved in (
Section 6, Example 1):
The existence and uniqueness results of fractional-order differential equations (FDEs) have been investigated extensively in the literature. Some of them are presented as follows: Fazli and Nieto [
11] investigated the existence and uniqueness of the solution of FDEs of Bagley–Torvik type by considering the existence of coupled lower and upper solutions. Pang et al. [
12] investigated the existence and uniqueness of the solution of the generalized FDEs with initial conditions by proposing a novel max-metric containing a Caputo derivative. Abbas [
13] studied the existence and uniqueness of the solution of FDEs by using Banach’s contraction principle together with Krasnoselskii fixed point theorem. For more works on existence and uniqueness results, we refer the reader to [
14,
15,
16,
17].
As most FDEs do not have closed-form solutions, different numerical techniques, including the finite difference method, variational iteration method and spectral methods, are preferably used. Among them, spectral methods have received considerable attention from the fractional community for solving FDEs, both ordinary and partial. Spectral methods are classified into three types, known as the collocation, tau and Galerkin methods. The basic idea of the spectral methods is to write the solution as a linear combination of basis vectors of global polynomials, typically Legendre, Jacobi and Chebyshev. The speed of convergence is considered the best advantage of the spectral methods, as the rate of convergence is exponential in these methods, which gives a high level of accuracy. Many efficient spectral techniques are obtained in the literature using the various global polynomials [
18,
19,
20,
21].
The construction of the operational matrices of fractional derivative operators defined with singular or nonsingular kernels has played a key role in the development of spectral methods. Many researchers have worked on the construction of the operational matrices of fractional derivatives using different types of global polynomials. For example, Benattia et al. [
22] introduced the operational matrix of the fractional derivatives to develop a numerical method that is based on the Chebyshev wavelet for solving FDEs. Saadatmandi et al. [
23] derived an operational matrix of derivatives of fractional order using the fractional-order Chebyshev functions. They also extended the results of [
23] for solving the coupled system of FDEs with variable coefficients [
24]. Additionally, Bharway et al. [
25] introduced a new shifted Chebyshev operational matrix of fractional integration for solving linear FDEs. Moreover, Talib et al. [
26] developed a new operational matrix based on the orthogonal shifted Legendre polynomials to numerically solve the fractional partial differential equations. Meanwhile, Rahimkhani et al. [
27] introduced a Bernoulli wavelet operational matrix of fractional integration for obtaining the approximate solution of a fractional delay differential equation. Kazem et al. [
18] derived an operational matrix that generalized the results presented in [
19]. Recently, Dehastani et al. [
28] calculated modified operational matrices of integration and pseudo-operational of fractional derivatives for the Lucas wavelet functions to compute the numerical solution of fractional Fredholm–Volterra integro-differential equations. Moreover, Dehastani et al. [
29,
30] also derived operational matrices of fractional-order derivatives and integration for fractional-order Bessel functions and fractional-order hybrid Bessel functions. In the derived numerical techniques, the operational matrices are applied to reduce the FDEs to a system of algebraic equations.
Dehestani et al. [
31] also presented a novel collocation method based on the Genocchi wavelet for the numerical solution of FDEs and time-fractional partial differential equations with delay.
Motivated by the aforementioned works, we extend the study of the spectral methods by constructing a numerical algorithm that is based on the fractional-order derivative operational matrix of VLPs in Caputo sense, together with the spectral tau method and spectral collocation method. The basis vectors of VLPs are used to approximate the solution of the problems. The derivative terms are approximated by using the fractional-order derivative operational matrix of VLPs. It is important to mention that the proposed algorithm is computer-oriented and is capable of reducing the FDEs to a system of algebraic equations, which greatly simplifies the problems. Subsequently, we use the operational matrices approach together with the spectral tau method in the case of linear FDEs, and the operational matrices approach together with the spectral collocation method in the case of nonlinear FDEs. Our proposed numerical algorithm produces highly efficient numerical results as obtained otherwise in the literature [
32,
33,
34].
The novel aspects of our proposed study are the development of the new fractional-order derivative operational matrix of VLPs in Caputo sense and the construction of the numerical algorithm that is based on this newly developed operational matrix. To the best of our knowledge, this is the first result where the numerical algorithm is presented by using the operational matrix of VLPs. Moreover, the proposed numerical algorithm is fit to solve both linear and nonlinear FDEs with initial conditions. In addition, our proposed method has advantages over other methods, such as the Homotopy perturbation method, because, in our case, the perturbation, linearization or discretization are not necessary to be implemented.
The structure of this paper is set in the following way. In
Section 2, we discuss the VLPs along with their properties. In
Section 3, the Vieta–Lucas operational matrix of fractional-order derivatives is derived. In
Section 4, the numerical method is developed by using the operational matrices of VLPs. In
Section 5, the error bound is determined. In
Section 6, the accuracy and the stability of the proposed method are analyzed by taking some numerical examples. In
Section 6, we conclude and give the summary of this paper.
6. Illustrative Examples
In this section, we give some numerical examples to show the accuracy of our proposed method.
Example 1. Consider the following fractional Bagley–Torvik equationsubject to the initial conditions with integer order The source term is as follows: The exact solution of the problem in Example 1 is: Now, we apply the technique that is described in Section 4.1 by choosing the first three terms of VLPs. We may write the approximation solution as Now, we have Vieta–Lucas polynomialsand The Vieta–Lucas operational matrices can be expressed as The residuals can be evaluated aswhere Now, using initial conditions, we have: Moreover, using the inner product of the residual with the Vieta–Lucas polynomials, we obtain a system of equations. If we take one equation from this system and two equations from Equation (65), then, by simultaneously solving these equations, we obtain hencewhich is the exact solution. Remark 1. The numerical results computed using our method are compared with the method of [33] by choosing various n. We observe that our proposed method produces efficient numerical results as compared to the numerical results obtained by using the method of [33] (see Table 1 and Table 2 and Figure 1). Moreover, for a small value of , the exact solution and the approximate solution computed by using our method coincide (see Table 1 and Figure 1). Remark 2. The numerical results of Example 1 computed at by using our method are compared with the results obtained by using the methods of [32,34] at . We observe that, for a small value of , the approximate solution obtained using our method coincides with the exact solution of Example 1 (see Table 3 and Table 4). However, the exact solution and the approximate solution computed by using the methods of [32,34] coincide at . This shows that our proposed method is numerically more efficient. Example 2. Consider the following linear initial value problem [38]: The exact solution of the problem is [39]. To solve the problem, we use the technique described in Section 4.1. The absolute error for and is shown in Table 5. An error plot is also shown in Figure 2 for these values. We can see in Table 5 that a good approximation has been achieved by using some initial terms of VLPs. Moreover, the numerical results for when and and 1 are plotted in Figure 3. The exact solution for is . It can be noted that the numerical solution converges to the analytical solution when α approaches 1. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state by also upon all of its historical states. Example 3. Consider the following initial value problem [10]:subject to the initial condition with integer order The source term is given as The exact solution at is given below: We test the behavior of our proposed method by solving Example (3) at various values of n. In Table 6, we list the and errors for different values of n. We compare the numerical results obtained by using our proposed method with the numerical results obtained in [10]. It can be observed that the errors computed by using our method are much smaller than those computed by using the method presented in [10]; see Table 6. This highlights the efficiency of our method for this problem. Note that the symbol means that the result for n is unavailable for the method [10]. Example 4. Consider the following nonlinear initial value problem [40]: The exact solution of the problem is [39]. We have solved the problem using the technique described in Section 4.2. The absolute error for and is shown in Table 5. An error plot is also shown in Figure 4 and Figure 5 for these values. We can see in Table 7 that a good approximation has been achieved. Numerical results for when and and 1 are plotted in Figure 6, along with the exact solutions at the given values of α. It can be noted that, as α approaches 1, the solution of the FDEs approaches that of the integer-order differential equations. We also analyze the nonlocal behavior of the fractional derivative by computing the results at various non-integer values of α, which highlights the advantage of using the fractional derivatives, as the next state of the system depends not only upon its current state but also upon all of its historical states (Table 8). Example 5. Consider the following nonlinear initial value problem: We solved this problem by using the same technique as described in
Section 4.2 with
.
The exact solution of the problem is , whereas the source term is .
The operational matrices can be expressed as
where
Now, using the initial condition, we obtain three equations:
Meanwhile, using the technique in
Section 4.2, we obtain the following equation:
Now, we collocate Equation (
74) at the first root of
and we obtain
.
Now, solving Equations (
73) and (
74), we obtain
which is the exact solution.
Example 6. Consider the following non-homogenous multi-order fractional problem:subject to the following initial conditions The source term is as below:
The exact solution corresponding to
,
is given below:
We can observe in Example 6 that a good approximation of the function has been achieved while using
as a scale level. The absolute error
Table 9 at different scale levels is given below. (see
Figure 7).
Moreover, a comparison between the exact solution and approximate solution at different scale levels for the values of
z is given in
Table 10.