# Fractional Leibniz integral rules for Riemann-Liouville and Caputo
fractional derivatives and their applications

Ismail T. Huseynov<EMAIL_ADDRESS>Arzu Ahmadova
<EMAIL_ADDRESS>Nazim I. Mahmudov<EMAIL_ADDRESS>Department of Mathematics, Faculty of Arts and Sciences, Eastern Mediterranean
University, Mersin 10, Gazimagusa, TRNC, Turkey

###### Abstract

In recent years, the theory for Leibniz integral rule in the fractional sense
has not been able to get substantial development. As an urgent problem to be
solved, we study a Leibniz integral rule for Riemann-Liouville and Caputo type
differentiation operators with general fractional-order of $n-1<\alpha\leq n$,
$n\in\mathbb{N}$ . A rule of fractional differentiation under integral sign
with general order is necessary and applicable tool for verification by
substitution for candidate solutions of inhomogeneous multi-term fractional
differential equations. We derive explicit analytical solutions of generalized
Bagley-Torvik equations in terms of recently defined bivariate Mittag-Leffler
type functions that based on fractional Green’s function method and verified
solutions by substitution in accordance by applying the fractional Leibniz
integral rule. Furthermore, we study an oscillator equation as a special case
of differential equations with multi-orders via the Leibniz integral rule.

###### keywords:

Caputo fractional derivative , Riemann-Liouville fractional derivative ,
Leibniz integral rule , Bagley-Torvik equation , oscillator equation ,
bivariate Mittag-Leffler function

##  1 Introduction

Fractional calculus is a generalization of the classical differential calculus
which has attracted growing attention due to the applications for many
problems in science and engineering such as reaction-diffusion systems [1],
viscoelasticity [2], electrical circuits [3, 4], control theory [5],
stochastic analysis [6] and time-delay systems [7].

One of the most frequently encountered tools in the theory of fractional
calculus is furnished by the Riemann-Liouvulle
$\prescript{RL}{t_{0}}{D}^{\alpha}_{t}$ and Caputo
$\prescript{C}{t_{0}}{D}^{\alpha}_{t}$ fractional differentiation operators. A
fractional analogue of Leibniz rule for differentiation is crucial and useful
properties of these operators. Podlubny in [8], Baleanu and Trujillo in [9]
give a proof the Leibniz rule for Riemann-Liouville and Caputo type
derivatives; these results are stated respectively below.

Assume that $0<\alpha<1$ and $f,g:[t_{0},T]\subset\mathbb{R}\to\mathbb{R}$
with all their derivatives are continuous. Then

$\displaystyle\prescript{RL}{t_{0}}{D}^{\alpha}_{t}\left\\{f(t)g(t)\right\\}=\sum_{k=0}^{\infty}\binom{\alpha}{k}f^{(k)}(t)\prescript{RL}{t_{0}}{D}^{\alpha-k}_{t}g(t),\quad
t\in(t_{0},T),$ (1.1)
$\displaystyle\prescript{C}{t_{0}}{D}^{\alpha}_{t}\left\\{f(t)g(t)\right\\}=\sum_{k=0}^{\infty}\binom{\alpha}{k}f^{(k)}(t)\prescript{RL}{t_{0}}{D}^{\alpha-k}_{t}g(t)-\frac{(t-t_{0})^{-\alpha}}{\Gamma(1-\alpha)}f(t_{0})g(t_{0}),\quad
t\in(t_{0},T),$ (1.2)

where binomial coefficients satisfy the identity:

$\binom{\alpha}{0}=1,\quad\binom{\alpha}{k}=\frac{\alpha(\alpha-1)\ldots(\alpha-k+1)}{k!}.$

Another essential property of Riemann-Liouville fractional differentiation
operator is obtained by Podlubny in [8] which is called fractional Leibniz
integral rule stated as below:

$\prescript{RL}{t_{0}}{D}^{\alpha}_{t}\int\limits_{0}^{t}K(t,\tau)\mathrm{d}\tau=\lim\limits_{\tau\to
t-0}\prescript{RL}{\tau}{D}^{\alpha-1}_{t}K(t,\tau)+\int\limits_{0}^{t}\prescript{RL}{\tau}{D}^{\alpha}_{t}K(t,\tau)\mathrm{d}\tau,\quad\alpha\in(0,1),$
(1.3)

where lower terminal $t_{0}=0$.

The following important particular case for convolution operator whenever we
have $K(t-\tau)f(\tau)$ instead of $K(t,\tau)$, the relationship (1.3) takes
the form:

$\prescript{RL}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}K(t-\tau)f(\tau)\mathrm{d}\tau=\lim\limits_{\tau\to
t-0}f(t-\tau)\prescript{RL}{\tau}{D}^{\alpha-1}_{t}K(\tau)+\int\limits_{0}^{t}\prescript{RL}{0}{D}^{\alpha}_{\tau}K(\tau)f(t-\tau)\mathrm{d}\tau,\quad\alpha\in(0,1).$
(1.4)

It is important to note that the above tools are necessary for checking by
substitution method for fractional differential equations with variable and
constant coefficients.

Fractional differential equations (FDEs) are differential equations involving
derivatives of arbitrary (fractional) order. FDEs provide one of the most
accurate tools to describe hereditary properties of natural phenomenon. Using
fractional derivatives instead of integer-order derivatives allows us for the
modeling of a wider variety of behaviours. However, sometimes, FDEs involving
one fractional order of differentiation are not sufficient to demonstrate
physical processes. Therefore, recently, several authors have studied more
general types of fractional-order models, such as multi-term equations [10,
11, 12, 13, 14] and multi-dimensional systems [4, 6, 7, 15, 16, 17, 18].

Multi-term differential equations with fractional-order have been studied and
solved using various mathematical methods, of which we mention a few as
follows. Luchko and several collaborates [10, 11] have used the method of
operational calculus to solve multi-term FDEs with constant coefficients with
regard to various types of fractional derivatives. Bazhlekova [12] has
considered multi-term fractional relaxation equations with Caputo fractional
derivatives by using a Laplace transform technique, and studied the
fundamental and impulse-response solutions of the initial value problem (IVP).
Kaczorek and Idczak [17] have considered existence and uniqueness results and
a Cauchy formula for the analytical solution of the time-varying linear system
with Caputo fractional derivative. Pak et al. [13] has recently investigated
multi-term FDEs with variable coefficients using a new method to construct
analytical solutions.

As one of the important special cases of multi-term FDEs, Bagley-Torvik
equations have been discussed in terms of analytical [8, 19, 20, 21] and
numerical methods [22, 23]. Bagley-Torvik equations with $\frac{1}{2}$-order
or $\frac{3}{2}$-order derivative describe the motion of real physical systems
in a Newtonian fluid [19]. In 1984, Bagley and Torvik [19] have considered the
following Cauchy problem under the homogeneous initial conditions:

$\displaystyle
my^{\prime\prime}(r)+\frac{2S\sqrt{\mu\rho}}{m}\left(\prescript{C}{}{D}^{\alpha}_{0_{+}}y\right)(r)+ky(r)=g(r),\quad
r>0,$ (1.5) $\displaystyle y(0)=y^{\prime}(0)=0,$

where $\left(\prescript{C}{}{D}^{\alpha}_{0_{+}}y\right)(\cdot)$ is Caputo
fractional differential operator of order $\alpha=\frac{1}{2}$ or
$\alpha=\frac{3}{2}$, $S$\- an area of the rigid plate, $\mu$\- viscosity,
$\rho$-fluid density, $m$-mass, $k$\- spring of stiffness and $g(\cdot)$-an
external force. An analytical solution of (1.5) has introduced by Podlubny [8]
in the form:

$y(r)=\int\limits_{0}^{r}G(r-\tau)g(\tau)\mathrm{d}\tau,\quad r>0,$ (1.6)

with

$G(r)=\frac{1}{m}\sum_{l=0}^{\infty}\frac{(-1)^{l}}{l!}\left(\frac{k}{m}\right)^{l}r^{2l+1}\mathcal{E}^{(l)}_{\frac{1}{2},2+\frac{3l}{2}}\left(\frac{-2S\sqrt{\mu\rho}}{m}\sqrt{r}\right),$

where $\mathcal{E}^{(l)}_{\alpha,\beta}(\cdot)$ is the $l$th-derivative of
two-parameter Mittag-Leffler function. In [20], Mahmudov et al. have studied
explicit analytical solutions for several families of generalized
multidimensional Bagley-Torvik equations with permutable matrices. In [21],
Wang et al. have modified the following Bagley-Torvik equation

$y^{\prime\prime}(r)+\mu\left(\prescript{C}{}{D}^{\alpha}_{0_{+}}y\right)(r)+y(r)=0,\quad\mu,r>0,$
(1.7)

where $\alpha=\frac{1}{2}$ or $\alpha=\frac{3}{2}$, to the sequential FDEs and
introduced a general solution of (1.7) by using the technique related to
characteristic roots. The numerical point of view Diethelm and Ford in [22]
have used linear multi-steps , Srivastava et al. in [23] have applied wavelet
approach to obtain approximate solutions of the Bagley-Torvik equations.

Therefore, the plan of this paper is systematized as below. Section 2 is a
mathematical preliminary section where we recall main definitions and results
from fractional calculus, special functions and necessary lemmas from
fractional differential equations. Section 3 is devoted to formulating the
Leibniz integral rule for higher order derivatives of Lebesgue integration
which depends on parameter in classical sense. In Section 4, we have
introduced fractional differentiation under the integral sign in Riemann-
Liouville and Caputo sense. Moreover, we have considered the derivative of
convolution operator which has more importance for differential equations with
classical or fractional order. In Section 5, we have acquired explicit
analytical solutions of Bagley-Torvik equations with Riemann-Liouville and
Caputo type fractional derivatives in terms of recently defined bivariate
Mittag-Leffler type functions in accordance with fractional Green’s function
method and tested the candidate solutions by using our newly defined tools
which are natural generalization of well-known Leibniz integral rule. At the
end, in Section 6 we give the conclusions and future directions.

##  2 Mathematical preliminaries

We embark on this section by briefly introducing the essential structure of
fractional calculus, special functions and fractional differential operators
(for the more salient details on the matter, see the textbooks [8, 25, 26, 28,
42, 44]). We begin by defining some notations, Riemann-Liouville and Caputo
fractional differentiation operators which are fundamental for fractional
calculus and fractional differential equations.

Let $\mathbb{R}^{n}$ be Euclidean space and $\mathbb{J}$ be some interval of
the real line, i.e. $\mathbb{J}\subset\mathbb{R}$. We suppose that
$\mathbb{J}=[t_{0},T]$ for some $t\in\mathbb{\hat{J}}$ and denote
$\mathbb{\hat{J}}=(t_{0},T)$. Assume that $f:\mathbb{J}\to\mathbb{R}$ is an
absolutely continuous function.

######  Definition 2.1 ([28, 42, 44]).

The Riemann-Liouville derivative operator of fractional order $n-1<\alpha\leq
n$ for $n\in\mathbb{N}$ is defined by

$\left(\prescript{RL}{t_{0}}{D}^{\alpha}_{t}g\right)(t)=\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\left(\prescript{}{t_{0}}{I}^{n-\alpha}_{t}g\right)(t)=\frac{1}{\Gamma(n-\alpha)}\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}\int\limits^{t}_{t_{0}}(t-s)^{n-\alpha-1}g(s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}},$ (2.1)

where $\prescript{}{t_{0}}{I}^{\alpha}_{t}$ is the Riemann-Liouville integral
operator of order $\alpha>0$ which is defined by

