An Implicit Hybrid Method for Solving Fractional Bagley-Torvik Boundary Value Problem

Abstract

In this article, a modified implicit hybrid method for solving the fractional Bagley-Torvik boundary (BTB) value problem is investigated. This approach is of a higher order. We study the convergence, zero stability, consistency, and region of absolute stability of the modified implicit hybrid method. Three of our numerical examples are presented.

1. Introduction

Several engineering and physical phenomena are modeled mathematically using fractional derivatives. Fractional derivatives have many applications, such as diffusion problems, liquid crystals, proteins, mechanics structural control, and biosystems [1,2,3,4,5,6,7,8,9]. Several analytical and numerical methods are used to solve fractional boundary and initial value problems, such as generalized differential transform, the Adomian decomposition method, the homotopy perturbation technique, fractional multistep methods, the spline approximation method, and the collocation method [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26].

In this article, we consider the Bagley-Torvik boundary (BTB) value problem in the form:

$$y^{″}\left(x\right)=a{y}^{\left(\frac{3}{2}\right)}\left(x\right)+by\left(x\right)+r\left(x\right)=f\left(x,y,y^{\left(\frac{3}{2}\right)}\right), x\in \left[0,X\right]$$

(1)

subject to

$$y\left(0\right)=\alpha , y\left(X\right)=\beta ,$$

(2)

where $a,b,\alpha ,$ and $\beta$ are constants, and $y,f\in L_1\left[0,X\right].$

BTB was discussed analytically in [27]. Then, several numerical approaches were used to solve it such as the discrete spline method [28], the Hybridizable discontinuous Galerkin method [29], generalizing the Taylor collocation method [30], and the operational matrix of Haar wavelet method [30]. Special attention was given when $y\left(0\right)=y\left(1\right)=0$. In addition, analytical solutions for such cases are investigated using the modified spectral method and the Adomian decomposition method [31].

We use a local fractional derivative, which is presented in [32,33]. This definition of fraction derivative works efficiently with the proposed method since it has several properties such as the product rule, power rule, and chain rule. These properties are given in the next section.

We modify an implicit hybrid method to solve Equations (1) and (2). We find an explicit formula to solve such a problem. We investigate some analytical properties of the proposed method such as consistency, stability, convergence, and order of convergence.

We organize our paper as follow. In Section 2, we mention some definitions and results which we use in this article. A modified fractional implicit hybrid method and its analytical properties are presented in Section 3. Three of our numerical examples are presented in Section 4. We compare our results with the results in [19,29]. Finally, we draw some conclusions in Section 5.

2. Preliminaries

First, we mention the definition of fractional derivative which we will use.

Definition 1.

Let $u\left(x\right)\in C_{\alpha }\left(0,X\right)$, the fractional derivatives of order $1\ge \alpha >0$ at $x_0$ is defined as [33,34]:

