Probabilistic solutions of fractional differential and partial differential
equations and their Monte Carlo simulations

Tamer Oraby, Harrinson Arrubla, Erwin Suazo∗

School of Mathematical and Statistical Sciences, The University of Texas Rio
Grande Valley,

1201 W. University Dr., Edinburg, Texas, USA,

∗Corresponding author<EMAIL_ADDRESS>
## Abstract

The work in this paper is four-fold. Firstly, we introduce an alternative
approach to solve fractional ordinary differential equations as an expected
value of a random time process. Using the latter, we present an interesting
numerical approach based on Monte Carlo integration to simulate solutions of
fractional ordinary and partial differential equations. Thirdly, we show that
this approach allows us to find the fundamental solutions for fractional
partial differential equations (PDEs), in which the fractional derivative in
time is in the Caputo sense and the fractional in space one is in the Riesz-
Feller sense. Lastly, using Riccati equation, we study families of fractional
PDEs with variable coefficients which allow explicit solutions. Those
solutions connect Lie symmetries to fractional PDEs.

Keywords: Caputo fractional derivative, Riesz-Feller fractional derivative,
Riccati equation, Lie symmetries, Green functions, Monte Carlo Integration,
Mittag-Leffler function.

##  1 Introduction

Although it was started in the second half of the eighteenth century by
Leibniz, Newton and l’Hôpital, [1], fractional calculus has received great
attention in the last two decades. Many physical, biological and
epidemiological models have found that fractional order models could perform
at least as good as well as their integer counterparts, [2, 3]. Integer order
models are also appearing as special cases of the fractional order. That makes
the dynamical behavior of those models richer and in some cases flexible. As
advances are made in fractional calculus and fractional modeling,
understanding of the physical interpretation of fractional derivatives is
becoming clearer. A memory kernel with algebraic decay is the most common
interpretation for the change in the dynamical behavior of the system [4, 5,
6]. Today fractional calculus is widely used in physical modeling. Examples
include the nonexponential relaxation in dielectrics and ferromagnets [7],
[8], the diffusion processes [9], [10] and the Hamiltonian Chaos [11] and
[12].

Anomalous diffusion could be modeled using fractional order stochastic
processes and their Fokker-Planck equations [13, 14, 15]. Continuous-time
random walk (CTRW) is one approach used to model anomalous diffusion [16]. In
particular, a CTRW with infinite-mean time to jump exhibits sub-diffusion
behavior. The time to jump could be modeled by a heavy tail distribution with
index $\beta$ such that $0<\beta<1$ leading to a mean square displacement that
is of order $t^{\beta}$ depicting the short-range jump and so thus sub-
diffusion. A long-range jump leads to super-diffusion, i.e., $\beta>1$, [17].

A subordinator is a non-decreasing Lévy process, i.e. a non-negative process
with independent and stationary increments. It is used as a random time, or
operational time, in defining time-changed-processes. Its density decays
algebraically as $1/t^{\alpha+1}$ as $t\rightarrow\infty$, with
$\alpha\in(0,1)$. A different type of time-change can be made by using the
inverse subordinator (i.e the hitting time of a subordinator) which sometimes
leads to sud-diffusion processes [18, 19].

The Riccati equation has played an important role in finding explicit
solutions for Fisher and Burgers equations, (see [20] and [21] and references
therein). Also, similarity transformations and the solutions of Riccati and
Ermakov systems have been extensively applied thanks to Lie groups and Lie
algebras [22], [23], [24], [25] and [26]. In this work, we show that this
approach allows us to find the fundamental solutions for fractional PDEs; the
fractional derivative in time is in the Caputo sense and the fractional
derivative in space one is in the Riesz-Feller sense. We establish a
relationship between the coefficients through the Riccati equation; we study
families of fractional PDEs with variable coefficients which allow explicit
solutions and find those explicit solutions. Those solutions connect Lie
symmetries to fractional PDEs.

Formulating solutions of fractional differential equations and partial
differential equations as expected values with respect to heavy-tail or power
law distributions could be enabled using Monte-Carlo integration methods and
Sampling Importance Integration to evaluate them. The mean problem is in the
number of simulations or random number generations that need to be done to
guarantee convergence and small standard error.

In Section 1, we will review fundamental definitions and classical results
needed from the classical theory of fractional differential equations. In
Section 2, we present the first main result of this work, Lemma 1, which allow
us to see a fractional functions and their fractional derivatives as Wright
type transformations of some functions and their derivatives. They could be
also interpreted as expected values of functions in a random time process.
Also, in Section 2 we present Theorem 2.1 which allow us to solve fractional
ordinary differential equations through solutions of regular ordinary
differential equations. In Section 3, we derive fractional green functions for
some important fractional partial differential equations, like diffusion,
Schrödinger, and wave equation. In Section 4, we derive green functions of
fractional partial differential equations with variable coefficients with
application to Fokker-Planck equations. In Section 5, we use the integral
transform or the expected value interpretation of the solutions of fractional
equations to carry out Monte Carlo simulations of their solution.

###  1.1 Preliminaries

In this section we give the required background of Caputo and Riesz-Feller
fractional differentiation.

#### Caputo Derivative.

Let $D^{n}$ be the Leibniz integer-order differential operator given by

$D^{n}f=\dfrac{d^{n}f}{dt^{n}}=f^{(n)},$

and let $J^{n}$ be an integration operator of integer order given by

$J^{n}f(t)=\dfrac{1}{n-1!}\int_{0}^{t}(t-\tau)^{n-1}f(\tau)d\tau,$ (1.1)

where $n\in\mathbb{Z}^{+}$. Let us use $D=D^{1}$ for the first derivative. For
fraction-order integrals, we use

$J^{n-\beta}f(t)=\dfrac{1}{\Gamma(n-\beta)}\int_{0}^{t}(t-\tau)^{n-\beta-1}f(\tau)d\tau,$
(1.2)

where $n-1<\beta\leq n$. Now, define the Caputo fractional differential
operator $D_{C}^{\beta}$ to be

$D_{C}^{\beta}f(t)=J^{n-\beta}D^{n}f(t),$

where $n-1<\beta\leq n$, for $n\in\mathbb{N}$.

#### Riemann-Liouville Derivative.

The Riemann-Liouville fractional differential operator $D_{RL}^{\beta}$ is
defined to be

$D_{RL}^{\beta}f(t)=D^{n}J^{n-\beta}f(t),$

where $n-1<\beta\leq n$, for $n\in\mathbb{N}$. We will use
$\partial^{\beta}_{t}F:=\dfrac{\partial^{\beta}F}{\partial t^{\beta}}$ and use
$\partial_{t}F:=\dfrac{\partial F}{\partial t}$.

The Riemann-Liouville fractional is related to the Caputo fractional
derivative through [27]:

$D_{RL}^{\beta}f(t)=D_{C}^{\beta}f(t)+\sum_{k=0}^{n-1}f^{(k)}(0)\dfrac{t^{k}}{k!}.$

While we will not discuss the Riemann-Liouville fractional derivatives in the
paper, the results presented in this paper are valid for Riemann-Liouville
fractional derivatives, when $f^{(k)}(0)=0$ for $k=0,1,\ldots,n-1$.

We will consider $n=1$ in this work; that is $0<\beta\leq 1$. Some of the
results be extended through the remark that for $0<\beta\leq 1$,
$D_{C}^{n+\beta-1}f(t)=D_{C}^{\beta}f^{(n-1)}(t)$ for $n\geq 1$. Note also
that when $0<\beta\leq 1$, $D_{RL}^{\beta}f(t)=D_{C}^{\beta}f(t)+f(0)$.

#### Riesz-Feller Derivative

The Riesz-Feller fractional differential operator $D_{RF}^{\alpha,\theta}$ is
defined to be [28]

$\displaystyle D_{RF}^{\alpha,\theta}f(x)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(1+\alpha)}{\pi}\sin((\alpha+\theta)\frac{\pi}{2})\int_{0}^{\infty}\dfrac{f(x+y)-f(x)}{y^{1+\alpha}}dy$
$\displaystyle+\frac{\Gamma(1+\alpha)}{\pi}\sin((\alpha-\theta)\frac{\pi}{2})\int_{0}^{\infty}\dfrac{f(x-y)-f(x)}{y^{1+\alpha}}dy$

for fractional order $0<\alpha\leq 2$, and the skewness parameter
$\theta\leq\min(\alpha,2-\alpha)$. The symmetric Riesz-Feller differential
operator is defined at $\theta=0$ and is simply denoted by $D_{RF}^{\alpha}$.

#### Transformations.

The Laplace transform of a function $f(t)$ is defined as

$\mathcal{L}(f)(s)=\widetilde{f}(s)=\int_{0}^{\infty}e^{-st}f(t)dt.$

The inverse Laplace transform is defined by

$\mathcal{L}^{-1}\left(\widetilde{f}\right)(t)=\dfrac{1}{2\pi
i}\int_{\mathcal{C}}e^{st}\widetilde{f}(s)ds$

where $\mathcal{C}$ is a contour parallel to the imaginary axis and to the
right of the singularities of $\widetilde{f}$. The Laplace transform of the
Caputo fractional derivative of a function is given by

$\mathcal{L}\left(D_{C}^{\beta}f\right)(s)=s^{\beta}\widetilde{f}(s)-\sum_{k=0}^{n-1}s^{\beta-1-k}f^{(k)}(0).$
(1.3)

The Fourier transform of a function $f(x)$ is defined as

$\mathcal{F}(f)(y)=\widehat{f}(y)=\int_{-\infty}^{\infty}e^{ixy}f(x)dx.$

The inverse Fourier transform is defined by

$\mathcal{F}^{-1}\left(\widehat{f}\right)(x)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}e^{-ixy}\widehat{f}(y)dy.$

The Fourier transform of the Riesz-Feller fractional derivative of a function
is given by

$\mathcal{F}\left(D_{RF}^{\alpha,\theta}f\right)(y)=-\psi_{\alpha}^{\theta}(y)\widehat{f}(y),$
(1.4)

