#  NUMERICAL SIMULATIONS FOR FRACTIONAL DIFFERENTIAL EQUATIONS OF HIGHER ORDER
AND A WRIGHT-TYPE TRANSFORMATION

M. Nacianceno1, T. Oraby1, H. Rodrigo1, Y. Sepulveda1, J. Sifuentes1, E.
Suazo1, 4, T. Stuck3, J. Williams2

###### Abstract.

In this work, a new relationship is established between the solutions of
higher fractional differential equations and a Wright-type transformation.
Solutions could be interpreted as expected values of functions in a random
time process. As applications, we solve the fractional beam equation,
fractional electric circuits with special functions as external sources, and
derive d’Alembert’s formula for the fractional wave equation. Due to this
relationship, we present two methods for simulating solutions of fractional
differential equations. The two approaches use the interpretation of the
Caputo derivative of a function as a Wright-type transformation of the higher
derivative of the function. In the first approach, we use the Runge-Kutta
method of hybrid orders 4 and 5 to solve ordinary differential equations
combined with the Monte Carlo integration to conduct the Wright-type
transformation. The second method uses a feedforward neural network to
simulate the fractional differential equation.

###### Key words and phrases:

Caputo derivative, fractional differential equations, numerical methods,
neural networks

###### 1991 Mathematics Subject Classification:

Primary: 60H15, 35R60; Secondary: 65M06, 65M50.

  1. (1)

School of Statistical and Mathematical Sciences, The University of Texas, Rio
Grande Valley, 1201 W University Dr, Edinburg, TX 78539

  2. (2)

Department of Mathematics, The University of Minnesota, 511 Bruininks Hall 222
Pleasant Street S.E.. Minneapolis MN 55455.

  3. (3)

Department of Physics, Lamar University, Archer Physics A112 Lamar University
4400 MLK Parkway P.O. Box 10046 Beaumont, Texas 77710

  4. (4)

Corresponding Author, erwin<EMAIL_ADDRESS>
##  1\. Introduction

While initiated in the latter part of the 18th century by Leibniz, Newton, and
l’Hopital, fractional calculus has garnered significant attention in the past
decades. Fractional differential equations (FDEs) involve fractional
derivatives of functions. These equations extend the concept of ordinary
(ODEs) and partial differential equations (PDEs) by allowing the use of
fractional orders of differentiation. Just as differential equations describe
how a function’s derivative relates to the function itself, fractional
differential equations describe how fractional derivatives relate to the
function. Solving fractional differential equations can be challenging due to
the nonlocal and noninteger nature of the derivatives involved. Analytical
solutions are not always possible and numerical methods are often employed to
approximate solutions.

Biological systems often exhibit behaviors that involve long-range and memory
interactions [18]. Fractional differential equations can be employed to model
processes such as enzyme reactions, population dynamics, and the spread of
diseases in a more accurate manner. Fractional calculus can be used to model
long-range memory financial processes, such as volatility in financial
markets. Fractional differential equations can provide more accurate models
for predicting financial fluctuations; for these and more applications see
[7], [19], [18] and references therein.

In this work, we use an alternative approach to solve fractional ordinary and
partial differential equations of higher order as the expected value of a
random time process $\mathcal{T}_{\beta}(t)$. Let $g_{\beta}(x;t)$ be the
probability density function of the random-time process
$\mathcal{T}_{\beta}(t)$ for all $t>0$ [14, 6]. The random time process
$\mathcal{T}_{\beta}(t)$ can be depicted as the inverse of a $\beta$ stable
subordinator with density $g_{\beta}(x,t)$, see [4, 19, 20]. 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^{\beta+1}$ as $t\rightarrow\infty$, with $\beta\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 sub-diffusion
processes [19]. Using the latter, we apply a numerical approach based on Monte
Carlo integration to simulate solutions of fractional ordinary and fractional
partial differential equations.

In addition, we present an alternative version of Lemma 1 presented in [21],
please see Section 3. Lemma 2 tries to establish a connection between
fractional derivatives of higher order and expected values of derivatives of
the same function, see (3.2). However, condition (3.1) is trivial for $n=1$,
but not trivial for $n>1$. Therefore, Lemma 2 in Section 3 provides a weaker
condition. It was inspired by the work of Dr. Mark Meerschaert and his work on
D’Alembert’s formula. Sadly, Dr. Mark Meerschaert passed away recently, but he
contributed tremendously to the theory of fractional calculus.

This paper is organized as follows: In Section 2, we review important results
of fractional calculus and in particular properties for $g_{\beta}(x;t)$ and a
Wright-type transformation fundamental for our work. In Section 3, in the
Lemma 2 a new relationship is established between the solutions of higher
fractional differential equations and the Wright-type transformation. Lemma 2
allows us to solve fractional differential equations of higher orders with
certain initial conditions, see Theorem 1. In Theorem 1 solutions could be
interpreted as expected values of functions in a random-time process. As a
first application, we solve the fractional beam equation. In Section 4, we
continue the applications of Theorem 1, solving fractional electric circuits
with special functions as external sources such as the Mittag-Leffler function
and other special functions. Similar work is also presented for the
homogeneous case. In Section 5, we use the expected value interpretation
presented in Theorem 1 of the solutions of fractional equations to perform
Monte Carlo simulations of their solution. The idea is to generate a new but
random timeline $L(t)$ using $g_{\beta}(\cdot,t)$ for each $t$ and then use
any ODE numerical solver, such as Runge-Kutta, to simulate it on that new
timeline and finally find the average of the solution, see [21]. In Section 6,
we provide a new derivation of d’Alembert’s formula for fractional wave
equations as an application of Lemma 2. Finally, in Section 7, we show how to
use feedforward neural networks to simulate solutions of fractional
differential equations. The main difference between fractional and ordinary
measurements is in the number of past measurements (memory length) required to
predict the measurement of the next step. We found that ODEs require one past
step, while FDEs require more than one (longer memory).

##  2\. On a Wright-type transformation

