ABSTRACT

There are few numerical techniques available to solve the Bagley-Torvik equation which occurs considerably frequently in various offshoots of applied mathematics and mechanics. In this paper, we show that Chelyshkov-tau method is a very effective tool in numerically solving this equation. To show the accuracy and the efficiency of the method, several problems are implemented and the comparisons are given with other methods existing in the recent literature. The results of numerical tests confirm that Chelyshkov-tau method is superior to other existing ones and is highly accurate.

Keywords:

Chelyshkov, Tau Method, Bagley-Torvik, Caputo Derivative, Residual Functions

1. Introduction

Recently, the fractional derivative has attracted a lot of attentions due to widely applied in various fields of physics and engineering. Fractional derivative is an excellent tool to describe memory and genetic characteristics of various materials and processes. Many phenomena in various fields of science and engineering such as control, porous media, electrochemistry, viscoelasticity, electromagnetic, etc. can be described by Fractional differential Equation (FDE). The evolution of effective and perfectly appropriate method for numerically solving FDEs has received great attention over the past years.

In this paper, we develop a new approach called Chelyshkov-tau method for solving Bagley-Torvik equation of the form

$A_2{u}^{″}\left(x\right)+A_1{u}^{\left(3/2\right)}\left(x\right)+A_{0}u\left(x\right)=f\left(x\right),x\in \left[0,1\right]$ (1)

subject to the boundary conditions

${\displaystyle \sum _{j=0}^1}{\alpha }_{ij}{u}^{\left(j\right)}\left(0\right)+{\beta }_{ij}{u}^{\left(j\right)}\left(1\right)={\gamma }_i,i=0,1$ (2)

where $A_2\ne 0$ represents mass of the thin rigid plate immersed in a Newtonian fluid, $A_1$ is constant depending on area of the thin rigid, viscosity and density of fluid and $A_0$ represents stiffness of the spring. $f\left(x\right)$ is a given function. $u\left(x\right)$ represents motion of the rigid plate. The questions of existence and uniqueness of the solution to this initial value problem have been discussed in [1] [2] , so there is no need to go into these matters here.

Several numerical and analytical methods of Equation (1) were considered by many authors, such as finite difference method [3] , collocation method based on Müntz polynomials [4] , Tau approximate [5] , Adomian decomposition method [6] [7] , variational iteration method [8] , the generalized block pulse operational matrix [9] , homotopy perturbation method [10] [11] , generalized differential transform method [12] , Legendre-collocation method [13] , Laplace transforms [14] , Fourier transforms [15] , eigenvector expansion [16] , fractional differential transform method [17] [18] , the fractional iteration method [19] , power series method [20] , Bessel collocation method [21] , wavelet [22] and the Haar wavelet method [23] .

Orthogonal polynomials play an important role in mathematics as well as in applications in mathematical physics, engineering and computer science. Chelyshkov polynomials are the most recent set of orthogonal polynomials. Chelyshkov polynomials have become increasingly important in numerical analysis. The efficiency of the method has been officially established by many researchers [24] [25] . With these backgrounds, we extend Chelyshkov-tau method for solving Bagley-Torvik equation of Equation (1).

The efficiency and accuracy of the numerical scheme is assessed on specific test problems. The numerical outcomes indicate that the method yields highly accurate results. The numerical solutions are compared with analytical and other existing numerical solutions in the literature.

The paper is organized as follows. Section 2 preliminarily provides some definitions which are crucial to the following discussion. In Section 3 we apply Chelyshkov-tau method for solving the model equation. In Section 4, we present numerical examples to exhibit the accuracy and the efficiency of the present method. where the numerical results presented in this paper are computed by Matlab programming. The conclusion is presented in the final section.

2. Preliminaries

2.1. Basic Definitions of Fractional

In this section, we introduce the basic necessary definitions and primary facts of the fractional calculus theory which will be more used in this work [26] [27] .

