Numerical Solution of the Distributed-Order Fractional Bagley-Torvik Equation

Ⅰ. INTRODUCTION

The fractional order calculus establishes the branch of mathematics dealing with differentiation and integration under an arbitrary order of the operation, that is the order can be any real or even complex number, not only the integer one. Although the history of fractional calculus is more than three centuries old, it only has received much attention and interest in the past 20 years; the reader may refer to [1], [2] for the theory and applications of fractional calculus. The generalization of dynamical equations using fractional derivatives proved to be useful and more accurate in mathematical modeling related to many interdisciplinary areas. Applications of fractional order differential equations include: electrochemistry [3], porous media [4] etc.

The fractional differential operator of distributed order

is a generalization of the single order $D^{\alpha}=\frac{d^{\alpha}}{dt^{\alpha}}\; (\alpha\in[\gamma_{1}, \gamma_{2}])$ which by considering a continuous or discrete distribution of fractional derivative is obtained. In fact, one can find that the time distributed-order differential equation is a continuous generalization of the single-order time-fractional one and the multi-term time-fractional in time. If the nonnegative density function in the distributed integral takes the single Dirac-$\delta$ function, the distributed-order time-differential equation is reduced to the single-order one. And if the linear combination of several Dirac-$\delta$ functions is taken as the nonnegative density function, then the multi-term one is recovered. The idea of fractional derivative of distributed order is stated by Caputo [5] and later developed by Caputo himself [6], Bagley and Torvik [7], [8]. Other researchers used this idea, and interesting reviews appeared to describe the related mathematical models of fractional differential equation of distributed order. For example, distributed-order differential equations have been used in [9] to model the stress-strain behavior of an anelastic medium and in [5] to find the eigenfunctions of the torsional models of anelastic or dielectric spherical shells and infinite plates; in [10] further dielectric models and diffusion equations lead to distributed-order differential equations. In [7], [8] distributed-order differential equations provide models of the input-output relationship of a linear time-variant system based on frequency domain observations. Furthermore, the stability analysis of distributed order fractional differential equations with respect to the nonnegative density function was studied in [11], [12].

The numerical methods of the distributed-order equations have been a main goal in researches. Hence, there are some papers discussing numerical methods of the distributed-order equations. For example, Diethelm et al. [13], and Katsikadelis [14] introduced effective numerical methods for solving linear and nonlinear distributed-order ordinary differential equations.

In this paper two numerical methods are presented for solving distributed-order Bagley-Torvik equation of the form [15]:

where $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{r}{w(\alpha)}\; ^{c}D_{t}^{\alpha}y(t)d\alpha$ $(r=1\; {\rm or}\; 2)$, with initial conditions

where $a, b, c, \upsilon_{1}$ and $\upsilon_{2}$ are constant coefficients and the function $f(x)$ is continuous. In particular, if $w(\alpha)=\delta(\alpha- {1}/{2})$ or $w(\alpha)=\delta(\alpha- {3}/{2})$, then (2) can model the frequency-dependent damping materials quite satisfactorily. It can also describe motion of real physical systems, the modeling of the motion of a rigid plate immersed in a Newtonian fluid and a gas in a fluid, respectively [16], [17]. Since the distributed-order operator in (2) is represented by an integral of weighted constant-order operator, thus, the recipe for the distributed-order Bagley-Torvik equation consist of two steps. Firstly, we transform the distributed-order equation into a multi-term time-fractional differential equation by applying the composite Boole's rule [18]. Then we apply the Grunwald-Letnikov method (GLM) [19] and the fractional differential transform method (FDTM) [20], to approximate the multi-term time-fractional differential equation. Finally, we compare the results obtained using the Grunwald-Letnikov method and the fractional differential transform method, with the exact solution for some cases to show the effectiveness of our methods.

Ⅱ. PRELIMINARIES AND NOTATIONS

In this section, we consider the main definitions of fractional derivative operators of single and distributed order.

The Riemann-Liouville definition is given as

where $m$ is the first integer which is not less than $\alpha$, i.e., $m-1 < \alpha\leq m$ and $\Gamma(.)$ is a Gamma function.

The Grunwald-Letnikov definition is given by

where $[\cdot]$ means the integer part.

Because the Grunwald-Letnikov's definition is the most straightforward from the viewpoint of numerical implementation so we will use it for solving fractional differential equations. The approximate Grunwald-Letnikov's definition is given below, where the step size of $h$ is assumed to be very small [21], [22],

where $t_{k}=kh$ $(k =1, \ldots, \frac{t-a}{h})$ and $c_{i}^{(\alpha)}$ $(i =0, 1, \ldots, k)$ are binomial coefficients, which can be computed as [22]

Finally, the Caputo fractional derivative of $y(t)$ is defined as

for $m-1 < \alpha\leq m$, $m\in {\mathbb N}$. The Caputo's definition has the advantage of dealing properly with initial value problems in which the initial conditions are given in terms of the field variables and their integer order which is the case in most physical processes.

