A numerical comparative analysis of methods for solving fractional differential equations

Abstract

In this paper, we present two collocation methods utilizing subdivision schemes and Bernstein polynomials as their basis functions to solve Caputo-type fractional differential equations. These methods transform fractional differential equations into systems of linear equations, which can be solved using various suitable techniques. Notably, these methods are versatile, making them applicable to both boundary value problems and initial value problems. To demonstrate their effectiveness, we applied these methods to three test problems of fractional differential equations, including the Bagley-Torvik equation and fractional oscillation equation. Our results reveal that both methods provide a more accurate approximation to the exact solution when compared to various algorithms found in the existing literature, showcasing their efficiency, accuracy, and consistency.

Keywords:

• Boundary conditions
• collocation method
• fractional differential equation
• initial conditions
• subdivision schemes

1. Introduction

Fractional Differential Equations (FDEs) constitute a potent mathematical framework that extends classical differential equations by involving derivatives of non-integer order. Their significance spans across various scientific domains owing to their unique capacity to encapsulate complex systems and phenomena characterized by memory, hereditary properties, or non-local effects. In the realm of physics, FDEs serve as indispensable tools, enabling precise descriptions of anomalous diffusion, viscoelastic materials, and memory-laden systems. Notably, in fluid dynamics, they excel in modeling non-Newtonian fluids and long-range interactions. Furthermore, their application in modeling biological systems, encompassing neuron firing patterns and disease behavior, provides nuanced insights into memory-dependent cellular processes and intricate biological dynamics. Within engineering, FDEs hold practical utility, accurately modeling systems manifesting non-local or memory-dependent behaviors, evident in control systems, signal processing, and material science. Moreover, their role in economics and finance is pivotal, facilitating a deeper understanding of financial markets, price dynamics, risk management, and complex financial time series analysis. This paper aims to illuminate the profound impact and significance of FDEs in comprehending the complexities of real-world phenomena that defy traditional differential equations, offering enriched insights into systems governed by memory effects, hereditary properties, and non-local interactions.

The use of fractional calculus has had a profound impact on the modeling of physical phenomena that are challenging to represent using integer-order derivatives. In the literature, there are various types of fractional operators available for handling these equations. Consequently, specific problems often lack a unique solution that satisfies all differential operators. The most commonly employed and widely recognized operators include Riemann-Liouville, Grnwald-Letnikov, and Caputo. Additionally, the Atangana-Baleanu and Caputo-Fabrizio operators have also proven to be valuable. Within the realm of fractional differential equations, various types exist, each with its own approximate and analytical solutions. Some of the most frequently encountered types include Bagley-Torvik fractional differential equations, Riccati-type fractional differential equations, and Fractal-Fractional differential equations.

Bagley-Torvik fractional differential equations have been addressed using both Riemann-Liouville and Caputo fractional operators. Researchers (Ray & Bera, 2005) have tackled these equations using the Adomian decomposition method, employing the Riemann-Liouville fractional derivative. The use of Caputo derivatives has led to solutions via operational matrices (Ray, 2012), the Chebyshev wavelet operational matrix method (Mohammadi, 2014), hybridizable discontinuous Galerkin methods (Karaaslan, Celiker, & Kurulay, 2016), and the Tau method (Mokhtary, 2016). To solve Riccati-type fractional differential equations, several approaches have been explored, including the Enhanced Homotopy Perturbation Method (HosseinNia, Ranjbar, & Momani, 2008), the Adomian Decomposition Method (Momani & Shawagfeh, 2006), and the utilization of Bernstein polynomials (Yzbasi, 2013). Similarly, the solution of fractal fractional differential equations has been investigated using methods such as the B-spline collocation method (Shloof & Gewily, 2021), operational matrix methods (Shloof, Senu, Ahmadian, & Salahshour, 2021), Taylor series technique for fractal bratu-type equation was introduced by He, Shen, Ji, & He (2020). A numerical approach based on the two-scale fractal transformation and the global residue harmonic balance method was proposed by (Lu & Chen, 2022) for finding the approximated solutions of a fractal modification of the YaoCheng oscillator with He’s fractal derivative. Two numerical approaches, the homotopy perturbation transform method and the method based on the fractional complex transform were presented by Lu & Sun (2021) for finding approximate solutions to the time fractional Boussinesq-Burgers equations with He’ s fractional derivative. Furthermore (Ain et al., 2022), explored the transmission dynamics of MERS-CoV between humans and camels using a non-linear fractional order model. A numerical approach based on the fractional complex transform and the homotopy perturbation method was proposed by Lu & Ma (2023) to solve the space-time fractional Benjamin-Bona-Mahony equation with Caputo fractional derivative. An efficient technique based on fractional complex transform and global residue harmonic balance method was proposed by Lu & Ma (2022). An efficient technique, based on the fractional complex transform and the global residue harmonic balance method, was proposed by Lu & Sun (2023) for solving the fractional modification of the nonlinear oscillator with coordinate-dependent mass. Different types of approaches have been applied by Anjum, Ain, & Li (2021), Ain, Anjum, & He (2021), Tao, Anjum, & Yang (2023) to solve the fractional differential equations.

Let us consider the following fractional differential equation: (1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) with initial or boundary conditions: (1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) (1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) where a₁, a₂, a₃, $\mathbb{A},\mathbb{B},\mathbb{C},\mathbb{E},$ and $\mathbb{F}$ are real constants. Additionally, y(z) is the unknown function to be determined, $D^{{\alpha }_1}$ and $D^{{\alpha }_2}$ are fractional differential operators of orders α₁ and α₂ respectively, where ${\alpha }_1,{\alpha }_2\notin \mathbb{N}.$ The parameter w appears in the inequality constraint related to the fractional orders of α₁ and α₂. Finally, f(z) is the given function.

Caputo’s fractional derivative stands out among various fractional differential operators, such as those by He or Jumarie, due to its unique incorporation of initial and boundary conditions. Its distinct advantages, including seamless extension of classical equations to fractional order and alignment with physical interpretations, make it a preferred choice for solving fractional differential equations in mathematical and practical applications. The mathematical expression for the Caputo fractional derivative of order α is as follows: (1.4) $$D^{\alpha }y(z)=\frac{1}{\Gamma (m-\alpha )}\underset{0}{\overset{z}{\int }}{(z-\tau )}^{m-\alpha -1}{y}^m(\tau )d\tau \hspace{1em}\mathit{\text{where}}\hspace{1em}m-1<\alpha <m.$$(1.4)

In this paper, we solve (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) by two collocation methods using subdivision collocation method and Bernstein collocation method. The Subdivision Collocation Method (SCM) enhances image quality in computer graphics and plays a vital role in simulating wave propagation and solving differential equations for fluid dynamics applications (Ejaz, Baleanu, Mustafa, Malik, & Chu, 2020; Ejaz & Mustafa, 2016; Mustafa, Abbas, Ejaz, Ismail, & Khan, 2017; Mustafa, Ejaz, Kouser, Ali, & Aslam, 2021). However, until now, subdivision schemes have not been employed to find numerical solutions for fractional differential equations.

In this paper, we use the collocation method with subdivision schemes as basis functions to solve Equation (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) . On the other hand, Bernstein polynomials have proven to be successful tools for solutions in various fields. Researchers, including (Khan, Mustafa, Omar, & Komal, 2017) and (Abdul Karim, Khan, & Basit, 2022), have used Bernstein polynomials to solve different integral equations, while (Jasim & Ibraheem, 2023) applied them to fractional Bagley-Torvik equations. We modify this method for general fractional differential equations defined in (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) by adjusting the node points, introducing our new technique: the Bernstein Collocation Method. In light of the literature mentioned above, it is evident that the fractional differential Equation (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) have not been addressed using these collocation methods.

1.1. Methodology

Here, we define two methodologies that have been employed to obtain the numerical solution of the second-order fractional differential equation:

1.1.1. Step by step methodology using Bernstein polynomials

• Consider the mth degree Bernstein polynomials.

