Study of Fractional-Order Boundary Value Problems Using Nondiscretization Numerical Method

. Tis paper is devoted to present a numerical scheme based on operational matrices to compute approximate solutions to fractional-order boundary value problems. For the mentioned operational matrices, we utilize fractional-order Bernoulli polynomials. Since fractional-order problems are usually difcult to treat for their corresponding analytical or exact solutions, therefore, we need sophisticated methods to fnd their best numerical solutions. Te presented numerical scheme has the ability to reduce the proposed problem to the corresponding algebraic equations. Te obtained algebraic equations are then solved by using the computational software MATLAB for the corresponding numerical results. Te used method has the ability to save much more time and also is reliable to secure the proper amount of memory. Several examples are solved by using the considered method. Also, the solutions are compared with their exact solution graphically. In addition, the absolute errors for diferent scale values are presented graphically. In addition, we compare our results with the results of the shifted Legendre polynomials spectral method.


Introduction
Newton and Leibnitz in the seventeenth century established the theory of calculus.Te notations of derivatives and integrals we use today have been introduced by Leibnitz.Later on, the concept of derivative and integral was extended from integer to any real order.First time in 1819, Lacroix gave the concept of noninteger order derivative [1].Te very frst application was investigated by Abel in 1823 [1].Terefore, Fourier, Liouville, Riemann, Grünwald, and Letnikov [2,3] gave tremendous attention to the said area.Arbitrary order diferentiation and integration provide a generalization to the classical order but fractional-order diferentiation and integration do have not unique defnitions.Terefore, various researchers have introduced various defnitions of arbitrary order derivatives and integrations.Among all these, the defnitions of Caputo and Riemann-Liouville are extensively applicable.
Te researchers have investigated that almost every model of physics and biology consists of arbitrary order derivatives.Fractional-order diferential equations (FODEs) are widely used in various felds of science and technology [4][5][6] and [7].Furthermore, various properties of numerous materials such as hereditary properties are described by FODEs.Terefore, the signifcance of FODEs has appealed to researchers for further development in the theory.Tere are many areas in the theory of FODEs, one of them being the study of the numerical approach of fractional-order initial value problems (FOIVPs) and fractional-order boundary value problems (FOBVPs).It has been observed that the initial and boundary value problems with ordinary derivatives have been greatly worked on by the researchers but FOIVPs and FOBVPs are in their beginning stages and need proper further exploration.Here, we remark that the aforementioned area has signifcant applications in various other areas.For instance, the author of [8] studied nonlocal kinetics using fractional calculus applications.In addition, researchers [9] investigated a blood therapy model using the concept of new fractional diferential operators very recently.Te mentioned operators have been used to evaluate the difusion process in [10].Recently, two diferent dynamical systems of prey-predator and infectious diseases were studied using FODEs in [11].
Te analytical solutions of most of the FODEs cannot be found because of the complexity of fractional order.Terefore, approximate analytical or numerical solutions are obtained for the said area.For this purpose, various analytical and numerical methods have emerged, for instance, eigenvector procedure [12], perturbation tools [13], iteration techniques [14], and transform methods [15].Authors studied numerically discrete-time prey-predator model and SIR-type model by using concepts of fractional calculus in [16] and [17].Authors [18] used the diferential transform method to study some systems of FODEs.Authors used decomposition schemes to compute approximate solutions of some problems in [19].Te fnite diference method has been applied in [20].Also, the Tau method has been applied in [21] and [22].Collocation techniques have also been utilized in [23].Wavelet analysis [24] has been used for numerical results increasingly.Recently, authors studied the time-fractional model of generalized Couette fow of couple stress nanofuid with heat and mass transfer in [25].Authors [26] studied fractional-order Caudrey-Dodd-Gibbon equations by using an integral transform.In addition, researchers [27] computed the approximate solution to wavelike equations.Additionally, authors [28] computed an approximate solution of the Noyes-feld model for the timefractional Belousov-Zhabotinsky reaction.In the same way, authors [29] have used the transform method to study fractional-order Swift-Hohenberg equations.
Here, we remark that spectral methods based on the operational matrices (OMs) for various orthogonal polynomials such as Jacobi, Legendre, and Laguerre have been developed for solving FODEs numerically [30].Te said OMs have been constructed from the mentioned polynomials by the techniques of discretization which occupies extra memory and consumes much time.Furthermore, the aforesaid methods have been applied to FOIVPs in a variety of cases.However, FOBVPs are very rarely investigated.To deal with boundary value problems, some extra operations are needed.For this purpose, we establish a numerical algorithm based on OMs for Bernoulli functions which are not orthogonal.We construct the said OMs through Bernoulli functions without discretization to save memory and time.By using these OMs, we establish a numerical algorithm to fnd the numerical solution for the considered problems.A new operational matrix corresponding to boundary conditions has also been obtained.
Inspired by the above discussion, we consider the given class of FOBVPs [31] under the fractional order described by where k 1 , k 2 are real constants and f ∈ C(I, R) is a source function with I � [0, 1].
In addition, we also extend the aforementioned scheme for the coupled system of FOBVPs as described by with BCs given by where the fractional orders are described by ) are any real constants.Here, we frst establish the mentioned operational matrices by using Bernoulli polynomials.By using these OMs, we convert the considered system of FODEs to Sylvester-type equations given by where the notions matrices and X is an unknown matrix to be computed.We solve the obtained equation ( 4) for unknown matrix X by using computational software like MATLAB.
Here, we demonstrate some novelty of our work.Boundary value problems (BVPs) constitute a very important branch of applied analysis because many engineering and physical problems to understand real-world procedures or phenomena are devoted to BVPs.Te concerned branch when studied under the fractional-order derivatives and integrals sense further enhances this feld of study, as many real-world processes involve short or long nature memory efects which cannot be detected through applications of ordinary derivatives.Hence, with the help of fractionalorder derivatives, we can explain long-and short-memory efects more clearly for best use in the area of engineering and physical sciences.Furthermore, fnding the exact solutions to fractional-order boundary value problems is a very tedious job and cannot be achieved in a simple way.Terefore, various numerical techniques have been developed to treat such problems.Among the available numerical methods, spectral numerical schemes based on OMs are powerful techniques.However, the operational matrices when established by using the discretization technique 2 Journal of Mathematics utilize more memory and time-consuming process.Terefore, omitting the discretization and obtaining the required matrices are our goal.In this work, we have utilized Bernoulli polynomials a powerful tool to establish operational matrices of fractional-order integration and diferentiation.Te concerned matrices convert the considered FODE to an algebraic equation of Sylvester-type matrix equation.Te obtained matrix equation is then solved by the use of the Gauss elimination procedure using computational software like MATLAB or Mathematica.In this way, we obtain the numerical scheme which needs no discretization as earlier used in research.Furthermore, the proposed polynomials have the ability to produce more accurate results than wavelet and other diferent methods.Here, it is interesting that in the last few years, spectral methods have been used very well for diferent problems of FODEs.For some more frequent works, we refer to [32][33][34][35][36] and [37].Here, researchers have used optimal control and Bernstein and Legendre polynomials to study various problems of FODEs.
Result devoted to stability has also been deduced by following the methodology of [38].
Te rest of the work of this article is organized as follows.In Section 2, we provide some necessary defnitions from the fundamental theory of calculus and FOBPs required in the subsequent development.Te function approximation procedure and the OMs of fractional integration and differentiation are given in Section 3. Also, in the same Section 3, we give new OM corresponding to BCs.Section 4 is devoted to the establishment of numerical algorithms for the solutions of FOBVPS and coupled system FOBVPS.Part 5 is related to investigating the convergency of the proposed method.In Section 6, the proposed method is applied to various examples, and the concerned numerical results are presented through graphs.Also, the exact and numerical solutions are compared in the same Section 6. Te last Section 7 is devoted to the conclusion of this article 7.