Now, we generalize the above definition of the fractional derivative of distributed-order in the Caputo sense with respect to order-density function $w(\alpha)$ as follows

where $^{c}_{0}D_{t}^{\alpha}=^{c} \!\!D_{t}^{\alpha}$. With out restricting the generality, the solution procedure is presented for $0\leq \nu < \beta \leq 2$.

Ⅲ. NUMERICAL METHODS

In this section, we introduce the composite Boole's rule, the Grunwald-Letnikov method and the fractional differential transform method to obtain approximate solutions of the distributed-order Bagley-Torvik equation.

A The Composite Boole's Rule

Here we demonstrate composite Boole's rule to discretize the distribution integral in (2). Suppose that the interval $[0, r]$ is subdivided into $4n$ subintervals $[\alpha_{k}, \alpha_{k+1}]$ of width $H= {r}/{4n}$ by using the equally spaced nodes $\alpha_{k} = \alpha_{0} + kH$, for $k=0, 1, \ldots, 4n$. Approximating the integral term in (2) using the composite Boole's rule with $4n$ equal subintervals, we obtain

with $^{c}D_{t}^{\alpha_0}y(t)=y(t)$ and $^{c}D_{t}^{\alpha_4n}y(t)= {^{c}D_{t}^{r}y(t)}$.

Hence, (2) may be approximated via

with initial conditions

Remark 1: When we approximate the integral in the distributed-order equation by a finite sum then the error for the composite Boole's rule is of the order $O(H^{6})$ [18].

In 2004, Diethelm and Ford [23] introduced the first scheme for the solution of multi-order fractional differential equations of the form (11). They converted equations with commensurate multiple fractional derivatives into a very simple linear system of fractional differential equations of low order and then the application of any single-term equation solver solves the problem.

Let $M$ be the least common multiple of the denominators of $\alpha_{1}, \alpha_{2}, \ldots, \alpha_{4n}$, and set $\gamma = 1/M$ and $N = 2M$. We have the following theorem on equivalence of a nonlinear system

Theorem 1 [23]: Equation (11), equipped with the initial conditions (12), is equivalent to the system of equations

together with the initial conditions

in the following sense:

1) Whenever $Y = (y_{0}, \ldots, y_{N-1})^{T}$ with $y_{0}\in {\mathbb C}^{2} [0, b]$ for some $b > 0$ is the solution of the system (13), equipped with the corresponding initial conditions, the function $y = y_{0}$ solves the multi-order (11), and it satisfies the initial conditions (12).

2) Whenever $y_{0}\in {\mathbb C}^{2} [0, b]$ is a solution of the multi-order (11) satisfying the initial conditions (12), the vector-valued function

satisfies the system (13) and the initial conditions (14).

Remark 2: If $w(\alpha)=\sum_{j=1}^{n}c_{j}\delta(\alpha-j\gamma)$, $(0 < n\gamma\leq1~\mbox{or}\; 2)$, then the distributed-order Bagley-Torvik (2) converts to multi-order time-fractional differential equation, therefore the first step is not required.

B The Grunwald-Letnikov Method

Now, numerical method to solving the system (13) and the initial conditions (14) is presented, which is based on the Grunwald-Letnikov's approximation described in [19].

Suppose $P$ is integer numbers so that $T = hP$. We use a uniform grid

Let $\tilde{y}(t_{k})$ denote the approximation to $y(t_{k})$ and also we have already calculated approximation $\tilde{y}(t_{k})\approx y(t_{k})$, $(k=0, 1, \ldots, n-1 \leq P)$. Therefore, a general numerical solution of the system (13) and the initial conditions (14) can be expressed by using the relations (6) and (7) as

subject to initial conditions

In [19], the consistence, convergence and stability of the numerical method mentioned above are studied.

Remark 3: The error associated with the approximation of a multi-term equation using the Grunwald-Letnikov method, is of the order $O(h^{(1+\gamma)})$.

C The Fractional Differential Transform Method

The fractional differential transform method is applied to solve fractional differential equations. For applying this method, one should be familiar with the basic idea of the fractional differential transform method. Let us expand the analytical and continuous function $f(t)$ in terms of a fractional power series as follows

where $\alpha$ is the order of fraction and $F(k)$ is the fractional differential transform of $f(t)$.

In order to avoid fractional initial and boundary conditions, we define the fractional derivative in the Caputo sense. The relation between the Riemann-Liouville operator and Caputo operator is given by

where $m-1 < q\leq m$ and $t>t_0$.

For the sake of simplicity of the following transformation the initial condition will be implemented to the integer order derivatives

for $k=0, 1, 2, \ldots, (q \alpha -1)$, where $q$ is the order of fractional differential equation considered.

Using (4) and (18), the following theorems of fractional differential transform method are given below. A complete discussion on the proofs of these theorems can be found in [20], [21].

Theorem 2: If $f(t)=g(t)\pm h(t)$, then $F(k)=G(k)\pm H(k)$.

Theorem 3: If $f(t)=g(t) h(t)$, then $F(k)=\sum_{l=0}^{k}G(l)\times$\\ $H(k-l)$.