$\left(\prescript{}{t_{0}}{I}^{\alpha}_{t}g\right)(t)=\frac{1}{\Gamma(\alpha)}\int\limits_{t_{0}}^{t}(t-s)^{\alpha-1}g(s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$ (2.2)

Furthermore, the following equality holds true:

$\left(\prescript{RL}{t_{0}}{D}^{\alpha}_{t}\left(\prescript{}{t_{0}}{I}^{\alpha}_{t}g\right)\right)(t)=g(t),\quad\alpha>0,\quad
t\in\mathbb{\hat{J}}.$ (2.3)

######  Definition 2.2 ([8, 25, 27]).

The Caputo derivative operator of fractional order $n-1<\alpha\leq n$ for
$n\in\mathbb{N}$ is defined by

$\left(\prescript{C}{t_{0}}{D}^{\alpha}_{t}g\right)(t)=\prescript{}{t_{0}}{I}^{n-\alpha}_{t}\left(\frac{\mathrm{d}^{n}}{\mathrm{d}t^{n}}g\right)(t)=\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}(t-s)^{n-\alpha-1}\frac{\mathrm{d}^{n}}{\mathrm{d}s^{n}}g(s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$ (2.4)

Moreover, the next relation holds true:

$\left(\prescript{C}{t_{0}}{D}^{\alpha}_{t}\left(\prescript{}{t_{0}}{I}^{\alpha}_{t}g\right)\right)(t)=g(t),\quad\alpha>0,\quad
t\in\mathbb{\hat{J}}.$ (2.5)

The relationship between Riemann-Liouville and Caputo fractional derivatives
are as follows [8]:

$\left(\prescript{C}{t_{0}}{D}^{\alpha}_{t}g\right)(t)=\left(\prescript{RL}{t_{0}}{D}^{\alpha}_{t}g\right)(t)-\sum_{k=0}^{n-1}\frac{(t-t_{0})^{k-\alpha}f^{(k)}(t_{0})}{\Gamma(k-\alpha+1)},\quad
n-1<\alpha\leq n,\quad n\in\mathbb{N}.$ (2.6)

The following results are useful in solving fractional differential equations.

######  Definition 2.3 ([47, 48]).

A function $g$ is said to be exponentially bounded on $[0,\infty)$ if it
satisfies an inequality of the form:

$|g(t)|\leq Me^{\sigma t},\quad t\geq T,$

for some real constants $M>0$, $T>0$ and $\sigma\in\mathbb{R}$.

######  Definition 2.4 ([47, 48]).

If $g:[0,\infty)\to\mathbb{R}$ is exponentially bounded for $t\geq 0$, then
the Laplace integral transform $\mathscr{L}\left\\{g(t)\right\\}(s)$ defined
by

$G(s)=\mathscr{L}\left\\{g(t)\right\\}(s)=\int\limits_{0}^{\infty}e^{-st}g(t)\mathrm{d}t,$

exists for $s\in\mathbb{C}$ and is an analytic function of $s$ for $\Re(s)>0$
and Laplace inversion formula is defined as

$\mathscr{L}^{-1}\left\\{G(s)\right\\}(t)\coloneqq\frac{1}{2\pi
i}\int\limits_{L}e^{st}G(s)ds,$

where $g(t)=\mathscr{L}^{-1}\left\\{G(s)\right\\}(t)$, $t\geq 0$ and $L$ is a
closed contour which enclosing the poles (singularities) of $g$.

######  Definition 2.5 ([8]).

The Laplace integral transform of Riemann-Liouville fractional derivative of
order
$\alpha\in(n-1,n]$, $n\in\mathbb{N}$ is given by [8]:

$\mathscr{L}\left\\{\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-\sum_{k=1}^{n}s^{k-1}\left(\prescript{RL}{0}{D}^{\alpha-k}_{t}y\right)(0),$
(2.7)

where $Y(s)$ represents the Laplace transform of the function $y(t)$ .

######  Remark 2.1 ([8]).

In the special cases, the Laplace integral transform of the Riemann-Liouville
fractional differentiation is:

  * 1.

If $\alpha\in(0,1]$, then

$\mathscr{L}\left\\{\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-\left(\prescript{RL}{0}{D}^{\alpha-1}_{t}y\right)(0).$

  * 2.

If $\alpha\in(1,2]$, then

$\mathscr{L}\left\\{\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-\left(\prescript{RL}{0}{D}^{\alpha-1}_{t}y\right)(0)-s\left(\prescript{RL}{0}{D}^{\alpha-2}_{t}y\right)(0).$

######  Definition 2.6.

The Laplace integral transform of Caputo fractional derivative of order
$\alpha\in(n-1,n]$, $n\in\mathbb{N}$ is given by [8]:

$\mathscr{L}\left\\{\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-\sum_{k=0}^{n-1}s^{\alpha-k-1}\prescript{}{}{y}^{(k)}(0),$
(2.8)

where $Y(s)$ represents the Laplace transform of the function $y(t)$ .

######  Remark 2.2.

In the special cases, the Laplace integral transform of the Caputo fractional
differentiation is:

  * 1.

If $\alpha\in(0,1]$, then

$\mathscr{L}\left\\{\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-s^{\alpha-1}y_{0},\quad\text{where}\quad
y_{0}=y(0).$

  * 2.

If $\alpha\in(1,2]$, then

$\mathscr{L}\left\\{\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)\right\\}(s)=s^{\alpha}Y(s)-s^{\alpha-1}y_{0}-s^{\alpha-2}y^{\prime}_{0},\quad\text{where}\quad
y_{0}=y(0)\quad\text{and}\quad y^{\prime}_{0}=y^{\prime}(0).$

######  Definition 2.7 ([47, 48]).

Let $f$ and $g$ be both piece-wise continuous functions on $[0,\infty)$. Then
the integral in

$f\ast g\coloneqq(f\ast g)(t)=\int\limits_{0}^{t}f(t-s)g(s)ds,$

is called the convolution operator of two functions $f$ and $g$ which is well-
defined and finite for any $t\geq 0$ and it has the commutativity property:

$f\ast g=g\ast f.$

######  Theorem 2.1 ([47, 48]).

Suppose that $f$ and $g$ are piece-wise continuous and exponentially bounded
functions on $[0,\infty)$. Then the Laplace transform of convolution operator
of two functions $f$ and $g$, given on $[0,\infty)$, has the following
property :

$\mathscr{L}\left\\{\left(f\ast
g\right)(t)\right\\}(s)=\mathscr{L}\left\\{f(t)\right\\}(s)\mathscr{L}\left\\{g(t)\right\\}(s),\quad
s\in\mathbb{C}.$

The Mittag-Leffler function is a generalization of the exponential function,
first proposed in 1903 [29] as a single-parameter function of one variable,
defined using a convergent infinite series. Extensions to two, three and
multi-parameters are well known and thoroughly studied in textbooks such as
[26, 31] which are involving single power series in one variable [32, 33, 34].
Extensions to two, three, or more variables, involving correspondingly double,
triple, or multiple power series, have been studied more recently [18, 35, 36,
37].

######  Definition 2.8 ([29]).

The classical Mittag-Leffler function is defined by

$E_{\alpha}(t)=\sum_{i=0}^{\infty}\frac{t^{i}}{\Gamma(i\alpha+1)},\quad\alpha>0,t\in\mathbb{R}.$

######  Remark 2.3 ([25]).

The Mittag-Leffler functions are often used in a form where the variable
inside the brackets is not $t$ but a fractional power $t^{\alpha}$, or even a
constant multiple $\lambda t^{\alpha}$, as follows:

$E_{\alpha}(\lambda
t^{\alpha})=\sum_{i=0}^{\infty}\frac{\lambda^{i}t^{i\alpha}}{\Gamma(i\alpha+1)},\quad\alpha>0,t,\lambda\in\mathbb{R}.$

The two-parameter Mittag-Leffler function [26] is given by

$E_{\alpha,\beta}(t)=\sum_{i=0}^{\infty}\frac{t^{i}}{\Gamma(i\alpha+\beta)},\quad\alpha>0,\beta\in\mathbb{R},t\in\mathbb{R}.$

The $l$-th derivative of two-parameter Mittag-Leffler function [26] is defined
by

$\frac{d^{l}}{dt^{l}}E_{\alpha,\beta}(t)=E^{(l)}_{\alpha,\beta}(t)=\sum_{i=0}^{\infty}\frac{(i+l)!}{i!}\frac{t^{i}}{\Gamma(i\alpha+l\alpha+\beta)},\quad
l\in\mathbb{N},\alpha>0,\beta\in\mathbb{R},t\in\mathbb{R}.$

The three-parameter Mittag-Leffler function [30] is determined by

$E_{\alpha,\beta}^{\gamma}(t)=\sum_{i=0}^{\infty}\frac{(\gamma)_{i}}{\Gamma(i\alpha+\beta)}\frac{t^{i}}{i!},\quad\alpha>0,\beta,\gamma\in\mathbb{R},t\in\mathbb{R},$

where $(\gamma)_{i}$ is the Pochhammer symbol denoting
$\frac{\Gamma(\gamma+i)}{\Gamma(\gamma)}$. These series are convergent,
locally uniformly in $\tau$ , provided the $\alpha>0$ condition is satisfied.
Note that

$E_{\alpha,\beta}^{1}(t)=E_{\alpha,\beta}(t),\quad
E_{\alpha,1}(t)=E_{\alpha}(t),\quad E_{1}(t)=\exp(t).$

The next lemma includes Laplace integral transform of three-parameter Mittag-
Leffler function which will be used throughout the proof of Lemma 2.2.

######  Lemma 2.1.

For $\alpha>\beta>0$, $\lambda\in\mathbb{R}$,
$l\in\mathbb{N}_{0}=\left\\{0,1,2,\dots\right\\}$ and $\Re(s)>0$, we have:

$\displaystyle\mathscr{L}^{-1}\Bigl{\\{}\frac{1}{(s^{\alpha}-\lambda
s^{\beta})^{l+1}}\Bigr{\\}}(\tau)$
$\displaystyle=t^{(l+1)\alpha-1}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{k}t^{k(\alpha-\beta)}}{\Gamma(k(\alpha-\beta)+(l+1)\alpha)}$
$\displaystyle\coloneqq
t^{(l+1)\alpha-1}E^{l+1}_{\alpha-\beta,(l+1)\alpha}(\lambda
t^{\alpha-\beta}).$

###### Proof.

By using the Taylor series representation of
$\frac{1}{(1-t)^{l+1}},l\in\mathbb{N}_{0}$ of the form

$\frac{1}{(1-t)^{l+1}}=\sum_{k=0}^{\infty}\binom{l+k}{k}t^{k},\quad|t|<1,$

we achieve that

$\displaystyle\frac{1}{(s^{\alpha}-\lambda
s^{\beta})^{l+1}}=\frac{1}{(s^{\alpha})^{l+1}}\frac{1}{(1-\frac{\lambda}{s^{\alpha-\beta}})^{l+1}}$
$\displaystyle=\frac{1}{s^{(l+1)\alpha}}\sum_{k=0}^{\infty}\binom{l+k}{k}\Big{(}\frac{\lambda}{s^{\alpha-\beta}}\Big{)}^{k}$
$\displaystyle=\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{k}}{s^{k(\alpha-\beta)+(l+1)\alpha}}.$

Taking inverse Laplace transform of the above function, we get the desired
result:

$\displaystyle\mathscr{L}^{-1}\Bigl{\\{}\frac{1}{(s^{\alpha}-\lambda
s^{\beta})^{l+1}}\Bigr{\\}}(t)$
$\displaystyle=\sum_{k=0}^{\infty}\lambda^{k}\binom{l+k}{k}\mathscr{L}^{-1}\Bigr{\\{}\frac{1}{s^{k(\alpha-\beta)+(l+1)\alpha}}\Bigr{\\}}(t)$
$\displaystyle=\sum_{k=0}^{\infty}\lambda^{k}\binom{l+k}{k}\frac{t^{k(\alpha-\beta)+(l+1)\alpha-1}}{\Gamma(k(\alpha-\beta)+(l+1)\alpha)}$
$\displaystyle=t^{(l+1)\alpha-1}E^{l+1}_{\alpha-\beta,(l+1)\alpha}(\lambda
t^{\alpha-\beta}),$

which is the required result. We have required an extra condition on $s$ for
convergence of the binomial type series in the Laplace domain, namely that

$\displaystyle s^{\alpha-\beta}>|\lambda|.$

However, this condition can be removed at the end, by analytic continuation of
both sides of the identity, to give the desired result for all
$s\in\mathbb{C}$ satisfying $\Re(s)>0$. The proof is complete. ∎

######  Definition 2.9 ([35]).

We consider the bivariate Mittag-Leffler function defined by

$E_{\alpha,\beta,\gamma}^{\delta}(u,v)=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\frac{(\delta)_{l+k}}{\Gamma(l\alpha+k\beta+\gamma)}\frac{u^{l}v^{k}}{l!k!},\quad\alpha,\beta>0,\gamma,\delta\in\mathbb{R},u,v\in\mathbb{R}.$
(2.9)

If we write $u=\lambda t^{\alpha}$ and $v=\mu t^{\beta}$ for a single variable
$t$, and multiply by a power function $t^{\gamma-1}$, we derive the following
univariate version:

$t^{\gamma-1}E_{\alpha,\beta,\gamma}^{\delta}(\lambda t^{\alpha},\mu
t^{\beta})=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\frac{(\delta)_{l+k}}{\Gamma(l\alpha+k\beta+\gamma)}\frac{\lambda^{l}\mu^{k}}{l!k!}t^{l\alpha+k\beta+\gamma-1}.$
(2.10)

Note that when $\delta=1$,

$\displaystyle E_{\alpha,\beta,\gamma}^{1}(\lambda t^{\alpha},\mu t^{\beta})$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\frac{(1)_{l+k}}{\Gamma(l\alpha+k\beta+\gamma)}\frac{\lambda^{l}\mu^{k}}{l!k!}t^{l\alpha+k\beta+\gamma-1}$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\frac{(l+k)!}{l!k!}\frac{\lambda^{l}\mu^{k}}{\Gamma(l\alpha+k\beta+\gamma)}t^{l\alpha+k\beta+\gamma-1}$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}}{\Gamma(l\alpha+k\beta+\gamma)}t^{l\alpha+k\beta+\gamma-1}.$

For simplicity, we denote $E_{\alpha,\beta,\gamma}^{1}(\lambda t^{\alpha},\mu
t^{\beta})\coloneqq E_{\alpha,\beta,\gamma}(\lambda t^{\alpha},\mu t^{\beta})$
in our results for this paper.

######  Lemma 2.2.

For $\alpha>\beta$, $\alpha>\gamma$, $\lambda,\mu\in\mathbb{R}$ and
$\Re(s)>0$, the following result holds true:

$\displaystyle\mathscr{L}^{-1}\Bigl{\\{}\frac{s^{\gamma}}{s^{\alpha}-\mu
s^{\beta}-\lambda}\Bigr{\\}}(t)$
$\displaystyle=t^{\alpha-\gamma-1}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}t^{l\alpha+k(\alpha-\beta)}}{\Gamma(l\alpha+k(\alpha-\beta)+\alpha-\gamma)}$
$\displaystyle=t^{\alpha-\gamma-1}E_{\alpha,\alpha-\beta,\alpha-\gamma}(\lambda
t^{\alpha},\mu t^{\alpha-\beta}).$