Definition 1. The Riemann-Liouville fractional integral operator $J^{\alpha }$ of order $\alpha$ on a usual Lebesgue space $L_1\left[a\mathrm{,}b\right]$ is given by

$J^{\alpha }\psi \left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^t}{\left(t-s\right)}^{\alpha -1}\psi \left(s\right)\text{d}s$

Some characteristics of this operator are:

$\begin{array}{l}{J}^0\psi \left(t\right)=\psi \left(t\right),J^{\alpha }{J}^{\beta }\psi \left(t\right)=J^{\alpha +\beta }\psi \left(t\right),J^{\alpha }{J}^{\beta }\psi \left(t\right)=J^{\beta }{J}^{\alpha }\psi \left(t\right),\left(\alpha ,\beta \ge 0\right)\\ J^{\alpha }{\left(t-s\right)}^{\mu }=\frac{\Gamma \left(\mu +1\right)}{\Gamma \left(\alpha +\mu +1\right)}{\left(t-s\right)}^{\alpha +\mu },\left(\mu >-1\right)\end{array}$

The Riemann-Liouville fractional derivative operator $D^{\alpha }$ is given by

$D^{\alpha }\psi \left(t\right)={\left(\frac{\text{d}}{\text{d}t}\right)}^m(J^{m-\alpha }\psi (t))$

where m is an integer, provided that $m-1<\alpha \le m$ .

However, its derivative has Some drawbacks when we try to make a model for a real phenomenon using fractional differential equations. Therefore, we will provide a modified fractional differential operator $D^{\alpha }$ proposed by by Caputo [28] .

Definition 2. The Caputo derivative definition is defined as

$D^{\alpha }\psi \left(t\right)=\frac{1}{\Gamma \left(m-\alpha \right)}{\displaystyle {\int }_0^x}{\left(t-s\right)}^{m-\alpha -1}{\psi }^{\left(m\right)}\left(s\right)\text{d}s\mathrm{,}m-\mathrm{1\; <mo>\; <\; </mo>}\alpha \le m\mathrm{,}m\in \mathbb{N}$ (3)

Hence, $\alpha >0$ and m is the smallest integer greater than or equal to $\alpha$ . For the Caputo fractional derivative we have