In this Section, we review important definitions and classical results that we
will need. 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,$ (2.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,$
(2.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}$.

The Wright function is a special function of importance to fractional calculus
and is defined by [13]:

$W_{\beta,\alpha}(z)=\sum\limits_{k=0}^{\infty}\dfrac{z^{k}}{k!\Gamma(\beta
k+\alpha)}.$ (2.3)

$\beta>-1,\alpha\in\mathbb{C},z\in\mathbb{C}$. Also, we will need the
definition of the Mitagg-Leffler function, which is a generalization of the
exponential function.

######  Definition 1.

The Mitagg-Leffler function is defined by

$E_{\beta,\alpha}(z)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma(\beta{k}+\alpha)},\hskip
5.0ptz\in\mathbb{C},\hskip 5.0pt\alpha,\beta>0.$

Next, we introduce the main character of this work, the $g_{\beta}(x;t)$,
which is a Wright-type function.

######  Definition 2.

$g_{\beta}$ is defined as a Wright-type function

$g_{\beta}(x;t)=\frac{1}{t^{\beta}}W_{-\beta,1-\beta}\left(\frac{-x}{t^{\beta}}\right)=\frac{1}{t^{\beta}}\sum_{k=0}^{\infty}\frac{\left(-x\right)^{k}}{t^{\beta
k}k!\Gamma(-\beta(k+1)+1)}$

where $0<\beta\leq 1,\hskip 5.0ptt>0,\hskip 5.0pts\geq 0.$ $g_{\beta}(x;t)$ is
a probability density function of the random time process
$\mathcal{T}_{\beta}(t)$ for all $t>0$ [14, 6]. $\mathcal{T}_{\beta}(t)$ can
be seen as the inverse of a $\beta$ stable subordinator with density
$g_{\beta}(x,t)$, see [4, 19, 20].

We use $g_{\beta}(x;t)$ to define the following transform.

######  Definition 3.

We define a Wright-type transform for a function $f(t)$ as

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

where $t>0,\hskip 5.0pt0<\beta\leq 1.$ We will use the alternative notation
$f_{\beta}(t)=$ $(f(t))_{\beta}.$

The following results are properties of $g_{\beta}(x;t)$. Some of them are
well known, see [18], and others have been recently proved in [21]. We set
these properties as a proposition for the convenience of the reader, and they
will play a fundamental role in our numerical simulations.

######  Proposition 1.

The following properties for $g_{\beta}$, $0<\beta<1$, are valid:

  1. (1)

$g_{\beta}$ is a probability distribution

$\int_{0}^{\infty}g_{\beta}(s;t)ds=1.$

  2. (2)

$g_{\beta}$ has the following Laplace transforms

$\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)$
(2.4)

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$ [23, 10]. At the same time

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

  3. (3)

$g_{\beta}$ has the following moments:

$\int_{0}^{\infty}s^{n}g_{\beta}(s;t)ds=E[S^{n}]=n!\frac{t^{n\beta}}{\Gamma(n\beta+1)}.$

  4. (4)

The following Taylor Expansion for $f_{\beta}(t)$ holds:

$f_{\beta}(t)=\int_{0}^{\infty}f(s)g_{\beta}(s;t)ds=\int_{0}^{\infty}\sum_{n=0}^{\infty}f^{(n)}(0)\frac{s^{n}}{n!}g_{\beta}(s;t)ds=\sum_{n=0}^{\infty}f^{(n)}(0)\frac{t^{n\beta}}{\Gamma(n\beta+1)}.$

From the proposition above we can make a table of the Wright beta transform
for different functions. This table will be useful in Section 4 where we will
solve nonhomogeneous fractional ODEs.

$f(t)$ | $f_{\beta}(t)$
---|---
$t$ | $\frac{t^{\beta}}{\Gamma(\beta+1)}$
$t^{2}$ | $\frac{2t^{2\beta}}{\Gamma(2\beta+1)}$
$t^{n}$ | $\frac{n!t^{n\beta}}{\Gamma(n\beta+1)}$
$\cos(t)$ | $\cos_{\beta}(t)=\sum_{n=0}^{\infty}\frac{(-1)^{n}t^{2n\beta}}{\Gamma(2n\beta+1)}=E_{2\beta}(-t^{2\beta})$
$\sin(t)$ | $\sin_{\beta}(t)=\sum_{n=0}^{\infty}\frac{(-1)^{n}t^{\beta(2n+1)}}{\Gamma((2n+1)\beta+1)}=t^{\beta}E_{2\beta,\beta+1}(-t^{2\beta})$
$e^{t}$ | $\sum_{n=0}^{\infty}\frac{t^{n\beta}}{\Gamma(n\beta+1)}=E_{\beta}(t^{\beta})$
$\cosh(t)$ | $\cosh_{\beta}(t)=\sum_{n=0}^{\infty}\frac{t^{(2n)\beta}}{\Gamma(2n\beta+1)}=E_{2\beta}(t^{2\beta})$
$\sinh(t)$ | $\sinh_{\beta}(t)=\sum_{n=0}^{\infty}\frac{t^{(2n+1)\beta}}{\Gamma((2n+1)\beta+1)}=t^{\beta}E_{2\beta,\beta+1}(t^{2\beta})$
$e^{t}\cos(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k}{(-1)^{k}}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$e^{-t}\cos(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k}{(-1)^{n-k}}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$e^{t}\sin(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k+1}{(-1)^{k}}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$e^{-t}\sin(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k+1}{(-1)^{n-k-1}}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$g(t)\cos(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k}{(-1)^{k}}{g^{(n-2k)}(0)}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$g(t)\sin(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k+1}{(-1)^{k}}{g^{(n-2k-1)}(0)}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$g(t)\cosh(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k}{g^{(n-2k)}(0)}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
$g(t)\sinh(t)$ | $\sum_{n=0}^{\infty}[\sum_{k=0}^{n}\binom{n}{2k+1}{g^{(n-2k-1)}(0)}]\frac{t^{\beta n}}{\Gamma(n\beta+1)}$
Table 1.  A Wright beta transform for different functions. This table will be
useful in Section 4 where we will solve nonhomogeneous fractional ODEs

##  3\. Higher order fractional ordinary differential equations

In this Section, we consider a new relationship between solutions of higher
fractional differential equations and the Wright-type transformation, thus
improving the following Lemma presented previously in [21].