where
$\psi_{\alpha}^{\theta}(y)=|y|^{\alpha}e^{\frac{i\text{sign}(y)\theta\pi}{2}}$.

#### Mittag-Leffler Function

The Mittag-Leffler function, which generalizes the exponential function, can
be written as follows,

$\begin{split}E_{\beta}(z)=\sum\limits_{k=0}^{\infty}\dfrac{z^{k}}{\Gamma(\beta
k+1)},\end{split}\qquad\begin{split}\beta\in\mathbb{R}^{+},\enspace
z\in\mathbb{C},\end{split}$ (1.5)

and the more general Mittag-Leffler function with two-parameters is defined to
be

$\begin{split}E_{\beta,\alpha}(z)=\sum\limits_{k=0}^{\infty}\dfrac{z^{k}}{\Gamma(\beta
k+\alpha)},\end{split}\qquad\begin{split}\beta,\alpha\in\mathbb{R}^{+},\enspace
z\in\mathbb{C}.\end{split}$ (1.6)

#### Wright function.

The Wright function is another special function of importance to fractional
calculus and is defined by [29],

$\begin{split}W_{\beta,\alpha}(z)=\sum\limits_{k=0}^{\infty}\dfrac{z^{k}}{k!\Gamma(\beta
k+\alpha)},\end{split}\qquad\begin{split}\beta>-1,\alpha\in\mathbb{C},\enspace
z\in\mathbb{C}.\end{split}$ (1.7)

The following Wright type function will be fundamental in the rest of this
work

$g_{\beta}(x;t)=\dfrac{1}{t^{\beta}}W_{-\beta,1-\beta}\left(-\dfrac{x}{t^{\beta}}\right),$
(1.8)

which is a probability density function of the random time process
$\mathcal{T}_{\beta}(t)$ for all $t>0$ [30, 31]. $\mathcal{T}_{\beta}(t)$ can
be seen as the inverse of a $\beta$ stable subordinator with density
$g_{\beta}(x,t)$, see [34, 35, 36]. It has a Laplace transform

$\mathcal{L}\left(g_{\beta}(\cdot;t)\right)(s)=\int_{0}^{\infty}e^{-sx}g_{\beta}(x;t)dx=E_{\beta}\left(-st^{\beta}\right)$
(1.9)

for $\Re(s)>0$ and moments
$\mathbb{E}\left[(\mathcal{T}_{\beta}(t))^{k}\right]=\Gamma(k+1)\dfrac{t^{k\beta}}{\Gamma(k\beta+1)}$
for $k\geq 1$ [37, 38]. At the same time

$\int_{0}^{\infty}e^{-st}g_{\beta}(x;t)dt=s^{\beta-1}e^{-xs^{\beta}}.$ (1.10)

It is shown in [30] that for $0<\beta,\alpha<1$,

$g_{\beta\alpha}(x;t)=\int_{0}^{\infty}g_{\beta}(x;s)g_{\alpha}(s;t)ds$ (1.11)

for $x>0$. For more details about Wright function see [30, 31]. From (1.9) and
(1.11), we get

$\int_{0}^{\infty}E_{\beta}\left(-st^{\beta}\right)g_{\alpha}(t;r)dt=\int_{0}^{\infty}\int_{0}^{\infty}e^{-sx}g_{\beta}(x;t)g_{\alpha}(t;r)dxdt=E_{\beta\alpha}\left(-sr^{\beta\alpha}\right),$

and from (1.10) and (1.11), we get

$\int_{0}^{\infty}E_{\alpha}\left(-tx^{\alpha}\right)g_{\beta}(x;t)dt=\int_{0}^{\infty}s^{\beta-1}e^{-xs^{\beta}}g_{\alpha}(s;x)ds.$

#### Lévy $\alpha$-stable distribution.

The Lévy $\alpha$-stable distribution with stability index $0<\alpha\leq 2$,
$L_{\alpha}^{\theta}(x)$ has a Fourier transform given by

$\mathcal{F}\left(L_{\alpha}^{\theta}(\cdot)\right)(y):=\widehat{L}_{\alpha}^{\theta}(y)=e^{-\psi_{\alpha}^{\theta}(y)}.$
(1.12)

The density $L_{\alpha}^{\theta}(x)$ has a fat tail proportional to $\lvert
x\rvert^{-(1+\alpha)}$.

Define
$L_{\alpha}^{\theta}(y;x)=\dfrac{1}{x^{\frac{1}{\alpha}}}L_{\alpha}^{\theta}\left(\dfrac{y}{x^{\frac{1}{\alpha}}}\right)$
for $y\in\mathbb{R}$ and $x>0$ which is a probability density function of the
$\alpha$-stable random process with asymmetry parameter $\theta$, denoted by
$\mathcal{L}_{\alpha}^{\theta}(x)$ for $x>0$, [39].

$\int_{-\infty}^{\infty}e^{isy}L_{\alpha}^{\theta}(y;x)dy=e^{-\psi_{\alpha}^{\theta}(s)x}.$
(1.13)

For $0<\alpha<1$, and $t,x>0$,

$L_{\alpha}^{-\alpha}(t;x)=\dfrac{x\alpha}{t}g_{\alpha}(x;t).$ (1.14)

See [40] for more details. Also, for $0<\alpha\leq 1$

$L_{\beta\alpha}^{\theta\alpha}(x;t)=\int_{0}^{\infty}L_{\beta}^{\theta}(x;s)L_{\alpha}^{-\alpha}(s;t)ds.$
(1.15)

##  2 Fractional Derivative as Expected Value with Respect to a Lévy
Distribution

In this Section we establish the first main result of this paper.

###  2.1 Caputo Fractional Derivative

The following lemma is one of the main results of this work for its wide
applicability.

###### Lemma 1

  1. 1.

Let $f\in C([0,\infty))$ and $f_{\beta}$ be a function. Hence,

$\widetilde{f}_{\beta}(s)=s^{\beta-1}\widetilde{f}\left(s^{\beta}\right)$
(2.1)

for $s\in[0,\infty)$ if and only if

$f_{\beta}(t)=\int_{0}^{\infty}f(x)g_{\beta}(x;t)dx$

for $0<\beta\leq 1$, and if integrals exist.

  2. 2.

Let $f\in C^{1}([0,\infty))$ and $f_{\beta}$ be a function such that

$\widetilde{f}_{\beta}(s)=s^{\beta-1}\widetilde{f}\left(s^{\beta}\right)$

for $s\in[0,\infty)$ and $f_{\beta}(0)=f(0)$ then

$D_{C}^{\beta}f_{\beta}(t)=\int_{0}^{\infty}f^{(1)}(x)g_{\beta}(x;t)dx$

for $0<\beta\leq 1$.

  3. 3.

Let $f\in C^{n}([0,\infty))$ such that

$\mathcal{L}\left(f_{\beta}^{(n-1)}(\cdot)\right)(s)=s^{\beta-1}\mathcal{L}\left(f^{(n-1)}(\cdot)\right)\left(s^{\beta}\right)$
(2.2)

for $s\in[0,\infty)$ and $f_{\beta}^{(n-1)}(0)=f^{(n-1)}(0)$. For $0<\beta\leq
1$, the following holds

$D_{C}^{n+\beta-1}f_{\beta}(t)=\int_{0}^{\infty}f^{(n)}(x)g_{\beta}(x;t)dx$

for $n\geq 1$.

  4. 4.

For $0<\alpha,\beta<1$ and such that $f_{\beta\alpha}(0)=f_{\beta}(0)$

$D_{C}^{\beta\alpha}f_{\beta\alpha}(t)=D_{C}^{\alpha}\left(D_{C}^{\beta}f_{\beta}(s)\right).$

###### Remark 1

The condition in equation (2.2) of part (3), is equivalent to stating that
$f_{\beta}^{(n-1)}=\left(f^{(n-1)}\right)_{\beta}$ in the sense of equation
(2.1).

##### Proof

  1. 1.

To show sufficiency, we use the Laplace transform as follows:

$\displaystyle\begin{split}\mathcal{L}\left(\int_{0}^{\infty}f(x)g_{\beta}(x;\cdot)dx\right)(s)&=\int_{0}^{\infty}\int_{0}^{\infty}e^{-st}f(x)g_{\beta}(x;t)dxdt\\\
\text{by equation
\eqref{eq:lw2}}&=\int_{0}^{\infty}f(x)s^{\beta-1}e^{-xs^{\beta}}dx\\\
&=s^{\beta-1}\widetilde{f}\left(s^{\beta}\right)\\\
&=\widetilde{f}_{\beta}(s)=\mathcal{L}\left(f_{\beta}(\cdot)\right)(s).\end{split}$

Necessity follows from the same lines.

  2. 2.

We will show that the Laplace transform of the right-hand side of equation (2)
is equal to that of the left hand side which is given by (1.3). The Laplace
transform of the right hand side of equation (2) is given by

$\displaystyle\begin{split}\mathcal{L}\left(\int_{0}^{\infty}f^{(1)}(x)g_{\beta}(x;\cdot)dx\right)(s)&=\int_{0}^{\infty}\int_{0}^{\infty}e^{-st}f^{(1)}(x)g_{\beta}(x;t)dxdt\\\
\text{by equation
\eqref{eq:lw2}}&=\int_{0}^{\infty}f^{(1)}(x)s^{\beta-1}e^{-xs^{\beta}}dx\\\
&=s^{\beta-1}\left(s^{\beta}\widetilde{f}(s^{\beta})-f(0)\right)\\\
&=s^{\beta}\widetilde{f}_{\beta}(s)-s^{\beta-1}f_{\beta}(0)\\\
&=\mathcal{L}\left(D_{C}^{\beta}f_{\beta}(\cdot)\right)(s).\end{split}$

  3. 3.

It is easy to prove that condition (2.1) is equivalent to

$f^{(n-1)}_{\beta}(t)=(f^{(n-1)})_{\beta}(t)$ (2.3)

By properties of Caputo derivative, conditions (2.3) and
$f_{\beta}^{(n-1)}(0)=f^{(n-1)}(0)$ we obtain