$\begin{array}{l}{D}^{\alpha }c=0,\left(c\text{is}\text{constant}\right),\\ D^{\alpha }{J}^{\alpha }\psi \left(t\right)=\psi \left(t\right),\\ J^{\alpha }{D}^{\alpha }\psi \left(t\right)=\psi \left(t\right)-{\displaystyle \sum _{i=0}^{m-1}}{\psi }^{\left(i\right)}\left(0^{+}\right)\frac{{\left(t-s\right)}^i}{i!},\\ D^{\alpha }{t}^{\beta }=\{\begin{array}{l}0\text{for}\beta <\lceil \alpha \rceil ,\\ \frac{\Gamma \left(\beta +1\right)}{\Gamma \left(\beta +1-\alpha \right)}{t}^{\beta -\alpha },\text{for}\beta \ge \lceil \alpha \rceil .\end{array}\end{array}$ (4)

It can be said when $\alpha \in \mathbb{N}$ , the Caputo differential operator matches with the integer-order differential operator.

For more details on fractional derivative definitions, theorems and its properties, you can see [26] [27] .

2.2. Some Properties of Chelyshkov Polynomials

We first review some important concepts and basics of the Chelyshkov function and conclude useful results that are important to this paper. Recently, these polynomials have established by Chelyshkov in [29] [30] [31] [32] [33] , which are orthogonal over the interval $\mathrm{<mo\; stretchy="false">\; [\; </mo>\; 0,1\; <mo\; stretchy="false">\; ]\; </mo>}$ , and are explicitly defined by

$C_{Nj}\left(x\right)={\displaystyle \sum _{i=0}^{N-j}}{\left(-1\right)}^i\left(\begin{array}{c}N-j\\ i\end{array}\right)\left(\begin{array}{c}N+j+k+1\\ N-j\end{array}\right)x^{N+k},j=0,1,\cdots ,N$ (5)

This gives the Rodrigues formula

$C_{Nj}\left(x\right)=\frac{1}{\left(N-j\right)!}\frac{1}{x^{j+1}}\frac{{\text{d}}^{N-j}}{\text{d}{x}^{N-j}}\left[x^{N+j+1}{\left(1-x\right)}^{N-j}\right],j=0,1,\cdots ,N$

and the orthogonality condition of Chelyshkov polynomials [32] is

${\displaystyle {\int }_0^1}{C}_{Nj}\left(x\right)C_{Nk}\left(x\right)\text{d}x=\{\begin{array}{l}\frac{1}{j+k+1},\text{for}j=k,j,k=0,1,\cdots ,N,N+1\\ 0,\text{for}j\ne k\mathrm{.}\end{array}$ (6)

Also it follows from this relation that

${\displaystyle {\int }_0^1}{C}_{Nj}\left(x\right)\text{d}x={\displaystyle {\int }_0^1}{x}^j\text{d}x=\frac{1}{j+1}$

By using the Cauchy integral formula for derivative and the Rodrigues type representation, we can get the integral relation

$C_{Nj}\left(x\right)=\frac{1}{2\text{π}i}\frac{1}{x^{j+2}}{\displaystyle {\int }_{\Omega 1}}\frac{{\overline{s}}^{N+j+2}{\left(1-s\right)}^{N-j}}{{\left(s-x^{-1}\right)}^{N-j+1}}\text{d}s$

such that the point $s=x^{-1}$ . lies in closed curve $\Omega 1$ .

Chelyshkov polynomials $C_{Nj}\left(x\right)$ provide a natural way to solve, expand, and interpret solutions. Actually, these polynomials can be expressed in terms of the Jacobi polynomials $P_k^{\left(\alpha \mathrm{,}\beta \right)}$ by the following relation,

$C_{Nj}\left(x\right)=x^j{P}_{N-j}^{\left(2j,1\right)}\left(1-2x\right),j=0,1,\cdots ,N$

Let function $u\left(x\right)$ , square integrable in $\left[\mathrm{0,1}\right]$ , can be expressed in terms of Chelyshkov polynomials as

$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum _{j=0}^N}{a}_j{C}_{Nj}\left(x\right),$ (7)

where the coefficients $a_j$ are the unknown Chelyshkov coefficients and $C_{Nj},j=0,1,\cdots ,N$ are Chelyshkov orthogonal polynomials of the degree N such that $N\ge 2$ . Also,

$\begin{array}{l}{u^{\prime }}_N(x)={\displaystyle \sum _{j=0}^N}{a}_j{C^{\prime }}_{Nj}\left(x\right),\\ {u^{″}}_N\left(x\right)={\displaystyle \sum _{j=0}^N}{a}_j{C^{″}}_{Nj}\left(x\right)\end{array}$ (8)

Then we can convert the solution expressed by (7) and its derivative (8) to matrix form

$\begin{array}{l}\left[u\left(x\right)\right]=C\left(x\right)A\text{or}\left[u\left(x\right)\right]=XCA\\ \left[u^{\prime }\left(x\right)\right]=\frac{\text{d}C\left(x\right)}{\text{d}x}A\text{or}\left[u^{\prime }\left(x\right)\right]=XMCA\\ \left[u^{″}\left(x\right)\right]=\frac{{\text{d}}^{2}C\left(x\right)}{\text{d}{x}^2}A\text{or}\left[u^{″}\left(x\right)\right]=X{M}^{2}CA\end{array}$ (9)

where

$A={\left[a_0,\cdots ,a_N\right]}^{\tau }\text{and}C\left(x\right)=\left[\begin{array}{ccccc}{C}_{N0}\left(x\right)& C_{N1}\left(x\right)& C_{N2}\left(x\right)& \cdots & C_{NN}\left(x\right)\end{array}\right]$

and

$X=\left[\begin{array}{ccccc}1& x& x^2& \cdots & x^N\end{array}\right],M={\left[\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 2& \dots & 0\\ 0& 0& 0& \ddots & 0\\ 0& 0& 0& \dots & N\\ 0& 0& 0& \dots & 0\end{array}\right]}_{\left(N+1\right)\times \left(N+1\right)}$

if N is odd,

if N is even,

Theorem 1. For $C_{Nj}\left(x\right)$ defined in (5) then the finite series can be converted

$u^{\left(3/2\right)}\left(x\right)\cong u_N^{\left(3/2\right)}\left(x\right)={\displaystyle \sum _{j=0}^N}{a}_j{C}_{Nj}^{\left(3/2\right)}(x)$

into matrix form

$\left[u^{\left(3/2\right)}\left(x\right)\right]=X{D}^{\left(3/2\right)}CA$ (10)

where

$D^{\left(3/2\right)}={\left[\begin{array}{ccccccc}0& 0& 0& 0& 0& \cdots & 0\\ 0& 0& 0& 0& 0& \cdots & 0\\ 0& 0& \frac{\Gamma \left(3\right)}{\Gamma \left(\frac{3}{2}\right)}{x}^{-3/2}& 0& 0& \cdots & 0\\ 0& 0& 0& \frac{\Gamma \left(4\right)}{\Gamma \left(\frac{5}{2}\right)}{x}^{-3/2}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& 0& 0& \dots & \frac{\Gamma \left(N+1\right)}{\Gamma \left(\frac{2N-1}{2}\right)}{x}^{-3/2}\end{array}\right]}_{\left(N+1\right)\times \left(N+1\right)}$ (11)

Proof. The proof is straightforward using Equation (4).

3. The Description of Chelyshkov Scheme

Let us seek the solution of (1) expressed in terms of Chelyshkov polynomials as

$u\left(x\right)\cong u_N\left(x\right)={\displaystyle \sum _{n=0}^N}{a}_n{C}_{Nn}\left(x\right).$

Replacing each terms of (1) with the corresponding approximations defined in (7), (8) and (10) and we obtain the following theorem.

Theorem 2. If the assumed approximate solution of the boundary-value problem (1)-(2) is (3), the discrete Chelyshkov-tau system for the determination of the unknown coefficients ${\left\{a_n\right\}}_{n=0}^N$ is given by

$A_2{\displaystyle \sum _{n=0}^N}{a}_n{C}_{Nn}^{\left(2\right)}\left(x\right)+A_1{\displaystyle \sum _{n=0}^N}{a}_n{C}_{Nn}^{\left(3/2\right)}\left(x\right)+A_0{\displaystyle \sum _{n=0}^N}{a}_n{C}_{Nn}\left(x\right)={\displaystyle \sum _{n=0}^N}{f}_n{C}_{Nn}\left(x\right)$ (12)

Equation (12), which can be written in the matrix form

$A_{2}X{M}^{2}CA+A_{1}X{D}^{\left(3/2\right)}CA+A_{0}XCA=XCF$ (13)

where

$F={\left[f_0\mathrm{,}{f}_1\mathrm{,}\cdots \mathrm{,}{f}_N\right]}^{\tau }$

The residual $R_N\left(x\right)$ for Equation (13) can be written as

$R_N\left(x\right)=X\left[A_2{M}^{2}CA+A_1{D}^{\left(3/2\right)}CA+A_{0}CA-CF\right]$ (14)

As in a typical tau method [34] we generate N-1 linear equations by applying

$\langle R_N\left(x\right),C_{Nn}\left(x\right)\rangle ={\displaystyle {\int }_0^1}{R}_N\left(x\right)C_{Nn}\left(x\right)\text{d}x=0,n=1,2,\cdots ,N-1$ (15)

The boundary condition is derived from Equation (2) and matrices for conditions are

${\displaystyle \sum _{j=0}^1}{\alpha }_{ij}X\left(0\right)M^{j}CA+{\beta }_{ij}X\left(1\right)M^{j}CA={\gamma }_i,i=0,1$ (16)

Equations (15) and (16) generate $N+1$ set of linear equations, respectively. These linear equations can be solved for unknown coefficients of the vector $A$ . Consequently, $u\left(x\right)$ given in Equation (7) can be calculated..

4. Numerical Results

In this section, we apply the Chelyshkov-tau method to various problems which were collected from the open literature [13] [35] [36] [37] [38] . Our primary interest is to compare our method with other methods on the same problems. All computations were carried out using Matlab on a personal computer. In the examples, the maximum absolute error at points is taken as

$\Vert E_{\text{Chelyshkov}}\Vert =\left|u_{\text{Exact}}-u_{\text{Chelyshkov}}\right|$

Example 1: [21] Consider the linear BVP

$u^{″}\left(x\right)+u^{\left(3/2\right)}\left(x\right)+u\left(x\right)=x+1,0<x<1$

subject to the boundary conditions

$u\left(0\right)=1\text{and}u\left(1\right)=2$

whose exact solution is

$u=x+1$

the approximate solution $u\left(x\right)$ by the truncated Chelyshkov polynomial for $N=2$ is

$u\left(x\right)=a_0{C}_{20}\left(x\right)+a_1{C}_{21}\left(x\right)+a_2{C}_{22}(x)$

Here, we have

$C=\left[\begin{array}{ccc}3& 0& 0\\ -12& 4& 0\\ 10& -5& 1\end{array}\right],M=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 2\\ 0& 0& 0\end{array}\right],D^{\left(3/2\right)}=\frac{4}{\sqrt{\text{π}}}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& x^{-3/2}\end{array}\right]$