######  Lemma 1.

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)$
(3.1)

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$
(3.2)

for $n\geq 1$.

We point out the following lemma, which characterizes higher fractional
derivatives of different orders in terms of Wright-type transformations; see
(3.3), which improves the presentation given in Lemma 1, see also [21]. In
particular, this Lemma does not require the restricted condition (3.1).

######  Lemma 2.

Under the assumption of $f_{\beta}(0)=f(0),$
$D_{c}^{\beta}f_{\beta}(0)=f^{\prime}(0),$
$D_{c}^{\beta}D_{c}^{\beta}f_{\beta}(0)=f^{{}^{\prime\prime}}(0)...$ and
$\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n$ times}}f_{\beta}(0)=f^{(n)}(0),$ the
following holds:

$\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n$
times}}f_{\beta}=\int_{0}^{\infty}f^{(n)}(x)g_{\beta}(x,t)dx.$ (3.3)

Proof

We will prove that both parts of the equation have the same Laplace transform.
The case $n=1$ has been proved previously in [21]. Indeed
$D_{c}^{\beta}f_{\beta}(t)=\int_{0}^{\infty}f^{\prime}(x)g_{\beta}(x;t)dx,$
where $0<{\beta}<1.$ Taking the Laplace transform on both sides of (3.3), on
the left-hand side we have

$\displaystyle\mathcal{L}\\{\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n$ times}}f_{\beta}\\}$ $\displaystyle=$
$\displaystyle s^{\beta}\overline{\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n-1$
times}}f_{\beta}(t)}-s^{\beta-1}\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n-1$ times}}f_{\beta}(0)$ $\displaystyle=$
$\displaystyle
s^{\beta}(s^{n\beta-1}\overline{f}(s^{\beta})-s^{(n-1)\beta-1}f_{\beta}(0)-s^{(n-2)\beta-1}D_{c}^{\beta}f_{\beta}(0)-...)$
$\displaystyle-s^{\beta-1}\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n-1$ times}}f_{\beta}(0)$ $\displaystyle=$
$\displaystyle
s^{(n+1)\beta-1}\overline{f}(s^{\beta})-s^{(n)\beta-1}f_{\beta}(0)-s^{n\beta-\beta-1}D_{c}^{\beta}f(0)...-s^{\beta-1}\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n-1$ times}}f_{\beta}(0).$

On the other hand, taking the Laplace transform on the right-hand side and
using our favorite equation (2.4), we have

$\displaystyle\mathcal{L}\left\\{{\int_{0}^{\infty}f^{(n)}(x)g_{\beta}(x;t)dx}\right\\}$
$\displaystyle=$ $\displaystyle
s^{\beta-1}\mathcal{L}\\{{f}^{(n)}{(\cdot)}\\}(s^{\beta})$ $\displaystyle=$
$\displaystyle
s^{\beta-1}\left(s^{n\beta}\overline{f}(s^{\beta})-s^{(n-1)\beta}f(0)-s^{(n-2)\beta}f^{\prime}(0)...-f^{(n-1)}(0)\right)$
$\displaystyle=$ $\displaystyle
s^{(n+1)\beta-1}\overline{f}(s^{\beta})-s^{(n)\beta-1}f(0)-s^{n\beta-\beta-1}f^{\prime}(0)...-s^{\beta-1}f^{(n-1)}(0).$

Therefore under the assumptions of $f_{\beta}(0)=f(0),$
$D_{c}^{\beta}f_{\beta}(0)=f^{\prime}(0),$
$D_{c}^{\beta}D_{c}^{\beta}f_{\beta}(0)=f^{{}^{\prime\prime}}(0)...$ and
$\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n$ times}}f_{\beta}(0)=f^{(n)}(0)$, we
obtain the equality.

Thanks to Lemma 2, we obtain the following Theorem which is a significant
improvement to Lemma 1. In the following Theorem we don’t require condition
(3.1) which seems to be quite restricted and was presented in [21].

######  Theorem 1.

For non-homogeneous fractional differential equations of the form:

$\sum_{k=1}^{n-1}a_{k}\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{k\text{
times}}D_{C}^{\beta}y(t)+a_{n+1}y(t)=F_{\beta}(t)$

such that $f_{\beta}(0)=z(0),$ $D_{c}^{\beta}f_{\beta}(0)=z^{\prime}(0),$
$D_{c}^{\beta}D_{c}^{\beta}f_{\beta}(0)=z^{{}^{\prime\prime}}(0)...$ and
$\underbrace{D_{c}^{\beta}\cdot\cdot\cdot
D_{c}^{\beta}D_{c}^{\beta}}_{\text{$n$ times}}f_{\beta}(0)=z^{(n)}(0)$ for
$0<\beta\leq 1$, the solution is given by

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

where the process $\mathcal{T}_{\beta}(t)$ can be seen as the inverse of a
$\beta$ stable subordinator whose density is $g_{\beta}(x,t)$,
$F_{\beta}(t)=\int_{0}^{\infty}F(s)g_{\beta}(s;t)ds$, and $z(t)$ is the
solution of the linear ordinary differential equation with the same initial
conditions.

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

$z^{n-1}(0)=z_{n},...z(0)=z_{0}.$

The following corollary solves second-order homogeneous fractional ODEs.

######  Corollary 1.

Consider the equation

$a_{2}D_{c}^{\beta}D_{c}^{\beta}y_{\beta}(t)+a_{1}D_{c}^{\beta}y_{\beta}(t)+a_{0}y_{\beta}(t)=0,$

if $r_{1}$ and $r_{2}$ are the roots of the characteristic polynomial
associated to the analogous second order differential equation, then the
solution has the form:
Case 1: Distinct and Real Roots

$\sum_{n=0}^{\infty}{[c_{2}r_{2}^{n}+c_{1}r_{1}^{n}]\frac{t^{n\beta}}{\Gamma{(n\beta+1)}}}=c_{2}E_{\beta}{(r_{2}t^{\beta})}+c_{1}E_{\beta}{(r_{1}t^{\beta}).}$

Case 2: Repeating Roots