Preliminaries
Here, we give some important defnitions from the basic fractional calculus theory and Bernoulli polynomials necessary for the subsequent development.
where c, k 1 , k 2 are constants.

Te Fractional-Order Bernoulli Polynomials.
In this section, we defne FOBPs and provide some of their properties.

Function Approximation
In this section, we give the procedure for how to approximate a function in terms of FOBPs.Here, we defne I � y: Now, suppose that is the set of FOBPs and Here, Z m is fnite-dimensional and closed vector space.So, every element g of L 2 (I) has a unique best approximation g 0 out of Z m , such that which can be written as where 〈, 〉 stands for inner product.Now, as g 0 ∈ Z m , the unique coefcients c 0 , c 1 , c 2 , . . ., c m are exist such that where M � m + 1, C M×1 is the coefcients matrix and B α M×1 (x) is fractional-order Bernoulli functions vector given by For evaluating C M×1 , consider Using ( 11), we get 4 Journal of Mathematics where Tis implies that which gives where where 3.1.Derivation of OMs.Here, in this section, we derive Bernoulli-type OMs of fractional-order integration and diferentiation.Moreover, we construct a new OM corresponding to some boundary value.

Theorem 3. Let B α M×1 (x) be the fractional-order Bernoulli functions vector, then
where P (],α) M×M is given by Proof.Consider Using the properties of the Riemann-Liouville fractional integral operator, we have where Now, we approximate x αr+] as Using (30) in (27), we get where (31) implies that Hence, we have Journal of Mathematics □ Lemma 4. Let B α i (x) be the FOBP, then Proof.By applying the properties (iii − v) of Caputo arbitrary order diferential operator in (11), one can easily prove the lemma.