###### Proof.

$\frac{s^{\gamma}}{s^{\alpha}-\mu s^{\beta}-\lambda}$ can be written via a
series expansion as follows:

$\frac{s^{\gamma}}{s^{\alpha}-\mu
s^{\beta}-\lambda}=\frac{s^{\gamma}}{s^{\alpha}-\mu
s^{\beta}}\frac{1}{1-\frac{\lambda}{s^{\alpha}-\mu
s^{\beta}}}=\sum_{l=0}^{\infty}\frac{\lambda^{l}s^{\gamma}}{(s^{\alpha}-\mu
s^{\beta})^{l+1}}.$

Then applying Lemma 2.1 to the last expression, we acquire that

$\displaystyle\frac{s^{\gamma}}{s^{\alpha}-\mu s^{\beta}-\lambda}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{\lambda^{l}s^{\gamma}}{s^{(l+1)\alpha}}\frac{1}{(1-\frac{\mu}{s^{\alpha-\beta}})^{l+1}}$
$\displaystyle=\sum_{l=0}^{\infty}\frac{\lambda^{l}s^{\gamma}}{s^{(l+1)\alpha}}\sum_{k=0}^{\infty}\binom{l+k}{k}\Big{(}\frac{\mu}{s^{\alpha-\beta}}\Big{)}^{k}$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}}{s^{(l+1)\alpha+k(\alpha-\beta)-\gamma}}.$

Taking inverse Laplace transform of the aforementioned function, we attain:

$\displaystyle\mathscr{L}^{-1}\Bigl{\\{}\frac{s^{\gamma}}{s^{\alpha}-\mu
s^{\beta}-\lambda}\Bigr{\\}}(t)$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\mu^{l}\lambda^{k}\mathscr{L}^{-1}\Bigl{\\{}\frac{1}{s^{(l+1)\alpha+k(\alpha-\beta)-\gamma}}\Bigr{\\}}(t)$
$\displaystyle=t^{\alpha-\gamma-1}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}t^{l\alpha+k(\alpha-\beta)}}{\Gamma(l\alpha+k(\alpha-\beta)+\alpha-\gamma)}$
$\displaystyle=t^{\alpha-\gamma-1}E_{\alpha,\alpha-\beta,\alpha-\gamma}(\lambda
t^{\alpha},\mu t^{\alpha-\beta}),$

which is the desired result. We have required extra conditions on $s$ for
convergence of the binomial type series in the Laplace domain, namely that

$\displaystyle s^{\alpha-\beta}>|\mu|,$ $\displaystyle|s^{\alpha}-\mu
s^{\beta}|>|\lambda|.$

However, these conditions can be removed at the end, by analytic continuation
of both sides of the identity, to give the desired result for all
$s\in\mathbb{C}$ satisfying $\mathrm{\Re}(s)>0$. The proof is complete. ∎

######  Lemma 2.3.

For any parameters $\alpha,\beta,\gamma,\lambda,\mu\in\mathbb{R}$ satisfying
$\alpha,\beta>0$ and $\gamma-1>\lfloor\alpha\rfloor$, we have

$\prescript{C}{0}{D}^{\alpha}_{t}\Big{[}t^{\gamma-1}E_{\alpha,\beta,\gamma}(\lambda
t^{\alpha},\mu
t^{\beta})\Big{]}=t^{\gamma-\alpha-1}E_{\alpha,\beta,\gamma-\alpha}(\lambda
t^{\alpha},\mu t^{\beta}),\quad t>0.$ (2.11)

###### Proof.

We have the following formula for Caputo derivatives of power functions [8,
27]:

$\prescript{C}{0}{D}^{\nu}_{t}\left(\frac{t^{\eta}}{\Gamma(\eta+1)}\right)=\begin{cases}\frac{t^{\eta-\nu}}{\Gamma(\eta-\nu+1)},\qquad\eta>\lfloor\nu\rfloor,\\\
\qquad 0,\quad\qquad\eta=0,1,2,\dots,\lfloor\nu\rfloor,\\\
\text{undefined},\qquad\text{otherwise}.\end{cases}$ (2.12)

Therefore, the given condition $\gamma-1>\lfloor\alpha\rfloor$, from (2.12) we
can attain

$\displaystyle\prescript{C}{0}{D}^{\alpha}_{t}\left[t^{\gamma-1}E_{\alpha,\beta,\gamma}(\lambda
t^{\alpha},\mu t^{\beta})\right]$
$\displaystyle=\prescript{C}{0}{D}^{\alpha}_{t}\left[\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}t^{l\alpha+k\beta+\gamma-1}}{\Gamma(l\alpha+k\beta+\gamma)}\right]$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\lambda^{l}\mu^{k}\prescript{C}{0}{D}^{\alpha}_{t}\left(\frac{t^{l\alpha+k\beta+\gamma-1}}{\Gamma(l\alpha+k\beta+\gamma)}\right)$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}t^{l\alpha+k\beta+\gamma-\alpha-1}}{\Gamma(l\alpha+k\beta+\gamma-\alpha)}$
$\displaystyle=t^{\gamma-\alpha-1}E_{\alpha,\beta,\gamma-\alpha}(\lambda
t^{\alpha},\mu t^{\beta}),\quad t>0.$

The proof is complete. ∎

######  Lemma 2.4.

For any parameters $\alpha,\beta,\gamma,\lambda,\mu\in\mathbb{R}$ satisfying
$\alpha,\beta,\gamma>0$, we have

$\prescript{RL}{0}{D}^{\alpha}_{t}\Big{[}t^{\gamma-1}E_{\alpha,\beta,\gamma}(\lambda
t^{\alpha},\mu
t^{\beta})\Big{]}=t^{\gamma-\alpha-1}E_{\alpha,\beta,\gamma-\alpha}(\lambda
t^{\alpha},\mu t^{\beta}),\quad t>0.$ (2.13)

###### Proof.

We have the following formula for Riemann-Liouville derivatives of power
functions [8, 27]:

$\prescript{RL}{0}{D}^{\nu}_{t}\left(\frac{t^{\eta}}{\Gamma(\eta+1)}\right)=\frac{t^{\eta-\nu}}{\Gamma(\eta-\nu+1)},\quad\nu,\eta\in\mathbb{R},\quad\eta>-1.$
(2.14)

Therefore, given the condition $\gamma>0$, in accordance with (2.14) we will
get the same result with (2.11). The proof is complete. ∎

######  Definition 2.10.

Let $\lambda_{i},\mu_{j}\in\mathbb{R}$, $\alpha_{i},\beta_{j}\in\mathbb{R}$,
$i=1,2,\ldots,p$, $j=1,2,\ldots,q$. Generalized Wright function or Fox-Wright
function $\prescript{}{p}{\Psi_{q}}(\cdot):\mathbb{R}\to\mathbb{R}$ is defined
by

$\prescript{}{p}{\Psi_{q}}(t)=\prescript{}{p}{\Psi_{q}}\left[\begin{array}[]{ccc}(\lambda_{i},\alpha_{i})_{1,p}\\\
(\mu_{j},\beta_{j})_{1,q}\end{array}\Big{|}t\right]=\sum_{k=0}^{\infty}\frac{\prod\limits_{i=1}^{p}\Gamma(\lambda_{i}+\alpha_{i}k)}{\prod\limits_{j=1}^{q}\Gamma(\mu_{j}+\beta_{j}k)}\frac{t^{k}}{k!}.$
(2.15)

The Fox-Wright function was established by Fox [45] and Wright [46]. If the
following condition holds

$\sum_{j=1}^{q}\beta_{j}-\sum_{i=1}^{p}\alpha_{i}>-1,$

then the series in (2.15) is convergent for arbitrary $t\in\mathbb{R}$.

##  3 Leibniz integral rule

In this section, we formulate Leibniz integral rule for higher order
derivatives on Lebesgue integration. It is known that according to the
suitable conditions, we can differentiate under the integral sign for Lebesgue
integrals [38]. We begin with the first derivative of a Lebesgue integral on
$X\subseteq\mathbb{R}$.

######  Theorem 3.2 ([38]).

Assume that $X,Y\subseteq\mathbb{R}$ are intervals. Suppose also that the
function $f:X\times Y\to\mathbb{R}$ satisfies the following assumptions:

(a) For every fixed $y\in Y$, the function $f(\cdot,y)$ is measurable on $X$;

(b) The partial derivative $\frac{\partial}{\partial y}f(x,y)$ exists for
every interior point $(x,y)\in X\times Y$;

(c) There exists a non-negative integrable function $g$ such that
$\left|\frac{\partial}{\partial y}f(x,y)\right|\leq g(x)$ exists for every
interior point $(x,y)\in X\times Y$;

(d) There exists $y_{0}\in Y$ such that $f(x,y_{0})$ is integrable on $X$.

Then for every $y\in Y$, the Lebesgue integral

$\int\limits_{X}f(x,y)\mathrm{d}x$

exists. Furthermore, the function $F:Y\to\mathbb{R}$, defined by

$F(y)=\int\limits_{X}f(x,y)\mathrm{d}x$

for every $y\in Y,$ is differentiable at every interior point of $Y$, and the
derivative of $F(y)$ satisfies

$F^{\prime}(y)=\int\limits_{X}\frac{\partial}{\partial y}f(x,y)\mathrm{d}x.$

The well-known rule for the differentiation of an integral depending on a
parameter with the upper limit also depends on the same parameter, namely:

######  Corollary 3.1 ([39]).

If $X=(y_{0},y)$ and assumptions of Theorem 3.2 are fulfilled, then the
following relation holds true for all $y\in Y$:

$\frac{d}{dy}\int\limits_{y_{0}}^{y}f(x,y)\mathrm{d}x=\int\limits_{y_{0}}^{y}\frac{\partial}{\partial
y}f(x,y)\mathrm{d}x+\lim\limits_{x\to y-0}f(x,y),\quad y\in X.$ (3.1)

So, the formula of differentiation under the integral sign for $K(t,s)$ with
respect to $t$ is

$\frac{d}{dt}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\int\limits_{t_{0}}^{t}\frac{\partial}{\partial
t}K(t,s)\mathrm{d}s+\lim\limits_{s\to t-0}K(t,s),\quad t\in\mathbb{\hat{J}}.$
(3.2)

Using the formula (3.2), we define the second-order derivative of the integral
depending on $t$:

$\displaystyle\frac{d^{2}}{dt^{2}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\frac{d}{dt}\left(\frac{d}{dt}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s\right)=\frac{d}{dt}\left(\lim\limits_{s\to
t-0}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial}{\partial
t}K(t,s)\mathrm{d}s\right)$ $\displaystyle=\frac{d}{dt}\lim\limits_{s\to
t-0}K(t,s)+\frac{d}{dt}\int\limits_{t_{0}}^{t}\frac{\partial}{\partial
t}K(t,s)\mathrm{d}s$ $\displaystyle=\frac{d}{dt}\lim\limits_{s\to
t-0}K(t,s)+\lim\limits_{s\to t-0}\frac{\partial}{\partial
t}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{2}}{\partial
t^{2}}K(t,s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{2}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-{0}}\frac{\partial^{2-l}}{\partial
t^{2-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{2}}{\partial
t^{2}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$

Then, the third-order differentiation of the integral will be:

$\displaystyle\frac{d^{3}}{dt^{3}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\frac{d}{dt}\left(\frac{d^{2}}{dt^{2}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s\right)$
$\displaystyle=$ $\displaystyle\frac{d}{dt}\left(\frac{d}{dt}\lim\limits_{s\to
t-0}K(t,s)+\lim\limits_{s\to t-0}\frac{\partial}{\partial
t}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{2}}{\partial
t^{2}}K(t,s)\mathrm{d}s\right)$ $\displaystyle=$
$\displaystyle\frac{d^{2}}{dt^{2}}\lim\limits_{s\to
t-0}K(t,s)+\frac{d}{dt}\lim\limits_{s\to t-0}\frac{\partial}{\partial
t}K(t,s)+\frac{d}{dt}\int\limits_{t_{0}}^{t}\frac{\partial^{2}}{\partial
t^{2}}K(t,s)\mathrm{d}s$ $\displaystyle=$
$\displaystyle\frac{d^{2}}{dt^{2}}\lim\limits_{s\to
t-0}K(t,s)+\frac{d}{dt}\lim\limits_{s\to t-0}\frac{\partial}{\partial
t}K(t,s)$ $\displaystyle+$ $\displaystyle\lim\limits_{s\to
t-0}\frac{\partial^{2}}{\partial
t^{2}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{3}}{\partial
t^{3}}K(t,s)\mathrm{d}s$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{3}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-{0}}\frac{\partial^{3-l}}{\partial
t^{3-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{3}}{\partial
t^{3}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$

Thus, we establish $n$th derivative of the integral for $n\in\mathbb{N}$ which
depends on $t$ by recursively as follows:

$\displaystyle\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\frac{d}{dt}\left(\frac{d^{n-1}}{dt^{n-1}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s\right)$
$\displaystyle=$ $\displaystyle\frac{d^{n-1}}{dt^{n-1}}\lim\limits_{s\to
t-0}K(t,s)+\frac{d^{n-2}}{dt^{n-2}}\lim\limits_{s\to
t-0}\frac{\partial}{\partial t}K(t,s)+\cdots+\frac{d}{dt}\lim\limits_{s\to
t-0}\frac{\partial^{n-2}}{\partial t^{n-2}}K(t,s)$ $\displaystyle+$
$\displaystyle\lim\limits_{s\to t-0}\frac{\partial^{n-1}}{\partial
t^{n-1}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{n}}{\partial
t^{n}}K(t,s)\mathrm{d}s$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{n}}{\partial
t^{n}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$

In the next theorem, we state and prove the Leibniz rule for higher order
derivatives.

######  Theorem 3.3.

Let the function $K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$ be such that the
following assumptions are fulfilled:

(a) For every fixed $t\in\mathbb{J}$, the function
$\frac{\partial^{n-1}}{\partial t^{n-1}}K(t,s)$ is measurable on $\mathbb{J}$
and integrable on $\mathbb{J}$ with respect to for some $t^{*}\in\mathbb{J}$;

(b) The partial derivative $\frac{\partial^{n}}{\partial t^{n}}K(t,s)$ exists
for every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(c) There exists a non-negative integrable function $g$ such that
$\left|\frac{\partial^{n}}{\partial t^{n}}K(t,s)\right|\leq g(s)$ for every
interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(d) The derivative $\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-{0}}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)$, $l=1,2,\ldots,n$ exists
for every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$.

Then, the following relation holds true for $n$th derivative under Lebesgue
integration for $n\in\mathbb{N}$:

$\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-{0}}\frac{\partial^{n-l}}{\partial
t^{n-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{n}}{\partial
t^{n}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$ (3.3)

###### Proof.

Using mathematical induction principle, we prove above theorem. It is obvious
that the equation (3.3) is true for $n=1$ [8]:

$\frac{d}{dt}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\lim\limits_{s\to
t-0}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial}{\partial
t}K(t,s)\mathrm{d}s.$

We assume that (3.3) holds true for $n=k$:

$\frac{d^{k}}{dt^{k}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\sum_{l=1}^{k}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{k-l}}{\partial
t^{k-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{k}}{\partial
t^{k}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$

We prove that (3.3) is true for $n=k+1$:

$\displaystyle\frac{d^{k+1}}{dt^{k+1}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\frac{d}{dt}\left(\frac{d^{k}}{dt^{k}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s\right)$
$\displaystyle=\frac{d}{dt}\left(\sum_{l=1}^{k}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{k-l}}{\partial
t^{k-l}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{k}}{\partial
t^{k}}K(t,s)\mathrm{d}s\right)$
$\displaystyle=\frac{d}{dt}\Big{(}\lim\limits_{s\to
t-0}\frac{\partial^{k-1}}{\partial
t^{k-1}}K(t,s)+\frac{d}{dt}\lim\limits_{s\to
t-0}\frac{\partial^{k-2}}{\partial
t^{k-2}}K(t,s)+\cdots+\frac{d^{k-1}}{dt^{k-1}}\lim\limits_{s\to
t-0}K(t,s)\Big{)}$
$\displaystyle+\frac{d}{dt}\int\limits_{t_{0}}^{t}\frac{\partial^{k}}{\partial
t^{k}}K(t,s)\mathrm{d}s$ $\displaystyle=\frac{d}{dt}\lim\limits_{s\to
t-0}\frac{\partial^{k-1}}{\partial
t^{k-1}}K(t,s)+\frac{d^{2}}{dt^{2}}\lim\limits_{s\to
t-0}\frac{\partial^{k-2}}{\partial
t^{k-2}}K(t,s)+\cdots+\frac{d^{k}}{dt^{k}}\lim\limits_{s\to t-0}K(t,s)$
$\displaystyle+\lim\limits_{s\to t-0}\frac{\partial^{k}}{\partial
t^{k}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{k+1}}{\partial
t^{k+1}}K(t,s)\mathrm{d}s$
$\displaystyle=\sum_{l=2}^{k+1}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{k-l+1}}{\partial t^{k-l+1}}K(t,s)+\lim\limits_{s\to
t-0}\frac{\partial^{k}}{\partial
t^{k}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{k+1}}{\partial
t^{k+1}}K(t,s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{k+1}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{k-l+1}}{\partial
t^{k-l+1}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{k+1}}{\partial
t^{k+1}}K(t,s)\mathrm{d}s,\quad t\in\mathbb{\hat{J}}.$

Therefore, the formula (3.3) holds true for all $n\in\mathbb{N}$ and
$t\in\mathbb{\hat{J}}$. ∎

The following important particular case must be defined for convolution
operator of the functions $f$ and $g$.

######  Corollary 3.2.

If $K(t,s)=f(t-s)g(s)$ and $t_{0}=0$, and assumptions of Theorem 3.3 are
satisfied, then the following relation is true for any $n\in\mathbb{N}$:

$\displaystyle\frac{d^{n}}{dt^{n}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\frac{\partial^{n}}{\partial
t^{n}}f(t-s)g(s)\mathrm{d}s,\quad t>0.$ (3.4)

###### Proof.

If we write $f(t-s)g(s)$ instead of $K(t,s)$ in (3.3), then we obtain

$\displaystyle\frac{d^{n}}{dt^{n}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\left(\lim\limits_{s\to
t-0}\left(\frac{\partial^{n-l}}{\partial t^{n-l}}f(t-s)g(s)\right)\right)$
$\displaystyle+\int\limits_{0}^{t}\frac{\partial^{n}}{\partial
t^{n}}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\left(\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}f(t-s)\lim\limits_{s\to
t-0}g(s)\right)$
$\displaystyle+\int\limits_{0}^{t}\frac{\partial^{n}}{\partial
t^{n}}f(t-s)g(s)\mathrm{d}s$ $\displaystyle=\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\frac{\partial^{n}}{\partial
t^{n}}f(t-s)g(s)\mathrm{d}s,\quad t>0.$

Thus, the proof is complete. ∎

##  4  Fractional Leibniz integral rules

Now, we are starting to prove fractional Leibniz integral rule for Riemann-
Liouville fractional derivative of order $\alpha\in(n-1,n],n\in\mathbb{N}$.
For this, firstly, let us consider partial Riemann-Liouville fractional
differentiation operator of order $n-1<\alpha\leq n$, $n\in\mathbb{N}$ with
respect to $t$ of a function $K(t,s)$ of two variables
$(t,s)\in\mathbb{J}\times\mathbb{J}$,
$K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$, defined by [28, 42, 44]:

$\prescript{RL,t}{t_{0}}{D}^{\alpha}_{t}K(t,s)=\frac{\mathrm{\partial}^{n}}{\mathrm{\partial}t^{n}}\prescript{t}{t_{0}}{I}^{n-\alpha}_{t}K(t,s)=\frac{1}{\Gamma(n-\alpha)}\frac{\mathrm{\partial}^{n}}{\mathrm{\partial}t^{n}}\int\limits^{t}_{t_{0}}(t-s)^{n-\alpha-1}K(s,\tau)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}},$ (4.1)

where $\prescript{t}{t_{0}}{I}^{\alpha}_{t}$ is the partial Riemann-Liouville
integral operator of order $\alpha>0$ which is given by:

$\prescript{t}{t_{0}}{I}^{\alpha}_{t}K(t,s)=\frac{1}{\Gamma(\alpha)}\int\limits_{t_{0}}^{t}(t-s)^{\alpha-1}K(s,\tau)\mathrm{d}s\\\
,\quad\text{for}\quad t\in\mathbb{\hat{J}}.$

The following important result in the theory of fractional calculus was first
proposed by Podlubny [8] for $\alpha\in(0,1]$ in Riemann-Liouville sense as
follows:

$\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\lim\limits_{s\to
t-0}\prescript{t}{s}{I^{1-\alpha}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t>t_{0}.$

Now, we are going to state and prove the following theorem for more general
case where $\alpha\in(n-1,n],n\in\mathbb{N}$ which is more useful tool for the
testing particular solution of inhomogeneous linear multi-order differential
equations with variable coefficients. Note that Matychyn has proposed [41]
Leibniz integral rule for Riemann-Liouville derivative of order $0<\alpha\leq
1$ on Lebesgue integration.

######  Theorem 4.4.

Let the function $K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$ be such that the
following assumptions are fulfilled:

(a) For every fixed $t\in\mathbb{J}$, the function
$\hat{K}(t,s)=\prescript{RL,t}{s}{D^{\alpha-1}_{t}}K(t,s)$ is measurable on
$\mathbb{J}$ and integrable on $\mathbb{J}$ with respect to some
$t^{*}\in\mathbb{J}$;

(b) The partial derivative $\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)$ exists
for every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(c) There exists a non-negative integrable function $g$ such that
$\left|\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\right|\leq g(s)$ for every
interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(d) The derivative $\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\prescript{RL,t}{s}{D^{\alpha-l}_{t}}K(t,s)$, $l=1,2,\ldots,n$ exists for
every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

Then, the following relation holds true for fractional derivative in Riemann-
Liouville sense under Lebesgue integration:

$\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\prescript{RL,t}{s}{D^{\alpha-l}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$ (4.2)

###### Proof.

Using the Definition 2.1 and Fubini’s theorem [38], we have

$\displaystyle\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}(t-\tau)^{n-\alpha-1}\mathrm{d}\tau\int\limits_{t_{0}}^{\tau}K(\tau,s)\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}\int\limits_{t_{0}}^{\tau}(t-\tau)^{n-\alpha-1}K(\tau,s)\mathrm{d}s\mathrm{d}\tau$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}K(\tau,s)\mathrm{d}\tau\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}\mathrm{d}s\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}K(\tau,s)\mathrm{d}\tau$
$\displaystyle=\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}\left(\frac{1}{\Gamma(n-\alpha)}\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}K(\tau,s)\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=\frac{d^{n}}{dt^{n}}\int\limits_{t_{0}}^{t}\prescript{t}{s}{I^{n-\alpha}_{t}}K(t,s)\mathrm{d}s.$