Theorem 4: If $f(t)=g_1(t) g_2(x)\ldots g_{n-1}(t) g_n(t)$, then

Theorem 5: If $f(t)=(t-t_0)^p$, then $F(k)=\delta(k-\alpha p)$, where

Theorem 6: If $f(t)=D_{t_0}^q[g(t)]$, then

Now assume that $q=\gamma$ and $\alpha=M$, thus, according to fractional differential transform method, by taking differential transform on both sides of the system (13) and the initial conditions (14) are transformed as follows

subject to initial conditions

Therefore, according to differential transform method the $K$-term approximations for (13) can be expressed as

Ⅳ. EXAMPLE

In this section, we present some examples to illustrate the efficiency of the methods. In the following, five different cases of the weighting function of order are discussed respectively.

Case 1: $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{1}\; {\delta(\alpha-\frac{1}{2})}\; ^{c}D_{t}^{\alpha}y(t)\; d\alpha$ and $a=b=c=1$, with the initial conditions

where

and the exact solution is $y(t)=t^3$. The obtained numerical result by means of the proposed FDTM solution, GLM solution and exact solution are shown in Fig. 1, with the final time $T=10$, $K=8$ and $h=0.01$ when ${w(\alpha)}$ be Case 1.

Figure  1.  Comparison between the exact solution, FDTM solution and the GL solution of (2) for $T=10$, $K=8$ and $h=0.01$ when ${w(\alpha)}$ be Case 1.

DownLoad: Full-Size Img PowerPoint

Case 2: $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{2}\; {\delta(\alpha-\frac{3}{2})}\; ^{c}D_{t}^{\alpha}y(t)\; d\alpha$ and $a=b=c=1$, with the initial conditions

where

and the exact solution is $y(t)=t^2$. Fig. 2 presents a comparison of the numerical solutions FDTM, GLM and the exact solution by step size $h=0.01$, $K=8$ and $T = 10$ when ${w(\alpha)}$ be Case 2.

Figure  2.  Comparison between the exact solution, FDTM solution and the GLM solution of equation (2) for $T=10$, $K=8$ and $h=0.01$ when ${w(\alpha)}$ be Case 2.

DownLoad: Full-Size Img PowerPoint

Case 3: $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{2}\; {\delta(\alpha-\frac{3}{2})+\delta(\alpha-1)+\delta(\alpha-\frac{1}{2})}$ $^{c}D_{t}^{\alpha}y(t)\; d\alpha$ and $a=b=c=1$, with the initial conditions

where

and the exact solution is $y(t)=t^3$. Comparison of numerical results with the exact solution is shown in Fig. 3 by step size $h=0.01$, $K=8$ and $T = 10$ when ${w(\alpha)}$ be Case 3.

Figure  3.  Comparison between the exact solution, FDTM solution and the GLM solution of (2) for $T=10$, $K=8$ and $h=0.01$ when ${w(\alpha)}$ be Case 3.

DownLoad: Full-Size Img PowerPoint

Case 4: $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{1} {6\alpha(1-\alpha)} ^{c}D_{t}^{\alpha}y(t) d\alpha$ and $a=1, b=0.5, c=1.5$, with the initial conditions

where

We show the numerical solutions in Fig. 4 with step size $h = 0.01$, $H=0.25$, $T =10$ and $K= {T}/{h}$ when ${w(\alpha)}$ be Case 4.

Figure  4.  GLM solution and FDTM solution of the distributed-order equation (2) for $T=10$, $h=0.01$, $H=0.25$ and $K= {T}/{h}$ when ${w(\alpha)}$ be Case 4.

DownLoad: Full-Size Img PowerPoint

Case 5: $^{c}D_{t}^{w(\alpha)}y(t)=\int_{0}^{1}\; {\Gamma(4-\alpha)}\; ^{c}D_{t}^{\alpha}y(t)\; d\alpha$ and $a=1, b=2, c=0.5$, with the initial conditions

where

We show the numerical solutions in Fig. 5 with step size $h = 0.01$, $H=0.25$, $T =10$ and $K= {T}/{h}$ when ${w(\alpha)}$ be Case 4.

Figure  5.  GLM solution and FDTM solution of the distributed-order (2) for $T=10$, $h=0.01$, $H=0.25$ and $K= {T}/{h}$ when ${w(\alpha)}$ be Case 5.

DownLoad: Full-Size Img PowerPoint

Ⅴ. CONCLUSION

In this paper, two schemes for the distributed-order Bagley-Torvik equation with respect to the nonnegative density function have been described. First, by approximating the integral term in the distributed-order equation using the composite Boole's rule, we obtained a multi-term time-fractional equation. Afterwards, the multi-term time-fractional differential equation is solved by the Grunwald-Letnikov method and the fractional differential transform method. Eventually, five examples of distributed-order Bagley-Torvik equation are presented to illustrate the efficiency and reliability of the methods. All numerical results are obtained using MATLAB 7.11.