$$D^{\alpha }u\left(x\right)=\underset{x\to x_0}{\lim }\frac{{\mathsf{\Delta }}^{\alpha }(u\left(x\right)-u\left(x_0\right)}{{(x-x_0)}^{\alpha }},$$

(3)

where ${\mathsf{\Delta }}^{\alpha }(u\left(x\right)-u\left(x_0\right)\approx \Gamma \left(\alpha +1\right)\mathsf{\Delta }\left(u\left(x\right)-u\left(x_0\right)\right).$

The power rule of this local fractional derivative is given in the following theorem.

Theorem 1.

This fractional derivative satisfies the following power rule:

$$D^{\alpha }{x}^p=\left\{\begin{array}{cc}\frac{\Gamma \left(p+1\right)}{\Gamma \left(p-\alpha +1\right)}{x}^{p-\alpha },& p\ge \alpha \\ 0& otherwise\end{array}\right\}.$$

(4)

In addition, it is easy to see that:

• $D^{\alpha }c=0$ for constant $c.$
• $D^{\alpha }\left(l\left(\circ m\right)\left(x\right)\right)={\left(\frac{dm}{dx}\right)}^{\alpha }{D}^{\alpha }l\left(m\left(x\right)\right)$.
• $D^{k\alpha }l\left(x\right)=D^{\alpha }{D}^{\alpha }\dots D^{\alpha }l\left(x\right).$
• $D^{\alpha }\left(l\left(x\right)m\left(x\right)\right)=l\left(x\right)D^{\alpha }m\left(x\right)+m\left(x\right)D^{\alpha }l\left(x\right).$

In this paper, a modified fractional implicit hybrid multistep method will be presented. To the best of our knowledge, no work has been done to discuss this problem using the implicit hybrid multistep methods.

Two major techniques are used to solve

$$w^{\left(m\right)}=f\left(x,w,w^{\prime },\dots ,w^{\left(m-1\right)}\right),w\left(a\right)=w_0,w^{\prime }\left(a\right)=w_1,\dots ,w^{\left(m-1\right)}\left(a\right)=w_{m-1},$$

(5)

which are one step methods, such as Taylor and Runge-Kutta methods, and multistep methods such as Adams-Bashforth and Adams-Moulton methods. One-step methods are suitable only for the first order since they have a very low order of accuracy. If the higher order Runge-Kutta method is used, more function evaluations per step are required. Hence, solving Equation (5) using one-step methods requires the transformation of the problem into a system of first order differential equations which makes the dimension of the problem high and its scale also high. As a result, it will be time consuming for large scale problems with a low accuracy.

On the other hand, multistep step methods do not need to transform Equation (5) into a system of first order differential equations. These methods give higher order accuracy. However, they are not efficient in terms of function evaluations as are the one step methods and require more than one value to start the integration process.

In this paper, we look for a method that is a continuous implicit hybrid one step method. This method is as efficient as the one step methods and has as high an accuracy as the multistep methods. Next, we define the $k$-step hybrid formula. Let $\left\{x_0,x_1,\dots ,x_N\right\}$ be a uniform partition of $\left[a,b\right]$ with $x_i=a+ih$, $i=0,1,\dots ,N$, and $h=\frac{b-a}{N}$.

Definition 2.

A $\mathit{k}$-step hybrid formula is defined by:

$${\sum }_{i=0}^k{a}_i{y}_{n+i}+{\sum }_{i=0}^l{a}_{n+v_i}{y}_{n+v_i}=h{\sum }_{i=0}^k{b}_i{f}_{n+i}+h{\sum }_{i=0}^l{b}_{n+v_i}{f}_{n+v_i},$$

where $a_k=1,$ $a_0$ and $b_0$ are nonzero, $v\notin \left\{0,1,\dots ,k\right\},$ $y_{n+i}=y\left(x_n+ih\right)$ and $f_{n+v_i }=f\left(x_{n+v_i},y_{n+v_i}\right)$. For more details, see in [34].

Definition 3.

Let:

$$\begin{gathered} ℒ\left[y\left[x_n\right];h\right]={\displaystyle \sum }_{i=0}^k{a}_i{y}_{n+i}+{\displaystyle \sum }_{i=0}^l{a}_{n+v_i}{y}_{n+v_i}=h{\displaystyle \sum }_{i=0}^k{b}_i{f}_{n+i}+h{\displaystyle \sum }_{i=0}^l{b}_{n+v_i}{f}_{n+v_i} \\ =c_0{y}_n+c_1{y}_n^{\prime }+\dots \end{gathered}$$

If $c_0=0,c_1=0,\dots ,c_{p+1}=0,c_{p+2}\ne 0$, then the order of the method is $p$ and the error constant is $c_{p+2}.$

Definition 4.

If the first and second characteristic polynomials are $\rho \left(z\right)={\sum }_{i=0}^k{\alpha }_i z^i$ and $\sigma \left(z\right)={\sum }_{i=0}^l{\beta }_i z^i$ with:

• ${\alpha }_0^2+{\beta }_0^2\ne 0,$
• The order is greater than or equal 1,
• ${\sum }_{i=0}^k{\alpha }_i=0,$
• $\rho \left(1\right)={\rho }^{\prime }\left(1\right)=0,$
• ${\rho }^{″}\left(1\right)=2\sigma \left(1\right)$,

then, it and its block method are called consistent.

Definition 5.

If no zeros of the first characteristic polynomial have a modulus greater than one and every root of modulus one has multiplicity not greater than one, then it is called zero stable.

Definition 6.

If the method is consistent and zero stable, it is convergent.

3. Method of Solution

To derive the modified fractional implicit hybrid method, we approximate the solution of Equation (1) by:

$$y\left(x\right)={\sum }_{k=0}^4{a}_x x^k,$$

(6)

with second derivative given by:

$$y^{″}\left(x\right)={\sum }_{k=2}^{4}k\left(k-1\right)a_x x^{k-2}.$$

(7)

Let $\left\{x_0,x_1,x_2,\dots ,x_m\right\}$ be a uniform partition of $\left[0,X\right]$ with $h=\frac{X}{m}$. Collocate Equation (7) at $x_n,x_{n+\frac{1}{2}},x_{n+1}$ and interpolate Equation (6) at $x_n,x_{n+\frac{1}{2}}$ to get:

$$\left(\begin{array}{ccc}1& x_n& \begin{array}{ccc}{x}_n^2& x_n^3& x_n^4\end{array}\\ 1& x_{n+\frac{1}{2}}& \begin{array}{ccc}{x}_{n+\frac{1}{2}}^2& x_{n+\frac{1}{2}}^3& x_{n+\frac{1}{2}}^4\end{array}\\ \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\end{array}& \begin{array}{c}\begin{array}{ccc}2& 6x_n& 12x_n^2\end{array}\\ \begin{array}{ccc}2& 6x_{n+\frac{1}{2}}& 12x_{n+\frac{1}{2}}^2\end{array}\\ \begin{array}{ccc}2& 6x_{n+1}& 12x_{n+1}^2\end{array}\end{array}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_1\\ \begin{array}{c}{a}_2\\ a_3\\ a_4\end{array}\end{array}\right)=\left(\begin{array}{c}{y}_n\\ y_{n+\frac{1}{2}}\\ \begin{array}{c}{f}_n\\ f_{n+\frac{1}{2}}\\ f_{n+1}\end{array}\end{array}\right),$$

(8)

where:

$$f_n=f\left(x_n, y_n,y_n^{\left(\frac{3}{2}\right)}\right), f_{n+\frac{1}{2}}=f\left(x_{n+\frac{1}{2}}, y_{n+\frac{1}{2}},y_n^{\left(\frac{3}{2}\right)}\right), f_{n+1}=f\left(x_{n+1}, y_{n+1},y_{n+1}^{\left(\frac{3}{2}\right)}\right).$$

Let $x_n=x-h{t}^{\frac{3}{2}}$. Then, $x_{n+\frac{1}{2}}=x-h{t}^{\frac{3}{2}}+\frac{h}{2}$ and $x_{n+1}=x-h{t}^{\frac{3}{2}}+h$. Solving Equation (8), we get:

$$y\left(x\right)={\alpha }_0\left(t\right)y_n+{\alpha }_1\left(t\right)y_{n+\frac{1}{2}}+{\beta }_0\left(t\right)f_n+{\beta }_1\left(t\right)f_{n+\frac{1}{2}}+{\beta }_2\left(t\right)f_{n+1},$$

(9)

where:

$${\alpha }_0\left(t\right)=1-2t^{\frac{3}{2}},$$

$${\alpha }_1\left(t\right)=2t^{\frac{3}{2}},$$

$${\beta }_0\left(t\right)=h^2\left(-\frac{7}{48}{t}^{\frac{3}{2}}+\frac{1}{2}{t}^3-\frac{1}{2}{t}^{\frac{9}{2}}+\frac{t^6}{6}\right),$$

$${\beta }_1\left(t\right)=h^2\left(-\frac{1}{8}{t}^{\frac{3}{2}}+\frac{2}{3}{t}^{\frac{9}{2}}-\frac{t^6}{3}\right),$$

$${\beta }_2\left(t\right)=h^2\left(\frac{1}{48}{t}^{\frac{3}{2}}-\frac{1}{6}{t}^{\frac{9}{2}}+\frac{t^6}{6}\right).$$

For $x=x_{n+1},$ $t=1$ and

$$y_{n+1}=-y_n+2y_{n+\frac{1}{2}}+\frac{h^2}{48}\left(f_{n+1}+10f_{n+\frac{1}{2}}+f_n\right).$$

(10)

Using the Taylor series about $x=x_n$ for Equation (9), we get:

$$\begin{gathered} y_{n+1}+y_n-2y_{n+\frac{1}{2}}-\frac{h^2}{48}\left(f_{n+1}+10f_{n+\frac{1}{2}}+f_n\right) \\ =-\frac{h^6}{15360}{y}^{\left(6\right)}\left(x_n\right)-\frac{h^7}{30360}{y}^{\left(7\right)}\left(x_n\right)-\cdots \end{gathered}$$

which means that the order of Equation (10) is 4 and the error constant is $\mathcal{O}\left(0.000065h^6\right)$.

The first and the second characteristic functions are given by:

$$\rho \left(z\right)=z-2z^{\frac{1}{2}}+1={\left(\sqrt{z}-1\right)}^2,$$

and

$$\sigma \left(z\right)=\frac{1}{48}\left(z+10z^{\frac{1}{2}}+1\right).$$

Simple calculation implies that:

• The roots of $\rho$ for which $\left|z\right|=1$ are simple.
• Sum of coefficients of $\rho$ is zero.
• ${\rho }^{\prime }\left(1\right)=0=\rho \left(1\right).$
• ${\rho }^{″}\left(1\right)=2!\sigma \left(1\right)=\frac{1}{2}.$

This means that Equation (10) is consistent and zero stable which means that it is convergent. To find the region of absolute stability, let:

$$g\left(z\right)=\frac{\rho \left(z\right)}{\sigma \left(z\right)}=\frac{48\left(z-2z^{\frac{1}{2}}+1\right)}{z+10z^{\frac{1}{2}}+1}.$$

Let $z=e^{i\theta }$, then:

$$g\left(\theta \right)=-\frac{96 {\sin }^2\left(\frac{\theta }{4}\right)}{5+\cos \left(\frac{\theta }{2}\right)}.$$

Thus, the interval of absolute stability is (−9.6, 0) and the region of absolute stability is given in Figure 1.

Differentiate ${\alpha }_i$ and ${\beta }_i$ to get:

$$D^{\frac{3}{2}}{\alpha }_0\left(t\right)=-\frac{3}{2}\sqrt{\pi },$$

$$D^{\frac{3}{2}}{\alpha }_1\left(t\right)=\frac{3}{2}\sqrt{\pi },$$

$$D^{\frac{3}{2}}{\beta }_0\left(t\right)=h^2\left(-\frac{7\sqrt{\pi }}{64}+\frac{4}{\sqrt{\pi }}{t}^{\frac{3}{2}}-\frac{315\sqrt{\pi }}{128}{t}^3+\frac{256}{63\sqrt{\pi }}{t}^{\frac{9}{2}}\right),$$

$$D^{\frac{3}{2}}{\beta }_1\left(t\right)=h^2\left(-\frac{3\sqrt{\pi }}{32}+\frac{105\sqrt{\pi }}{32}{t}^3-\frac{512}{63\sqrt{\pi }}{t}^{\frac{9}{2}}\right),$$

$$D^{\frac{3}{2}}{\beta }_2\left(t\right)=h^2\left(\frac{\sqrt{\pi }}{64}-\frac{105\sqrt{\pi }}{128}{t}^3+\frac{256}{63\sqrt{\pi }}{t}^{\frac{9}{2}}\right).$$

Similarly

$$D_t^{\frac{3}{2}}y\left(x\right)=\frac{\sqrt{27}}{8}{\pi }^{\frac{3}{4}}{h}^{\frac{3}{2}}{D}_t^{\frac{3}{2}}y\left(t\right),$$

Thus

$$D_t^{\frac{3}{2}}y\left(t\right)=\frac{8}{\sqrt{27}{\pi }^{\frac{3}{4}}{h}^{\frac{3}{2}}}\left(\begin{array}{c}{D}_t^{\frac{3}{2}}{\alpha }_0\left(t\right)y_n+D_t^{\frac{3}{2}}{\alpha }_1\left(t\right)y_{n+\frac{1}{2}}+D_t^{\frac{3}{2}}{\beta }_0\left(t\right)f_n\\ +D_t^{\frac{3}{2}}{\beta }_1\left(t\right)f_{n+\frac{1}{2}}+D_t^{\frac{3}{2}}{\beta }_2\left(t\right)f_{n+1}\end{array}\right).$$

(11)

Then, at $x_n,x_{n+\frac{1}{2}},x_{n+1},$ $t=0,\frac{1}{\sqrt[3]{4}},1$, which imply that:

$${\mathrm{D}}^{\frac{3}{2}}{\mathrm{y}}_{\mathrm{n}}=\frac{4}{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}\left(-y_n+y_{n+\frac{1}{2}}\right)-\frac{8\sqrt{h}}{\sqrt[4]{\pi }\sqrt{27}}\left(7f_n+6f_{n+\frac{1}{2}}-f_{n+1}\right),$$

(12)

$${\mathrm{D}}^{\frac{3}{2}}{\mathrm{y}}_{\mathrm{n}+\frac{1}{2}}=\frac{4}{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}\left(-y_n+y_{n+\frac{1}{2}}\right) +\frac{8\sqrt{h}}{\sqrt{27}\sqrt[4]{{\pi }^5}}\left(\begin{array}{c}\frac{90896-23373\pi }{32256}{f}_n\\ +\frac{-9192+5859\pi }{8064}{f}_{n+\frac{1}{2}}\\ +\frac{16384-6611\pi }{326}{f}_{n+1}\end{array}\right),$$

(13)

$${\mathrm{D}}^{\frac{3}{2}}{\mathrm{y}}_{\mathrm{n}+1}=\frac{4}{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}\left(-y_n+y_{n+\frac{1}{2}}\right) +\frac{8\sqrt{h}}{\sqrt{27}\sqrt[4]{{\pi }^5}}\left(\begin{array}{c}\frac{65024-20727\pi }{8064}{f}_n\\ -\frac{8192-3213\pi }{1008}{f}_{n+\frac{1}{2}}\\ +\frac{32768-6489\pi }{8064}{f}_{n+1}\end{array}\right),$$

(14)

From Equation (12), we get:

$$y_{n+\frac{1}{2}}=y_n -\frac{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}{4}{\mathrm{D}}^{\frac{3}{2}}{\mathrm{y}}_{\mathrm{n}}-\frac{2h^2}{3\pi }\left(7f_n+6f_{n+\frac{1}{2}}-f_{n+1}\right).$$

(15)

Let $D^{\frac{3}{2}}{y}_0=\theta$. Then, $y_1,\dots ,y_m$ are functions of $\theta$. Using the shooting method, we find the value of $\theta$.

4. Numerical Results

In this section, we present three of our examples. Comparison with References [19] and [29] will be presented.

Example 1.

Consider the following problem:

$$D^{2}y\left(t\right)+D^{\frac{3}{2}}y\left(t\right)+y\left(t\right)=2+4\sqrt{\frac{t}{\pi }}+t^2+\alpha =g\left(t\right),$$

subject to

$$y\left(0\right)=\alpha ,y\left(5\right)=25+{\alpha }^2,$$

where the exact solution is $y\left(t\right)=t^2+{\alpha }^2.$ Let $h=0.01$ and $x_i=ih, i=0, 1, \dots , 500.$ Let $D^{\frac{3}{2}}y\left(0\right)=\theta .$ Using the modified fractional implicit hybrid method, we get the following system:

$$y_{n+1}=-y_n+2y_{n+\frac{1}{2}}-\frac{{10}^{-4}}{48}\left(\begin{array}{c}{D}^{\frac{3}{2}}\left(y_{n+1}+10y_{n+\frac{1}{2}}+y_n-\right)\\ +y_{n+1}+10y_{n+\frac{1}{2}}+y_n\end{array}\right)+\frac{{10}^{-4}}{4}g\left(t_n\right),$$

$$y_{n+\frac{1}{2}}=y_n -\frac{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}} }{4}{D}^{\frac{3}{2}}{y}_n-\frac{2\ast {10}^{-4}}{3\pi }\left(\begin{array}{c}{D}^{\frac{3}{2}}\left(7y_{n+1}+6y_{n+\frac{1}{2}}-y_n\left(t\right)\right)\\ +7y_{n+1}+6y_{n+\frac{1}{2}}-y_n\end{array}\right)+\frac{24\ast {10}^{-4}}{3\pi }g\left(t_n\right),$$