□
Theorem 5. Let B α M×1 (x) be the FOBPs vector, then where D (],α) M×M is given by Proof.Consider By applying the properties of Caputo arbitrary order diferential operator, we have Now, we approximate x αr− ] by Using (41) in (38), we have Journal of Mathematics where (44) Also, using Lemma 4, we get from (44) Hence, from (44) and (45), we have □ Theorem 6.Let B α M×1 (x) be a FOBPs matrix and let ϕ(x) be any function given ϕ(x) � lx n , l � 0, 1, 2, . . ., l ∈ R and where M×M is given by where Φ can be calculated by using equation (38). Proof.
Using the well-known property of the Beta function given by in equation (47), we have which can be approximated in FOBPs as Hence, we have the required result. □

Numerical Algorithms
Here, we show the fundamental importance of OMs developed in the previous section by applying them to various FOBVPs.For this purpose, we consider the following two steps.

Numerical Scheme for Scaler
Problem (1).Here, we develop a numerical scheme for a single problem given in (1).
To obtain the solution in terms of FOBPs, we assume that Applying I ] and property (ii) for the Caputo diferential operator, we get from (54) Using the given boundary conditions, we have c 0 � u 0 and Putting the values of c 0 and c 1 in (55), one can get Using the afore-established OMs, and after simplifcation, we get from (57) where x and ϕ � x.Now by applying D ] 1 to (58), we have Putting (54), (58), and (59) in (1) yields where After simplifcation, we obtain where Tis is a Sylvester-type equation.Calculating C M×1 from ( 23) and substituting back in (20), one can have the desired approximate solution of the given FOBVP.

Numerical Scheme for Coupled System (2) with BCs (3).
Here, we develop a numerical scheme for the computation of an approximate solution for a coupled system (2) with BCs (3).In this regard, we assume that Applying I ] , I ω and property (ii) for the Caputo differential operator, we get Using the given BCs, we have Inserting ( 65) in (64), we obtain After simplifcation of (66), we get where and D ω 2 to (67), we have Putting ( 63), (66), and (68) in coupled system (2) yields where