Using the formula (3.3) for the last part of above expression, we get a
desired result:

$\displaystyle\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}\prescript{t}{s}{I^{n-\alpha}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\frac{\partial^{n}}{\partial
t^{n}}\prescript{t}{s}{I^{n-\alpha}_{t}}K(t,s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\prescript{RL,t}{s}{D^{\alpha-l}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$

∎

######  Corollary 4.1.

If we have $K(t,s)=f(t-s)g(s)$, $t_{0}=0$ and assumptions of Theorem 4.4 are
fulfilled, then following equality holds true for convolution operator in
Riemann-Liouville sense for any $n\in\mathbb{N}$:

$\displaystyle\prescript{RL}{0}{D^{\alpha}_{t}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\prescript{RL,t}{s}{D^{\alpha-l}_{t}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s,\quad
t>0.$ (4.3)

###### Proof.

If we write $f(t-s)g(s)$ instead of $K(t,s)$ in (4.2), then we obtain

$\displaystyle\prescript{RL}{0}{D^{\alpha}_{t}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\left(\lim\limits_{s\to
t-0}\left(\frac{\partial^{n-l}}{\partial
t^{n-l}}\prescript{t}{s}{I^{n-\alpha}_{t}}f(t-s)g(s)\right)\right)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}\prescript{t}{s}{I^{n-\alpha}_{t}}f(t-s)\lim\limits_{s\to t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\prescript{RL,t}{s}{D^{\alpha-l}_{t}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s,\quad
t>0.$

Thus, the proof is complete. ∎

Then, we are going to introduce fractional differentiation under the integral
sign in Caputo sense which will be useful for checking the candidate solutions
of fractional differential equations with multi-orders. For this, firstly, let
us consider partial Caputo fractional differentiation operator of order
$n-1<\alpha\leq n$, $n\in\mathbb{N}$ with respect to $t$ of a function
$K(t,s)$ of two variables $(t,s)\in\mathbb{J}\times\mathbb{J}$,
$K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$, defined by [28, 42, 44]:

$\prescript{C,t}{t_{0}}{D}^{\alpha}_{t}K(t,s)=\prescript{t}{t_{0}}{I}^{n-\alpha}_{t}\frac{\mathrm{\partial}^{n}}{\mathrm{\partial}t^{n}}K(t,s)=\frac{1}{\Gamma(n-\alpha)}\int\limits^{t}_{t_{0}}(t-s)^{n-\alpha-1}\frac{\mathrm{\partial}^{n}}{\mathrm{\partial}s^{n}}K(s,\tau)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}},$ (4.4)

Matychyn and Onyshchenko [40] showed that the fractional Leibniz integral rule
for Caputo fractional derivative coincide with Riemann-Liouville one when
$\alpha\in(0,1]$:

$\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\lim\limits_{s\to
t-0}\prescript{t}{s}{I^{1-\alpha}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$

More generally, the fractional Leibniz integral rule for fractional derivative
of order $\alpha\in(n-1,n]$, $n\geq 2$ in Caputo sense is stated and proved in
the following theorem.

######  Theorem 4.5.

Let the function $K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$ be such that the
following assumptions are fulfilled.

(a) For every fixed $t\in\mathbb{J}$, the function
$\hat{K}(t,s)=\prescript{C,t}{s}{D^{\alpha-1}_{t}}K(t,s)$ is measurable and
integrable on $\mathbb{J}$ with respect to some $t^{*}\in\mathbb{J}$;

(b) The partial derivative $\prescript{C,t}{s}{D^{\alpha}_{t}}K(t,s)$ exists
for every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(c) There exists a non-negative integrable function $g$ such that
$\left|\prescript{C,t}{s}{D^{\alpha}_{t}}K(t,s)\right|\leq g(s)$ for every
interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(d) The integral
$\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$ ,
$l=1,2,\ldots,n$, $n\in\mathbb{N}$ exists for every interior point
$(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

Then, the following relation holds true for fractional derivative in Caputo
sense under Lebesgue integration:

$\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}K(t,s)\right\\}+\int\limits_{t_{0}}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$ (4.5)

###### Proof.

Using the Definition 2.2, Fubini’s theorem [38], and the formula (3.3) we have

$\displaystyle\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}(t-\tau)^{n-\alpha-1}\mathrm{d}\tau\frac{d^{n}}{d\tau^{n}}\int\limits_{t_{0}}^{\tau}K(\tau,s)\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}(t-\tau)^{n-\alpha-1}\sum_{l=1}^{n}\frac{d^{l-1}}{d\tau^{l-1}}\lim\limits_{s\to{\tau-0}}\frac{\partial^{n-l}}{\partial\tau^{n-l}}K(\tau,s)\mathrm{d}\tau$
$\displaystyle+\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}(t-\tau)^{n-\alpha-1}\mathrm{d}\tau\int\limits_{t_{0}}^{\tau}\frac{\partial^{n}}{\partial\tau^{n}}K(\tau,s)\mathrm{d}s$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}\int\limits_{t_{0}}^{\tau}(t-\tau)^{n-\alpha-1}\frac{\partial^{n}}{\partial\tau^{n}}K(\tau,s)\mathrm{d}s\mathrm{d}\tau$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}\frac{\partial^{n}}{\partial\tau^{n}}K(\tau,s)\mathrm{d}\tau\mathrm{d}s$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\frac{1}{\Gamma(n-\alpha)}\int\limits_{t_{0}}^{t}\mathrm{d}s\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}\frac{\partial^{n}}{\partial\tau^{n}}K(\tau,s)\mathrm{d}\tau$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\int\limits_{t_{0}}^{t}\left(\frac{1}{\Gamma(n-\alpha)}\int\limits_{s}^{t}(t-\tau)^{n-\alpha-1}\frac{\partial^{n}}{\partial\tau^{n}}K(\tau,s)\mathrm{d}\tau\right)\mathrm{d}s$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\int\limits_{t_{0}}^{t}\prescript{t}{s}{I^{n-\alpha}_{t}}\frac{\partial^{n}}{\partial
t^{n}}K(t,s)\mathrm{d}s$
$\displaystyle=\prescript{t}{t_{0}}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}K(t,s)\right\\}$
$\displaystyle+\int\limits_{t_{0}}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$

Therefore, the proof is compete. ∎

However, the fractional Leibniz integral rule for Caputo derivative of order
$0<\alpha\leq 1$ is different from the general case which is given in the
relation (4.5).

######  Theorem 4.6.

Let the function $K:\mathbb{J}\times\mathbb{J}\to\mathbb{R}$ be such that the
following assumptions are fulfilled.

(a) For every fixed $t\in\mathbb{J}$, the function
$\hat{K}(t,s)=\prescript{t}{s}{I^{1-\alpha}_{t}}K(t,s)$ is measurable on
$\mathbb{J}$ and integrable on $\mathbb{J}$ with respect to some
$t^{*}\in\mathbb{J}$;

(b) The partial derivative $\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)$ exists
for every interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

(c) There exists a non-negative integrable function $g$ such that
$\left|\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\right|\leq g(s)$ for every
interior point $(t,s)\in\mathbb{\hat{J}}\times\mathbb{\hat{J}}$;

Then, the Caputo fractional derivative under Lebesgue integration coincides
with the fractional differentiation of an integral in Riemann-Liouville sense
for $0<\alpha\leq 1$:

$\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\lim\limits_{s\to
t-0}\prescript{t}{s}{I^{1-\alpha}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$ (4.6)

###### Proof.

In accordance the formula (2.6), it is obvious for $\alpha\in(0,1]$:

$\left(\prescript{C}{t_{0}}{D^{\alpha}_{t}}g\right)(t)=\left(\prescript{RL}{t_{0}}{D^{\alpha}_{t}}g\right)(t)-\frac{(t-t_{0})^{-\alpha}}{\Gamma(1-\alpha)}g(t_{0}),\quad
t\in\mathbb{\hat{J}}.$ (4.7)

Since the formula (4.7), we can attain that

$\displaystyle\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s-\left[\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s\right]_{t=t_{0}}\times\frac{(t-t_{0})^{-\alpha}}{\Gamma(1-\alpha)}$
$\displaystyle=\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$

Therefore, fractional Leibniz integral rule for Caputo derivative is identical
with the Riemann-Liouville one whenever $0<\alpha\leq 1$:

$\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int\limits_{t_{0}}^{t}K(t,s)\mathrm{d}s=\lim\limits_{s\to
t-0}\prescript{t}{s}{I^{1-\alpha}_{t}}K(t,s)+\int\limits_{t_{0}}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}}.$

∎

It is important to introduce Caputo fractional derivative of convolution
operator in general sense which is so accurate tool for testing particular
solution of Caputo type multi-term FDEs.

######  Corollary 4.2.

If we have $K(t,s)=f(t-s)g(s)$, $t_{0}=0$, and assumptions of Theorem 4.6 are
fulfilled, then following equality holds true for convolution operator in
Caputo’s sense of order $\alpha\in(n-1,n]$, where $n\geq 2$:

$\prescript{C}{0}{D^{\alpha}_{t}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s=\prescript{t}{0}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)\right\\}+\int\limits_{0}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}f(\tau-s)g(s)\mathrm{d}s,\quad
t>0.$ (4.8)

###### Proof.

If we write $f(t-s)g(s)$ instead of $K(t,s)$ in (4.8), then we acquire

$\displaystyle\prescript{C}{0}{D^{\alpha}_{t}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}f(t-s)g(s)\right\\}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\frac{d^{l-1}}{dt^{l-1}}\left(\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial t^{n-l}}f(t-s)\lim\limits_{s\to
t-0}g(s)\right)\right\\}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{n-\alpha}_{t}}\left\\{\sum_{l=1}^{n}\lim\limits_{s\to
t-0}\frac{\partial^{n-l}}{\partial
t^{n-l}}f(t-s)\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to t-0}g(s)\right\\}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s,\quad
t>0.$

Thus, the proof is complete. ∎

######  Corollary 4.3.

If we have $K(t,s)=f(t-s)g(s)$, $t_{0}=0$, and assumptions of Theorem 4.6 are
fulfilled, then following equality holds true for convolution operator in
Caputo’s sense for $\alpha\in(0,1]$:

$\prescript{C}{0}{D^{\alpha}_{t}}\int\limits_{0}^{t}f(t-s)g(s)\mathrm{d}s=\lim\limits_{s\to
t-0}\prescript{t}{s}{I^{1-\alpha}_{t}}f(t-s)\lim\limits_{s\to
t-0}g(s)+\int\limits_{0}^{t}\prescript{RL,t}{s}{D^{\alpha}_{t}}f(t-s)g(s)\mathrm{d}s,\quad
t>0.$ (4.9)

###### Proof.

If we make use of the substitution $K(t,s)=f(t-s)g(s)$ in the relation (4.6),
the proof is straightforward. So, we omit it here. ∎

######  Theorem 4.7.

The relationship between Leibniz integral rule for Riemann-Liouville and
Caputo fractional differentiation operators of order $n-1<\alpha\leq n$,
$n\geq 2$ holds true:

$\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s=\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s+\sum_{i=1}^{n-1}\sum_{l=1}^{i}\left[\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{i-l}}{\partial
t^{i-l}}K(t,s)\right]_{t=t_{0}}\times\frac{(t-t_{0})^{i-\alpha}}{\Gamma(i-\alpha+1)},t\in\mathbb{\hat{J}}.$
(4.10)

###### Proof.

Using the relationship between Riemann-Liouville and Caputo fractional
derivatives (2.6), we get

$\displaystyle\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s$
$\displaystyle=\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s+\sum_{i=0}^{n-1}\frac{(t-t_{0})^{i-\alpha}}{\Gamma(i-\alpha+1)}\Big{[}\frac{d^{i}}{dt^{i}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s\Big{]}_{t=t_{0}}$
$\displaystyle=\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s+\sum_{i=1}^{n-1}\frac{(t-t_{0})^{i-\alpha}}{\Gamma(i-\alpha+1)}\Big{[}\sum_{l=1}^{i}\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{i-l}}{\partial t^{i-l}}K(t,s)\Big{]}_{t=t_{0}}$
$\displaystyle+\frac{(t-t_{0})^{-\alpha}}{\Gamma(1-\alpha)}\left[\int_{t_{0}}^{t}K(t,s)\mathrm{d}s\right]_{t=t_{0}}+\sum_{i=0}^{n-1}\frac{(t-t_{0})^{i-\alpha}}{\Gamma(i-\alpha+1)}\Big{[}\int_{t_{0}}^{t}\frac{\partial^{i}}{\partial
t^{i}}K(t,s)\mathrm{d}s\Big{]}_{t=t_{0}}$
$\displaystyle=\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s+\sum_{i=1}^{n-1}\sum_{l=1}^{i}\left[\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}\frac{\partial^{i-l}}{\partial
t^{i-l}}K(t,s)\right]_{t=t_{0}}\times\frac{(t-t_{0})^{i-\alpha}}{\Gamma(i-\alpha+1)},t\in\mathbb{\hat{J}}.$
(4.11)

∎

######  Corollary 4.4.

If we replace $K(t,s)$ with $f(t-s)g(s)$ and consider $t_{0}=0$ for lower
bound of the integral in Theorem 4.7, the relationship between Riemann-
Liouville and Caputo type Leibniz integral rules for convolution operator of
the functions $f$ and $g$ holds true for $n-1<\alpha\leq n,n\geq 2$:

$\displaystyle\prescript{RL}{0}{D^{\alpha}_{t}}\int_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle=\prescript{C}{0}{D^{\alpha}_{t}}\int_{0}^{t}f(t-s)g(s)\mathrm{d}s$
$\displaystyle+\sum_{i=1}^{n-1}\sum_{l=1}^{i}\lim\limits_{s\to
t-0}\frac{\partial^{i-l}}{\partial
t^{i-l}}f(t-s)\left[\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)\right]_{t=0}\frac{t^{-\alpha}}{\Gamma(1-\alpha)},\quad t>0.$ (4.12)

######  Corollary 4.5.

The fractional Leibniz rule for Riemann-Liouville and Caputo type fractional
differential operators coincides for $0<\alpha\leq 1$:

$\prescript{RL}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s=\prescript{C}{t_{0}}{D^{\alpha}_{t}}\int_{t_{0}}^{t}K(t,s)\mathrm{d}s,\quad
t\in\mathbb{\hat{J}},$ (4.13)
$\prescript{RL}{0}{D^{\alpha}_{t}}\int_{0}^{t}f(t-s)g(s)\mathrm{d}s=\prescript{C}{0}{D^{\alpha}_{t}}\int_{0}^{t}f(t-s)g(s)\mathrm{d}s,\quad
t>0.$ (4.14)

##  5 Fractional Green’s function method

The Laplace transform is a convenient technique for solving the Cauchy problem
associated with multi-term FDEs with constant coefficients. For instance, let
us consider linear in-homogeneous FDE with multi-orders in Caputo’s sense and
constant coefficients:

$\left\\{\prescript{C}{0}{D}^{\alpha_{n}}_{t}+\lambda_{1}\prescript{C}{0}{D}^{\alpha_{n-1}}_{t}+\lambda_{2}\prescript{C}{0}{D}^{\alpha_{n-2}}_{t}+\ldots+\lambda_{n-1}\prescript{C}{0}{D}^{\alpha_{1}}_{t}+\lambda_{n}\right\\}y(t)=g(t),\quad
t>0,$ (5.1)

under the homogeneous initial conditions:

$y^{(k)}(0)=0,\quad k=0,1,\ldots,n-1,$ (5.2)

where $\prescript{C}{0}{D}^{\alpha_{i}}_{(\cdot)}y(\cdot),i=1,2,\ldots,n$, are
the Caputo fractional differentiation operators of orders
$i-1\leq\alpha_{i}\leq i$, $\lambda_{i}\in\mathbb{R}$ for $i=1,2,\ldots,n$
denote constants and $g\in C([0,\infty),\mathbb{R})$ is the continuous force
or input function.

The classical analogue of the same problem is considered by Miller in [24].
Let us consider IVP for $n$-th order linear differential equation with
constant coefficients:

$\left\\{D^{n}+\lambda_{1}D^{n-1}+\lambda_{2}D^{n-2}+\ldots+\lambda_{n-1}D+\lambda_{n}\right\\}y(t)=0,\quad
t>0,$ (5.3)

with zero initial conditions

$D^{k}y(0)=0,\quad 0\leq k\leq n-1.$ (5.4)

The fractional Green’s function is a very useful and applicable practical as
well as theoretical tool for solving the IVP (5.1)-(5.2) for multi-order FDE.

If we let

$P(x)=x^{n}+\lambda_{1}x^{n-1}+\lambda_{2}x^{n-2}+\ldots\lambda_{n-1}x+\lambda_{n}$
(5.5)

be a polynomial which is related to the equation (5.3), then according to the
Laplace transform method, the unique solution of the following differential
system:

$\begin{cases}P(D)y(t)=g(t),\quad t>0,\\\ D^{k}y(0)=0,\quad
k=0,1,\ldots,n-1,\end{cases}$

can be represented in terms of a convolution integral

$y(t)=\int\limits_{0}^{t}H(t-s)g(s)\mathrm{d}s,\quad t>0,$ (5.6)

where $H(\cdot)$ is the Green or weight function associated with the
differential operator $P(D)$ that is evaluated by taking inverse Laplace
transform of the transfer function.

###  5.1 Applications of fractional Leibniz rules

In this subsection, we study applications of the fractional Leibniz integral
rule in Riemann-Liouville and Caputo sense using the generalized Bagley-Torvik
equations. Moreover, we have used Leibniz integral rule for checking candidate
solution of the oscillator equation in classical sense.

In the following cases, to obtain analytical representation of solutions for
the Cauchy problem we will apply fractional Green’s function method as we
mentioned in Section 5.

Case 1: We consider the IVP for generalized Bagley-Torvik equations with
Riemann-Liouville fractional derivatives of order $1<\alpha\leq 2$ and
$0<\beta\leq 1$ in the form of:

$\displaystyle\begin{cases}\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)-\mu\left(\prescript{RL}{0}{D}^{\beta}_{t}y\right)(t)-\lambda
y(t)=g(t),\quad t>0,\\\
\prescript{}{0}{I^{2-\alpha}_{t}}y(t)|_{t=0}=\prescript{}{0}{I^{1-\alpha}_{t}}y(t)|_{t=0}=0,\quad\quad\lambda,\mu\in\mathbb{R}.\end{cases}$
(5.7)

######  Theorem 5.8.

A unique solution $y\in C^{2}([0,\infty),\mathbb{R})$ of the Cauchy problem
(5.7) has the following formula:

$y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.8)

###### Proof.

We assume that (5.7) has a unique solution $y(t)$ and $g(t)$ is continuous on
$\left[0,\infty\right)$ and exponentially bounded, then $y(t)$,
$\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)$, and
$\left(\prescript{RL}{0}{D}^{\beta}_{t}y\right)(t)$ are exponentially bounded,
thus their Laplace transform exist.

Applying Laplace integral transform for Riemann-Liouville fractional
derivative using the formula (2.7) to the both sides of (5.7) yields:

$\displaystyle\left(s^{\alpha}-\mu s^{\beta}-\lambda\right)Y(s)=G(s).$ (5.9)

Then we solve (5.9) with respect to $Y(s)$,

$\displaystyle Y(s)=\frac{G(s)}{s^{\alpha}-\mu s^{\beta}-\lambda}.$ (5.10)

Taking inverse Laplace transform of (5.10) and applying Lemma 2.2, we find an
explicit representation of solution to (5.7):

$\displaystyle
y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.11)

∎

Verification by substitution. Having found explicit form for $y(t)$, it
remains to confirm that $y(t)$ is an analytical solution of (5.7) indeed.
Firstly, for make the use of checking by substitution, we apply fractional
Leibniz integral rule in Riemann-Liouville sense for the first and second
terms of (5.7). Then the first term will be as follows:

$\displaystyle\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)$
$\displaystyle=\prescript{RL}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{2}\lim\limits_{s\to
t-0}\prescript{RL,t}{0}{D}^{\alpha-l}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to
t-0}\prescript{RL,t}{0}{D}^{\alpha-1}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$ $\displaystyle+\lim\limits_{s\to
t-0}\prescript{RL,t}{0}{D}^{\alpha-2}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d}{dt}\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to
t-0}E_{\alpha,\alpha-\beta,1}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$ $\displaystyle+\lim\limits_{s\to
t-0}(t-s)E_{\alpha,\alpha-\beta,2}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d}{dt}\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$

From now on, we apply Pascal’s rule for binomial coefficients to the first
term of above expression and the limit of the second term is equal to zero as
$s\to t-0$, we obtain

$\displaystyle\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)$
$\displaystyle=\prescript{RL}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=g(t)+\lambda\lim\limits_{s\to
t-0}(t-s)^{\alpha}E_{\alpha,\alpha-\beta,\alpha+1}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$ $\displaystyle+\mu\lim\limits_{s\to
t-0}(t-s)^{\alpha-\beta}E_{\alpha,\alpha-\beta,\alpha-\beta+1}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=g(t)+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$

Now, using Lemma 2.4 and again applying Pascal’s rule for the expression under
above integral, we attain

$\displaystyle\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\alpha}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{t}(t-s)^{-1}E_{\alpha,\alpha-\beta,0}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}(t-s)^{l\alpha+k(\alpha-\beta)-1}}{\Gamma(l\alpha+k(\alpha-\beta))}g(s)\mathrm{d}s$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{t}\frac{(t-s)^{-1}}{\Gamma(0)}g(s)\mathrm{d}s+\int\limits_{0}^{t}\sum_{l=1}^{\infty}\sum_{k=0}^{\infty}\binom{l+k-1}{k}\frac{\lambda^{l}\mu^{k}(t-s)^{l\alpha+k(\alpha-\beta)-1}}{\Gamma(l\alpha+k(\alpha-\beta))}g(s)\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=1}^{\infty}\binom{l+k-1}{k-1}\frac{\lambda^{l}\mu^{k}(t-s)^{l\alpha+k(\alpha-\beta)-1}}{\Gamma(l\alpha+k(\alpha-\beta))}g(s)\mathrm{d}s$
$\displaystyle=$
$\displaystyle\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l+1}\mu^{k}(t-s)^{(l+1)\alpha+k(\alpha-\beta)-1}}{\Gamma((l+1)\alpha+k(\alpha-\beta))}g(s)\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k+1}(t-s)^{l\alpha+(k+1)(\alpha-\beta)-1}}{\Gamma(l\alpha+(k+1)(\alpha-\beta))}g(s)\mathrm{d}s$
$\displaystyle=$
$\displaystyle\lambda\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle+$
$\displaystyle\mu\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$