$X=\left[\begin{array}{ccc}1& x& x^2\end{array}\right],F=\left[\begin{array}{c}\frac{1}{3}\\ \frac{5}{4}\\ \frac{35}{12}\end{array}\right]$ (17)

By applying Equation (15) We obtain

$\left(\frac{23}{3}+\frac{16}{7\sqrt{\text{π}}}\right)a_0-\left(\frac{10}{3}+\frac{8}{7\sqrt{\text{π}}}\right)a_1+\left(\frac{2}{3}+\frac{8}{35\sqrt{\text{π}}}\right)a_2-\frac{1}{3}=0$ (18)

By applying Equation (16), we have

$3a_0-1=0$ (19)

$a_0-a_1+a_2-2=0$ (20)

By solving Equations (18)-(20), we get

$a_0=\frac{1}{3},a_1=\frac{5}{4},a_2=\frac{35}{12}$

Thus we can write

$y\left(x\right)=\left[\begin{array}{ccc}\frac{1}{3}& \frac{5}{4}& \frac{35}{12}\end{array}\right]\left[\begin{array}{c}10x^2-12x+3\\ -5x^2+4x\\ x^2\end{array}\right]=1+x$

which is the exact solution.

Example 2: [13] [35] [36] Consider the linear BVP

$A_2{u}^{″}\left(x\right)+A_1{u}^{\left(3/2\right)}\left(x\right)+A_{0}u\left(x\right)=f\left(x\right),0<x<1$