$\sum_{n=0}^{\infty}{[c_{2}r^{n}+c_{1}nr^{n-1}]\frac{t^{n\beta}}{\Gamma{(n\beta+1)}}}=c_{2}E_{\beta}(rt^{\beta})+c_{1}\frac{\partial}{\partial
r}E_{\beta}(rt^{\beta}).$

Case 3: Complex Roots Let $r_{2}=\lambda+\mu i,r_{1}=\lambda-\mu i$

$\sum_{n=0}^{\infty}{[c_{2}r_{2}^{n}+c_{1}r_{1}^{n}}]\frac{t^{n\beta}}{\Gamma{(n\beta+1)}}=c_{2}E_{\beta}{(r_{2}t^{\beta})}+c_{1}E_{\beta}{(r_{1}t^{\beta}).}$

###  3.1. Application to Fractional Embedded Beam

Euler studied as an eigenvalue problem how a thin elastic column buckles under
a compressive axial. The differential equation governing the deflection
$\phi(t)$ of a thin elastic column subject to a constant compressive axial
force (or load) $P$ when $P$ is applied to its top is given by

$\frac{d^{2}}{dt^{2}}\left(EI\frac{d^{2}\phi(t)}{dt^{2}}\right)+P\frac{d^{2}\phi(t)}{dt^{2}}=0$
(3.4)

where $E$ is Young’s modulus of elasticity and $I$ is the moment of inertia of
a cross-section about a vertical line through a centroid, see [28]. If we
assume that the column is uniform, $EI$ is a constant. In the next example, we
solve the equivalent fractional high-order differential equation.

Consider a beam of length $L$. As a first application, we solve the fractional
embedded beam equation

$EI\,D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}\phi(t)=w_{0},$

subject to
$\phi(0)=0,\phi(L)=0,D_{c}^{\beta}\phi(0)=0,D_{c}^{\beta}\phi(L)=0.$ Using
Theorem 1 we obtain

$\displaystyle\phi_{\beta}(t)$ $\displaystyle=$
$\displaystyle\frac{w_{0}}{24EI}\int_{0}^{\infty}(x^{4}-2x^{3}L+L^{2}x^{2})g_{\beta}(x;t)dx$
$\displaystyle=$
$\displaystyle\frac{w_{0}}{24EI}\left(4!\frac{t^{4\beta}}{\Gamma(4\beta+1)}-12\frac{t^{3\beta}L}{\Gamma(3\beta+1)}+L^{2}\frac{t^{2\beta}L}{\Gamma(2\beta+1)}\right).$

Figure 1. Monte Carlo simulation of the solution of
$EI\,D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}\phi(t)=w_{0}$ with
$L=1$ and $EI=w_{0}$.

Figure 1 shows a Monte Carlo simulation of the solution, for more details see
Section 5. As a second application we solve, the following equivalent
fractional equation corresponding to (3.4)

$EID_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}\phi(t)+PD_{c}^{\beta}D_{c}^{\beta}\phi(t)=0,$

subject to $\phi(0)=0,\phi(L)=0,D_{c}^{\beta}D_{c}^{\beta}\phi(0)=0$ and
$D_{c}^{\beta}D_{c}^{\beta}\phi(L)=0.$ By Theorem 1, the solution is given by

$\displaystyle\phi_{\beta}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}(c_{1}\cos(n\pi x/L)+c_{2}\sin(\alpha
x)+c_{3}x+c_{4})g_{\beta}(x;t)dx$ $\displaystyle=$ $\displaystyle
c_{1}cos_{\beta}((n\pi/L)^{\frac{1}{\beta}}t)+c_{2}sin_{\beta}((n\pi/L)^{\frac{1}{\beta}}t)+\frac{t^{\beta}L}{\Gamma(\beta+1)}+c_{4}$
$\displaystyle=$ $\displaystyle
c_{1}E_{2\beta}(-((n\pi/L)^{\frac{1}{\beta}}t)^{2\beta})+c_{2}(n\pi/L)t^{\beta}E_{2\beta,\beta+1}(-(n\pi/L)^{2}t^{2\beta})+\frac{t^{\beta}L}{\Gamma(\beta+1)}+c_{4}.$

##  4\. Solving Fractional Electric Circuits

In this Section, we apply Theorem 1 to solve extensions of fractional electric
circuits with an external source, including an external source coming from the
Mittagg-Leffler function and other special functions; similar work was
presented for the homogeneous case, [1]. The basic equations of electric
circuits involving resistors, capacitors, and inductors are [28]:

$\displaystyle I^{\prime\prime}(t)+\frac{1}{LC}I(t)$ $\displaystyle=$
$\displaystyle 0,$ (4.1) $\displaystyle V^{\prime}(t)+\frac{1}{RC}V(t)$
$\displaystyle=$ $\displaystyle 0,$ (4.2) $\displaystyle
I^{\prime}(t)+\frac{R}{L}I(t)$ $\displaystyle=$ $\displaystyle\frac{V}{L},$
(4.3)

where the $I(t)$ is the current in the circuit at time $t,$ $L$ is the
inductance, $C$ is the capacitance, $R$ is the resistance and $V$ is the
voltage drop across the circuit. Equation (4.1) represents the LC (inductor
-capacitor) circuit, equation (4.2) represents RC (resistor- capacitor)
circuit.

######  Example 1 (RC circuit).

By Theorem 1, the solution of the fractional circuit
$D_{c}^{\beta}V(t)+V(t)/RC=0$ with $V(0)=V_{0}$ is given by

$V_{\beta}(t)=\int_{0}^{\infty}V_{0}e^{-x/RC}g_{\beta}(x;t)dx=V_{0}E_{\beta}\left(-t^{\beta}/RC\right).$

This solution agrees with the literature [7]. Furthermore by Theorem 1,
$V_{\beta}(t)=\mathbb{E}[y_{0}e^{\lambda\mathcal{T}_{\beta}(t)}]$.

######  Example 2 (LC Circuits).

The solution of the fractional LC circuit
$D_{c}^{\beta}D_{c}^{\beta}V(t)+V(t)/LC=0$ with $V(0)=0$ and
$V^{\prime}(0)=1$, where $\omega=1/LC$ is a real-valued constant is given by