• Calculating the integer and fractional order derivatives of Bernstein polynomials.

• Begin by constructing the inconsistent system of equations for the second-order fractional differential equation using the Bernstein polynomials as basis functions of the chosen Bernstein polynomials. Subsequently, refine the system by incorporating initial or boundary conditions.

• The solution to the refined system of equations can be obtained using any numerical method. Specifically, we employed the Gaussian elimination method to determine the solution.

• The solution to the system ultimately serves as the solution to the second-order fractional differential equation.

1.1.2. Step by step methodology using subdivision scheme

• Identify a suitable linear-binary interpolating subdivision scheme that ensures the generation of at least a C ²-continuous curve. In simpler terms, the basis functions of the scheme should be at least twice continuously differentiable. In addition, the scheme itself should have at least an approximation order of six.

• Approximating the fractional order derivatives, Caputo’s derivative has been utilized, where the integral term has been approximated using the midpoint rule.

• Begin by constructing the inconsistent system of equations for the second-order fractional differential equation using the basis functions of the chosen scheme. Subsequently, refine the system by incorporating initial or boundary conditions and applying the extrapolation method utilizing a quintic polynomial.

• The solution to the refined system of equations can be obtained using any numerical method. Specifically, we employed the Gaussian elimination method to determine the solution.

• The solution to the system ultimately serves as the solution to the second-order fractional differential equation.

This research work is organized as follows: Section 2, discusses some of the key characteristics of Bernstein polynomials and the binary subdivision scheme. Section 3, develops the Bernstein collocation method based on Bernstein polynomials for the solution of a specific problem (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) . Section 4, introduces a method based on subdivision collocation for the solution of specific problem (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) . In Section 5, we present error estimations for both methods. Section 6 provides numerical findings based on the suggested algorithm and includes comparisons with other existing approaches. Finally, Section 7 presents conclusions, limitations, and future work.

2. Preliminaries

In this section, a brief introduction to Bernstein polynomials and subdivision schemes is provided below.

2.1. Bernstein polynomials

The Bernstein polynomials of degree m are defined as (2.1) $${\mathcal{B}}_{l,m}(s)=\left(\begin{array}{c}m\\ l\end{array}\right)s^l{(1-s)}^{m-l},l=0,1,2,\dots ,m \text{with} \left(\begin{array}{c}m\\ l\end{array}\right)=\frac{m!}{l!(m-l)!}.$$(2.1)

For mth degree, we have m + 1 Bernstein polynomials and for simplicity it was supposed that $B_{l,m}=0\mathit{\text{\hspace{0.5em}if\hspace{0.5em}l}}<0\mathit{\text{\hspace{0.5em}or\hspace{0.5em}l}}>m.$

2.2. Subdivision scheme and derivatives

Consider a 6-point binary interpolatory subdivision scheme as presented by Qu & Agarwal (1996) (2.2) $$\begin{array}{l}{P}_{2i}^{l+1}=P_i^l\\ P_{2i+1}^{l+1}=\left(\frac{9}{16}+2k\right)(P_i^l+P_{i+1}^l)\\ -\left(\frac{1}{16}+3k\right)(P_{i-1}^l+P_{i+2}^l)+k(P_{i-2}^l+P_{i+3}^l)\end{array}$$(2.2) where k is the tension parameter and this scheme produces C² interpolatory curve for $0<k\le \frac{3}{256}.$ The scheme defined in (2.2)(2.2) $$\begin{array}{l}{P}_{2i}^{l+1}=P_i^l\\ P_{2i+1}^{l+1}=\left(\frac{9}{16}+2k\right)(P_i^l+P_{i+1}^l)\\ -\left(\frac{1}{16}+3k\right)(P_{i-1}^l+P_{i+2}^l)+k(P_{i-2}^l+P_{i+3}^l)\end{array}$$(2.2) satisfies the following two scale relation $$\begin{array}{l}\Psi (z)=\Psi (2z)+\frac{1}{256}\{150(\Psi (2z-1)\\ +\psi (2z+1))-25(\Psi (2z-3)+\psi (2z+3))\\ +3(\Psi (2z-5)+\Psi (2z+5))\},\end{array}$$