$$D^{\frac{3}{2}}{y}_{n+\frac{1}{2}}=\frac{4}{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}\left(-y_n+y_{n+\frac{1}{2}}\right)-\frac{0.8}{\sqrt{27}\sqrt[4]{{\pi }^5}}\left(\begin{array}{c}\frac{90896-23373\pi }{ 32256}\left(D^{\frac{3}{2}}{y}_n+y_n-g\left(t_n\right)\right)\\ +\frac{-9192+5859\pi }{ 8064}(D^{\frac{3}{2}}{y}_{n+\frac{1}{2}}+y_{n+\frac{1}{2}}-g\left(t_n\right)\\ +\frac{16384-6611\pi }{ 326}\left(D^{\frac{3}{2}}{y}_{n+1}+y_{n+1}-g\left(t_n\right)\right)\end{array})\right),$$

$$D^{\frac{3}{2}}{y}_{n+1}=\frac{4}{\sqrt{3}\sqrt[4]{\pi }{h}^{\frac{3}{2}}}\left(-y_n+y_{n+\frac{1}{2}}\right) +\frac{0.8}{\sqrt{27}\sqrt[4]{{\pi }^5}}\left(\begin{array}{c}\frac{65024-20727\pi }{ 8064}\left(D^{\frac{3}{2}}{y}_n+y_n-g\left(t_n\right)\right)\\ -\frac{8192-3213\pi }{1008}(D^{\frac{3}{2}}{y}_{n+\frac{1}{2}}+y_{n+\frac{1}{2}}-g\left(t_n\right))\\ +\frac{32768-6489\pi }{ 8064}\left(D^{\frac{3}{2}}{y}_{n+1}+y_{n+1}-g\left(t_n\right)\right)\end{array}\right),$$