Therefore, we have

$\displaystyle\left(\prescript{RL}{0}{D}^{\alpha}_{t}y\right)(t)$
$\displaystyle=g(t)+\lambda\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle+\mu\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.12)

Similarly, the second term of (5.7) will be

$\displaystyle\left(\prescript{RL}{0}{D}^{\beta}_{t}y\right)(t)=\prescript{RL}{0}{D}^{\beta}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.13)

Taking linear combination of and (5.1) and (5.1) together with (5.8), we get
the desired result.

Case 2: We consider the Cauchy problem for generalized Bagley-Torvik equations
which is the special case of Caputo type fractional multi-term differential
equations with constant coefficients of order $1<\alpha\leq 2$ and
$0<\beta\leq 1$ in the form of:

$\displaystyle\begin{cases}\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)-\mu\left(\prescript{C}{0}{D}^{\beta}_{t}y\right)(t)-\lambda
y(t)=g(t),\quad t>0,\\\
y(0)=y^{\prime}(0)=0,\quad\quad\lambda,\mu\in\mathbb{R}.\end{cases}$ (5.14)

######  Theorem 5.9.

A unique solution $y\in C^{2}([0,\infty),\mathbb{R})$ of the Cauchy problem
(5.14) has the following formula:

$y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.15)

###### Proof.

We assume that (5.14) has a unique solution $y(t)$ and $g(t)$ is continuous on
$\left[0,\infty\right)$ and exponentially bounded, then $y(t)$,
$\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)$, and
$\left(\prescript{C}{0}{D}^{\beta}_{t}y\right)(t)$ are exponentially bounded,
thus their Laplace transform exist.

Applying the formula of Laplace transform for Caputo fractional derivative
(2.8) to the both sides of (5.14) yields:

$\displaystyle\left(s^{\alpha}-\mu s^{\beta}-\lambda\right)Y(s)=G(s).$ (5.16)

Then we solve (5.16) with respect to $Y(s)$,

$\displaystyle Y(s)=\frac{G(s)}{s^{\alpha}-\mu s^{\beta}-\lambda}.$ (5.17)

Taking inverse Laplace transform of (5.17) and applying Lemma 2.2, we find an
explicit representation of solution to (5.14):

$\displaystyle
y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.18)

∎

######  Remark 5.1.

Since initial conditions equal to zero, in accordance with the formula (2.6),
analytical solutions should be coincide with each other for the Cauchy
problems in Riemann-Liouville (5.7) and Caputo (5.14) senses.

Verification by substitution. Having found explicit form for $y(t)$, it
remains to confirm that $y(t)$ is an analytical solution of (5.14) indeed.

Now, we again make use of checking by substitution via Caputo fractional
Leibniz integral rule. In this case, we apply first Pascal’s rule before
applying fractional Leibniz rule since
$\prescript{C}{s}{D}^{\alpha}_{t}\left(\frac{(t-s)^{\alpha-1}}{\Gamma(\alpha)}\right)$
is undefined in accordance with (2.12). Then according to the formula (2.3),
we obtain

$\displaystyle\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)=\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=g(t)+\lambda\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle+\mu\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{2\alpha-\beta-1}E_{\alpha,\alpha-\beta,2\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.19)

Then using the formula for fractional Leibniz integral rule in Caputo sense
(4.8) of order $1<\alpha\leq 2$, we have

$\displaystyle\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{2-\alpha}_{t}}\Big{[}\lim\limits_{s\to
t-0}\sum_{l=1}^{2}\frac{\partial^{2-l}}{\partial
t^{2-l}}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)\Big{]}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{2-\alpha}_{t}}\Big{[}\lim\limits_{s\to
t-0}\frac{\partial}{\partial
t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)\Big{]}$
$\displaystyle+\prescript{t}{0}{I^{2-\alpha}_{t}}\Big{[}\lim\limits_{s\to
t-0}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d}{dt}\lim\limits_{s\to
t-0}g(s)\Big{]}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\prescript{t}{0}{I^{2-\alpha}_{t}}\Big{[}\lim\limits_{s\to
t-0}(t-s)^{2\alpha-2}E_{\alpha,\alpha-\beta,2\alpha-1}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)\Big{]}$
$\displaystyle+\prescript{t}{0}{I^{2-\alpha}_{t}}\Big{[}\lim\limits_{s\to
t-0}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\frac{d}{dt}\lim\limits_{s\to
t-0}g(s)\Big{]}$
$\displaystyle+\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$

Thus, by Lemma 2.3, we get

$\displaystyle\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-1}E_{\alpha,\alpha-\beta,2\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.20)

Similarly, by applying the formula (4.8), we also have

$\displaystyle\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{2\alpha-\beta-1}E_{\alpha,\alpha-\beta,2\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\prescript{C,t}{0}{D}^{\alpha}_{t}(t-s)^{2\alpha-\beta-1}E_{\alpha,\alpha-\beta,2\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.21)

On the other hand, since $0<\beta\leq 1$, we will apply the formula (4.9) for
second term of the equation (5.14):

$\displaystyle\prescript{C}{0}{D}^{\beta}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to
t-0}\prescript{t}{0}{I}^{1-\beta}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\beta}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to
t-0}(t-s)^{\alpha-\beta}E_{\alpha,\alpha-\beta,\alpha-\beta+1}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})\lim\limits_{s\to
t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\beta}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\prescript{RL,t}{0}{D}^{\beta}_{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$
(5.22)

where

$\prescript{t}{0}{I}^{1-\beta}_{t}\left\\{\frac{(t-t_{0})^{\alpha-1}}{\Gamma(\alpha)}\right\\}=\frac{(t-t_{0})^{\alpha-\beta}}{\Gamma(\alpha-\beta+1)},\quad
t>0.$

Next we plug (5.1) and (5.1) into (5.1), we therefore get

$\displaystyle\prescript{C}{0}{D}^{\alpha}_{t}\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=g(t)+\lambda\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle+\mu\int\limits_{0}^{t}(t-s)^{\alpha-\beta-1}E_{\alpha,\alpha-\beta,\alpha-\beta}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s.$

Taking linear combination of above equation together with (5.15) and (5.1), we
arrive at

$\displaystyle\left(\prescript{C}{0}{D}^{\alpha}_{t}y\right)(t)-\mu\left(\prescript{C}{0}{D}^{\beta}_{t}y\right)(t)-\lambda
y(t)=g(t),\quad t>0.$

######  Remark 5.2.

Kilbas et al. [27] have obtained the analytical solutions of the Cauchy
problems (5.7) and (5.14) in terms of Fox-Wright functions below:

$y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}H_{\alpha,\beta;\lambda,\mu}(t-s)g(s)\mathrm{d}s,\quad
t>0,$ (5.23)

where

$H_{\alpha,\beta;\lambda,\mu}(t)\coloneqq\sum_{l=0}^{\infty}\frac{\lambda^{l}t^{l\alpha}}{l!}\prescript{}{1}{\Psi_{1}}\left[\begin{array}[]{ccc}(l+1,1)\\\
(l\alpha+\alpha,\alpha-\beta)\end{array}\Big{|}\mu t^{\alpha-\beta}\right].$

###### Proof.

Using the definition of Fox-Wright function [45, 46], we arrive at

$\displaystyle y(t)$
$\displaystyle=\int\limits_{0}^{t}\sum_{l=0}^{\infty}\frac{\lambda^{l}(t-s)^{l\alpha+\alpha-1}}{l!}\prescript{}{1}{\Psi_{1}}\left[\begin{array}[]{ccc}(l+1,1)\\\
(l\alpha+\alpha,\alpha-\beta)\end{array}\Big{|}\mu(t-s)^{\alpha-\beta}\right]g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}(t-s)^{l\alpha+k(\alpha-\beta)+\alpha-1}}{\Gamma(l\alpha+k(\alpha-\beta)+\alpha)}g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s,\quad
t>0.$

Therefore, our solution in terms of univariate version of bivariate Mittag-
Leffler type functions coincide with the solution by means of Fox-Wright type
functions shown in [27]. ∎

######  Remark 5.3.

Podlubny [8] have attained the analytical solutions of the Cauchy problems
(5.7) and (5.14) in terms of $l$-th derivative of two-parameter Mittag-Leffler
functions below:

$y(t)=\int\limits_{0}^{t}(t-s)^{\alpha-1}H_{\alpha,\beta;\lambda,\mu}(t-s)g(s)\mathrm{d}s,\quad
t>0,$ (5.24)

where

$H_{\alpha,\beta;\lambda,\mu}(t)\coloneqq\sum_{l=0}^{\infty}\frac{\lambda^{l}t^{l\alpha}}{l!}\mathcal{E}^{(l)}_{\alpha-\beta,\alpha+l\beta}(\mu
t^{\alpha-\beta}).$

###### Proof.

Using the definition of $l$-th derivative of two-parameter Mittag-Leffler
function, we arrive at

$\displaystyle y(t)$
$\displaystyle=\int\limits_{0}^{t}\sum_{l=0}^{\infty}\frac{\lambda^{l}(t-s)^{l\alpha+\alpha-1}}{l!}\mathcal{E}^{(l)}_{\alpha-\beta,\alpha+l\beta}(\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\sum_{l=0}^{\infty}\frac{\lambda^{l}}{l!}(t-s)^{(l+1)\alpha-1}\sum_{k=0}^{\infty}\frac{(l+k)!}{k!}\frac{\mu^{k}(t-s)^{k(\alpha-\beta)}}{\Gamma(k(\alpha-\beta)+l(\alpha-\beta)+l\beta+\alpha)}g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\sum_{l=0}^{\infty}\sum_{k=0}^{\infty}\binom{l+k}{k}\frac{\lambda^{l}\mu^{k}(t-s)^{l\alpha+k(\alpha-\beta)+\alpha-1}}{\Gamma(l\alpha+k(\alpha-\beta)+\alpha)}g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha-\beta,\alpha}(\lambda(t-s)^{\alpha},\mu(t-s)^{\alpha-\beta})g(s)\mathrm{d}s,\quad
t>0.$

Therefore, our solution in terms of univariate version of bivariate Mittag-
Leffler type functions coincide with the solution by means of $l$-th
derivative of two-parameter Mittag-Leffler type functions shown in [8]. ∎

Case 3: In special case, we substitute $\alpha=2$ and $\beta=1$ in (5.7) and
(5.14), then we get the following Cauchy problem for the classical second
order linear differential equation - the oscillator equation with constant
coefficients:

$\begin{cases}y^{\prime\prime}(t)-\mu y^{\prime}(t)-\lambda y(t)=g(t),\quad
t>0,\\\ y^{\prime}(0)=y(0)=0,\quad\lambda,\mu\in\mathbb{R},\end{cases}$ (5.25)

######  Theorem 5.10.

A unique solution $y\in C^{2}([0,\infty),\mathbb{R})$ of the Cauchy problem
(5.25) has the following formula:

$y(t)=\int\limits_{0}^{t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s,\quad
t>0.$ (5.26)

###### Proof.

By using verification by substitution, we have

$\displaystyle\frac{d^{2}}{dt^{2}}\int\limits_{0}^{t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\sum_{l=1}^{2}\lim\limits_{s\to
t-0}\frac{\partial^{2-l}}{\partial
t^{2-l}}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))\frac{d^{l-1}}{dt^{l-1}}\lim\limits_{s\to
t-0}g(s)$ $\displaystyle+\int\limits_{0}^{t}\frac{\partial^{2}}{\partial
t^{2}}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to t-0}\frac{\partial}{\partial
t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))\lim\limits_{s\to t-0}g(s)$
$\displaystyle+\lim\limits_{s\to
t-0}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))\frac{d}{dt}\lim\limits_{s\to
t-0}g(s)$ $\displaystyle+\int\limits_{0}^{t}\frac{\partial^{2}}{\partial
t^{2}}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\lim\limits_{s\to
t-0}E_{2,1,1}(\lambda(t-s)^{2},\mu(t-s))\lim\limits_{s\to t-0}g(s)$
$\displaystyle+\int\limits_{0}^{t}(t-s)^{-1}E_{2,1,0}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=g(t)+\lambda\lim\limits_{s\to
t-0}(t-s)^{2}E_{2,1,3}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle+\mu\lim\limits_{s\to
t-0}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle+\int\limits_{0}^{t}(t-s)^{-1}E_{2,1,0}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=g(t)+\int\limits_{0}^{t}(t-s)^{-1}E_{2,1,0}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s.$
(5.27)