Which satisfies $$\Psi (x)=\{\begin{array}{c}1,\hspace{1em}\mathit{\text{if}}z=0,\hfill \\ 0,\hspace{1em}\mathit{\text{if}}z\ne 0.\hfill \end{array}$$

The first and second order derivatives of the scheme (2.2)(2.2) $$\begin{array}{l}{P}_{2i}^{l+1}=P_i^l\\ P_{2i+1}^{l+1}=\left(\frac{9}{16}+2k\right)(P_i^l+P_{i+1}^l)\\ -\left(\frac{1}{16}+3k\right)(P_{i-1}^l+P_{i+2}^l)+k(P_{i-2}^l+P_{i+3}^l)\end{array}$$(2.2) was presented by Qu & Agarwal (1996), given below in Table 1.

Table 1. First and second order derivatives of the scheme (2.2)(2.2) $$\begin{array}{l}{P}_{2i}^{l+1}=P_i^l\\ P_{2i+1}^{l+1}=\left(\frac{9}{16}+2k\right)(P_i^l+P_{i+1}^l)\\ -\left(\frac{1}{16}+3k\right)(P_{i-1}^l+P_{i+2}^l)+k(P_{i-2}^l+P_{i+3}^l)\end{array}$$(2.2) .

i	0	±1	±2	±3	±4
$\Psi \prime (i)$	0	$\mp \frac{272}{365}$	$\pm \frac{53}{365}$	$\mp \frac{16}{1095}$	$\mp \frac{1}{2920}$
$\Psi ″(i)$	$-\frac{295}{56}$	$\frac{356}{105}$	$-\frac{92}{105}$	$\frac{4}{35}$	$\frac{3}{560}$

3. Bernstein collocation method

The unknown function y(z), as it appears in (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) , can be approximated by v(z) using a linear combination with Bernstein polynomials as the basis, defined as follows. (3.1) $$V(z)=\sum _{i=0}^n{v}_i{\mathcal{B}}_{i,n}(z),\hspace{1em}zϵ[0,1]$$(3.1) where v_i are unknown to be determined. $B_{i,n}(z)$ is the Bernstein polynomial as defined in (2.1)(2.1) $${\mathcal{B}}_{l,m}(s)=\left(\begin{array}{c}m\\ l\end{array}\right)s^l{(1-s)}^{m-l},l=0,1,2,\dots ,m \text{with} \left(\begin{array}{c}m\\ l\end{array}\right)=\frac{m!}{l!(m-l)!}.$$(2.1) . Equation (3.1)(3.1) $$V(z)=\sum _{i=0}^n{v}_i{\mathcal{B}}_{i,n}(z),\hspace{1em}zϵ[0,1]$$(3.1) can be expressed as (3.2) $$V(z)={(\mathcal{B}(z))}^T\mathbb{V},$$(3.2) where (3.3) $$\mathbb{V}={\left[\begin{array}{cccc}{v}_0& v_1& ,\dots ,& v_n\end{array}\right]}^T,\text{and} \mathcal{B}{(z)}^T=\left[\begin{array}{ccccc}{\mathcal{B}}_{0,n}(z)& {\mathcal{B}}_{1,n}(z)& {\mathcal{B}}_{2,n}(z)& \dots & {\mathcal{B}}_{n,n}(z)\end{array}\right].$$(3.3)

Theorem 3.1.

The vector B(z) defined in (3.3)(3.3) $$\mathbb{V}={\left[\begin{array}{cccc}{v}_0& v_1& ,\dots ,& v_n\end{array}\right]}^T,\text{and} \mathcal{B}{(z)}^T=\left[\begin{array}{ccccc}{\mathcal{B}}_{0,n}(z)& {\mathcal{B}}_{1,n}(z)& {\mathcal{B}}_{2,n}(z)& \dots & {\mathcal{B}}_{n,n}(z)\end{array}\right].$$(3.3) has an equivalent form in terms of an invertible matrix M of size $(N+1)\times (N+1)$ and vector X(z) as follows(3.4) $$\mathcal{B}(z)=\mathit{\text{MX}}(z),$$(3.4) with(3.5) $$\mathcal{B}(z)=\left[\begin{array}{c}{\mathcal{B}}_{0,n}(z)\hfill \\ {\mathcal{B}}_{1,n}(z)\hfill \\ {\mathcal{B}}_{2,n}(z)\hfill \\ ⋮\hfill \\ {\mathcal{B}}_{n,n}(z)\hfill \end{array}\right],X(z)={\left[\begin{array}{cccccc}1& z& z^2& z^3& \dots & z^n\end{array}\right]}^T,$$(3.5) and(3.6) $$M=\left[\begin{array}{ccccc}{(-1)}^0\left(\begin{array}{c}n\\ 0\end{array}\right)& {(-1)}^1\left(\begin{array}{c}n\\ 1\end{array}\right)& {(-1)}^2\left(\begin{array}{c}n\\ 2\end{array}\right)& \dots & {(-1)}^n\left(\begin{array}{c}n\\ n\end{array}\right)\\ 0& {(-1)}^0\left(\begin{array}{c}n\\ 1\end{array}\right)\left(\begin{array}{c}1\\ 1\end{array}\right)& {(-1)}^1\left(\begin{array}{c}n\\ 2\end{array}\right)\left(\begin{array}{c}2\\ 1\end{array}\right)& \dots & {(-1)}^{n-1}\left(\begin{array}{c}n\\ n\end{array}\right)\left(\begin{array}{c}n\\ 1\end{array}\right)\\ 0& 0& {(-1)}^0\left(\begin{array}{c}n\\ 2\end{array}\right)\left(\begin{array}{c}2\\ 2\end{array}\right)& \dots & {(-1)}^{n-2}\left(\begin{array}{c}n\\ n\end{array}\right)\left(\begin{array}{c}n\\ 2\end{array}\right)\\ ⋮& ⋮& ⋮& \dots & ⋮\\ 0& 0& 0& \dots & {(-1)}^{n-n}\left(\begin{array}{c}n\\ n\end{array}\right)\left(\begin{array}{c}n\\ n\end{array}\right)\end{array}\right].$$(3.6)

3.1. Integer order derivatives of Bernstein polynomials

First and second derivative of $\mathcal{B}(x)$ are calculated by using (3.7) $$\frac{d^k}{d{x}^k}(\mathcal{B}(x))=M{\phi }_{k}X(x),\hspace{1em}k=1,2$$(3.7) where k is the order of the derivative, M & X are defined in (3.4)(3.4) $$\mathcal{B}(z)=\mathit{\text{MX}}(z),$$(3.4) , and ${\phi }_1$ & ${\phi }_2$ are defined in (3.8)(3.8) $${\phi }_1=\left[\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ 1& 0& 0& \dots & 0& 0\\ 0& 2& 0& \dots & 0& 0\\ 0& 0& 3& \dots & 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & n& 0\end{array}\right]\text{and} {\phi }_2=\left[\begin{array}{ccccccc}0& 0& 0& \dots & 0& 0& 0\\ 0& 0& 0& \dots & 0& 0& 0\\ (1)(2)& 0& 0& \dots & 0& 0& 0\\ 0& (2)(3)& 0& \dots & 0& 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & (n)(n-1)& 0& 0\end{array}\right]$$(3.8) . (3.8) $${\phi }_1=\left[\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ 1& 0& 0& \dots & 0& 0\\ 0& 2& 0& \dots & 0& 0\\ 0& 0& 3& \dots & 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮\\ 0& 0& 0& \dots & n& 0\end{array}\right]\text{and} {\phi }_2=\left[\begin{array}{ccccccc}0& 0& 0& \dots & 0& 0& 0\\ 0& 0& 0& \dots & 0& 0& 0\\ (1)(2)& 0& 0& \dots & 0& 0& 0\\ 0& (2)(3)& 0& \dots & 0& 0& 0\\ ⋮& ⋮& ⋮& \dots & ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & (n)(n-1)& 0& 0\end{array}\right]$$(3.8)

3.2. Fractional order derivative of Bernstein polynomials

In this section, we construct the matrix of the fractional derivative within the framework of the Caputo fractional derivative. As demonstrated by Jasim & Ibraheem (2023), the α-order fractional derivative of a polynomial basis is (3.9) $$D^{\alpha }{z}^k=\{\begin{array}{c}\frac{\Gamma (k+1)}{\Gamma (k+1-\alpha )}{z}^{k-\alpha },\hspace{1em}\mathrm{kϵN},k>\alpha \hfill \\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{kϵN},k<\alpha .\hfill \end{array}$$(3.9)

The α₁, α₂ derivative of (3.4)(3.4) $$\mathcal{B}(z)=\mathit{\text{MX}}(z),$$(3.4) in the consequence of (3.9)(3.9) $$D^{\alpha }{z}^k=\{\begin{array}{c}\frac{\Gamma (k+1)}{\Gamma (k+1-\alpha )}{z}^{k-\alpha },\hspace{1em}\mathrm{kϵN},k>\alpha \hfill \\ 0,\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{kϵN},k<\alpha .\hfill \end{array}$$(3.9) are (3.10) $$D^{{\alpha }_k}(\mathcal{B}(z))=M{\mu }_k(z),k=1,2$$(3.10) where M is defined in (3.4)(3.4) $$\mathcal{B}(z)=\mathit{\text{MX}}(z),$$(3.4) , and $${\mu }_1(z)=\left[\begin{array}{c}0\hfill \\ \frac{\Gamma (2)}{\Gamma (2-{\alpha }_1)}{z}^{1-{\alpha }_1}\hfill \\ \frac{\Gamma (3)}{\Gamma (3-{\alpha }_1)}{z}^{2-{\alpha }_1}\hfill \\ ⋮\hfill \\ \frac{\Gamma (n+1)}{\Gamma (n+1-{\alpha }_1)}{z}^{n-{\alpha }_1}\hfill \end{array}\right]\text{and} {\mu }_2(z)=\left[\begin{array}{c}0\hfill \\ 0\hfill \\ \frac{\Gamma (3)}{\Gamma (3-{\alpha }_2)}{z}^{2-{\alpha }_2}\hfill \\ ⋮\hfill \\ \frac{\Gamma (n+1)}{\Gamma (n+1-{\alpha }_2)}{z}^{n-{\alpha }_2}\hfill \end{array}\right],$$

3.3. Numerical approximation

Now using (3.1)(3.1) $$V(z)=\sum _{i=0}^n{v}_i{\mathcal{B}}_{i,n}(z),\hspace{1em}zϵ[0,1]$$(3.1) , (3.7)(3.7) $$\frac{d^k}{d{x}^k}(\mathcal{B}(x))=M{\phi }_{k}X(x),\hspace{1em}k=1,2$$(3.7) , (3.10)(3.10) $$D^{{\alpha }_k}(\mathcal{B}(z))=M{\mu }_k(z),k=1,2$$(3.10) in (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) , we get the following system (3.11) $$\mathbb{A}{(M{\phi }_{2}X(z))}^T\mathbb{V}+\mathbb{B}{(M{\mu }_2(z))}^T\mathbb{V}+\mathbb{C}{(M{\phi }_1(z))}^T\mathbb{V}+\mathbb{E}{(M{\mu }_1(z))}^T\mathbb{V}+\mathbb{F}{(B(z))}^T\mathbb{V}=f(z),$$(3.11) this implies (3.12) $$[\mathbb{A}{(M{\varphi }_{2}X(z_i))}^T+\mathbb{B}{(M{\mu }_2(z_i))}^T+\mathbb{C}{(M{\phi }_1(z_i))}^T+\mathbb{E}{(M{\mu }_1(z_i))}^T+\mathbb{F}{(B(z_i))}^T]\mathbb{V}=f(z_i),$$(3.12) with the collocation nodes are defined as $z_j=\frac{j}{n}+ϵ,\hspace{1em}j=1,2,\dots ,n-1,$ where ϵ is very small positive number. Thus, from (3.12)(3.12) $$[\mathbb{A}{(M{\varphi }_{2}X(z_i))}^T+\mathbb{B}{(M{\mu }_2(z_i))}^T+\mathbb{C}{(M{\phi }_1(z_i))}^T+\mathbb{E}{(M{\mu }_1(z_i))}^T+\mathbb{F}{(B(z_i))}^T]\mathbb{V}=f(z_i),$$(3.12) , system of $(n-1)$ equations with $(n+1)$ unknown is obtained and given in (3.13)(3.13) $$J_1(z)\mathbb{V}=F_1(z),$$(3.13) . (3.13) $$J_1(z)\mathbb{V}=F_1(z),$$(3.13) where (3.14) $$J_1(z_i)=\mathbb{A}{(M{\varphi }_{2}X(z_i))}^T+\mathbb{B}{(M{\mu }_2(z_i))}^T+\mathbb{C}{(M{\phi }_1(z_i))}^T+\mathbb{E}{(M{\mu }_1(z_i))}^T+\mathbb{F}(B(z_i),$$(3.14) and (3.15) $$F_1={(f(z_1),f(z_2),\dots ,f(z_{n-1}))}^T.$$(3.15)

To ensure a unique solution of the system (3.12)(3.12) $$[\mathbb{A}{(M{\varphi }_{2}X(z_i))}^T+\mathbb{B}{(M{\mu }_2(z_i))}^T+\mathbb{C}{(M{\phi }_1(z_i))}^T+\mathbb{E}{(M{\mu }_1(z_i))}^T+\mathbb{F}{(B(z_i))}^T]\mathbb{V}=f(z_i),$$(3.12) , we incorporate initial or boundary conditions to obtain a consistent system of linear equations. Thus, the requried consistant system of equations are defined in (3.16)(3.16) $$J\mathbb{V}=F,$$(3.16) . (3.16) $$J\mathbb{V}=F,$$(3.16) with J and F are defined as $$J={({(y(0))}^T,{(J_1(Z))}^T,{(y\prime (0))}^T)}^T\text{\hspace{0.5em}\&\hspace{0.5em}}F={(a_1,{(F_1(z))}^T,a_2)}^T,$$ if (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) is given with the conditions (1.2)(1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) or $$J={({(y(0))}^T,{(J_1(Z))}^T,{(y(1))}^T)}^T\text{\hspace{0.5em}\&\hspace{0.5em}}F={(a_1,{(F_1(z))}^T,a_3)}^T,$$ if (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) is given with the conditions (1.3)(1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) , where $J_1(z),F_1(z),$ y(0), $y\prime (0)$ and y(1) are defined in (3.14)(3.14) $$J_1(z_i)=\mathbb{A}{(M{\varphi }_{2}X(z_i))}^T+\mathbb{B}{(M{\mu }_2(z_i))}^T+\mathbb{C}{(M{\phi }_1(z_i))}^T+\mathbb{E}{(M{\mu }_1(z_i))}^T+\mathbb{F}(B(z_i),$$(3.14) , (3.15)(3.15) $$F_1={(f(z_1),f(z_2),\dots ,f(z_{n-1}))}^T.$$(3.15) , (1.2)(1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) and (1.3)(1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) respectively. Thus by solving (3.16)(3.16) $$J\mathbb{V}=F,$$(3.16) we obtained the unique solution of (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) with the conditions (1.2)(1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) or (1.3)(1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) .

3.4. Convergence analysis of bernsitein collocation method

Theorem 3.2.

(Jasim & Ibraheem,2023). Suppose the function $V:[0,1)\to R$ is n + 1 times continuously differentiable, $y\in [0,1)$ and $V_1=\mathit{\text{Span}}\{B_{0,n}(z),B_{1,n}(z),\dots ,B_{n,n}(z)\}$, if ${(B(z))}^T\mathbb{V}$ is the best approximation to y from V then the error bound is displayed as follows $\left|\right|y-{(B(z))}^T|{|}_2\le \frac{l}{(n+1)!\sqrt{(2n+3)}}$, where $k=|y^{n+1}(z),z\in [0,1).$

Proof.

(See (Jasim & Ibraheem, 2023)) □

3.5. Error estimation

Let y(z) is exact and V(z) is the approximate solution of (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) then $$\mathit{\text{AbsoluteError}}=\left|\right|y(z)-V(z)\left|\right|$$

4. Subdivision collocation method

Here we use subdivision schemes to formulate this method, consider N be a positive integer $(N\ge 4),h=\frac{1}{N}$ and z_i = ih for $i=-4,-2,0,\dots ,N+4$ Let (4.1) $$u(z)=\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)\hspace{1em},0\le z\le 1$$(4.1) be the approximate solution of (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) . The first and second derivative of (4.1)(4.1) $$u(z)=\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)\hspace{1em},0\le z\le 1$$(4.1) are (4.2) $$u\prime (z)=\frac{1}{h}\sum _{i=-4}^{N+4}{c}_i\Psi \prime \left(\frac{z-z_i}{h}\right)\text{and} u″(z)=\frac{1}{h^2}\sum _{i=-4}^{N+4}{c}_i\Psi ″\left(\frac{z-z_i}{h}\right).$$(4.2)

Midpoint rule of integration is defined as: (4.3) $$\underset{c}{\overset{d}{\int }}f(z)\mathit{\text{dz}}=(d-c)f\left(\frac{d+c}{2}\right)$$(4.3) where d, c are upper and lower limits of integration.

When we approximate the integral appeared in (1.4)(1.4) $$D^{\alpha }y(z)=\frac{1}{\Gamma (m-\alpha )}\underset{0}{\overset{z}{\int }}{(z-\tau )}^{m-\alpha -1}{y}^m(\tau )d\tau \hspace{1em}\mathit{\text{where}}\hspace{1em}m-1<\alpha <m.$$(1.4) by using midpoint rule of integration, the caputo fractional derivative takes the following form (4.4) $$D^{\alpha }y(z)=\frac{z}{\Gamma (m-\alpha )}{\left(\frac{z}{2}\right)}^{m-\alpha -1}{y}^m\left(\frac{z}{2}\right)\hspace{1em}\mathit{\text{where}}\hspace{1em}m-1<\alpha <m$$(4.4)

Putting (4.1)(4.1) $$u(z)=\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)\hspace{1em},0\le z\le 1$$(4.1) , (4.2)(4.2) $$u\prime (z)=\frac{1}{h}\sum _{i=-4}^{N+4}{c}_i\Psi \prime \left(\frac{z-z_i}{h}\right)\text{and} u″(z)=\frac{1}{h^2}\sum _{i=-4}^{N+4}{c}_i\Psi ″\left(\frac{z-z_i}{h}\right).$$(4.2) and (4.4)(4.4) $$D^{\alpha }y(z)=\frac{z}{\Gamma (m-\alpha )}{\left(\frac{z}{2}\right)}^{m-\alpha -1}{y}^m\left(\frac{z}{2}\right)\hspace{1em}\mathit{\text{where}}\hspace{1em}m-1<\alpha <m$$(4.4) in (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) , we get $$\begin{array}{l}\frac{\mathbb{A}}{h^2}\sum _{i=-4}^{N+4}{c}_i\Psi ″\left(\frac{z-z_i}{h}\right)\\ +\frac{\mathbb{B}}{\Gamma (2-{\alpha }_2)}{z}^{2-{\alpha }_2}\frac{2^{{\alpha }_2-1}}{h^2}\sum _{i=-4}^{N+4}{c}_i\Psi ″\left(\frac{\frac{z}{2}-z_i}{h}\right)\\ +\frac{\mathbb{C}}{h}\sum _{i=-4}^{N+4}{c}_i\Psi \prime \left(\frac{z-z_i}{h}\right)\\ +\frac{\mathbb{E}}{\Gamma (1-{\alpha }_1)}{z}^{1-{\alpha }_1}\frac{2^{{\alpha }_1}}{h}\sum _{i=-4}^{N+4}{c}_i\Psi \prime \left(\frac{\frac{z}{2}-z_i}{h}\right)\\ +\mathbb{F}\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)=f(z).\end{array}$$

Taking summation common and putting $z=z_j=\mathit{\text{jh}}$ and $z_i=\mathit{\text{ih}},$ we get (4.5) $$\begin{array}{l}\sum _{i=-4}^{N+4}{c}_i[\frac{\mathbb{A}}{h^2}\Psi ″(j-i)+\frac{\mathbb{B}}{\Gamma (2-{\alpha }_2)}{(\mathit{\text{jh}})}^{2-{\alpha }_2}\frac{2^{{\alpha }_2-1}}{h^2}\Psi ″\left(\frac{j}{2}-i\right)\\ +\frac{\mathbb{C}}{h}\Psi \prime (j-i)+\frac{\mathbb{E}}{\Gamma (1-{\alpha }_1)}{(\mathit{\text{jh}})}^{1-{\alpha }_1}\frac{2^{{\alpha }_1}}{h}\Psi \prime \left(\frac{j}{2}-i\right)\\ +\mathbb{F}\Psi (j-i)]=f(\mathit{\text{jh}}),j=0,1,\dots ,N.\end{array}$$(4.5)

Making substitution of coefficients of $\Psi ,\Psi \prime ,\Psi ″,$ such as $${\mathbb{A}}_1=\frac{\mathbb{A}}{h^2},\hspace{1em}{\mathbb{B}}_j=\frac{\mathbb{B}}{\Gamma (2-{\alpha }_2)}{2}^{{\alpha }_2-1}{(\mathit{\text{jh}})}^{2-{\alpha }_2}\frac{1}{h^2},{\mathbb{C}}_1=\frac{\mathbb{C}}{h},\hspace{1em}{\mathbb{E}}_j=\frac{\mathbb{E}{2}^{{\alpha }_1}}{\Gamma (2-{\alpha }_2)}\frac{{(\mathit{\text{jh}})}^{1-\alpha 1}}{h},$$ for $j=0,1,2,\dots ,N$ and putting all in (4.5)(4.5) $$\begin{array}{l}\sum _{i=-4}^{N+4}{c}_i[\frac{\mathbb{A}}{h^2}\Psi ″(j-i)+\frac{\mathbb{B}}{\Gamma (2-{\alpha }_2)}{(\mathit{\text{jh}})}^{2-{\alpha }_2}\frac{2^{{\alpha }_2-1}}{h^2}\Psi ″\left(\frac{j}{2}-i\right)\\ +\frac{\mathbb{C}}{h}\Psi \prime (j-i)+\frac{\mathbb{E}}{\Gamma (1-{\alpha }_1)}{(\mathit{\text{jh}})}^{1-{\alpha }_1}\frac{2^{{\alpha }_1}}{h}\Psi \prime \left(\frac{j}{2}-i\right)\\ +\mathbb{F}\Psi (j-i)]=f(\mathit{\text{jh}}),j=0,1,\dots ,N.\end{array}$$(4.5) , we get (4.6) $$\begin{array}{c}\sum _{i=-4}^{N+4}{c}_i[{\mathbb{A}}_1\Psi ″(j-i)+{\mathbb{B}}_j\Psi ″\left(\frac{j}{2}-i\right)+{\mathbb{C}}_1\Psi \prime (j-i)\\ +{\mathbb{E}}_j\Psi \left(\frac{i}{2}-i\right)+\mathbb{F}\Psi (j-i)]=f(z_j),\end{array}$$(4.6) for $j=0,1,2,\dots ,N,$ this implies (4.7) $$\sum _{i=-4}^{N+4}{c}_i[{\mathbb{G}}_{j,i}+{\mathcal{R}}_{j,i}]=f(\mathit{\text{jh}}),$$(4.7) for $j=0,1,2,\dots ,N,$ where (4.8) $${\mathbb{G}}_{j,i}={\mathbb{A}}_1\Psi ″(j-i)+{\mathbb{C}}_1\Psi \prime (j-i)+\mathbb{F}\Psi (j-i),$$(4.8) and (4.9) $$R_{j,i}={\mathbb{B}}_j\Psi ″\left(\frac{j}{2}-i\right)+{\mathbb{E}}_j\Psi \prime \left(\frac{j}{2}-i\right).$$(4.9)

Equation (4.7)(4.7) $$\sum _{i=-4}^{N+4}{c}_i[{\mathbb{G}}_{j,i}+{\mathcal{R}}_{j,i}]=f(\mathit{\text{jh}}),$$(4.7) can be expressed as (4.10) $$\begin{array}{l}{c}_{-4}[{\mathbb{G}}_{j,-4}+{\mathcal{R}}_{j,-4}]+c_{-3}[{\mathbb{G}}_{j,-3}+{\mathcal{R}}_{j,-3}]+c_{-2}[{\mathbb{G}}_{j,-2}+{\mathcal{R}}_{j,-2}]\\ \hspace{1em}+\dots \dots .+c_{N+2}[{\mathbb{G}}_{j,N+2}+{\mathcal{R}}_{j,N+2}]+c_{N+3}[{\mathbb{G}}_{j,N+3}+{\mathcal{R}}_{j,N+3}]\\ \hspace{1em}+c_{N+4}[{\mathbb{G}}_{j,N+4}+{\mathcal{R}}_{j,N+4}]=f(\mathit{\text{jh}}),\end{array}$$(4.10) $$j=0,1,2,\dots ,N.$$

The functions $\mathbb{G}\hspace{0.5em}\mathit{\text{and}}\hspace{0.5em}\mathcal{R}$ gives us non-zero values at some fixed ranges of i and j, which are in the theorems described below

Theorem 4.1.

For $i=-4,-3,-2,\dots ,N+4,j=0,1,2,\dots ,N$ and using (4.8)(4.8) $${\mathbb{G}}_{j,i}={\mathbb{A}}_1\Psi ″(j-i)+{\mathbb{C}}_1\Psi \prime (j-i)+\mathbb{F}\Psi (j-i),$$(4.8) $${\mathbb{G}}_{j,i}=\{\begin{array}{c}{\mathbb{A}}_1\Psi ″(0)+{\mathbb{C}}_1\Psi \prime (0)+\mathbb{F},\hspace{1em}\mathit{\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}if\hspace{0.5em}j}}=i\hfill \\ {\mathbb{A}}_1\Psi ″(j-i)+{\mathbb{C}}_1\Psi \prime (j-i),\hspace{1em}\mathit{\text{if}}-4\le j-i\le 4,j-i\ne 0\hfill \\ 0,\hspace{1em}\mathit{\text{otherwise}}.\hfill \end{array}$$