$\displaystyle V_{s}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}g_{\beta}(x,t)\sin(\omega
x)\,dx=\sum_{n=0}^{\infty}\frac{(-1)^{n}(\omega)^{2n+1}}{(2n+1)!}\int_{0}^{\infty}g_{\beta}(x,t)\frac{x^{2n+1}}{(2n+1)!}\,dx$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{(-1)^{n}(\omega^{\frac{1}{\beta}}t)^{\beta(2n+1)}}{\Gamma((2n+1)\beta+1)}=\omega
t^{\beta}E_{2\beta,\beta+1}(-\omega^{2}t^{2\beta})=sin_{\beta}(\omega^{\frac{1}{\beta}}t).$

We have applied Theorem 1 to $z(x)=\sin(\omega x)$ which solves
$D^{2}z(x)=-\omega^{2}z(x)$, with $z(0)=0$. That is,
$y(t)=\mathbb{E}[\sin(\omega\mathcal{T}_{\beta}(t))]$. This agrees with A.A.
Stanislasky’s finding, see equation [25]. Furthermore, the solution of the
frational LC circuit $D_{C}^{\beta}D_{C}^{\beta}V(t)+V(t)/LC=0$ with $y(0)=1$
and $y^{\prime}(0)=0$ where $\omega=1/LC$ is a real-valued constant is given
by

$\displaystyle V_{c}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}g_{\beta}(x,t)\cos(\omega
x)\,dx=\sum_{k=0}^{\infty}\frac{(-1)^{k}(\omega)^{2k}}{(2k)!}\int_{0}^{\infty}g_{\beta}(x,t)\frac{x^{2k}}{(2k)!}\,dx$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\frac{(-1)^{n}\omega^{2k}}{(2k)!}\frac{t^{\beta(2k)}(2k)!}{\Gamma(2k\beta+1)}=\sum_{k=0}^{\infty}\frac{(-1)^{k}(\omega^{\frac{1}{\beta}}t)^{\beta
2k}}{\Gamma(2k\beta+1)}$ $\displaystyle=$ $\displaystyle
E_{2\beta}(-(\omega^{\frac{1}{\beta}}t)^{2\beta})=cos_{\beta}(\omega^{\frac{1}{\beta}}t).$

This also agrees with A.A. Stanislasky’s finding
$\mathbf{E}_{2\beta,1}(\omega^{2}t^{2\beta})$ [25].

Next, we consider examples of fractional electric circuits with an external
source.

######  Example 3.

Solve
$D_{c}^{\beta}y_{\beta}(t)+ay_{\beta}(t)=F_{\beta}(t)=E_{\beta}(t^{\beta}),$
$0<\beta<1$ for the case that $F(t)=e^{t}$. Combining the solution for the
regular ODE $y^{\prime}+ay=F(t)$ and by Theorem 1, we have

$y_{\beta}(t)=\int_{0}^{\infty}\left(\frac{1}{a+1}(e^{s}-e^{-as})+e^{-as}\right)g_{\beta}(s;t)ds,$

and by propostion 1, we obtain

$y_{\beta}(t)=\frac{1}{a+1}E_{\beta}(t^{\beta})+\frac{1}{a+1}E_{\beta}(-at^{\beta})+y_{0}E_{\beta}(-at^{\beta}).$

Figure 2. A Monte Carlo simulation of the solution of
$D_{c}^{\beta}y_{\beta}(t)+ay_{\beta}(t)=E_{\beta}(t^{\beta}),$ with $a=1$.

Figure 2 shows a Monte Carlo simulation of the solution; for more details see
Section 5.

######  Example 4.

Let’s consider the following non-homogeneous fractional equation

$D_{c}^{\beta}D_{c}^{\beta}y_{\beta}(t)+\omega^{2}y_{\beta}(t)=\left(t^{2}\right)_{\beta}=\frac{2t^{2\beta}}{\Gamma(2\beta+1)},0<\beta<1.$

Using Theorem 1, we apply the Wright-type transform to the solution of the
standard ODE $y^{{}^{\prime\prime}}(t)+\omega^{2}y(t)=t^{2},$ obtaining

$\displaystyle y_{\beta}(t)$
$\displaystyle=\int_{0}^{\infty}\left[c_{1}\cos(\omega s)+c_{2}\sin(\omega
s)+\left(\frac{s^{2}}{\omega^{2}}-\frac{2}{\omega^{4}}\right)\right]g_{\beta}(s;t)ds$
$\displaystyle=c_{1}E_{2\beta}(-(\omega t)^{2\beta})+c_{2}(\omega
t)^{\beta}E_{2\beta,\beta+1}(-(\omega
t)^{2\beta})+\frac{1}{\omega^{2}}\left[\frac{2t^{2\beta}}{\Gamma(2\beta+1)}\right]-\frac{2}{\omega^{4}}.$

Figure 3. A Monte Carlo simulation of the solution of
$D_{c}^{\beta}D_{c}^{\beta}y_{\beta}(t)+\omega^{2}y_{\beta}(t)=\frac{2t^{2\beta}}{\Gamma(2\beta+1)},0<\beta<1$,
and $c_{1}=c_{2}=1$.

Figure 3 shows a Monte Carlo simulation of the solution; for more details see
Section 5.

######  Example 5.

$D_{c}^{\beta}D_{c}^{\beta}$$y(t)+\omega^{2}y(t)=t^{\beta}E_{2\beta,\beta+1}(-t^{2\beta})=sin_{\beta}(t),$
$0<\beta<1$

Using Theorem 1, we apply the Wright-type transform to the solution of the
standard ODE $y^{{}^{\prime\prime}}(t)+\omega^{2}y(t)=\sin{(t)},$ obtaining