Applying Pascal’s rule for binomial coefficients for the last term of above
equality, we get

$\displaystyle\int\limits_{0}^{t}(t-s)^{-1}E_{2,1,0}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\frac{(t-s)^{-1}}{\Gamma(0)}g(s)\mathrm{d}s+\lambda\int\limits_{0}^{t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle+\mu\int\limits_{0}^{t}E_{2,1,1}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\lambda\int\limits_{0}^{t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle+\mu\int\limits_{0}^{t}E_{2,1,1}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s.$
(5.28)

Next, we have

$\displaystyle\frac{d}{dt}\int\limits_{0}^{t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}\frac{\partial}{\partial
t}(t-s)E_{2,1,2}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s$
$\displaystyle=\int\limits_{0}^{t}E_{2,1,1}(\lambda(t-s)^{2},\mu(t-s))g(s)\mathrm{d}s.$
(5.29)

Again taking linear combination of above equations (5.1) and (5.1) together
with (5.26), we prove the desired result. ∎

##  6 Conclusions and future work

The theory of Leibniz integral rule allows us to study particular solutions of
classical and fractional multi-term differential equations. To the best of our
knowledge, we derive explicit analytical solutions of well-known Bagley-Torvik
and oscillator equations in terms of bivariate Mittag-Leffler functions via
the technique of fractional Green’s function, since this theory has not been
presented in recent literature.

The major contributions of our research work are as below:

  * 1.

we have proposed a Leibniz rule for higher order derivatives in classical
sense which is more productive tool for testing solutions of multi-order
differential equation;

  * 2.

we have introduced fractional Leibniz rule for Riemann-Liouville and Caputo
type fractional differentiation operators;

  * 3.

we have investigated differentiation of convolution operator which is more
crucial in theory of differential equations with constant coefficients of
classical and fractional-order derivatives;

  * 4.

analytical explicit solutions of the generalized Bagley-Torvik and oscillator
equations are derived in terms of univariate version of bivariate Mittag-
Leffler type functions in accordance with the method of Laplace integral
transform;

  * 5.

We have showed that our analytical solutions are coincide with Fox-Wright type
and $l$-th derivative of two-parameter Mittag-Leffler type functions;

  * 6.

we tested the candidate solutions of Cauchy problems for Bagley-Torvik
equations with fractional-order sense and oscillator equation with classical-
order one via our new fractional Leibniz integral rules.

There are a number of potential directions in which the results acquired here
can be extended. Our future work will proceed to study the Leibniz integral
rule results for $\psi$-Hilfer and Hadamard type fractional derivatives and
the analytical explicit solutions of multi-term fractional differential
equations in terms of natural extensions of Mittag-Leffler type functions.

## References

  * [1] B. Datsko, V. Gafiychuk, Complex spatio-temporal solutions in fractional reaction-diffusion systems near a bifurcation point. _Fract. Calc. Appl. Anal._ 21(2018), 237-253.
  * [2] H. Schissel, R. Metzler, A. Blumen, T. F. Nonnenmacher, Generalized viscoelastic models: their fractional equations with solutions. _J. Phys. A: Math. Gen._ 28(1995), 6567-6584.
  * [3] A. Ahmadova, N. I. Mahmudov, Langevin differential equations with general fractional orders and their applications to electric circuit theory, _J. Comput. Appl. Math._ (2020), accepted.
  * [4] T. Kaczorek, Positive linear systems with different fractional orders._Bull. Pol. Acad. Sci. Tech. Sci._ 58(2010), 453–458.
  * [5] N. I. Mahmudov, Approximate controllability of semilinear deterministic and stochastic evolution equations in abstract spaces. _SIAM J. Control. Optim._ 42(2003), 1604-1622.
  * [6] A. Ahmadova, N. I. Mahmudov, Existence and uniqueness results for a class of fractional stochastic neutral differential equations, _Chaos Solitons &Fract._ 139(2020), https://doi.org/10.1016/j.chaos.2020.110253.
  * [7] I. T. Huseynov, N. I. Mahmudov, Delayed analogue of three-parameter Mittag-Leffler functions and their applications to Caputo type fractional time-delay differential equations. _Math. Meth. Appl. Sci._ 2020, https://doi.org/10.1002/mma.6761.
  * [8] I. Podlubny, Fractional Differential Equations. New York: Academic Press, 1999.
  * [9] D. Baleanu, J. I. Trujillo, A new method of finding the fractional Euler–Lagrange and Hamilton equations within Caputo fractional derivatives. _Commun. Nonlinear Sci. Numer. Simulat._ 15(2010), 1111–1115.
  * [10] R. Hilfer, Y. F. Luchko, Z. Tomovski, Operational method for the solution of fractional differential equations with generalized Riemann-Liouville fractional derivatives, _Fract. Calc. Appl. Anal._ 12(2009), 299–318.
  * [11] Y. F. Luchko, R. Gorenflo, An operational method for solving fractional differential equations with Caputo derivatives. _Acta. Math. Vietnam._ 24(1999), 207-233.
  * [12] E. Bazhlekova, Properties of the fundamental and the impulse-response solutions of multi-term fractional differential equations. _Complex Analysis and Applications’13 (Proc. Intern. Conf., Sofia), Bulg. Acad. Sci. Sofia_ (2013) 55–64.
  * [13] S. Pak, H. K. Choi, K. R. Sin, Analytical solutions of linear inhomogeneous fractional differential equation with continuous variable coefficients, _Adv. Differ. Equ._ (2019), 256.
  * [14] K. Diethelm, N. J. Ford, Multi-order fractional differential equations and their numerical solution, _Appl. Math. Comp._ 154(2004), 621-640.
  * [15] B. Bonilla, M. Rivero, J. J. Trujillo, On systems of linear fractional differential equations with constant coefficients,_Appl. Math. Comp._ 187(2007), 68-78.
  * [16] K. Diethelm, S. Siegmund, H. T. Tuan, Asymptotic behavior of solutions of linear multi-order fractional differential equations systems, _Fract. Cal. Appl. Anal._ 20(2017), 1165-1195.
  * [17] T. Kaczorek, D. Idczak, Cauchy formula for the time-varying linear systems with Caputo derivative, _Fract. Cal. Appl. Anal._ 20(2017), 494-505.
  * [18] I. T. Huseynov, A. Ahmadova, A. Fernandez, N. I. Mahmudov, Explicit analytic solutions of incommensurate fractional differential equation systems. _Appl. Math. Comp._ (2020), 125590, https://doi.org/10.1016/j.amc.2020.125590.
  * [19] P. J. Torvik, R. L. Bagley, On the appearance of the fractional derivative in the behavior of real materials, _J. Appl. Mech._ 51(1984), 294-298.
  * [20] N. I. Mahmudov, I. T. Huseynov, N. A. Aliev, F. A. Aliev, Analytical approach to a class of Bagley-Torvik equations. _TWMS J. Pure Appl. Math._ 11(2020), 238-258.
  * [21] Z. Wang, X. Wang, General solution of the Bagley-Torvik equation with fractional-order derivative. _Commun. Nonlinear Sci. Numer. Simul._ 15(2010), 1279-1285.
  * [22] K. Diethelm, J. Ford, Numerical Solution of the Bagley-Torvik Equation._BIT Numer. Math._ 42(2002), 490–507, https://doi.org/10.1023/A:1021973025166.
  * [23] H. M. Srivastava, F. A. Shah, R. Abbas, An application of the Gegenbauer wavelet method for the numerical solution of the fractional Bagley-Torvik equation. _Russian J. Math. Phys._ 26(2019), 77–93.
  * [24] K. S. Miller, Linear differential equation in the real domain, W. W. Norton and Co., New York, 1963.
  * [25] K. Diethelm, The analysis of fractional differential equations. An application-oriented exposition using differential operators of Caputo type, Lecture notes in Mathematics. Berlin: Springer-Verlag, 2010.
  * [26] R. Gorenflo, A. A. Kilbas, F. Mainardi, S. V. Rogosin, Mittag-Leffler Functions, Related Topics and Applications. Berlin: Springer-Verlag, 2014.
  * [27] A. A. Kilbas, H. M. Srivastava, J. J. Trujillo, Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier, 2006.
  * [28] K. S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations. New York: Wiley, 1993.
  * [29] G. Mittag-Leffler, Sur la nouvelle fonction $E_{\alpha}(x)$. _C. R. Acad. Sci Paris_ , 137(1903), 554–558.
  * [30] T. R. Prabhakar, A singular integral equation with a generalized Mittag-Leffler function in the kernel. _Yokohama Math. J._ 19(1971), 7-15.
  * [31] V. S. Kiryakova, Generalized Fractional Calculus and Applications, Longman & J. Wiley, Harlow - N. York, 1994.
  * [32] V. Kiryakova, Multiple (multiindex) Mittag-Leffler functions and relations to generalized fractional calculus, J. Comput. Appl. Math. 118 (2000), 241–259.
  * [33] V. S. Kiryakova, The multi-index Mittag-Leffler function as an important class of special functions of fractional calculus, _Comput. Math. Appl._ 59(2010), 1885–1895.
  * [34] V. S. Kiryakova, The special functions of fractional calculus as generalized fractional calculus operators of some basic functions, _Comput. Math. Appl._ 59(2010), 1128–1141
  * [35] A. Fernandez, C Kürt, M. A. Özarslan, A naturally emerging bivariate Mittag-Leffler function and associated fractional-calculus operators, _Comp. Appl. Math._ 39(2020), https://doi.org/10.1007/s40314-020-01224-5.
  * [36] I. T. Huseynov, A. Ahmadova, G. O. Ojo, N. I. Mahmudov, A natural extension of Mittag-Leffler function associated with a triple infinite series, Mathematics, under review. Preprint available at arXiv:2011.03999.
  * [37] R. K. Saxena, S. L. Kalla, R.Saxena, Multivariate analogue of generalised Mittag-Leffler function, _Integr. Transf. Spec. F._ 22(2011), 533-548.
  * [38] W. Chen, Introduction to Lebesgue Integration. Imperial College, London, 2003.
  * [39] G. I. Fikhtengoltz, Course of Differential and Integral Calculus, vol. 2, Nauka, Moscow, 1969.
  * [40] I. Matychyn, V. Onyshchenko, Solution of linear fractional order systems with variable coefficients, _Fract. Cal. Appl. Anal._ 23(2020), 753-763.
  * [41] I. Matychyn, Analytical solution of linear fractional systems with variable coefficients involving Riemann-Liouville and Caputo derivatives, _Symmetry_ 11(2019), 1366.
  * [42] K. B. Oldham, J. Spanier, The Fractional Calculus. San Diego: Academic Press, 1974.
  * [43] E. D. Rainville, Special functions. New York, Macmillan, 1960.
  * [44] S. G. Samko, A. A. Kilbas, O. I. Marichev, Fractional Integrals and Derivatives: Theory and Applications. London: Taylor & Francis , 2002.
  * [45] C. Fox, The asymptotic expansion of generalized hypergeometric functions, _Proc. London Math. Soc._ 27(1928), 389-400.
  * [46] E. M. Wright, The asymptotic expansion of the generalized hypergeometric function, _J. London Math. Soc_. 10(1935), 286-293.
  * [47] I. N. Sneddon, The use of integral transforms. New Delhi: Tata McGraw-Hill, 1979.
  * [48] E. T. Whittaker, G. N. Watson, A course of modern analysis, fourth ed. Cambridge: Cambridge University Press, 1927.