$D_{C}^{\beta+n-1}f_{\beta}(t)=D_{C}^{\beta}f^{(n-1)}_{\beta}(t)=D_{C}^{\beta}(f^{(n-1)})_{\beta}(t)=\int_{0}^{\infty}f^{(n)}(x)g_{\beta}(x;t)dx.$

as we wanted.

  4. 4.

Using equation (1.11),

$D_{C}^{\beta\alpha}f_{\beta\alpha}(t)=\int_{0}^{\infty}f^{\prime}(x)g_{\beta\alpha}(x;t)dx=\int_{0}^{\infty}f^{\prime}(x)\int_{0}^{\infty}g_{\beta}(x;s)g_{\alpha}(s;t)dsdx,$

and using part 2, we obtain

$\displaystyle\begin{split}D_{C}^{\beta\alpha}f_{\beta\alpha}(t)&=\int_{0}^{\infty}g_{\alpha}(s;t)\int_{0}^{\infty}f^{\prime}(x)g_{\beta}(x;s)dxds\\\
&=\int_{0}^{\infty}D_{C}^{\beta}f_{\beta}(s)g_{\alpha}(s;t)ds=D_{C}^{\alpha}\left(D_{C}^{\beta}f_{\beta}(t)\right).\end{split}$

$\square$

###### Theorem 1

The linear fractional differential equation

$\sum_{k=1}^{n}a_{k}D_{C}^{k+\beta-1}y(t)+a_{n+1}y(t)=F_{\beta}(t)$ (2.4)

such that $y(0)=y_{0}$ and $D_{C}^{\beta+i-1}y(0)=y_{i}$ for $i=1,\ldots,n-1$
and $0<\beta<1$, has a solution given by

$y_{\beta}(t)=\int_{0}^{\infty}z(x)g_{\beta}(x;t)dx=\mathbb{E}\left[z(\mathcal{T}_{\beta}(t))\right].$

Where $z(t)$ is the solution of the linear ordinary differential equation

$\sum_{k=1}^{n}a_{k}z^{(k)}(t)+a_{n+1}z(t)=F(t)$ (2.5)

such that $z^{(k)}(0)=y_{i}$ for $k=0,\ldots,n-1$ and both solutions satisfy
the following condition

$\mathcal{L}\left(y^{(n-1)}(\cdot)\right)(s)=s^{\beta-1}\mathcal{L}\left(z^{(n-1)}(\cdot)\right)\left(s^{\beta}\right)$

for $s\in[0,\infty)$. Where
$F_{\beta}(t)=\int_{0}^{\infty}F(x)g_{\beta}(x;t)dx$ and the process
$\mathcal{T}_{\beta}(t)$ can be seen as the inverse of a $\beta$ stable
subordinator whose density is $g_{\beta}(x,t)$.

Remark. The function $F(x)$ could be found using Laplace transform and Lemma 1
part a.

$\widetilde{F}_{\beta}(s)=s^{\beta-1}\widetilde{F}\left(s^{\beta}\right).$

##### Proof

Taking a $\beta-$Wright type transformation on both sides of equation (2.5),
we obtain

$\sum_{k=1}^{n}a_{k}\int_{0}^{\infty}z^{(k)}(x)g_{\beta}(x;t)dx+a_{n+1}\int_{0}^{\infty}z(x)g_{\beta}(x;t)dx=\int_{0}^{\infty}F(x)g_{\beta}(x;t)dx,$

and applying Lemma 1 we obtain

$\sum_{k=1}^{n}a_{k}D_{C}^{k+\beta-1}y_{\beta}(t)+a_{n+1}D_{C}^{\beta-1}y_{\beta}(t)=F_{\beta}(t)$

as we wanted. $\square$

In the following, let $0<\beta\leq 1$. The best method to find

$f_{\beta}(t)=\int_{0}^{\infty}f(x)g_{\beta}(x;t)dx$

is by using Taylor’s expansion of $f(t)$ about $0$ giving

$f_{\beta}(t)=\sum_{n=0}^{\infty}f^{(n)}(0)\dfrac{t^{\beta
n}}{\Gamma(n\beta+1)}.$ (2.6)

That relationship can be used to find many fractional analogue functions that
can be found in the literature. It will be also seen through the following
examples.

######  Example 1 (Fractional velocity)

The solution of the FDE $D_{C}^{\beta}y_{\beta}(t)=c$ with
$y_{\beta}(0)=y_{0}$, where $c$ is a real-valued constant is given by

$y_{\beta}(t)=\int_{0}^{\infty}(y_{0}+cx)g_{\beta}(x;t)dx=y_{0}+\dfrac{ct^{\beta}}{\Gamma(\beta+1)}$

since $y(t)=y_{0}+ct$ solves $Dy(t)=c$, with $y(0)=y_{0}$, this solution
agrees with the literature [41]. That is,
$y_{\beta}(t)=\mathbb{E}[y_{0}+c\mathcal{T}_{\beta}(t)]$.

######  Example 2 (Fractional growth/decay models)

The solution of the FDE $D_{C}^{\beta}y(t)=\lambda y(t)$ with $y(0)=y_{0}$
where $\lambda$ is a real-valued constant is given by

$y(t)=\int_{0}^{\infty}y_{0}e^{\lambda
x}g_{\beta}(x;t)dx=y_{0}E_{\beta}\left(\lambda t^{\beta}\right)$

since $z(x)=y_{0}e^{\lambda x}$ solves $Dz(x)=\lambda z(x)$, with
$z(0)=y_{0}$, this solution agrees with the literature [41]. That is,
$y_{\beta}(t)=\mathbb{E}[y_{0}e^{\lambda\mathcal{T}_{\beta}(t)}]$. See section
5 for graphical representation.

######  Example 3 (Fractional oscillations)

The solution of the FDE $D_{C}^{\beta+1}y(t)=-\omega^{2}y(t)$ with $y(0)=0$
and $y^{\prime}(0)=1$, where $\omega$ is a real-valued constant is given by

$y(t)=\int_{0}^{\infty}\sin(\omega x)g_{\beta}(x;t)dx=\sin_{\beta}\left(\omega
t\right)$

since $z(x)=\sin(\omega x)$ solves $D^{2}z(x)=-\omega^{2}z(x)$, with $z(0)=0$.
That is, $y_{\beta}(t)=\mathbb{E}[\sin(\omega\mathcal{T}_{\beta}(t))]$.

Using (2.6) the fractional analogue of sine and cosine functions are given by

$\sin_{\beta}\left(t\right)=\sum_{n=0}^{\infty}(-1)^{n}\dfrac{t^{(2n+1)\beta}}{\Gamma((2n+1)\beta+1)}=t^{\beta}E_{2\beta,\beta+1}\left(-t^{2\beta}\right)$

and,

$\cos_{\beta}\left(t\right)=\sum_{n=0}^{\infty}(-1)^{n}\dfrac{t^{2n\beta}}{\Gamma(2n\beta+1)}=E_{2\beta}\left(-t^{2\beta}\right),$

which agrees with [42]. It is not hard then to see that
$D_{C}^{\beta}\cos_{\beta}\left(t\right)=-\sin_{\beta}\left(t\right)$ and
$D_{C}^{\beta}\sin_{\beta}\left(t\right)=\cos_{\beta}\left(t\right)$.

##  3 Green functions for fractional partial differential equations

In the following we will use the notation ${}_{C}D_{t}^{\beta}$ and
${}_{RF}D_{x}^{\alpha}$ for Caputo and Riesz-Feller derivatives, respectively,
to identify the variable with respect to which the derivatives are calculated.

###  3.1 The fractional diffusion equation

In this example, using the approach of the previous section, we find the Green
function for the fractional diffusion equation