with initial conditions

$u\left(0\right)=0\text{and}{u}^{\prime }\left(0\right)=0$

which is known to have analytical solution as

$\begin{array}{c}u\left(x\right)={\displaystyle {\int }_0^x}{G}_3\left(x-\tau \right)f\left(\tau \right)\text{d}\tau G_3\left(x\right)\\ =\frac{1}{A_2}{\displaystyle \sum _{k=0}^{\infty }}\frac{{\left(-1\right)}^k}{k!}{\left(\frac{A_0}{A_2}\right)}^k{x}^{2k+1}{E}_{\frac{1}{2},2+\frac{3k}{2}}^{\left(k\right)}\left(-\frac{A_1}{A_2}\sqrt{x}\right)\end{array}$

where $E_{\lambda \mathrm{,}\mu }^{\left(k\right)}\left(u\right)$ is the kth derivative of the Mittag-Leffler function with parameters $\lambda$ and $\mu$ given by

$E_{\lambda ,\mu }^{\left(k\right)}\left(y\right)={\displaystyle \sum _{j=0}^{\infty }}\frac{\left(j+k\right)!y^j}{j!\Gamma \left(\lambda j+\lambda k+\mu \right)},k=0,1,2,\cdots$

and the $G_3\left(x\right)$ three-term Green’s equation. Let $A_2=1$ , $A_1=A_0=0.5$ and $f\left(x\right)=8$ .