where $y_{500}=25+{\alpha }^2-0.0000790783\theta$. To find $\theta$, we set:

$$y_{500}=25+{\alpha }^2$$

to get $\theta =0.$ The effect of $\mathsf{\alpha }$ on the solution is given in

Table 1. Comparison between our results and the results in Reference [20] for Example 1.

α	$\mathit{e}\left(\mathit{\alpha }\right)$ in [19]	$\mathit{e}\left(\mathit{\alpha }\right)$ Using Our Method
0	$1.71\times {10}^{-10}$	$1.12\times {10}^{-16}$
0.001	$2.10\times {10}^{-10}$	$2.21\times {10}^{-16}$
0.01	$3.84\times {10}^{-10}$	$2.27\times {10}^{-16}$
0.1	$7.14\times {10}^{-10}$	$3.01\times {10}^{-16}$
0.5	$2.92\times {10}^{-10}$	$3.02\times {10}^{-16}$
1	$8.73\times {10}^{-10}$	$3.07\times {10}^{-16}$

where:

$$e\left(\alpha \right)=max\left\{\left|y\left(t_n\right)-y_n\right|,n=0, 1,\dots , 100\right\}.$$

We compare our results with the results in Reference [19].

From Table 1, we see that there is no significant effect for the initial condition. In addition, the proposed method gives more accurate results than Reference [19]. The approximate and exact solutions are given in Figure 2.