$\displaystyle y_{\beta}(t)$ $\displaystyle=\int_{0}^{\infty}C_{1}\cos{(\omega
x)}+C_{2}\sin{(\omega x)}+\frac{\sin{(x)}}{\omega^{2}-1}g_{\beta}(x;t)dx$
$\displaystyle=\int_{0}^{\infty}C_{1}\cos{(\omega
x)}g_{\beta}(x;t)dx+\int_{0}^{\infty}C_{2}\sin{(\omega
x)}g_{\beta}(x;t)dx+\int_{0}^{\infty}\frac{\sin{(x)}}{\omega^{2}-1}g_{\beta}(x;t)dx$
$\displaystyle=C_{1}\cos_{\beta}(\omega t)+C_{2}\sin_{\beta}(\omega
t)+\frac{\sin_{\beta}(t)}{\omega^{2}-1}$
$\displaystyle=C_{1}E_{2\beta}(-\omega^{2}t^{2\beta})+C_{2}\omega
t^{\beta}E_{2\beta,\beta+1}(-\omega^{2}t^{2\beta})+\frac{t^{\beta}E_{2\beta,\beta+1}(-t^{2\beta})}{\omega^{2}-1}.$

Figure 4. A Monte Carlo simulation of the solution of
$D_{c}^{\beta}D_{c}^{\beta}$$y(t)+\omega^{2}y(t)=t^{\beta}E_{2\beta,\beta+1}(-t^{2\beta})=sin_{\beta}(t),$
$0<\beta<1$ with $\omega=2$, and $c_{1}=c_{2}=1$.

Figure 4 shows a Monte Carlo simulation of the solution; for more details see
Section 5

##  5\. Monte Carlo Simulations of Fractional Differential Equations

In this section, we will show how to simulate solutions of fractional
differential equations using Monte Carlo methods. The following algorithm was
introduced for first derivatives in [21]. In this work, we extend these ideas
to higher derivatives. The idea is to generate a new but random timeline
$L(t)$ using $g_{\beta}(\cdot,t)$ for each $t$ and then use any ODE numerical
solver, like Runge-Kutta, to simulate it over that new timeline and finally
find the average of the solution, see [21].

Monte-Carlo Algorithm

  1. (1)

Step 1: For each $t_{i}$, where $i=1,2,\ldots,n$, generate random numbers
$s_{1},\ldots,s_{m}$ from $g_{\beta}(\cdot,t_{i})$.

  2. (2)

Step 2: Discretize the timeline $L(t)$ from $\min(s_{1},\ldots,s_{m})$ to
$\max(s_{1},\ldots,s_{m})$.

  3. (3)

Step 3: Solve the ODE using Runge-Kutta from $\min(s_{1},\ldots,s_{m})$ to
$\max(s_{1},\ldots,s_{m})$ and interpolate the solution at the rest of
$s_{1},\ldots,s_{m}$.

  4. (4)

Step 4: Find the average of the solution over those $s_{1},\ldots,s_{m}$.

In the following three examples, we simulate some fractional differential
equations using the Monte Carlo algorithm.

######  Example 6.

Let $0<\beta<1$, and consider

$D_{c}^{\beta}y(t)+y(t)=\sin_{\beta}{(t)}$

such that $y(0)=0$. It is simulated by solving $z^{\prime}(t)+z(t)=\sin{(t)}$
such that $z(0)=0$ using Runge-Kutta of 4 and 5 over simulated timelines
$L(t)$ for each $t$. See Figure 5.

Figure 5. A Monte Carlo simulation of the solution of
$D_{C}^{\beta}y(t)+y(t)=\sin_{\beta}{(t)}$ with $y(0)=0$ evaluated at
$\beta=0.5$ (left) and $\beta=0.9$ (right).

######  Example 7.

Let $0<\beta<1$, and consider

$D_{c}^{\beta}D_{c}^{\beta}y(t)+\omega^{2}y(t)=\sin_{\beta}{(t)}$

such that $y(0)=0$ and $y^{\prime}(0)=0$. It is simulated by solving
$z^{\prime\prime}(t)+w^{2}z(t)=\sin{(t)}$ such that $z(0)=0$ and
$z^{\prime}(0)=0$ using Runge-Kutta of 4 and 5 over simulated timelines $L(t)$
for each $t$. See Figure 6.

Figure 6. A Monte Carlo simulation of the solution of
$D_{c}^{\beta}D_{c}^{\beta}y(t)+\omega^{2}y(t)=\sin_{\beta}{(t)}$ with
$y(0)=0$ and $y^{\prime}(0)=0$.

######  Example 8.

Let $0<\beta<1$, and consider

$D_{c}^{\beta}D_{c}^{\beta}y(t)+\omega^{2}y(t)=\frac{2t^{2\beta}}{\Gamma(2\beta+1)}$

such that $y(0)=0$ and $y^{\prime}(0)=1$. It is simulated by solving
$z^{\prime\prime}(t)+w^{2}z(t)=t^{2}$ such that $z(0)=0$ and $z^{\prime}(0)=0$
using Runge-Kutta of 4 and 5 over simulated timelines $L(t)$ for each $t$.
Figure 7 shows a Monte Carlo simulation of the solution; for more details see
Section 5.

Figure 7. A Monte Carlo simulation of the solution of
$D_{c}^{\beta}D_{c}^{\beta}y(t)+\omega^{2}y(t)=\frac{2t^{2\beta}}{\Gamma(2\beta+1)}$
with $y(0)=0$ and $y^{\prime}(0)=1$ at $\beta=0.2$ (top left), $\beta=0.5$
(top right), $\beta=0.7$ (bottom left), and $\beta=0.9$ (bottom right).

##  6\. d’Alembert fractional formula for a fractional wave equation

In this Section we present the last application of Theorem 1 and Lemma 2, the
d’Alembert’s formula for a fractional wave equation. We also use the
combination of Monte Carlo simulations with $g_{\beta}(\cdot,t)$ to show how
the solutions change on time for different values of fractional derivatives,
see Figure 9.

######  Corollary 2.

(d’Alembert’s formula) The fractional wave equation

$\displaystyle D_{c}^{\beta}D_{c}^{\beta}u_{\beta}(x,t)$ $\displaystyle=$
$\displaystyle c^{2}\frac{\partial^{2}u_{\beta}(x,t)}{\partial x^{2}},\text{
}-\infty<x<\infty,\text{ }t>0$ (6.1) $\displaystyle u_{\beta}(x,0)$
$\displaystyle=$ $\displaystyle f_{\beta}(x)$ (6.2)
$\displaystyle\frac{\partial u_{\beta}(x,0)}{\partial t}$ $\displaystyle=$
$\displaystyle 0$ (6.3)