Table 1 exhibits a comparison between the exact, the results obtained by using Chelyshkov tau for $N=14$ with analogous results of Çenesiz et al. [35] for underlying the generalized Taylor collocation method (GTCM) and Setia [36] , who used second kind Chebyshev wavelet method (CWM) and with analogous results of El-Gamel and Abd El-Hady [13] for underlying Legendre-collocation method.

Figure 1 displays the estimated absolute error function for $N=14$ with the present method.

Example 3: [13] [19] Now we turn to IVP

$u^{″}\left(x\right)+u^{\left(3/2\right)}\left(x\right)+u\left(x\right)=7x+\frac{8}{\sqrt{\text{π}}}{x}^{3/2}+x^3+1,0<x<1$

subject to the boundary conditions

$u\left(0\right)=1\text{and}{u}^{\prime }\left(0\right)=1$

whose exact solution is

$u\left(x\right)=x^3+x+1$

Table 2 exhibits a comparison between the results obtained by using Chelyshkov tau for $N=8$ with analogous results of Mekkaoui and Hammouch [19] for underlying the variational iteration method (VIM), the fractional iteration method (FIM) and with analogous results of El-Gamel and Abd El-Hady [13] for underlying Legendre-collocation method.

Figure 2 displays the estimated absolute error function for $N=8$ with the present method.

Example 4: [38] Consider the linear BVP

$u^{″}\left(x\right)+\frac{8}{17}{u}^{\left(3/2\right)}\left(x\right)+\frac{13}{51}u\left(x\right)=f\left(x\right),0<x<1$

where

$f\left(x\right)=\frac{1}{89250\sqrt{\text{π}x}}\left[48p\left(x\right)+7\sqrt{x}q\left(x\right)\right]$

$p\left(x\right)=16000x^4-32480x^3+21280x^2-4746x+149$

and

$q\left(x\right)=3250x^5-9425x^4+264880x^3-448107x^2+233262x-34578$

and subject to the boundary conditions

$y\left(0\right)=0\text{and}y\left(1\right)=0$

whose exact solution is

$y\left(x\right)=x^5-\frac{29}{10}{x}^4+\frac{76}{25}{x}^3-\frac{339}{250}{x}^2+\frac{27}{125}x$

Table 3 exhibits a comparison between the absolute errors obtained by using Chelyshkov tau for $N=8$ with analogous errors of Rehman and Ali Khan [38] for underlying Haar wavelets method.

Figure 3 displays the estimated absolute error function for $N=8$ with the present method.

5. Conclusion

In this paper, Chelyshkov operational matrix of fractional derivative has been derived. Our approach was based on the tau method. The proposed technique is easy to implement efficiently and yield accurate results. Moreover, only a small number of Chelyshkov polynomials is needed to obtain a satisfactory result. In addition, an interesting feature of this method is to find the analytical solution if the equation has an exact solution that is polynomial functions. Numerical examples are included and a comparison is made with an existing result.

Cite this paper

El-Gamel, M., Abd-El-Hady, M. and El-Azab, M. (2017) Chelyshkov-Tau Approach for Solving Bagley-Torvik Equation. Applied Mathematics, 8, 1795-1807. https://doi.org/10.4236/am.2017.812128

References