Example 2.

Consider the following problem:

$$D^{2}y\left(t\right)+D^{\frac{3}{2}}y\left(t\right)+y\left(t\right)=1+\frac{8}{\pi }{t}^{\frac{3}{2}}+6t+t^3=g\left(t\right)$$

subject to

$$y\left(0\right)=1,y\left(1\right)=2$$

then, the exact solution is $y\left(t\right)=t^3+1.$ Let $h=0.01$ and $x_i=ih, i=0,1,\dots ,100.$ Let $D^{\frac{3}{2}}y\left(0\right)=\theta .$ Following the procedure described in Example 1, we find that $\theta =1.2\times {10}^{-14}$.

The errors for Example 2 are given in

Table 2. The errors for Example 2.

$\mathit{t}$	$\left|\mathit{y}\left({\mathit{t}}_{\mathit{n}}\right)-{\mathit{y}}_{\mathit{n}}\right|$
0	$1.01\times {10}^{-16}$
0.1	$1.21\times {10}^{-16}$
0.2	$1.35\times {10}^{-16}$
0.3	$1.39\times {10}^{-16}$
0.4	$2.62\times {10}^{-16}$
0.5	$3.53\times {10}^{-16}$
0.6	$2.89\times {10}^{-16}$
0.7	$2.21\times {10}^{-16}$
0.8	$1.98\times {10}^{-16}$
0.9	$1.45\times {10}^{-16}$
1.0	$1.11\times {10}^{-16}$