has the following solution

$u_{\beta}(x,t)=\frac{1}{2}[f_{\beta}(x+ct)+f_{\beta}(x-ct)].$

Proof: Taking the Fourier transform on both sides of (6.1)-(6.3) and using the
representation of Caputo derivative presented in Lemma 2 , we obtain

$\displaystyle D_{c}^{\beta}\text{\ }D_{c}^{\beta}\mathcal{U}(k,t)$
$\displaystyle=$ $\displaystyle-c^{2}k^{2}\mathcal{U}(k,t)$ (6.4)
$\displaystyle\mathcal{U}(k,0)$ $\displaystyle=$ $\displaystyle F(k)$
$\displaystyle\text{\ }\frac{d\mathcal{U}}{dt}(k,0)$ $\displaystyle=$
$\displaystyle 0.$

(where $\mathcal{U}(k,t)$ is the Fourier transform of $u(x,t)$) in order to
solve equation 6.4 using Theorem 1. We first solve the regular ODE

$\displaystyle\text{\ }\frac{d^{2}\mathcal{U}}{dt^{2}}(k,t)$ $\displaystyle=$
$\displaystyle-c^{2}k^{2}\mathcal{U}(k,t)$ $\displaystyle\mathcal{U}(k,0)$
$\displaystyle=$ $\displaystyle F(k)$ $\displaystyle\text{\
}\frac{d\mathcal{U}}{dt}(k,0)$ $\displaystyle=$ $\displaystyle 0$

The solution is given by

$\mathcal{U}(k,t)=\frac{F(k)}{2}\left(e^{ickt}+e^{-ickt}\right).$

Taking the beta transformation, we obtain the solution given by

$\displaystyle\mathcal{U}_{\beta}(k,t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}g_{\beta}(z;t)\,\mathcal{U}(k,z)dz$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}g_{\beta}(z;t)\frac{F(k)}{2}\left(e^{ickz}+e^{-ickz}\right)dz.$

Now, taking the Fourier inverse

$\frac{1}{\sqrt{2\pi}}\int\limits_{\mathbb{R}}e^{i\kappa
x}\mathcal{U}_{\beta}(k,t)dk=\frac{1}{2}\frac{1}{\sqrt{2\pi}}\int\limits_{\mathbb{R}}e^{i\kappa
x}\int_{0}^{\infty}g_{\beta}(z;t)\frac{F(k)}{2}\left(e^{ickz}+e^{-ickz}\right)dzdk$

and interchanging the order of integration

$\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}g_{\beta}(z;t)\int\limits_{\mathbb{R}}e^{i\kappa
x}\mathcal{U}(k,z)dkdz=\frac{1}{2}\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}g_{\beta}(z;t)\int\limits_{\mathbb{R}}e^{i\kappa
x}\frac{F(k)}{2}\left(e^{ickz}+e^{-ickz}\right)dkdz,$

we get

$\displaystyle\int_{0}^{\infty}g_{\beta}(z;t)u(x,z)dz$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}g_{\beta}(z;t)\int\limits_{\mathbb{R}}e^{i(\kappa
x+ckz)}\frac{F(k)}{2}dkdz$
$\displaystyle+\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}g_{\beta}(z;t)\int\limits_{\mathbb{R}}e^{i(\kappa
x-ckz)}\frac{F(k)}{2}dkdz.$

Simplifying, we obtain

$u_{\beta}(x,t)=\frac{1}{2}\int_{0}^{\infty}g_{\beta}(z;t)f(x+cz)dz+\frac{1}{2}\int_{0}^{\infty}g_{\beta}(z;t)f(x+cz)dz$

Thus, we obtain the d’Alembert solution of the wave equation of the form

$u_{\beta}(x,t)=\frac{1}{2}[f_{\beta}(x+ct)+f_{\beta}(x-ct)]=\frac{1}{2}\mathbb{E}\left[f(x+\mathcal{T}_{\beta}(t))+f(x-\mathcal{T}_{\beta}(t))\right]$

where the process $\mathcal{T}_{\beta}(t)$ can be seen as the inverse of a
$\beta$ stable subordinator whose density is $g_{\beta}(x,t)$; for a detailed
discussion of the d’Alembert solution of the fractional wave equation in
general spaces, see [17] and references therein.

Figure 8. A Monte Carlo simulation of the solution for fractional wave
equation
$D_{c}^{\beta}D_{c}^{\beta}u_{\beta}(x,t)=c^{2}\frac{\partial^{2}u_{\beta}(x,t)}{\partial
x^{2}},u_{\beta}(x,0)=f_{\beta}(x),\frac{\partial u_{\beta}(x,0)}{\partial
t}=0$ evaluated at $\beta=0.1$ no diffusion (top), $\beta=0.5$ showing
diffusion (middle), and $\beta=0.9$ d’Alembert’s formula (bottom) with $c=.5$
and $f(x)=\exp(-10x^{2})$.

##  7\. Deep Learning to Simulate Fractional Differential Equations

In this section, we will show how to use Feedforward Neural Networks to
simulate solutions of ordinary and fractional differential equations. The main
difference between the fractional and ordinary is in the number of past
measurements (memory length) required to predict the measurement of the next
step. We found that ODEs require one past step while FDEs require more than
one (longer memory).

The following algorithm is used to simulate a solution for the fractional
differential equations mentioned above. The algorithm uses Runge-Kutta with
hybrid orders 4 and 5 to solve the differential equations and then uses the
Monte Carlo interpolation to conduct the Wright-type transformation. The data
from the transform is used as training data for a feedforward neural network
(FFNN). The neural network then provides predictions on future solutions to
the fractional differential equation. Due to the definition of derivatives, in
fractional calculus being different than traditional derivatives we approach
the construction of the FFNN differently. For normal differential equations we
would be able to use a single moment from the past to predict a moment in the
future. For our fractional differential equation example we needed to take
three moments from the past to predict the future. The neural network was
constructed out of three layers with 10 neurons in each layer.

Algorithm For Feedforward Neural Network

  1. (1)

Step 1: Simulate training, validation, and test Monte-Carlo simulations of the
fractional differential equations.

  2. (2)