Theorem 4.2.

For $i=-4,-3,\dots ,N+3,N+4,$ $j=0,1,2,\dots ,N$ and using (4.9)(4.9) $$R_{j,i}={\mathbb{B}}_j\Psi ″\left(\frac{j}{2}-i\right)+{\mathbb{E}}_j\Psi \prime \left(\frac{j}{2}-i\right).$$(4.9) , we get: $${\mathcal{R}}_{j,i}=\{\begin{array}{c}{\mathbb{B}}_j\Psi ″\left(\frac{j}{2}-i\right)+{\mathbb{E}}_j\Psi \prime \left(\frac{j}{2}-i\right),\hspace{0.5em}-4\le \frac{j}{2}-i\le 4\hspace{0.5em}\mathit{\text{and}}\hspace{0.5em}\frac{j}{2}-\mathit{\text{i\hspace{0.5em}be\hspace{0.5em}integer}}\hfill \\ 0,\hspace{1em}\mathit{\text{otherwise}}\hfill \end{array}$$

Theorem 4.3.

On the sight of theorem (4.1) and theorem (4.2), the system (4.7)(4.7) $$\sum _{i=-4}^{N+4}{c}_i[{\mathbb{G}}_{j,i}+{\mathcal{R}}_{j,i}]=f(\mathit{\text{jh}}),$$(4.7) in terms of matrix is(4.11) $$\mathit{\text{PC}}=(P_1+P_2)C=F,$$(4.11)

where$$\begin{array}{c}\begin{array}{c}C={(c_{-4},c_{-3},\dots ,c_{N+3},c_{N+4})}^T,\hfill \\ F={(f(0),f(h),\dots ,f((N-1)h),f(\mathit{\text{Nh}})=f(1))}^T,\hfill \end{array}\\ P_1=\left(\begin{array}{cccccccccc}{\mathbb{G}}_{0,-4}& {\mathbb{G}}_{0,-3}& {\mathbb{G}}_{0,-2}& {\mathbb{G}}_{0,-1}& {\mathbb{G}}_{0,0}& {\mathbb{G}}_{0,1}& \dots & 0& 0& 0\\ 0& {\mathbb{G}}_{1,-3}& {\mathbb{G}}_{1,-2}& {\mathbb{G}}_{1,-1}& {\mathbb{G}}_{1,0}& {\mathbb{G}}_{1,1}& \dots & 0& 0& 0\\ 0& 0& {\mathbb{G}}_{2,-2}& {\mathbb{G}}_{2,-1}& {\mathbb{G}}_{2,0}& {\mathbb{G}}_{2,1}& \dots & 0& 0& 0\\ 0& 0& 0& {\mathbb{G}}_{3,-1}& {\mathbb{G}}_{3,0}& {\mathbb{G}}_{3,1}& \dots & 0& 0& 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& 0& 0& 0& \dots & {\mathbb{G}}_{N-1,N+2}& {\mathbb{G}}_{N-1,N+3}& 0\\ 0& 0& 0& 0& 0& 0& \dots & {\mathbb{G}}_{N,N+2}& {\mathbb{G}}_{N,N+3}& {\mathbb{G}}_{N,N+4}\end{array}\right),\\ P_2=\left(\begin{array}{cccccccccc}{\mathcal{R}}_{0,-4}& {\mathcal{R}}_{0,-3}& {\mathcal{R}}_{0,-2}& {\mathcal{R}}_{0,-1}& {\mathcal{R}}_{0,0}& {\mathcal{R}}_{0,1}& \dots & 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \dots & 0& 0& 0\\ 0& 0& {\mathcal{R}}_{2,-2}& {\mathcal{R}}_{2,-1}& {\mathcal{R}}_{2,0}& {\mathcal{R}}_{2,1}& \dots & 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \dots & 0& 0& 0\\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& 0& 0& 0& \dots & {\mathcal{R}}_{N,N+2}& {\mathcal{R}}_{N,N+3}& {\mathcal{R}}_{N,N+4}\end{array}\right).\end{array}$$

4.1. Approximation of derivative condition

It is very difficult to impose derivative condition as it is in matrix form. So, we first approximate it by using finite difference method. As we have approximate the unknown function with u(z). So, we approximate the derivative of u(z). For h as defined in (4.1)(4.1) $$u(z)=\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)\hspace{1em},0\le z\le 1$$(4.1) , the first order derivative satisfies the following (4.13) $$u\prime (z)=\frac{1}{h}\sum _{i=i_{\mathrm{\text{min}}}}^{i_{\mathrm{\text{max}}}}{b}_{i}c(z+\mathit{\text{ih}})+O(h^p),\hspace{1em}p>0$$(4.13) where $O(h^p)$ denotes the order of error and p is the integer. Set i_min = 0 and $i_{\mathrm{\text{max}}}=p$ for forward difference. The vector $b=(b_{i_{\mathrm{\text{min}}}},\dots .,b_{i_{\mathrm{\text{max}}}})$ is the unknowns called convolution mask for approximation. To find these unknowns, it is necessary that (4.14) $$\sum _{i_{\mathrm{\text{min}}}}^{i_{\mathrm{\text{max}}}}{i}^n{b}_i=\{\begin{array}{c}0,\hspace{1em}0\le n\le p\hspace{1em}\mathit{\text{and}}\hspace{1em}n\ne 1\hfill \\ 1,\hspace{1em}n=1\hfill \end{array}$$(4.14)

For our method we consider the order of error $O(h^6)$ i.e. p = 6 then i_min = 0 and i_max = 6. (4.14)(4.14) $$\sum _{i_{\mathrm{\text{min}}}}^{i_{\mathrm{\text{max}}}}{i}^n{b}_i=\{\begin{array}{c}0,\hspace{1em}0\le n\le p\hspace{1em}\mathit{\text{and}}\hspace{1em}n\ne 1\hfill \\ 1,\hspace{1em}n=1\hfill \end{array}$$(4.14) implies $$\sum _{i_{\mathrm{\text{min}}}}^{i_{\mathrm{\text{max}}}}{i}^n{b}_i=\{\begin{array}{c}0,\hspace{1em}0\le n\le 6\hspace{1em}\mathit{\text{and}}\hspace{1em}n\ne 1\hfill \\ 1,\hspace{1em}n=1\hfill \end{array}$$

The solution of above equations gives us the values of b_i then, putting it in (4.13)(4.13) $$u\prime (z)=\frac{1}{h}\sum _{i=i_{\mathrm{\text{min}}}}^{i_{\mathrm{\text{max}}}}{b}_{i}c(z+\mathit{\text{ih}})+O(h^p),\hspace{1em}p>0$$(4.13) , we get $$u\prime (0)=\left(\frac{N}{60}\right)\left(-147c_0+360c_1-450c_2+400c_3-225c_4+72c_5-10c_6\right)$$ this is the approximation of first derivative condition.

4.2. Essential conditions

The initial/boundary conditions along with (4.11)(4.11) $$\mathit{\text{PC}}=(P_1+P_2)C=F,$$(4.11) makes a system of order $(N+3)\times (N+9).$ Six more conditions are still needed for unique solution, that can be calculated using quintic polynomial interpolating $(x_i,c_i),0\le i\le 5.$ Precisely, we have (4.15) $$S(-z_i)=c_{-i}\hspace{1em}i=1,2$$(4.15) where (4.16) $$S(z_i)=\sum _{j=1}^6\left(\begin{array}{c}6\\ j\end{array}\right){(-1)}^{j+1}c(z_{i-j})$$(4.16)

The above polynomial along with the conditions (1.2)(1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) or (1.3)(1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) gives us the following left and right end conditions respectively. (For details see (Ejaz & Mustafa, 2016)). (4.17) $$U_{\mathit{\text{left}}}={\left(\begin{array}{cccccccccccc}0& 1& -6& 15& -20& 15& -6& 1& 0& 0& \dots & 0\\ 0& 0& 1& -6& 15& -20& 15& -6& 1& 0& \dots & 0\\ 0& 0& 0& 1& -6& 15& -20& 15& -6& 1& \dots & 0\\ 0& 0& 0& 0& -147& 360& -450& 400& -255& 72& \dots & 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& \dots & 0\end{array}\right)}_{5\times (N+9)},$$(4.17) and (4.18) $$U_{\mathit{\text{right}}}={\left(\begin{array}{cccccccccccc}0& \dots & 1& -6& 15& -20& 15& -6& 1& 0& 0& 0\\ 0& \dots & 0& 1& -6& 15& -20& 15& -6& 1& 0& 0\\ 0& \dots & 0& 0& 1& -6& 15& -20& 15& -6& 1& 0\end{array}\right)}_{3\times (N+9)},$$(4.18) or (4.19) $${\mathbb{U}}_{\mathit{\text{left}}}={\left(\begin{array}{cccccccccccc}0& 1& -6& 15& -20& 15& -6& 1& 0& 0& \dots & 0\\ 0& 0& 1& -6& 15& -20& 15& -6& 1& 0& \dots & 0\\ 0& 0& 0& 1& -6& 15& -20& 15& -6& 1& \dots & 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& \dots & 0\end{array}\right)}_{4\times (N+9)},$$(4.19) and (4.20) $${\mathbb{U}}_{\mathit{\text{right}}}={\left(\begin{array}{cccccccccccc}0& \dots & 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& \dots & 1& -6& 15& -20& 15& -6& 1& 0& 0& 0\\ 0& \dots & 0& 1& -6& 15& -20& 15& -6& 1& 0& 0\\ 0& \dots & 0& 0& 1& -6& 15& -20& 15& -6& 1& 0\end{array}\right)}_{4\times (N+9)}$$(4.20)

Thus we get a stable system defined as (4.21) $$\mathit{\text{HC}}=F,$$(4.21) where $H={(U_{\mathit{\text{left}}}^T,P^T,U_{\mathit{\text{right}}}^T)}^T$ or $H=({\mathbb{U}}_{\mathit{\text{left}}}^T,P^T,{\mathbb{U}}_{\mathit{\text{right}}}^T)$ depending on the conditions (1.2)(1.2) $$\mathit{\text{IC}}:y(0)=a_1,\hspace{1em}y\prime (0)=a_2,$$(1.2) or (1.3)(1.3) $$\mathit{\text{BC}}:y(0)=a_1,\hspace{1em}y(1)=a_3,$$(1.3) respectively, P, U_left, U_right, ${\mathbb{U}}_{\mathit{\text{left}}}$ & ${\mathbb{U}}_{\mathit{\text{right}}}$ are defined in (4.11)(4.11) $$\mathit{\text{PC}}=(P_1+P_2)C=F,$$(4.11) , (4.17)(4.17) $$U_{\mathit{\text{left}}}={\left(\begin{array}{cccccccccccc}0& 1& -6& 15& -20& 15& -6& 1& 0& 0& \dots & 0\\ 0& 0& 1& -6& 15& -20& 15& -6& 1& 0& \dots & 0\\ 0& 0& 0& 1& -6& 15& -20& 15& -6& 1& \dots & 0\\ 0& 0& 0& 0& -147& 360& -450& 400& -255& 72& \dots & 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& \dots & 0\end{array}\right)}_{5\times (N+9)},$$(4.17) , (4.18)(4.18) $$U_{\mathit{\text{right}}}={\left(\begin{array}{cccccccccccc}0& \dots & 1& -6& 15& -20& 15& -6& 1& 0& 0& 0\\ 0& \dots & 0& 1& -6& 15& -20& 15& -6& 1& 0& 0\\ 0& \dots & 0& 0& 1& -6& 15& -20& 15& -6& 1& 0\end{array}\right)}_{3\times (N+9)},$$(4.18) , (4.19)(4.19) $${\mathbb{U}}_{\mathit{\text{left}}}={\left(\begin{array}{cccccccccccc}0& 1& -6& 15& -20& 15& -6& 1& 0& 0& \dots & 0\\ 0& 0& 1& -6& 15& -20& 15& -6& 1& 0& \dots & 0\\ 0& 0& 0& 1& -6& 15& -20& 15& -6& 1& \dots & 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& \dots & 0\end{array}\right)}_{4\times (N+9)},$$(4.19) & (4.20)(4.20) $${\mathbb{U}}_{\mathit{\text{right}}}={\left(\begin{array}{cccccccccccc}0& \dots & 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& \dots & 1& -6& 15& -20& 15& -6& 1& 0& 0& 0\\ 0& \dots & 0& 1& -6& 15& -20& 15& -6& 1& 0& 0\\ 0& \dots & 0& 0& 1& -6& 15& -20& 15& -6& 1& 0\end{array}\right)}_{4\times (N+9)}$$(4.20) respectively and the vectors C & F are defined in (4.12). Hence unique solution of (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) is obtained by solving (4.21)(4.21) $$\mathit{\text{HC}}=F,$$(4.21) .

4.3. Convergence analysis of subdivision collocation method

Theorem 4.4.

(Qu & Agarwal,1996). Suppose the function y(z) is 6 times continuously differentiable, $y\in C^6[0,1)$ and $\{u(z_i)\}$ is the best approximation to y obtained from (4.21)(4.21) $$\mathit{\text{HC}}=F,$$(4.21) with the sixth order boundary treatment at the end points then $\left|\right|y^l-u^l(z)|{|}_{\infty }=O(h^{6-j}),j=0,1,2.$

Proof.

(See (Qu & Agarwal, 1996)) □

4.4. Error estimation

Let y(z) is exact and u(z) is the approximate solution of (1.1)(1.1) $$\mathbb{A}y″(z)+\mathbb{B}{D}^{{\alpha }_2}y(z)+\mathbb{C}y\prime (z)+\mathbb{E}{D}^{{\alpha }_1}y(z)+\mathbb{F}y(z)=f(z),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}w-1<{\alpha }_w(w=1,2)<w$$(1.1) then $$\mathit{\text{AbsoluteError}}=\left|\right|y(z)-u(z)\left|\right|$$

5. Numerical implementation of proposed methods and discussions

In this section, we evaluate the accuracy of the proposed methods by solving three test problems and comparing the results with existing methods from the literature. In Example 1, numerical results are presented in Tables 2 and 3, indicating a significant agreement between our computed solutions and the exact solutions. Similarly for Example 2, we compare exact solutions with both the approximate solutions and the results from Momani & Odibat (2007) in Tables 4 and 5, showcasing the superiority of our proposed method. Finally, for Example 3, numerical results are compared with those from Zahra & Elkholy (2013) in Tables 6 and 7, highlighting the accuracy of our proposed approach.

Table 2. Exact solution, numerical solutions u(z) and V(z) of Example 5.1 for $N=n=8.$

z	Exact Solutions y(z)	Numerical Solutions V(z)	Numerical Solutions u(z)
		Cpu time = 0.078s	Cpu time = 0.080s
0	1	1.000000000000000	0.999999999999957
0.125	1.125	1.125000000000000	1.12499999999996
0.25	1.25	1.250000000000000	1.24999999999996
0.375	1.375	1.375000000000000	1.37499999999997
0.5	1.5	1.500000000000002	1.49999999999997
0.625	1.625	1.625000000000004	1.62499999999998
0.75	1.75	1.750000000000004	1.74999999999999
0.875	1.875	1.875000000000009	1.87499999999999

Table 3. Absolute errors of Example 5.1 for $N=n=8.$

z	$|y(z)-V(z)|$	$|y(z)-u(z)|$
0	0	0
0.125	4.4408e-16	4.3299e-14
0.25	2.2204ew-16	3.97459e-14
0.375	4.4408e-16	3.30846e-14
0.5	2.4424e-15	2.59792e-14
0.625	4.2188e-15	1.88737e-14
0.75	4.440e-15	1.24344e-14
0.875	9.325e-15	7.99360e-15

Table 4. Exact solution, numerical solution u(z), V(z) and (Momani & Odibat, 2007) of Example 5.2 for $N=n=10.$

z		Numerical	Numerical
		Solution V(z)	Solution u(z)
	Exact Solution y	Cpu time = 0.036s	Cpu time = 0.0785s	by Momani & Odibat (2007)
0	0.0000	0.0000	0.0000	0.0000
0.1	0.039750	0.039744	0.033661	0.039874
0.2	0.157036	0.157031	0.143874	0.158512
0.3	0.347370	0.347366	0.331165	0.353625
0.4	0.604695	0.604692	0.5890570	0.622083
0.5	0.921768	0.921765	0.910457	0.960047
0.6	1.290457	1.290454	1.285579	1.363093
0.7	1.702008	1.702007	1.705004	1.826257
0.8	2.147287	2.147286	2.156361	2.344224
0.9	2.617001	2.617001	2.628194	2.911278

Table 5. Absolute errors of Example 5.2 for $N=n=8.$

z	$|y(x)-V(z)|$	$|y(x)-u(z)|$	by Momani & Odibat (2007)
0	0.000	0	0.0000
0.1	5.65e-6	6.089e-3	1.24e-4
0.2	5.185e-6	1.3162e-2	1.476e-3
0.3	4.373e-6	1.6205e-2	6.255e-3
0.4	3.468e-6	1.5638e-2	1.7388e-2
0.5	3.202e-6	1.1311e-2	3.8279e-2
0.6	2.490e-6	4.878e-3	7.2636e-2
0.7	1.280e-6	2.996e-3	1.24249e-1
0.8	5.370e-7	9.074e-3	1.96937e-1
0.9	3.207e-7	1.1193e-2	2.94277e-1

Table 6. Exact solution, numerical solution u(z), V(z) and (Zahra & Elkholy, 2013) of Example 5.3 for $N=n=8.$

z		Numerical	Numerical
	Exact	Solution V(z)	Solution u(z)
	Solution y	Cpu time = 0.046s	Cpu time = 0.077s	by Zahra & Elkholy (2013)
0	0	2.9526e-16	0	0
0.125	−0.00021362	−0.000213623046875	−0.00052683	−0.00221
0.25	−0.0029297	−0.002929687500000	−0.0036692	−0.00701
0.375	−0.01236	−0.012359619140625	−0.013414	−0.01819
0.5	−0.03125	−0.031250000000000	−0.032556	−0.03810
0.625	−0.05722	−0.057220458984375	−0.058702	−0.06403
0.75	−0.079102	−0.079101562500000	−0.08055	−0.08467
0.875	−0.073273	−0.073272705078125	−0.074281	−0.07654

Table 7. Absolute errors of Example 5.3 for $N=n=8.$

z	$|y(z)-V(z)|$	$|y(z)-u(z)|$	by Zahra & Elkholy (2013)
0	0	0	0
0.125	1.08e-18	3.1321e-4	2.00e-3
0.25	4.33e-19	7.3955e-4	4.08e-3
0.375	0	1.0541e-3	5.83e-3
0.5	0	1.3057e-3	6.85e-3
0.625	0	1.4818e-3	6.81e-3
0.75	0	1.448e-3	5.57e-3
0.875	1.38e-17	1.008e-3	3.26e-3
1	0	0	0

Example 5.1.

Consider special type of fractional differential equation known as Bagley-Torvik equation(5.1) $$\mathit{\text{Ay}}″(z)+B{D}^{3/2}y(z)+\mathit{\text{Cy}}(z)=f(z),\hspace{1em}0\le z\le 1$$(5.1) subject to initial conditions(5.2) $$y(0)=1,y\prime (0)=1$$(5.2)

The exact solution of (5.1)(5.1) $$\mathit{\text{Ay}}″(z)+B{D}^{3/2}y(z)+\mathit{\text{Cy}}(z)=f(z),\hspace{1em}0\le z\le 1$$(5.1) is(5.3) $$y(z)=1+z$$(5.3)

The Example 5.1 has been solved by Arikoglu & Ozkol (2007) for $f(z)=1+z,$ A = 1, B = 1, and C = 1, and we consider similar values. Table 2 shows the exact solution y(z) and the numerical solutions obtained by the Bernstein collocation method V(z) and the Subdivision collocation method u(z) for Example 5.1 at different node points z, while Table 5 shows their absolute errors. Figure 1 displays their graphical representation. where V(z) and u(z) are the numerical solutions of Berstein collocation method and Subdivision collocation method defined in Equations (3.1)(3.1) $$V(z)=\sum _{i=0}^n{v}_i{\mathcal{B}}_{i,n}(z),\hspace{1em}zϵ[0,1]$$(3.1) and (4.1)(4.1) $$u(z)=\sum _{i=-4}^{N+4}{c}_i\Psi \left(\frac{z-z_i}{h}\right)\hspace{1em},0\le z\le 1$$(4.1) respectively.

Figure 1. Graphical representation of example 5.1 for $N=n=8.$

Example 5.2.

Consider the following fractional differential equation known as fractional oscillation equation(5.4) $$y″(z)-a{D}^{\alpha }y(z)-\mathit{\text{by}}=8,\hspace{1em}z>0,0<\alpha <2$$(5.4) subject to initial conditions:(5.5) $$y(0)=0,y\prime (0)=0$$(5.5)

Following (Momani & Odibat, 2007), we consider $a=b=-1$ and $\alpha =0.5.$ Table 4 shows the exact solution y(z) and numerical solutions for both methods of Example (5.2) while Table 5 shows their absolute errors. And Figure 2 displays their graphical representation. It can be easily seen that both of our methods are in good agreement with exact solution.

Figure 2. Graphical representation of example 5.2 for $N=n=10.$

Example 5.3.

Consider the following fractional differential equation:(5.6) $$y″(z)+\eta D^{\alpha }y(z)+\mu y(z)=f(z)$$(5.6) with boundary conditions(5.7) $$y(0)=0,y(1)=0$$(5.7)

The exact solution of (5.6)(5.6) $$y″(z)+\eta D^{\alpha }y(z)+\mu y(z)=f(z)$$(5.6) is:(5.8) $$y(z)=z^4(z-1)$$(5.8)

Following (Zahra & Elkholy, 2013), we also take $\alpha =0.3,\eta =0.5$ and μ = 1 in (5.6)(5.6) $$y″(z)+\eta D^{\alpha }y(z)+\mu y(z)=f(z)$$(5.6) . Table 5 shows the exact solution y(z) and numerical solutions for both methods of Example 5.3 while Table 7 shows their absolute errors. And Figure 3 displays their graphical representation. It can be easily seen that both of our methods are in good agreement with exact solution.

Figure 3. Graphical representation of example 5.3 for $N=n=8.$

6. Concluding remarks

Two collocation methods, based on Bernstein polynomials and subdivision schemes, are utilized in this paper for solving fractional differential equations using the Caputo operator. The Midpoint rule of integration is used to derive the Subdivision collocation method. Both methods do not distinguish between whether the problems are initial value or boundary value. We have presented applications for both methods by considering various types of fractional differential equations. These methods give us better results, even for problems for which a closed-form exact solution is not available, as demonstrated in Example 5.2. In this article, we initially select an appropriate subdivision scheme and Bernstein polynomial for a specified interval. Next, we compute derivatives of the subdivision schemes. Then, we transform the problem into a system of linear equations using both techniques. Following that, we compute a suitable solution for the system using a numerical technique. We compare the obtained results with existing methods in the literature, demonstrating that the proposed methods yield superior outcomes. Finally, we assess the convergence rate and present the results through graphs and tables, providing detailed explanations for each.

• It is evident from Table 2 and Table 5 that both methods gives us approximation of high accuracy.

Table 4 and Table 5 demonstrates the comparison of these methods with the method described by Momani & Odibat (2007).

Table 5 and Table 7 tells us reliability of both methods than the method defined by Zahra & Elkholy (2013).

• Figures 1–3 gives us visual confirmation.

Finally, one can easily conquer with graphs and tables, Bernstein collocation method and subdivision collocation method are more encouraging.

6.1. Limitations of methods

The proposed methods are designed for fractional differential equation in a specified interval. We apply midpoint rule of integration to approximate the definite integral appears in caputo definition. This would be replaced with better approximation to enhance the accuracy of subdivision based method.

6.2. Future work

The extension of proposed methods for the solution of high order fractional differential equations on an arbitrary interval is another direction for future work.

Authors’ contributions

Conceptualization, Syeda Tehmina Ejaz; Formal analysis, Namra Zulqarnain, Jihad Younis and Saima Bibi; Methodology, Syeda Tehmina Ejaz and Namra Zulqarnain; Supervision, Syeda Tehmina Ejaz; Writing original draft, Namra Zulqarnain and Syeda Tehmina Ejaz; Writing, review and editing, Syeda Tehmina Ejaz Jihad Younis and Saima bibi.

Disclosure statement

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Data availability statement

The data used to support the findings of the study are available within this paper.