Journal of Mathematics
By taking transpose and rearranging the terms of (69), we have where O M×M and O M are zero matrix and zero vector, respectively.By simplifcation, (71) can be written as Tis is a system of Sylvester-type equations.Calculating the matrices C T M×1 and E T M×1 from (72) and substituting back in (67), one can obtain the desired solution.

Convergence Analysis
Here, we discuss the convergence of the proposed method developed above.For this purpose, we provide the following theorems to show that the error matrices E (],α) I and E (],α) D getting smaller and smaller as we increase the number of FOBPs.
Theorem 7. Let V ⊆ H be a closed subspace of H.space H with fnite dimension and S � v 1 , v 2 , v 3 , . . ., v n   form basis for V. Let h ∈ H and h have the unique best approximation v * in V. Ten, where (75) Ten, we have lim Proof.Te proof is given in [35].

I
of OM P (],α) defned by is bounded by the following inequality: where Proof.It has been proved in [35].
given by is bounded by the following inequality: where Journal of Mathematics (82) Proof.It has been proved in [35].
From the above theorems, we come to the conclusion that as we increase the number of FOBPs, the error matrices E (],α) I and E (],α) D approach to zero.Furthermore, here, we deduce some results regarding stability following [38].Let U(x) ∈ C[0, 1] be exact solution, and U M (x) be the approximate solution of the proposed problem, then the error bounded is computed as where B M , U M are the maximum values of B α M (x), U M (x) in [0, 1], respectively.Also, if M ⟶ ∞, then ε ⟶ 0. Furthermore, the given remark holds.

□
Remark 11.For K > 0, the relation Theorem 1 .Let U(x) be the exact solution to the problem under our consideration, and U M (x) �  U(B α M (x)) T be the approximate solution of the proposed problem, where  U is the Bernoulli coefcient matrix that is determined by solving the algebraic equation (62), then using Remark 11, the method is stable if K < 1.
Proof.Consider the exact solution as U(x) and  U(x) be any approximate solution to the proposed problem, we have using the above remark Taking maxima over [0, 1] of both sides and rearranging the terms, (85) implies that Hence, the solution is stable.Furthermore, if there exists a nondecreasing function say φ: (0, 1) ⟶ (0, ∞), such that φ(ε) � ε, where φ(0)� 0, then from relation (87), we conclude that which further yields that the solution of the proposed linear problem of fractional order is generalized stable.□

Numerical Examples
Here, we solve some examples through the scheme we developed in Section 4 and present the concerned numerical results through graphs.Also, we compare the exact and numerical solutions of the said examples.
Let us assume that the given problem has an exact solution at α� 2 and β� 1 which is given by Additionally, f(x) can be approximated as In Figure 1, we compare the actual and approximate solutions at diferent fractional orders, and the corresponding absolute errors are also presented graphically.We see that as the fractional order approaches the corresponding integer value 2, the concerned curve tends to the corresponding integer order curve.Te concerned error also reduces as the order increases.In Figure 2, we provide a graphical presentation of numerical solutions at various scale levels.Here, we see that as the value of M is increasing, the corresponding absolute error also decreases.Hence, the mentioned method is also scale-oriented.
Example 14.Here, we take another FOBVP given by We solve this problem for various fractional orders.Let us assume that the given problem has an exact solution at ]� 2 and ] 1 � 1 which is given by with the source term given by In Figure 3, we compare the exact and approximate solutions at diferent fractional orders, and the corresponding absolute errors are also presented graphically.We see that as the fractional order approaches the corresponding integer value 2, the concerned curve tends to the corresponding integer order curve.Te concerned error also reduces as the order increases.In Figure 4, we provide 12 Journal of Mathematics a graphical presentation of numerical solutions at various scale levels.We provide a graphical presentation of numerical solutions at various scale levels.Here, we see that as the value of M is increasing, the corresponding absolute error also decreases.Hence, the mentioned method is also scale-oriented.
Example 15.Consider the system of FOBVPs [36] with BCs given by We solve this problem under the set of parameters Let us assume that this problem has the exact solution given below under the given parameters Journal of Mathematics 13 Ten, f(x) and g(x) are given by  14