Step 2: Train the FFNN on the training data set and use the validation data
set.

  3. (3)

Step 3: Use three data points or more to predict one future point.

######  Example 9.

We show the FFNN solution of the FDE given by

$D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}y(t)+2D_{c}^{\beta}D_{c}^{\beta}y(t)+D_{c}^{\beta}y(t)+5y(t)=E_{\beta}{(-t)}$

with $0<\beta<1$ with random $y(0)$.

Figure 9. The graph shows neural network and Monte Carlo integration with
Runge-Kutta solutions for
$D_{c}^{\beta}D_{c}^{\beta}D_{c}^{\beta}y(t)+2D_{c}^{\beta}D_{c}^{\beta}y(t)+D_{c}^{\beta}y(t)+5y(t)=E_{\beta}{(-t^{\beta})}$
for $\beta=0.5$ (left) and $\beta=0.9$ (right).

The two graphs produced agree well. When $\beta=0.5,$ the mean squared error
between the neural network and the Monte Carlo simulations is $3.9833e-10$.
When $\beta=0.9$ the mean squared error between the neural network and the
Monte Carlo simulations is $2.0942e-9$. See Figure 9.

##  8\. conclusion

A new relationship (Lemma 2) is established between the solutions of higher
fractional differential equations and a Wright-type transformation. Theorem 1
allows us to solve fractional differential equations of higher orders with
certain initial conditions. Furthermore, solutions could be interpreted as
expected values of functions in a random-time process. Applications include
the fractional beam equation, fractional electric circuits with special
functions as external sources such as the Mittag-Leffler function and other
special functions, and we provide a new derivation of d’Alembert’s formula for
fractional wave equations.

We have shown that formulating solutions of fractional differential equations
and partial differential equations as expected values with respect to power-
law distributions could be enabled using Monte Carlo integration methods along
with other numerical methods such as Runge-Kutta to numerically solve
fractional differential equations. The results of combining Monte Carlo
integration methods along with Runge-Kutta methods showed an excellent
approximation to the exact solution combined with Monte Carlo integration.
Using a feedforward neural network along with a number of steps to make a
prediction was also very successful in making similar simulation results like
those done by Monte Carlo integration.

######  Acknowledgement 1.

All the authors are honored and grateful to report that M. Nacianceno, Y.
Sepulveda, T. Stuck and J. Williams were fully funded during their research
experience for undergraduates (REU) by National Science Foundation with Award
# 2150478 as undergraduate researchers. Also, the faculty E. Suazo, T. Oraby,
H. Rodrigo and J. Sifuentes were partially funded by the same grant.

## References

  * [1] A. Alsaedi, J. J. Nieto and V. Venktesh, Fractional electric circuits, Advances Mechanical Engineering, vol. 7 (12) 1-7, 2015.
  * [2] B. Mandelbrot and H. M. . Taylor, On the Distribution of Stock Price Differences, _Operations Research_ , vol. 15, no. 6, pp. 1057–1062, 1967.
  * [3] O. E. Barndorff-Nielsen and A. Shiryaev, _Change of Time and Change of Measure_.  Advanced Series on Statistical Science & Applied Probability, World Scientific, 2010.
  * [4] N.H. Bingham, Limit theorems for occupation times of Markov processes, Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete vol. 17, pages1–22 (1971)_._
  * [5] Y. Fujita, Integrodifferential equation which interpolates the heat equation and the wave equation, Osaka J. Math. 27 (1990), 309-321.
  * [6] R. Gorenflo, Mittag-leffler waiting time, power laws, rarefaction, continuous time random walk, diffusion limit, arXiv preprint arXiv:1004.4413, 2010.
  * [7] Bangti Jin, Fractional Differential Equations: An Approach via Fractional Derivatives, Applied Mathematical Sciences, vol 206. Springer, 2021.
  * [8] M. Hahn, K. Kobayashi and S. Umarov, SDEs Driven by a Time-Changed Levy Process and Their Associated Time-Fractional Order Pseudo-Differential Equations, _J Theor Probab 25_ , 262-279 (2012).
  * [9] 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.
  * [10] 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.
  * [11] F. Mainardi, Fractional relaxation-oscillation and fractional diffusion-wave phenomena F Mainardi Chaos, Solitons Fractals 7 (9), 1461-1477.
  * [12] 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.
  * [13] 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.
  * [14] F. Mainardi and A. Consiglio, The Wright functions of the second kind in mathematical physics, _Mathematics 2020_ , vol. 8 No 6, p. 884, 2020.
  * [15] F. Mainardi and G. Pagnini, Mellin-Barnes integrals for stable distributions and their convolutions, _Fractional Calculus and Applied Calculus_ , vol. 11, no. 4, 2008.
  * [16] 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.
  * [17] Cheng-Gang Li, Miao Li, Sergey Piskarev, Mark M. Meerschaert, The fractional d’Alembert’s formulas, Journal of Functional Analysis, vol. 277, issue 12, 2019.
  * [18] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity, An Introduction to Mathematical Models, Imperial College Press, London, World Scientific, Singapore, May 2010.
  * [19] M. M. Meerschaert and A. Sikorskii, Stochastic Models for Fractional Calculus, De Gruyter, Berlin, 2012.
  * [20] M. M. Meerschaert and P. Straka, Inverse Stable Subordinators, Math. Model. Nat. Phenom. Vol. 8, No. 2, 2013, pp. 1-16, 2013.
  * [21] T. Oraby, E. Suazo, H. Arrubla, Probabilistic solutions of fractional differential and partial differential, Chaos, Solitons & Fractals 166, 112901.
  * [22] 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.
  * [23] A. Piryatinska, Inference for the levy models and their applications in medicine and statistical physics, _Department of Statistics, Case Western Reserve University_ , 2005.
  * [24] 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
  * [25] A. A. Stanislavsky, Fractional oscillator, Physical Review E, 70, 051103 (2004).
  * [26] A. A. Stanislavsky, Probability interpretation of the integral of fractional order, Theoretical and Mathematical Physics, 138 (3): 418-431 (2004).
  * [27] 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.
  * [28] D. G. Zill, Differential Equations with Boundary-Value Problems, Cengage Learning; 9th edition, 2017.