. The approximate and exact solutions are given in Figure 3.

Example 3.

Consider the following problem:

$$D^{2}y\left(t\right)+D^{\frac{3}{2}}y\left(t\right)+y\left(t\right)=\frac{8}{\sqrt{\pi }}{t}^{\frac{3}{2}}+t^3+7t+1=g\left(t\right),$$

subject to

$$y\left(0\right)=1,y\left(1\right)=3.$$

then, the exact solution is $y\left(t\right)=t^3+t+1.$ Let $h=0.01$ and $x_i=ih, i=0, 1,\dots , 300.$ Let $D^{\frac{3}{2}}y\left(0\right)=\theta .$ Following the procedure described in Example 1, we find that $\theta =2.6\times {10}^{-15}$.

We compare our results with the results in Reference [29]. Let

$$error=\max \left\{\left|y\left(t_n\right)-y_n\right|:n=0,1,\dots ,300\right\}.$$

In

Table 3. Comparison between our results and the results in Reference [29] for Example 3.

Method	Error
Our method	$2.89\times {10}^{-16}$
Results in Reference [29]	$1.98\times {10}^{-13}$

, we present the comparison between our results and the results in Reference [29] for Example 3. The approximate and exact solutions are given in Figure 4.

5. Conclusions

In this paper, the modified fractional implicit hybrid method is presented for solving a class of fractional BTB. The proposed method is based on an implicit hybrid multistep method. We study the convergence, zero stability, consistency, and region of absolute stability. Three numerical examples are presented. We notice the following:

• The modified implicit hybrid method is consistent, zero stable, and the order of convergence is 4.
• The interval of convergence is (−9.6, 0) and the region of absolute stability is given in Figure 1.
• The order of convergence is high without the need to refer to more initial conditions.
• As seen in Table 1 and Table 3, the modified fractional implicit hybrid method gives more accurate results than other methods.
• From Table 1, Table 2 and Table 3 and Figure 2, Figure 3 and Figure 4, we see that the approximate solutions are very accurate and very close to the exact solutions.
• The modified fractional implicit hybrid method can be applied to more physical and engineering applications.