Journal of Mathematics
Here, we give a plot at M � 6 and diferent values of fractional order ω, ].In Figure 5, we compare the actual and approximate solutions at diferent fractional orders, and the corresponding absolute errors are also presented graphically.We see that as the fractional order approaches the corresponding integer value 2, the concerned curves tend to the corresponding integer order curves.Te concerned error also reduces as the order increases.In Figures 6 and 7, we provide plots of approximate solutions at various scale levels.Also, we provide a graphical presentation of numerical solutions at various scale levels.Here, we see that as the value of M is increasing, the corresponding absolute error also decreases.Hence, the mentioned method is also scaleoriented.In Figure 7, we provide plots of absolute errors at diferent scale levels.Here, we compare the numerical results for Example 15 with those given in [36] in Table 1.We see that the existing method produces slightly good results than the numerical results obtained in [36] for the given problem by using the shifted Legendre polynomials spectral method.In addition, here in Table 2, we compare the CPU time of our proposed method with that of [36].

Journal of Mathematics
Let us assume that this problem has the following exact solution under the set of parameters.(102) Here, we compare the actual and numerical solutions and the corresponding absolute errors in Figures 8 and 9, respectively.We compare the actual and approximate solutions at diferent fractional orders, and the corresponding absolute errors are also presented graphically.We see that as the fractional orders approach the corresponding integer value 2, the concerned curves tend to the corresponding integer order curves.Te concerned error also reduces as the order increases.Furthermore, we provide a graphical presentation of approximate solutions at various scale levels and the corresponding absolute errors in Figures 10 and 11, respectively.In addition, we provide a graphical presentation of numerical solutions at various scale levels.Here, we see that as the value of M is increasing, the corresponding absolute error also decreases.Hence, the mentioned method is also scale-oriented.

Conclusion and Discussion
In this article, we have investigated diferent classes of FODEs under boundary conditions for numerical solutions.We have used FOBPs along with some important properties to construct some operational matrices corresponding to fractional-order derivative and integration.Based on these matrices, we have converted the proposed problems to a system of Sylvester-type operational matrices.Upon using MATLAB, we solved several examples to demonstrate the procedure.From the numerical investigation, we see that the method is powerful and can be used to investigate various FODEs for numerical solutions.Te efciency of the method can be further improved by enlarging the scale level.Te greater the scale level, the greater the accuracy, and vice versa.Also, we have computed the maximum absolute error by using diferent scale levels.In some examples, we have also investigated the problem by using various fractional orders to simulate the results.From the numerical results, we see that the operational method based on FOBPs can also be used as a powerful spectral method to handle FODEs for numerical purposes.We have compared our numerical results with those given in [36] by using the shifted Legendre polynomials spectral method.Our results are good than those given in the mentioned reference.Te CPU time of the proposed method has been compared with the CPU time of 18 Journal of Mathematics the shifted Legendre spectral method.Te CPU time of the proposed method is better than the spectral numerical method based on Legendre polynomials.In the future, we will extend the mentioned OMs method to variable order problems.Also, we will study FODEs involving nonsingular type derivatives by using our proposed method.

Absolute error at M = 6
Absolute error at M = 8 Absolute error at M = 10

6 Figure 11 :
Figure 11: Absolute errors at diferent scale levels of M.

Table 1 :
[36]arison between numerical results of Example 15 by our proposed method and that used in[36].

Table 2 :
[36]arison of CPU time to compute solution in Example 15 by our proposed method and that used in[36]using Cori7 generation HP machine of 1024 GHz.