$\left\\{\begin{split}{}_{C}D_{t}^{\beta}u(x,t)&=\,_{RF}D_{x}^{\alpha,\theta}u(x,t)\\\
u(x,0)&=f(x)\end{split}\qquad\begin{split}(x,t)&\in\mathbb{R}\times[0,\infty)\\\
x&\in\mathbb{R}\end{split}\right.$ (3.1)

with $0<\beta\leq 1$, $0<\alpha\leq 2$ and $u(\pm\infty,t)=0$, $t>0$. Applying
the Fourier transform to $x$, we obtain

${}_{C}D_{t}^{\beta}\widehat{u}(k,t)$
$\displaystyle=-\psi_{\alpha}^{\theta}(k)\widehat{u}(k,t)$
$\displaystyle\widehat{u}(k,0)$ $\displaystyle=\widehat{f}(k).$

Solving this ordinary differential equation using Lemma 1, we obtain

$\widehat{u}(k,t)=\int_{0}^{\infty}\widehat{f}(k)e^{-\psi_{\alpha}^{\theta}(k)x}g_{\beta}(x;t)dx=\widehat{f}(k)E_{\beta}\left(-\psi_{\alpha}^{\theta}(k)t^{\beta}\right),$

see example 2.1. Using the Inverse Fourier Theorem and the convolution theorem
we get

$u(x,t)=\int_{-\infty}^{\infty}G_{\alpha,\beta}^{\theta}(x-y,t)f(y)dy$

where the Green function is given by

$G_{\alpha,\beta}^{\theta}(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikx}E_{\beta}\left(-\psi_{\alpha}^{\theta}(k)t^{\beta}\right)dk.$

In the particular case of $\theta=0$,

$G_{\alpha,\beta}^{0}(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikx}E_{\beta}\left(-\lvert
k\rvert^{\alpha}t^{\beta}\right)dk.$ (3.2)

For the case $\beta=1,\alpha=2$, we obtain the classical heat equation

$\left\\{\begin{split}\partial_{t}u(x,t)&=\partial_{x}^{2}u(x,t)\\\
u(x,0)&=f(x).\end{split}\qquad\begin{split}(x,t)&\in\mathbb{R}\times[0,\infty)\\\
x&\in\mathbb{R}\end{split}\right.$

Since $E_{1}(\lambda t)=e^{\lambda t}$, we obtain for $x\in\mathbb{R}$ and
$t>0$

$G_{1,2}^{0}(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}e^{ikx-
tk^{2}}dk=\dfrac{1}{\sqrt{4\pi t}}e^{-\frac{x^{2}}{4t}}$

or $L_{2}^{0}(x,t)$, by using standard formulas from Fourier transform. This
result agrees with the Green function of the standard heat equation, see for
example [43].

###### Corollary 1

The solution of the fractional heat equation (3.1) could be written as

$u_{\alpha,\beta}(x,t)=\int_{0}^{\infty}\int_{0}^{\infty}L_{2}^{0}(x,\tau)L_{\alpha/2}^{-\alpha/2}(\tau,s)g_{\beta}(s,t)d\tau
ds,$ (3.3)

where $L_{2}^{0}(x,\tau)$ is the solution of the classical heat equation.

##### Proof

As was shown in [28], it could be shown that the solution in (3.2) could be
rewritten as

$u_{\alpha,\beta}(x,t)=\int_{0}^{\infty}L_{\alpha}^{0}(x,s)g_{\beta}(s,t)ds,$

but for $0<\alpha/2\leq 1$

$L_{\alpha}^{0}(x,s)=\int_{0}^{\infty}L_{2}^{0}(x,\tau)L_{\alpha/2}^{-\alpha/2}(\tau,s)d\tau,$
(3.4)

by equation (1.15). Then the result follows. $\square$

Based on this corollary,
$u_{\alpha,\beta}(x,t)=\mathbb{E}\left[L_{2}^{0}(x,\mathcal{L}_{\alpha/2}^{-\alpha/2}(\mathcal{T}_{\beta}(t)))\right]$
where $\mathcal{L}_{\alpha/2}^{-\alpha/2}(t)$ is the $\alpha/2-$Lévy process
with $\theta=-\alpha/2$ and $\mathcal{T}_{\beta}(t)$ is the random
$\beta-$time process. See section 5 for graphical representation.

The moments of $X_{\alpha,\beta}(t)\sim u_{\alpha,\beta}(\cdot,t)$ are given
by [39],

$\mathbb{E}(|X_{\alpha,\beta}(t)|^{s})=t^{\frac{s\beta}{\alpha}}\frac{\Gamma(1-\frac{s}{\alpha})\Gamma(1+\frac{s}{\alpha})\Gamma(1+s)}{\Gamma(1-\frac{s}{2})\Gamma(1+\frac{s}{2})\Gamma(1+\frac{s\beta}{\alpha})}$
(3.5)

for $-\min(1,\alpha)<\mathcal{R}(s)<\alpha$. That formula shows different
regimes of diffusion, subdiffusion and superdiffusion based on whether
$\frac{\alpha}{2}=\beta$, $\frac{\alpha}{2}>\beta$, or
$\frac{\alpha}{2}<\beta$, respectively.

Moreover, Equation (3.3) shows that the $\alpha-\beta$ fractional process
$\\{X_{\alpha,\beta}(t),t\geq 0\\}$ is equivalent in distribution to the
subordinated process
$\\{B(\mathcal{L}_{\alpha/2}^{-\alpha/2}(\mathcal{T}_{\beta}(t))),t\geq 0\\}$,
where $\\{B(t),t\geq 0\\}$ is a Brownian motion. That relationship postulates
that Lévy flights are random dilation, with probability distribution
$L_{\alpha/2}^{-\alpha/2}$, of the standard deviation or the time parameter in
the Brownian motion. That dilation results in an expand in the range of
possible displacement by magnitude beyond the regular tails of the standard
Gaussian distribution. The process
$\\{B(\mathcal{L}_{\alpha/2}^{-\alpha/2}(t)),t\geq 0\\}$ was introduced in
[44], and used in [45] to model stock price differences. See also more about
those subordinated processes in [46].

###  3.2 Fractional wave equation

If we consider d’Alembert solution $u(x,t)=(f(x-kt)+f(x+kt))/2$ of the wave
equation

$\frac{\partial^{2}u(x,t)}{\partial
t^{2}}=k^{2}\frac{\partial^{2}u(x,t)}{\partial x^{2}},$ (3.6)

where $f$ is a smooth enough function. And if we take the $\beta-$ Wright type
transformation on both sides of equation (3.6), we obtain

$\int_{0}^{\infty}\dfrac{\partial^{2}u(x,s)}{\partial
s^{2}}g_{\beta}(s;t)ds=k^{2}\int_{0}^{\infty}\dfrac{\partial^{2}}{\partial
x^{2}}u(x,s)g_{\beta}(s;t)ds=k^{2}\dfrac{\partial^{2}}{\partial
x^{2}}\int_{0}^{\infty}u(x,s)g_{\beta}(s;t)ds$

by Lemma (1) we see that

$u_{\beta}(x,t)=1/2\int_{0}^{\infty}u(x,s)g_{\beta}(s;t)ds=\int_{0}^{\infty}[f(x-ks)+f(x+ks)]g_{\beta}(s;t)ds=(f_{\beta}(x-kt)+f_{\beta}(x+kt))/2$
(3.7)

is a solution of the fractional wave equation

${}_{C}D_{t}^{\beta+1}u_{\beta}(x,t)=k^{2}\dfrac{\partial^{2}u_{\beta}}{\partial
x^{2}}(x,t).$

We refer to the reader to [47], [48] and [49], for previous work on fractional
d’Alambert solution.

###  3.3 Examples of Green functions for fractional Schrödinger equation.

We show in this section that we can apply Lemma (1) to Schrödinger equation.
Since the work of [34] the connection with inverse stable subordinators and
Wright type functions was made for difussion equations. These examples might
raise the interest of specialists inverse stable subordinators to find similar
results for dispersion equations.

One of the classical one dimensional fractional Schrödinger equations with
Caputo derivative in time considered in the literature is given by

$i{}_{0}D_{t}^{\beta}\phi(x,t)=k\frac{\partial^{2}\phi(x,t)}{\partial
x^{2}}+c\phi(x,t),$ (3.8)

Given that when the skewness $\theta=0$, we have
${}_{x}D_{0}^{\alpha}\phi=\frac{\partial^{2}\phi(x,t)}{\partial x^{2}}$, we
consider a more general version

$\left\\{\begin{aligned}
\begin{split}i{}_{0}D_{t}^{\beta}\phi_{\beta,\alpha}(x,t)&=-h{}_{x}D_{\theta}^{\alpha+1}\phi_{\beta,\alpha}(x,t)+c\phi_{\beta,\alpha}\\\
\phi_{\beta,\alpha}(x,0)&=f(x)\end{split}\qquad\begin{split}(x,t)&\in\mathbb{R}\times[0,\infty)\\\
x&\in\mathbb{R}\end{split}\end{aligned}\right.$

Let $0<\beta,\alpha\leq 1$. The Green function of the time and space-
fractional linear Schrödinger equation is as follows. The linear Schrödinger
Fourier transform is

$\left\\{\begin{aligned}
{}_{0}D_{t}^{\beta}\widehat{\phi_{\beta,\alpha}}(x,t)&=-i\widehat{\phi_{\beta,\alpha}}(k,t)\left[h\psi^{\theta}_{\alpha+1}(k)+ic\right]\\\
\widehat{\phi_{\beta,\alpha}}(k,0)&=F(k)\end{aligned}\right.$

Then, applying lemma 3, the solution is given by
$\widehat{\phi_{\beta,\alpha}}(k,t)=\int_{0}^{\infty}z(x)g_{\beta}(x;t)dx$
where the function
$z(t)=F(k)\exp\left(-it\left[h\psi_{\alpha+1}^{\theta}(k)+c\right]\right)$ and
solves the last Fractional PDE initial problem. Now, setting the new solution
we have that

$\widehat{\phi_{\beta,\alpha}}(k,t)=F(k)E_{\beta}\left(-it^{\beta}\left[h\psi_{\alpha+1}^{\theta}(k)+c\right]\right).$

Next, applying the Fourier transform inverse and convolution theorems lead us
to

$\phi(x,t)=\mathcal{F}^{-1}\left\\{F(k)E_{\beta}\left(-it^{\beta}\left[h\psi^{\theta}_{\alpha+1}(k)+c\right]\right)\right\\}=\int_{-\infty}^{\infty}G(x-y,t)f(y)dy$

where the Green function is given by

$G(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}\exp(ikx)E_{\beta}\left(-it^{\beta}\left[h\psi^{\theta}_{\alpha+1}(k)+c\right]\right)dk$

In the particular case $\theta=0$, and $\beta=\alpha=1$ it results to

$G(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}\exp(-ithk^{2}+ixk-itc).$

And by complex integration see [32] we obtain

$\displaystyle G(x,t)=\sqrt{\dfrac{1}{4i\pi
th}}\exp\left(\dfrac{-x^{2}}{4ith}-itc\right).$

Similarly, we can consider the fractional Schrödinger equation given by

$i{}_{0}D_{t}^{\beta}\phi(x,t)=k\frac{\partial^{2}\phi(x,t)}{\partial
x^{2}}+\frac{\partial\phi(x,t)}{\partial x},$ (3.9)

we find the Green function for the more general fractional Schrödinger
equation

$\left\\{\begin{aligned}
\begin{split}{}_{0}D_{t}^{\beta}\phi_{\beta,\alpha}(x,t)&=ih{}_{x}D_{\theta}^{\alpha+1}\phi_{\beta,\alpha}(x,t)-ic{}_{x}D_{\theta}^{\alpha}\phi_{\beta,\alpha}\\\
\phi_{\beta,\alpha}(x,0)&=f(x)\end{split}\qquad\begin{split}(x,t)&\in\mathbb{R}\times[0,\infty)\\\
x&\in\mathbb{R}\end{split}\end{aligned}\right.$

After taking the Fourier transform we obtain

$\left\\{\begin{aligned}
{}_{0}D_{t}^{\beta}\widehat{\phi_{\beta,\alpha}}(x,t)&=-i\widehat{\phi_{\beta,\alpha}}(k,t)\left[h\psi^{\theta}_{\alpha+1}(k)-c\psi^{\theta}_{\alpha}(k)\right]\\\
\widehat{\phi_{\beta,\alpha}}(k,0)&=F(k)\end{aligned}\right.$

Then, applying lemma 3, the solution is given by
$\widehat{\phi_{\beta,\alpha}}(k,t)=\int_{0}^{\infty}z(x)g_{\beta}(x;t)dx,$where
the function
$z(x)=F(k)\exp\left(-ix\left[h\psi_{\alpha+1}^{\theta}(k)-c\psi_{\alpha}^{\theta}(k)\right]\right)$
and solves the last Fractional PDE initial problem. Now, setting the new
solution we have that

$\widehat{\phi_{\beta,\alpha}}(k,t)=F(k)E_{\beta}\left(-it^{\beta}\left[h\psi_{\alpha+1}^{\theta}(k)-c\psi_{\alpha}^{\theta}(k)\right]\right).$

Next, applying the inverse Fourier transform and convolution theorems lead us
to

$\phi(x,t)=\mathcal{F}^{-1}\left\\{F(k)E_{\beta}\left(-it^{\beta}\left[h\psi^{\theta}_{\alpha+1}(k)-c\psi^{\theta}_{\alpha}(k)\right]\right)\right\\}=\int_{-\infty}^{\infty}G(x-y,t)f(y)dy$

where the Green function is given by

$G(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}\exp(ikx)E_{\beta}\left(-it^{\beta}\left[h\psi^{\theta}_{\alpha+1}(k)-c\psi^{\theta}_{\alpha}(k)\right]\right)dk$

In the particular case $\theta=0$, and $\beta=\alpha=1$ it results to

$G(x,t)=\dfrac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(-ithk^{2}+itc|k|+ixk\right)dk.$

For the Green functions for the regular versions of Schrödinger equations we
refer the reader to [32].

##  4 Green functions for fractional PDEs with variable coefficients

In this section we use Lie symmetry group methods and Lemma 1 to introduce
families of fractional partial differential equations with variabe
coefficients exhibiting explicit solutions. More specifically, we show that
solutions for FDE with variable coefficients of the form

${}_{C}D_{t}^{\beta}u_{\beta}(x,t)=\sigma(x)\dfrac{\partial
u_{\beta}}{\partial x^{2}}(x,t)+\mu(x)\dfrac{\partial u_{\beta}}{\partial
x}(x,t),0<\beta\leq 1.$ (4.1)

can be expressed as a Wright type transformation for PDEs that we define as

$u_{\beta}(x,t)=\int_{0}^{\infty}u(x,s)g_{\beta}(s,t)ds,$ (4.2)

where where $u(x,t)$ is the solution of the associated standard PDE

$u_{t}(x,t)=\sigma(x)u_{xx}(x,t)+\mu(x)u_{x}(x,t).$ (4.3)

###### Theorem 2

Let’s consider the fractional partial differential equation with a variable
coefficient of the form (4.1). The solution is given by the Wright-type
transformation define by (4.2) of the analogous PDE of the form (4.3) which
admits a Green function. Furthermore,

$\lim_{\beta\rightarrow 1^{-1}}u_{\beta}(x,t)=u(x,t).$

##### Proof

Taking the $\beta-$ Wright type transformation on both sides of (4.1), we get

$\displaystyle\int_{0}^{\infty}\dfrac{\partial u(x,s)}{\partial
s}g_{\beta}(s;t)ds$
$\displaystyle=\sigma(x)\int_{0}^{\infty}u_{xx}(x,s)g_{\beta}(s;t)ds$
$\displaystyle+\mu(x)\int_{0}^{\infty}u_{x}(x,s)g_{\beta}(s;t)ds.$

By Lemma 1, we obtain

${}_{C}D_{t}^{\beta}u_{\beta}(x,t)$
$\displaystyle=\sigma(x)\dfrac{\partial^{2}}{\partial
x^{2}}\int_{0}^{\infty}u(x,s)g_{\beta}(s;t)ds$
$\displaystyle+\mu(x)\dfrac{\partial}{\partial
x}\int_{0}^{\infty}u(x,s)g_{\beta}(s;t)ds.$

Finally, we obtain

${}_{C}D_{t}^{\beta}u_{\beta}(x,t)=\sigma(x)\dfrac{\partial^{2}u_{\beta}(x,t)}{\partial
x^{2}}+\mu(x)\dfrac{\partial u_{\beta}(x,t)}{\partial x}.$

For $0<v<1$, we recall that

$\lim_{v\rightarrow 1^{-}}W_{-v,1-v}(-z)=\delta(x-1),$

therefore, we get

$\lim_{v\rightarrow 1^{-}}g_{v}(s;t)=\lim_{v\rightarrow
1^{-}}t^{-v}W_{-v,1-v}(-st^{-v})=\frac{1}{t}\delta(st^{-1}-1)=\delta(s-t)$

Therefore, if we assume $u(x,t)=\int_{0}^{\infty}G(x,y,t)f(y)dy,$ $u$ has a
Green function, using Fubini-Tonelli’s theorem and assuming $f$ is a suitable
function, then

$\displaystyle\lim_{\beta\rightarrow 1^{-}}u_{\beta}(x,t)$ $\displaystyle=$
$\displaystyle\lim_{\beta\rightarrow
1^{-}}\int_{0}^{\infty}g_{\beta}(s;t)u(x,s)ds=\lim_{\beta\rightarrow
1^{-}}\int_{0}^{\infty}g_{\beta}(s;t)\int_{0}^{\infty}G(x,y,s)f(y)dyds$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\delta(s-t)G(x,y,s)f(y)dy=\int_{0}^{\infty}G(x,y,t)f(y)dy.$

$\square$ Next, we will present the following families, exhibiting explicit
solutions thanks to the Theorem above.

###  4.1 Solutions for a fractional Fokker-Planck equation with a forcing
function.

Let $\mathcal{T}_{\beta}(t)$ for all $t>0$ be a random time process that is an
inverse of a $\beta$ stable subordinator. Let it also be independent of a
standard Brownian motion $W(t)$. The Langevin stochastic differential equation
with continuous time change using $\mathcal{T}_{\beta}(t)$ is given by

$dX(t)=-f(X(t))d\mathcal{T}_{\beta}(t)+\sqrt{2D}dW(\mathcal{T}_{\beta}(t))$
(4.4)

such that $X(0)=x_{0}$ and $f$ is a Lipschitz function, [33]. The solution
process $X(t)$ of the time-changed Langevin equation (4.4) is a non-stationary
stochastic process and is completely specified by finding the probability
density $p(x,t|x_{0})\geq 0$, which satisfies the fractional Fokker-Planck
equation

${}_{C}D_{t}^{\beta}p$ $\displaystyle=$ $\displaystyle
D\frac{\partial^{2}p}{\partial x^{2}}+\frac{\partial}{\partial
x}\left[f(x)p\right]$ $\displaystyle p(x,0|x_{0})$ $\displaystyle=$
$\displaystyle\delta(x-x_{0}).$

Applying the Theorem of this Section, we obtain the following Corollaries.

###### Corollary 2

If $f$ satisfies the Riccati equation of the form

$2f^{\prime}(x)-f^{2}(x)+\beta^{2}x^{2}-\gamma+\dfrac{16v^{2}-1}{x^{2}}=0$

and $\beta,$ $\gamma$ and $v$ are constants, then the family of the fractional
Fokker-Planck equation

${}_{C}D_{t}^{\beta}u(x,t)=\dfrac{\partial^{2}u}{\partial
x^{2}}+\dfrac{\partial}{\partial x}\left[f(x)u(x,t)\right]$ (4.5)

where $f(x)=-f(-x)$ ($f$ is an odd function) admits an explicit solution of
the form (4.2) where $u$ is given by

$u(x,t)=\int_{0}^{\infty}p(x,t|y)\varphi(y)dy,$ (4.6)

$p$ is given by

$\displaystyle p(x,t|x_{0})$ $\displaystyle=$ $\displaystyle
F\left(\frac{x}{\sqrt{4\sinh^{2}(\beta t)}}\right)\frac{e^{\gamma
t/4}}{\left(4\sinh^{2}(\beta t)\right)^{1/4}}$
$\displaystyle\times\exp\left[-\left(\frac{\beta}{4}\coth(\beta
t)x^{2}+\dfrac{1}{2}\int f(x)dx\right)+\frac{x_{0}^{2}\beta}{2(1-e^{2\beta
t})}\right],$

and

$F(z)=\left\\{\begin{array}[]{cc}z^{1/2}\left[A_{1}I_{2v}(kz)+A_{2}I_{-2v}(kz)\right],&\text{for}\enspace
z>0\\\ \lvert z\rvert^{1/2}\left[B_{1}K_{2v}(k\lvert
z\rvert)+B_{2}I_{2v}(k\lvert z\rvert)\right],&\text{for}\enspace
x<0\end{array}\right.$

where $k=\beta x_{0}$ and $A_{1},$ $A_{2},$ $B_{1}$ and $B_{2}$ are arbitrary
constants be determined by boundary and continuity conditions and if $\alpha$
is a real number

$\displaystyle I_{\alpha}(z)$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{\infty}\frac{\left(z/2\right)^{2m+\alpha}}{m!\Gamma(m+\alpha+1)},$
$\displaystyle K_{\alpha}(z)$ $\displaystyle=$
$\displaystyle\frac{\pi}{2}\frac{I_{-\alpha}(z)-I_{\alpha}(z)}{\sin(\alpha\pi)}.$

##### Proof

This Corollary is a direct consequence of the Theorem of this section and the
similarity solution presented in [22] using the method of Lie group symmetries
for the equation

$\displaystyle\dfrac{\partial p}{\partial t}$
$\displaystyle=\dfrac{\partial^{2}p}{\partial x^{2}}+\dfrac{\partial}{\partial
x}\left[f(x)p\right]$ (4.7) $\displaystyle p(x,0|x_{0})$
$\displaystyle=\phi(x-x_{0})$ (4.8)

when $f$ is odd. $\square$

###### Corollary 3

A particular case of interest is obtained if

$f(x)=ax+\frac{b}{x},\text{ such that }a>0\text{ and }-\infty<b<1.$

The transition probability density is

$\displaystyle p(x,t|x_{0})$ $\displaystyle=$ $\displaystyle
ax_{0}^{1/2}\left(\frac{x}{x_{0}}\right)^{-\frac{b}{2}}z^{\frac{1}{2}}I_{-\left(\frac{1}{2}+\frac{b}{2}\right)}(kz)$
$\displaystyle\times\frac{e^{\gamma
t/4}}{\left(4\sinh^{2}(at)\right)^{1/4}}\exp\left(-\left(\frac{a}{4}\coth(at)\right)^{2}-\dfrac{ax^{2}}{4}\right)$

for $x\geq 0$.

###### Example 4

Let’s consider the standard fractional Fokker-Planck equation of the form

${}_{C}D_{t}^{\beta}u(x,t)$ $\displaystyle=\frac{\partial^{2}u}{\partial
x^{2}}+x\frac{\partial u}{\partial x}+u$ (4.9) $\displaystyle u(x,0)$
$\displaystyle=f(x).$ (4.10)

By the Theorem of this Section, we obtain (4.2) where $p(x,t)$ satisfies the
classical Fokker-Planck equation (FPE)

$\dfrac{\partial p}{\partial t}=\dfrac{\partial^{2}p}{\partial
x^{2}}+x\dfrac{\partial p}{\partial x}+p.$

It is possible to find the Green function of FPE, see for example [43].
Therefore we obtain explicit expression for the solution

$u(x,t)=\int_{0}^{\infty}\int_{0}^{\infty}\frac{\exp\left[-\frac{\left(x-e^{-s}y\right)^{2}}{2\left(1-e^{-2s}\right)}\right]}{\sqrt{2\pi\left(1-e^{-2s}\right)}}g_{\beta}(s;t)dsf(y)dy.$

###### Corollary 4

The family of fractional PDEs with space variable coefficients of the form

${}_{C}D_{t}^{\beta}u(x,t)$ $\displaystyle=xu_{xx}(x,t)+f(x)u_{x}(x,t)$ (4.11)
$\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$ (4.12)

where $f$ satisfies the Riccati equation of the form

$xf^{\prime}-f+\frac{1}{2}f^{2}=Ax^{\frac{3}{2}}+Cx-\frac{3}{8}$

admits an explicit solution of the form (4.2) where $u$ is given by (4.6) and
$p(x,y,t)$ is the inverse Laplace transform of

$\displaystyle U_{\lambda}(x,t)$
$\displaystyle=\sqrt{\frac{\sqrt{x}\left(1+\lambda
t\right)}{\sqrt{x}\left(1+\lambda
t\right)-A\lambda\frac{t^{3}}{12}}}\exp[S(x,y,t)]$
$\displaystyle\times\exp\left[-\frac{1}{2}\left(F(x)-F\left(\frac{(12(1+\lambda
t)\sqrt{x}-A\lambda t^{3})^{2}}{144\left(1+\lambda
t\right)^{4}}\right)\right)\right]$

where $F^{\prime}(x)=\frac{f(x)}{x}$,

$S(\lambda,x,t)=-\dfrac{\lambda(x+Ct^{2}/2)}{1+\lambda
t}-\dfrac{2At^{2}\sqrt{x}(3+\lambda t)}{3(1+\lambda
t)^{2}}+\dfrac{A^{2}t^{4}(2\lambda t(3+\lambda t/2)-3)}{108(1+\lambda t)^{3}}$

for $\lambda\geq 0.$

##### Proof

This Corollary is a direct consequence of the Theorem of this section and the
similarity solution presented in Theorem 6.1 on [54] for

$\displaystyle\dfrac{\partial u}{\partial t}$
$\displaystyle=x\dfrac{\partial^{2}u}{\partial x^{2}}+f(x)\dfrac{\partial
u}{\partial x}$ (4.13) $\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$
(4.14)

when $f$ is an odd function. $\square$

###### Corollary 5

The family of fractional PDEs with space variable coefficients of the form

${}_{C}D_{t}^{\beta}u(x,t)$ $\displaystyle=xu_{xx}(x,t)+f(x)u_{x}(x,t)$ (4.15)
$\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$ (4.16)

where $f$ satisfies the Riccati equation of the form

$xf^{\prime}-f+\dfrac{1}{2}f^{2}=Ax+B$

admits an explicit solution of the form (4.2) where $u$ is given by (4.6) and
$p(x,y.t)$ is the inverse Laplace transform of

$U_{\lambda}(x,t)=\exp\left[-\dfrac{\lambda(x+At^{2}/2)}{1+\lambda
t}-\dfrac{1}{2}\left(F(x)-F\left(\dfrac{x}{\left(1+\lambda
t\right)^{2}}\right)\right)\right]$

where $F^{\prime}(x)=\frac{f(x)}{x}$ and for $\lambda\geq 0$.

##### Proof

This Corollary is a direct consequence of the Theorem of this section and the
similarity solution presented in Theorem 4.1 on [54] for

$\displaystyle\dfrac{\partial u}{\partial t}$
$\displaystyle=x\dfrac{\partial^{2}u}{\partial x^{2}}+f(x)\dfrac{\partial
u}{\partial x}$ (4.17) $\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$
(4.18)

when $f$ is odd. $\square$

###### Example 5

${}_{C}D_{t}^{\beta}u(x,t)$
$\displaystyle=xu_{xx}(x,t)+\left(\dfrac{1+3\sqrt{x}}{2\left(1+\sqrt{x}\right)}\right)u_{x}(x,t)$
(4.19) $\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$ (4.20)

has a solution of the form (4.2) where
$u(x,t)=\int_{0}^{\infty}G(x,y,t)\varphi(y)dy$ and

$G(x,y,t)=\dfrac{\cos\left(\dfrac{2\sqrt{xy}}{t}\right)}{\sqrt{\pi
yt}\left(1+\sqrt{x}\right)}\left(1+\sqrt{y}\tanh\left(\dfrac{2\sqrt{xy}}{t}\right)\right)\exp\left[-\dfrac{x+y}{t}\right].$

###### Example 6

${}_{C}D_{t}^{\beta}u(x,t)$
$\displaystyle=xu_{xx}(x,t)+\left(\dfrac{1}{2}+\sqrt{x}\coth\left(\sqrt{x}\right)\right)u_{x}(x,t)$
(4.21) $\displaystyle u(x,0)$ $\displaystyle=\varphi(x)$ (4.22)

has a solution of the form (4.2) where
$u(x,t)=\int_{0}^{\infty}G(x,y,t)\varphi(y)dy$ and

$G(x,y,t)=\dfrac{\sinh\left(\dfrac{2\sqrt{xy}}{t}\right)}{\sqrt{\pi
yt}}\dfrac{\sinh\left(\sqrt{y}\right)}{\sinh\left(\sqrt{x}\right)}\exp\left[-\dfrac{x+y}{t}-\dfrac{1}{4}t\right].$

###### Theorem 3

A solution for

$_{C}D_{t}^{\beta\alpha}u_{\beta\alpha}(x,t)=\sigma(x)\dfrac{\partial^{2}u_{\beta\alpha}(x,t)}{\partial
x^{2}}+\mu(x)\dfrac{\partial u_{\beta\alpha}(x,t)}{\partial x}$ (4.23)

can be obtained as an $\alpha-$Wright type transformation of $u_{\beta}$,
where $u_{\beta}$ is a solution for

${}_{C}D_{t}^{\beta}u_{\beta}(x,t)=\sigma(x)\dfrac{\partial^{2}u_{\beta}(x,t)}{\partial
x^{2}}+\mu(x)\dfrac{\partial u_{\beta}(x,t)}{\partial x}.$

##### Proof

Let’s consider

${}_{C}D_{t}^{\beta}u_{\beta}(x,t)=\sigma(x)\dfrac{\partial^{2}u_{\beta}(x,t)}{\partial
x^{2}}+\mu(x)\dfrac{\partial u_{\beta}(x,t)}{\partial x}.$

By Lemma 1, it can be written as

$\displaystyle\int_{0}^{\infty}\dfrac{\partial u(x,w)}{\partial
w}g_{\beta}(w;s)dw$ $\displaystyle=\sigma(x)\dfrac{\partial^{2}}{\partial
x^{2}}\int_{0}^{\infty}u(x,w)g_{\beta}(w;s)dw$
$\displaystyle+\mu(x)\dfrac{\partial}{\partial
x}\int_{0}^{\infty}u(x,w)g_{\beta}(w;s)dw.$

Taking a $\alpha$-Wright type transformation on both sides of the equation we
obtain,

$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\dfrac{\partial
u(x,w)}{\partial w}g_{\beta}(w;s)dwg_{\alpha}(s;t)ds$
$\displaystyle=\int_{0}^{\infty}\sigma(x)\dfrac{\partial^{2}}{\partial
x^{2}}\int_{0}^{\infty}u(x,w)g_{\beta}(w;s)dwg_{\alpha}(s;t)ds$
$\displaystyle+\int_{0}^{\infty}\mu(x)\dfrac{\partial}{\partial
x}\int_{0}^{\infty}u(x,w)g_{\beta}(w;s)dwg_{\alpha}(s;t)ds$

and using (1.11) we obtain

$\displaystyle\int_{0}^{\infty}\dfrac{\partial u(x,w)}{\partial
w}g_{\beta\alpha}(w;t)dw$
$\displaystyle=\sigma(x)\dfrac{\partial^{2}}{\partial
x^{2}}\int_{0}^{\infty}u(x,w)g_{\beta\alpha}(w;t)dw$
$\displaystyle+\mu(x)\dfrac{\partial}{\partial
x}\int_{0}^{\infty}u(x,w)g_{\beta\alpha}(w;t)dw,$

which completes the proof. $\square$

##  5 Numerical Simulations

We use Monte Carlo integration to simulate the solutions of factional
differential equations and partial differential equations given their
solutions. As expected, simulations of heavy tail distributions require the
use of a large amount of random numbers for adequate coverage. We first
present the lemma by M. Kanter [55] to generate random numbers distributed as
the Lévy $\alpha$-stable distribution with stability index $0<\alpha\leq 2$,
$L_{\alpha}^{\theta}(x)$, defined by equation (1.12). See also [56].

######  Lemma 2 (Lemma 4.1 [55])

Let $\alpha\in(0,1)$ and let $L_{\alpha}^{-\alpha}(x)$ as defined in (1.12).
Then for $x\geq 0$

$L_{\alpha}^{-\alpha}(x)=\dfrac{1}{\pi}\left(\dfrac{\alpha}{1-\alpha}\right)\left(\dfrac{1}{x}\right)^{(1-\alpha)^{-1}}\int_{0}^{\pi}a(\varphi)\exp\left(-\left(\dfrac{1}{x}\right)^{\alpha/(1-\alpha)}\right)d\varphi$
(5.1)

where

$a(\varphi)=\left(\dfrac{\sin(\alpha\varphi)}{\sin(\varphi)}\right)^{(1-\alpha)^{-1}}\left(\dfrac{\sin((1-\alpha)\varphi)}{\sin(\alpha\varphi)}\right).$
(5.2)

######  Corollary 6 (Corollary 4.1 [55])

Let $U_{1}$ and $U_{2}$ be independent random variables where $U_{1}$ is
uniformly distributed on $[0,1]$, and $U_{2}$ is uniformly distributed on
$[0,\pi]$. Then for $\alpha\in(0,1)$, $L_{\alpha}^{-\alpha}(x)$ is the density
of $(-a(U_{2})/\log(U_{1}))^{(1-\alpha)/\alpha}$ where $a$ is given by
equation (5.2).

###### Corollary 7

If $X\sim{L}_{\beta}^{-\beta}(\cdot)$, then $g_{\beta}(\cdot,t)$ is the
probability density function of $Y=\dfrac{t^{\beta}}{X^{\beta}}$.

##### Proof

For $t>0$, the probability density function of $Y$ is given by

$\displaystyle\begin{split}f_{Y}(y)&={L}_{\beta}^{-\beta}\left(\dfrac{t}{y^{1/\beta}}\right)\left|\dfrac{dx}{dy}\right|\\\
&=\dfrac{t}{\beta
y}\dfrac{1}{y^{1/\beta}}{L}_{\beta}^{-\beta}\left(\dfrac{t}{y^{1/\beta}}\right)\\\
&=\dfrac{t}{\beta y}{L}_{\beta}^{-\beta}(t,y)\\\ &=g_{\beta}(y,t).\end{split}$
(5.3)

$\square$

###### Corollary 8

If $X\sim{L}_{\alpha}^{-\alpha}(\cdot)$, then
${L}_{\alpha}^{-\alpha}(\cdot,t)$ is the probability density function of
$Y=Xt^{1/\alpha}$.

##### Proof

Similar proof as above. $\square$ ${L}_{\alpha}^{-\alpha}(\cdot,t)$ can be
seen as the density of a $\alpha$ stable subordinator, here the skewness
parameter equals -$\alpha$. By using the above corollaries, a powerful
computational approach is performed using $g_{\beta}$ instead of simulating
the Wright function. The codes are done on Python using numba and matplotlib
libraries.

######  Example 7 (Fractional growth/decay models)

The solution of the Fractional differential equation

$\displaystyle\left\\{\begin{split}D_{C}^{\beta}y(t)&=-y(t)\\\
y(0)&=1\end{split}\qquad\right.$ (5.4)

has a solution given by
$y(t)=\int_{0}^{\infty}e^{-x}g_{\beta}(x;t)dx=E_{\beta}(-t^{\beta})$.

Figure 5.1: The Monte Carlo integration repeated 100 times of simulations of
the solution of equation (5.4) for each $\beta=0.5$ compared to the exact
solution.

As a conclusion, the numerical solution converges to the actual solution when
the random number values of $g_{\beta}$ are considerably high. Convergence is,
however, good starting from $10^{4}$ randomly generated values from
$g_{\beta}$.

###### Example 8

Let’s consider the standard fractional Fokker-Planck equation of the form

${}_{C}D_{t}^{\beta}p_{\beta}(x,t)$ $\displaystyle=$
$\displaystyle\frac{\partial^{2}p_{\beta}}{\partial x^{2}}+x\frac{\partial
p_{\beta}}{\partial x}+p_{\beta}$ (5.5) $\displaystyle p(x,0)$
$\displaystyle=$ $\displaystyle f(x).$ (5.6)

We obtained a solution in the form

$p_{\beta}(x,t)=\int_{0}^{\infty}\int_{0}^{\infty}\frac{\exp\left[-\frac{\left(x-e^{-s}y\right)^{2}}{2\left(1-e^{-2s}\right)}\right]}{\sqrt{2\pi\left(1-e^{-2s}\right)}}g_{\beta}(s;t)dsf(y)dy.$
(5.7)

See Figure 5.2 for numerical simulation of the solution and see Figure 5.3 for
the mean absolute error for $N=10^{6}$ with respect to $N=10^{7}$.

|  |
---|---|---
(a) | (b) | (c)
Figure 5.2: Fokker Plank Monte Carlo integration simulation of solution in equation (5.7) using a number of $10^{7}$ random values for $\beta=0.1$, $0.5$ and $0.9$, respectively. |  |
---|---|---
(a) | (b) | (c)
Figure 5.3: Mean error in Fokker Plank Monte Carlo integration simulation of
solution in equation (5.7) over time $t\in[1,10]$. The mean absolute error
between using a number of $10^{7}$ and $10^{6}$ random values from $g_{\beta}$
for $\beta=0.1$, $0.5$ and $0.9$ respectively.

###### Example 9

The solution of the fractional heat/diffusion equation

$u_{\alpha,\beta}(x,t)=\int_{0}^{\infty}\int_{0}^{\infty}L_{2}^{0}(x,\tau)L_{\alpha/2}^{-\alpha/2}(\tau,s)g_{\beta}(s,t)d\tau
ds$ (5.8)

could be numerically computed for each $x$ and $t$ using a large number of
simulations from $g_{\beta}(s,t)$, and for each $s$ simulate a large number of
$\tau$ are simulated from $L_{2}^{0}(x,\tau)$
$L_{\alpha/2}^{-\alpha/2}(\tau,s)$. Then the values of $L_{2}^{0}(x,\tau)$ are
averaged up over all values of $\tau$ and then over all values of $s$. See
Figure 5.4.

A similar procedure could be done to simulate random numbers from
$u_{\alpha,\beta}(x,t)$, but after simulating one $s$ from $g_{\beta}(s,t)$,
and then one $\tau$ from $L_{\alpha/2}^{-\alpha/2}(\tau,s)$, then we simulate
one $x$ from $L_{2}^{0}(x,\tau)$.

|  |
---|---|---
(a) | (b) | (c)
|  |
(d) | (e) | (f)
|  |
(g) | (h) | (i)
Figure 5.4: Simulation of the solution of fractional heat equation
$u_{\alpha,\beta}$ of solution in equation (5.8) for $\beta=0.1$ in (a), (d)
and (g), $\beta=0.5$ in (b), (e), and (h) and $\beta=0.9$ in (c), (f), and (i)
and $\alpha=0.2$ in (a), (b) and (c), $0.6$ in (d), (e) and (f) and $1$ in
(g), (h) and (i) for $(x,t)\in[-5,5]\times[1,10]$. The Monte Carlo integration
was performed using a number of $10^{7}$ random values from $g_{\beta}$.

Linear ordinary differential equations and partial differential equations
could be also solved numerically. Solutions then get interpolated at the
values generated by $g_{\beta}(\cdot,t)$ and averaged at each value of $t$.

###### Example 10

The solution of the fractional differential equation

$D_{C}^{\beta+1}y(t)=-y(t)$ (5.9)

with $y(0)=0$ and $y^{\prime}(0)=1$, could be numerically computed for each
$t$ by numerically solving

$D^{2}z(x)=-z(x),$ (5.10)

with $z(0)=0$, and $z^{\prime}(0)=1$ at $N$ randomly generated $s$ values from
$g_{\beta}(s,t)$ and then average the values up for each $t$. See Figure 5.5.
We used the Runge-Kutta method of hybrid order 4 and 5 from minimum generated
$s$ to maximum generated $s$ and then interpolated the solution at the rest of
the randomly generated $s$ values. The solution is also given by is given by

$\sin_{\beta}\left(t\right)=\int_{0}^{\infty}\sin(x)g_{\beta}(x;t)dx=t^{\beta}E_{2\beta,\beta+1}\left(-t^{2\beta}\right).$

Figure 5.5: Monte Carlo simulation of the solution of fractional differential
equation (5.9), for $\beta=0.5$. The Monte Carlo integration was performed
using a number of $N=10^{7}$ random values from $g_{\beta}$. The figures also
show direct Monte Carlo integration of the sine function and the Mittag-
Leffler function definition of the fractional sine function.

##  6 Conclusion

In this work, we have shown an alternative way to find the fundamental
solutions for fractional partial differential equations (PDEs). Indeed,
Riccati equation have allow us to study families of fractional diffusion PDEs
with variable coefficients which allow explicit solutions. Those solutions
connect Lie symmetries to fractional PDEs. We expect similar results for
fracational dispersive equations. These results will appear in another work.

In our approach, we have taken advantage of Lie symmetries applied to
fractional diffusion PDEs with variable coefficients. We predict that Feynman
path integrals can play a similar role with fractional dispersive equations.

We conjecture that a general solution similar to that in equation (3.3) can be
shown to hold for a larger class of fractional partial differential equations.
A conjecture for which an interesting consequence that the Lévy flight could
be a result of the wide expanse allowed for the normal diffusion to make a
jump over.

Monte-Carlo integration of solutions of ordinary differential equations and
partial differential equations with respect to heavy-tailed distributions is a
new approach that can prove valuable for other general fractional equations.
Evaluating such solutions take small amount of time thanks to the availability
of fast computing devices. A general formula like equation (3.3) can make
numerical solutions easier.

## References

  * [1] A. Consiglio and F. Mainardi, On the evolution of fractional diffusive waves, Ricerche di Matematica, 2019.
  * [2] M. R. Islam, A. Peace, D. Medina, and T. Oraby, Integer versus fractional order SEIR deterministic and stochastic models of measles, International Journal of Environmental Research and Public Health, 2020.
  * [3] R. Almeida, A. B. da Cruz, N. Martins, and M. T. T. Monteiro, An epidemiological MSEIR model described by the Caputo fractional derivative, International Journal of Dynamics and Control, pp. 1–14, 2018.
  * [4] M. Stynes, Fractional-order derivatives defined by continuous kernels are too restrictive, Appl. Math. Lett., vol. 85, pp. 22–26, nov 2018.
  * [5] K. Diethelm, R. Garrappa, A. Giusti, and M. Stynes, Why Fractional Derivatives with Nonsingular Kernels Should Not Be Used, Fractional Calculus and Applied Analysis 2020 23:3, vol. 23, pp. 610–634, jul 2020.
  * [6] Dokuyucu and A. Mustafa, On the Fractional Derivative and Integral Operators, Fractional Order Analysis, pp. 1–41, sep 2020.
  * [7] R. R. Nigmatullin and Y. E. Ryabov, Renewal processes of Mittag-Leffler and wright type, Fiz. Tverd. Tela., vol. 39, p. 101, 1997.
  * [8] K. Weron and A. Klauser, Probabilistic basis for the Cole-Cole relaxation law, Current Developments in Mathematics, vol. 236, p. 59-69, 2000.
  * [9] R. Metzler, E. Barkai, and J. Klafter, Anomalous diffusion and relaxation close to thermal equilibrium: A fractional Fokker-Planck equation approach, Phys. Rev. Lett., vol. 82, p. 3563., 1999.
  * [10] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep., vol. 339, p. 1., 2000.
  * [11] M. F. Shlesinger and G. M. Zaslavsky and J. Klafter, Strange kinetics, Nature, vol. 363, p. 31., 1993.
  * [12] A. I. Saichev and G. M. Zaslavsky, Fractional kinetic equations: solutions and applications, Chaos, vol. 7, p. 753, 1997.
  * [13] J. Klafter, A. Blumen, and M. Shlesinger, Stochastic pathway to anomalous diffusion, Physical Review A, vol. 35, pp. 3081–3085, apr 1987.
  * [14] I. Sokolov and J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos, vol. 15, no. 2, pp. 1–18, 2005\.
  * [15] V.Yanovsky, A.Chechkin, D. Schertzer, and A.Tur, Lévy anomalous diffusion and fractional Fokker-Planck equation, Physica A: Statistical Mechanics and its Applications, vol. 282, pp. 13–34, Jul 2000.
  * [16] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Physics Reports, vol. 339, pp. 1–77, dec 2000.
  * [17] R. Metzler, J. H. Jeon, A. Cherstvy, and E. Barkai, Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking, Physical Chemistry Chemical Physics, vol. 16, pp. 24128–24164, Oct 2014.
  * [18] E. Orsingher and F. Polito, Fractional pure birth processes, Bernoulli, vol. 16, no. 3, pp. 858–881, 2010.
  * [19] E. Orsingher, C. Ricciuti, and B. Toaldo, Population models at stochastic times, Advances in Applied Probability, vol. 48, no. 2, pp. 481–498, 2016.
  * [20] S. Eule and R. Friedrich, A note on the forced Burgers equation, Physics Letters A: General, Atomic and Solid State Physics, vol. 351, pp. 234 – 241., 2006.
  * [21] Z. Feng, Traveling wave behavior for a generalized Fisher equation, Chaos Solitons Fractals, vol. 38, p. 481â€“488., 2008.
  * [22] G. W. Bluman, Similarity solutions of the one dimensional Fokker-Planck equation, Int. J. Non-Lin. Mech., vol. 6, pp. 143–153., 1971.
  * [23] G. Bluman and S. Kumei, Symmetries and differential equations, Springer-Verlag, p. New York., 1989.
  * [24] G. Bluman and S. Anco, Symmetry and integration methods for differential equations, Applied Mathematical Sciences, Springer-Verlag, p. New York., 2002.
  * [25] G. Bluman and J. D. Cole, The general similarity solution of the heat equation, J. Math. Mech., vol. 18, p. 1025-1042, 1969.
  * [26] P. Olver, Application of Lie groups to differential equations, Springer-Verlag, p. New York., 1986.
  * [27] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications.  Academic Press, 1999.
  * [28] F. Mainardi, Y. Luchko, and G. Pagnini, The fundamental solution of the space-time fractional diffusion equation, An International Journal for Theory and Applications, vol. 4 No 2, pp. 153 – 192., 2001.
  * [29] F. Mainardi, R. Glorenflo, and A. Vivoli, Renewal processes of Mittag-Leffler and wright type, Fractional Calculus and Applied Analysis, vol. 8 No 1, pp. 7 – 38., 2001.
  * [30] F. Mainardi and A. Consiglio, The Wright functions of the second kind in mathematical physics, Mathematics 2020, vol. 8 No 6, p. 884., 2020.
  * [31] R. Gorenflo, Mittag-leffler waiting time, power laws, rarefaction, continuous time random walk, diffusion limit, 2010.
  * [32] R Cordero-Soto, RM Lopez, E Suazo, SK Suslov, Propagator of a charged particle with a spin in uniform magnetic and perpendicular electric fields, Letters in Mathematical Physics, 84 (2), 159-178
  * [33] M. Hahn, K. Kobayashi and S. Umarov, SDEs Driven by a Time-Changed Lévy Process and Their Associated Time-Fractional Order Pseudo-Differential Equations, J Theor Probab 25, 262–279 (2012).
  * [34] N.H. Bingham, Limit theorems for occupation times of Markov processes, Z. Warsch. verw. Geb. 17, 1–22, 1971.
  * [35] M. M. Meerschaert and A. Sikorskii, Stochastic Models for Fractional Calculus, De Gruyter, Berlin, 2012.
  * [36] M. M. Meerschaert and P. Straka, Inverse Stable Subordinators, Math Model Nat Phenom, 2013.
  * [37] A. Piryatinska, Inference for the levy models and their applications in medicine and statistical physics, Department of Statistics, Case Western Reserve University, 2005.
  * [38] K. K. Kataria and P. Vellaisamy, On the convolution of Mittag -Leffler distributions and its applications to fractional point processes, Stochastic Analysis and Applications, vol. 37, no. 1, pp. 115–122, 2019.
  * [39] F. Mainardi and G. Pagnini, Mellin-Barnes integrals for stable distributions and their convolutions, Fractional Calculus and Applied Calculus, vol. 11, no. 4, 2008.
  * [40] F. Mainardi, A. Mura, and G. Pagnini, The M-Wright function in time-fractional diffusion processes: A tutorial survey, International Journal of Differential Equations, vol. 2010, 2009.
  * [41] Bangti Jin, Fractional Differential Equations: An Approach via Fractional Derivatives, Applied Mathematical Sciences, 2021
  * [42] A. A. Stanislavsky, Fractional oscillator, Phys. Rev. E, vol. 70, 5, p. 051103, 2004.
  * [43] E. Suazo, S. K. Suslov, and J. M. Vega, The Riccati differential equation and a diffusion-type equation, New York J. Math, vol. 17a, pp. 225 – 224., 2011.
  * [44] B. W. Huff, The Strict Subordination of Differential Processes, Sankhy Ä: The Indian Journal of Statistics, Series A, vol. 31, no. 4, pp. 403–412, 1969.
  * [45] B. Mandelbrot and H. M. . Taylor, On the Distribution of Stock Price Differences, Operations Research, vol. 15, no. 6, pp. 1057–1062, 1967\.
  * [46] O. E. Barndorff-Nielsen and A. Shiryaev, Change of Time and Change of Measure.  Advanced Series on Statistical Science & Applied Probability, WORLD SCIENTIFIC, 2010.
  * [47] Y. Fujita: Integrodifferential equation which interpolates the heat equation and the wave equation, Osaka J. Math. 27 (1990), 309-321.
  * [48] Cheng-Gang Li, Miao Li, Sergey Piskarev, Mark M. Meerschaert, The fractional d’Alembert’s formulas, Journal of Functional Analysis, Volume 277, Issue 12, 2019\.
  * [49] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity. Imperial College Press, London, 2010.
  * [50] E. Orsingher and L. Beghin, Time-fractional telegraph equations and telegraph processes with Brownian time, Probability Theory and Related Fields, vol. 128, pp. 141–160, Jan 2004.
  * [51] F. Huang, Analytical solution for the time-fractional telegraph equation, Journal of Applied Mathematics, vol. 2009, 2009.
  * [52] J. Chen, F. Liu, and V. Anh, Analytical solution for the time-fractional telegraph equation by the method of separating variables, Journal of Mathematical Analysis and Applications, vol. 338, pp. 1364–1377, feb 2008.
  * [53] R. C. Cascaval, E. C. Eckstein, C. Frota, and J. Goldstein, Fractional telegraph equations, Journal of Mathematical Analysis and Applications, vol. 276, pp. 145–159, dec 2002.
  * [54] M. Craddock, Symmetry group methods for fundamental solutions, J. Differential Equations, vol. 207, pp. 285–302., 2004.
  * [55] M. Kanter, Stable densities under change of scale and total variation inequalities, The annals of Probability, vol. 3, No. 4,, pp. 697 – 707., 1975.
  * [56] D. Cahoy, Estimation and Simulation for the M-Wright Function, Communications in Statistics-Theory and Methods, vol. 41, no. 8, pp. 1466–1477, 2012.