Sage Journals: Discover world-class research

Abstract

This article describes an efficient algorithm based on residual power series to approximate the solution of a class of partial differential equations of time-fractional Fokker–Planck model. The fractional derivative is assumed in the Caputo sense. The proposed algorithm gives the solution in a form of rapidly convergent fractional power series with easily computable coefficients. It does not require linearization, discretization, or small perturbation. To test simplicity, potentiality, and practical usefulness of the proposed algorithm, illustrative examples are provided. The approximate solutions of time-fractional Fokker–Planck equations are obtained by the residual power series method are compared with those obtained by other existing methods. The present results and graphics reveal the ability of residual power series method to deal with a wide range of partial fractional differential equations emerging in the modeling of physical phenomena of science and engineering.

Keywords

Caputo fractional derivatives Fokker–Planck equation residual power series numerical algorithms

Introduction

Recent decades have witnessed great attention toward fractional calculus, which can be considered as a generalization of classical integer-order integration and differentiation. Many definitions have been suggested for fractional derivatives such as Riesz, Riemann-Liouville, Grunwald-Letnikov, Caputo, and conformable fractional definitions.^1–4 Such differential and integral operators of non-integer order include all historical states of the function in a weighted form called the memory effect. Anyhow, a large number of physical systems are modeled using fractional differential equations (FDEs), particularly fractional partial differential equations (FPDEs). The FPDEs have achieved significance and publicity due to their tremendous use in different fields such as electrochemistry, electrical circuits, theoretical biology, and quantum mechanics.^5–8 Furthermore, the most significant feature for using the FPDEs in such and other applications is the non-local property, while the differential operator of integer order is local. In this light, since the next state of a fractional system depends not only on the current state but also on its entire historical states. This leads to deep consistency of the mathematical model components in dynamic systems and physical processes. Nevertheless, solving those FDEs is a challenge, especially for numerical calculations. Thus, an effective, reliable, and appropriate numerical methods are needed in handling the partial differential equations (PDEs) with a fractional order of physical interest.

The Fokker–Planck equation is one of the classical widely used equations of statistical physics, which was first presented by Fokker and Planck for describing the Brownian motion for particles and the change of probability of a random function in space and time.⁹ In addition, chemical Fokker–Planck equation can be derived from an uncontrolled, second-order truncation of the Kramers–Moyal expansion of the chemical master equation. This equation turns out to be more accurate than the linear-noise approximation of the chemical master equation. Anyhow, Fokker–Planck equation arises in the modeling of many natural science phenomena, including quantum optics, electron relaxation, polymer dynamics, solid-state system, probability flux, and other theoretical and practical models.¹⁰

The concern of this analysis is to consider the numerical approximate solutions of the Fokker–Planck PDE with time-fractional derivative of the following form

\begin{matrix} D_{t}^{β} v (x, t) = [- \frac{\partial}{\partial x} f_{1} (x, t) + \frac{\partial^{2}}{\partial x^{2}} f_{2} (t, x)] v (x, t), \\ x \in R, t > 0, 0 < β \leq 1 \end{matrix}

(1)

along with the initial condition

v (x, 0) = v_{0} (x)

(2)

where $D_{t}^{β}$ is the Caputo time-fractional derivatives, $v_{0} (x)$ is the given analytic function on $R$ , and $v (x, t)$ is suggested to be a causal function of space and time, which means that $v (x, t)$ is vanishing for $x < 0$ and $t < 0$ . The smooth functions $f_{1} (x, t)$ and $f_{2} (x, t)$ are the diffusion and drift coefficients, respectively. Under such assumptions, we assume that the initial value problem (IVP) (1) and (2) has a unique analytical solution on the domain of interest.

The fractional Fokker–Planck equation (F-FPE) has been successfully used in biological molecules, chemical physics, energy consumption, and engineering. Indeed, fractional diﬀusion, a specific type of F-FPE, has been also applied to several situations such as frequency-dependent damping behavior of materials, viscoelasticity, and diﬀusion processes.⁹ Unfortunately, it is not easy to obtain the exact solution for FDEs in general. So, many numerical and analytical techniques are employed to approximate these solutions. The multistep reduced differential transform method,¹⁰ the predictor–corrector approach,¹¹ the Laplace transform method,¹² the variational iteration method (VIM),¹³ and Adomain decomposition method (ADM)¹³ are some of the advanced numerical and approximate methods that have been applied for F-FPEs.

In this article, the residual power series (RPS) method is implemented for solving IVPs (1) and (2). The RPS technique was developed to approximate solutions for certain class of fuzzy differential equations.¹⁴ Later, it was applied in solving different types of differential equations due to its simplicity, accuracy, and efficiency.^15–19 The residual power series method (RPSM) has many advantages; first, it is easy to construct a power series solution for handling both linear and nonlinear equations without the terms of linearization, discretization, or perturbation. Second, the present method provides the solutions in Taylor expansions; therefore, the exact solutions will be available when the solutions are polynomials.^20–24 This technique is a direct way to ensure the rate of convergence for series solution, as it depends on minimizing the residual error related. Third, the solutions along with their derivatives can be applied for each arbitrary point in the given interval. Fourth, the RPSM does not require modifications while converting from lower to higher order. Consequently, it has to be easily applied to the proposed system by selecting an appropriate value for the initial guesses approximations. Fifth, the RPS technique needs minor computational requirements with less time and more accuracy. In addition, the presented method is not affected by round-off errors, since it gives the solution in a form of fractional power series (FPS), and substituting values for the solution variables happens as a final step. Finally, it is of global nature in terms of the solutions obtained as well as its ability to solve various types of mathematical, physical, and engineering problems.^25–29

The rest of this work is organized as follows. In section “Preliminaries and notations,” some essential definitions and primary results relating to fractional calculus are given. In sections “Description of the RPS algorithm,” the RPS algorithm is presented for solving time-fractional Fokker–Planck equations (TF-FPEs). Numerical and analytical results for some illustrative examples using the RPS algorithm are introduced in section “Numerical experiments.” Meanwhile, numeric comparison between the proposed method and those available in the literature is discussed. Concluding remarks are given in the last section.

Preliminaries and notations

In this section, we revisit some essential definitions and basic properties of popular fractional operators, Riemann–Liouville fractional integral and Caputo fractional derivative. Then, we survey the most important results of the FPS representation. Throughout this analysis, the set of real numbers and the set of natural numbers are denoted by $R$ and $N$ , respectively, while the gamma function is denoted by Γ.

Definition 2.1

The integral operator for Riemann–Liouville of order $β \geq 0$ is given by²

J_{t}^{β} v (x, t) = {\begin{matrix} \frac{1}{Γ (β)} \int_{0}^{t} \frac{v (x, η)}{{(t - η)}^{1 - β}} d η, & x \in I, 0 < η 〈 t, β 〉 0 \\ v (x, t) & β = 0 \end{matrix}

(3)

where I is the domain of interest for x.

Next, we present the Caputo fractional derivative³ of order $β$ , which is an alternative operator to the Riemann–Liouville fractional operator as follows

\begin{matrix} D_{t}^{β} u (x) = \frac{1}{Γ (n - β)} \\ \int_{0}^{t} \frac{u^{(n)} (s)}{{(x - s)}^{β - n + 1}} ds, n - 1 < β < n for n \in N \end{matrix}

(4)

In fact, Caputo fractional derivative allows us to include the classical initial and boundary conditions in the formulation of the model, whereas the derivative of a constant is 0. For such reasons, the Caputo sense is considered in this analysis to handle the Fokker–Planck equation.

Definition 2.2

For $n - 1 < β < n$ , $n \in N$ . The Caputo time-fractional derivative operator of order $β$ is defined by¹

D_{t}^{β} v (x, t) = {\begin{matrix} \frac{1}{Γ (n - β)} \int_{0}^{t} \frac{\partial^{n} v (x, ζ)}{\partial ζ^{n}} {(t - ζ)}^{n - β - 1} d ζ, n - 1 < β < n \\ \frac{\partial^{n} v (x, ζ)}{\partial ζ^{n}}, β = n \in N \end{matrix}

(5)

Similarly, the Caputo space-fractional derivative operator of order $β$ is defined by

\begin{matrix} D_{x}^{β} v (x, t) = \\ {\begin{matrix} \frac{1}{Γ (n - β)} \int_{0}^{t} \frac{\partial^{n} v (z, t)}{\partial z^{n}} {(x - z)}^{n - β - 1} d ζ, n - 1 < β < n \\ \frac{\partial^{n} v (z, t)}{\partial z^{n}}, β = n \in N \end{matrix} \end{matrix}

(6)

Theorem 2.1

If $n - 1 < α \leq n$ and $n \in N$ , then $D_{t}^{β} J_{t}^{β} v (x, t) = v (x, t)$ and $J_{t}^{β} D_{t}^{β} v (x, t) = v (x, t) - \sum_{j = 0}^{n - 1} ((\partial^{j} v (x, 0^{+})) / \partial t^{j}) (t^{j} / j!)$ , where $x \in I, t > 0$ .¹

The following are some properties of the operators $D_{a}^{β} and J_{a}^{β}$ for $β > 0$ , $q > - 1$ , that can be found in Podlubny¹

D_{t}^{β} t^{q} = \frac{Γ (q + 1)}{Γ (q - β + 1)} t^{q - β}

(7)

J_{t}^{β} t^{q} = \frac{Γ (q + 1)}{Γ (q + β + 1)} t^{q + β}

(8)

Definition 2.3

An FPS representation at $t_{0}$ has the following form¹⁹

\begin{matrix} \sum_{m = 0}^{\infty} v_{m} (x) {(t - t_{0})}^{m β} = v_{0} (x) + v_{1} (x) {(t - t_{0})}^{β} \\ + v_{2} (x) {(t - t_{0})}^{2 β} + \dots \end{matrix}

(9)

where $0 \leq n - 1 < β \leq n$ , $x \in I$ and $t \geq t_{0}$ is called multiple fractional power series (MFPS) about $t_{0}$ .

Theorem 2.2

Suppose that v(x, t) has the following MFPS representation at t = t₀.¹⁹

\begin{matrix} v (x, t) = \sum_{m = 0}^{\infty} v_{m} (x) {(t - t_{0})}^{m β}, \\ 0 \leq n - 1 < β \leq n, x \in I, t_{0} \leq t < t_{0} + ρ . \end{matrix}

(10)

If $D_{t}^{m β} v (x, t), m = 0, 1, 2, \dots$ are continuous on $I \times (t_{0}, t_{0} + ρ)$ , then the coefficients $v_{m} (x)$ will be in the form $v_{m} (x) = ((D_{t}^{m β} v (x, t_{0})) / (Γ (m β + 1)))$ , such that $D_{t}^{m β} = D_{t}^{β} \cdot D_{t}^{β} \cdot \cdot \cdot D_{t}^{β}$ , where $ρ$ is a positive real number, called the radius of convergence for the MFPS.

Description of the RPS algorithm

The main goal of this section is to present the methodology of the RPS technique in obtaining the MFPS approximation of the time-fractional Fokker–Planck model based on the formula of generalized Taylor in Caputo sense by providing a fractional recursion formula to obtain the coefficients of the MFPS depending on minimizing the residual function. To do this, let us assume the solution $v (x, t)$ of TF-FPEs (1) and (2) has the following MFPS expansion about $t_{0} = 0$

v (x, t) = \sum_{n = 0}^{\infty} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)}, 0 < β \leq 1, t \geq 0, x \in R

(11)

By starting with the initial guess approximation $v (x, 0) = v_{0} (x)$ , the series solution of equation (11) can be rewritten as

v (x, t) = v_{0} (x) + \sum_{n = 1}^{\infty} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)}

(12)

To obtain the MFPS approximate solution, let $v_{k} (x, t)$ indicate the $k th$ -truncated MFPS of $v (x, t)$ such that

v_{k} (x, t) = v_{0} (x) + \sum_{n = 1}^{k} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)}

(13)

Now, define the $k th$ -residual function as follows

{Res}_{v}^{k} (x, t) = D_{t}^{β} v_{k} (x, t) + [\frac{\partial}{\partial x} f_{1} (x, t) - \frac{\partial^{2}}{\partial x^{2}} f_{2} (t, x)] v (x, t)

(14)

where the residual function can be given in the form

Re s_{v} (x, t) = lim_{k \to \infty} {Res}_{v}^{k} (x, t)

(15)

Evidently, $Re s_{v} (x, t) = 0$ and $lim_{k \to \infty} {Res}_{v}^{k} (x, t) = Re s_{v} (x, t) = 0$ for each $x \in R$ and $0 < t < ρ$ , where $ρ$ is the radius of convergence for the MFPS (11). According to the RPSM,^25–27,30 it can be noted that $D_{t}^{j β} Re s_{v} (x, t) = 0$ . In addition, $D_{t}^{(j - 1) β} Re s_{v} (x, t)_{t = 0} = D_{t}^{(j - 1) β} {Res}_{v}^{j} (x, t)_{t = 0}$ for each $j = 1, 2, \dots, k$ . Therefore, the following relations can hold

D_{t}^{(k - 1) β} {Res}_{v}^{k} {(x, t)}_{t = 0} = 0, k = 1, 2, 3, \dots

(16)

These relations help us to determine the values of the coefficients $v_{n} (x), n = 1, 2, \dots, k$ . And so, the approximate solution for TF-FPEs (1) and (2) has been completely constructed.

Anyhow, the next algorithm clarifies the procedure in obtaining the unknown coefficients of equation (13).

Algorithm 3.1

To determine the required coefficients of $v_{k} (x, t)$ , do the following steps:

Step 1: The initial condition $v (x, 0) = v_{0} (x)$ , which is the zeroth FPS approximate solution of $v (x, t)$ .

Step 2: The $k th$ -truncated MFPS $v_{k} (x, t) = v_{0} (x) + \sum_{n = 1}^{k} v_{n} (x) (t^{n β} / (Γ (n β + 1)))$ into the $k th$ -residual function ${Res}_{v}^{k} (x, t)$ of equation (14).

Step 3: $D_{t}^{(k - 1) β} {Res}_{v}^{k} (x, t)$ for $k = 1, 2, \dots, N$ .

Step 4: The resulting fractional equations $D_{t}^{(k - 1) β} {Res}_{v}^{k} (x, 0) = 0$ for $v_{k} (x)$ .

Step 5: The obtained coefficients, for $k = 1, 2, \dots, N$ , back into equation (13).

Numerical experiments

The purpose of this section is to show the high degree of accuracy, efficiency, and applicability of this algorithm. The approximate analytical solutions of TF-FPEs are constructed in a rapidly convergent FPS form. Numeric comparisons of the results obtained by the proposed method, ADM¹³ and VIM¹³ are provided. The tabular and graphical results reveal that the RPS approach is easy to implement and accurate when applied to the TF-FPEs, as well as it introduces a promising tool for solving many fractional PDEs. The present computations are performed using Mathematica $10$ (Wolfram Mathematica) software package.

Example 4.1

Consider the following TF-FPE

\begin{matrix} D_{t}^{β} v (x, t) = [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] v (x, t), \\ x \in [0, 2], t \geq 0, 0 < β \leq 1 \end{matrix}

(17)

with the initial conditions

v (x, 0) = x

(18)

The exact solution of IVPs (17) and (18) for standard case at $β = 1$ is given by $v (x, t) = x e^{t}$

Using the last description of RPS algorithm, the solution of IVPs (17) and (18) is

v (x, t) = x + \sum_{n = 1}^{\infty} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)}

and $k th$ -truncated MFPS is $v_{k} (x, t) = x + \sum_{n = 1}^{k} v_{n} (x) (t^{n β} / (Γ (n β + 1)))$ . The $k th$ -residual function for equation (17) can be given by

{Res}_{v}^{k} (x, t) = D_{t}^{β} v_{k} (x, t) + \frac{\partial}{\partial x} (x v_{k} (x, t)) - \frac{\partial^{2}}{\partial x^{2}} (\frac{x^{2}}{2} v_{k} (x, t))

(19)

For $k = 1$ , the first residual function is

\begin{matrix} {Res}_{v}^{1} (x, t) = D_{t}^{β} (v_{0} (x) + \sum_{n = 1}^{1} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)}) \\ + \frac{\partial}{\partial x} (x (v_{0} (x) + \sum_{n = 1}^{1} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)})) \\ - \frac{\partial^{2}}{\partial x^{2}} (\frac{x^{2}}{2} (v_{0} (x) + \sum_{n = 1}^{1} v_{n} (x) \frac{t^{n β}}{Γ (n β + 1)})) \\ = D_{t}^{β} (x + v_{1} (x) \frac{t^{β}}{Γ (β + 1)}) \\ + \frac{\partial}{\partial x} (x (x + v_{1} (x) \frac{t^{β}}{Γ (β + 1)})) \\ - \frac{\partial^{2}}{\partial x^{2}} (\frac{x^{2}}{2} (x + v_{1} (x) \frac{t^{β}}{Γ (β + 1)})) \\ = v_{1} (x) - x - x \frac{t^{β}}{Γ (β + 1)} \end{matrix}

(20)

Depending on equation (16), the first unknown coefficient of MFPS expansion is $v_{1} (x) = x$ . Hence, the first FPS approximate solution is

v_{1} (x, t) = x + x \frac{t^{β}}{Γ (β + 1)}

(21)

As the former, to determine the second coefficient $v_{2} (x)$ , consider $k = 2$ in the MFPS (13) and then substitute $v_{2} (x, t)$ into the second residual function ${Res}_{v}^{2} (x, t)$ of equation (19) to get

R e s_{v}^{2} (x, t) = (2 v_{2} (x) - 2 x) \frac{t^{β}}{Γ (β + 1)} - x (2 {v^{'}}_{2} (x) + x {v^{″}}_{2} (x)) \frac{t^{2 β}}{Γ (2 β + 1)}

(22)

By applying $D_{t}^{β}$ on both sides of equation (22), it follows that

D_{t}^{β} R e s_{v}^{2} (x, t) = (2 v_{2} (x) - 2 x) - (2 {v^{'}}_{2} (x) + x {v^{″}}_{2} (x)) \frac{t^{β}}{Γ (β + 1)}

(23)

Using the fact that $D_{t}^{β} {Res}_{v}^{2} (x, 0) = 0$ , the second unknown coefficient in the MFPS is $v_{2} (x) = x$ . Therefore, the second FPS approximate solution is

v_{2} (x, t) = x + x \frac{t^{β}}{Γ (β + 1)} + x \frac{t^{2 β}}{Γ (2 β + 1)}

(24)

Applying similar argument for $k = 3$ , the third unknown coefficient in the MFPS (13) will be $v_{3} (x) = x$ . Moreover, the third FPS approximate solution can be written as

v_{3} (x, t) = x + x \frac{t^{β}}{Γ (β + 1)} + x \frac{t^{2 β}}{Γ (2 β + 1)} + x \frac{t^{3 β}}{Γ (3 β + 1)}

(25)

In the same manner, the process can be repeated till the arbitrary order and then the coefficients of the MFPS solution (13) can be obtained. Consequently, we have $v_{k} (x) = x$ for $k \geq 1$ . Furthermore, if we collect all the last results, the solution $v (x, t)$ can be given as follows

\begin{matrix} v (x, t) = x (1 + \frac{t^{β}}{Γ (β + 1)} + \frac{t^{2 β}}{Γ (2 β + 1)} + \frac{t^{3 β}}{Γ (3 β + 1)} \\ + \dots + \frac{t^{k β}}{Γ (k β + 1)} + \dots) \end{matrix}

(26)

where $E_{γ} (t) = \sum_{j = 0}^{\infty} (t^{j} / (Γ (γ j + 1)))$ is the Mittag–Leffler function.¹ By setting $β = 1$ , equation (26) can be reduced to $v (x, t) = x e^{t}$ , which is the exact solution of the classical form of IVPs (17) and (18).

In view of the obtained previous results and without loss of generality, the geometric behavior of the 10th FPS approximate solution of IVPs (17) and (18) has been studied by drawing the three-dimensional (3D) space figures at different values of $β$ for $x \in [0, 2]$ and $t \in [0, 10]$ . Figure 1 represents a comparison between the 10th FPS approximate solution for $β \in {1.0, 0.95, 0.75, 0.5}$ and the exact solution. Regarding these figures, it is worth mentioning that running the RPS algorithm via Mathematica 10.0 to get such graphs took $0.66, 80.77,$ $88.41,$ and 90.22 times in seconds for $β \in {1.0, 0.95, 0.75, 0.5}$ , respectively. Moreover, to show the efficiency and accuracy of the present algorithm, the numerical results at $x = 1$ with some selected grid points t with step size 0.1 on $[0, 1]$ are given in Table 1. The results in Table 1 show that the FPS approximate solutions are in good agreement with exact solutions.

Figure 1.

Geometric behavior of the 10th FPS solution of Example 4.1 for $x \in [0, 2]$ and $t \in [0, 10]$ : (a) the FPS approximate solution at $β = 1$ , and (b) the FPS approximate solution at $β = 0.95$ , (c) the FPS approximate solution at $β = 0.75$ , and (d) the FPS approximate solution at $β = 0.5 .$

Table 1.

Numerical results of the 10th FPS approximate solution for Example 4.1 at $β = 1$ .

t	Exact solution	Approximation solution	Absolute error	Relative error
0.1	1.1051709180756477	1.1051709180756475	$2.22045 \times 10^{- 16}$	$2.00914 \times 10^{- 16}$
0.2	1.2214027581601699	1.2214027581601692	$6.66134 \times 10^{- 16}$	$5.45384 \times 10^{- 16}$
0.3	1.3498588075760032	1.3498588075759577	$4.55191 \times 10^{- 14}$	$3.37214 \times 10^{- 14}$
0.4	1.4918246976412703	1.4918246976401834	$1.08691 \times 10^{- 12}$	$7.28576 \times 10^{- 12}$
0.5	1.6487212707001282	1.6487212706873657	$1.76250 \times 10^{- 11}$	$7.74082 \times 10^{- 11}$
0.6	1.8221188003905090	1.8221188002948574	$9.56517 \times 10^{- 11}$	$5.24948 \times 10^{- 11}$
0.7	2.0137527074704766	2.0137527069445813	$5.25895 \times 10^{- 10}$	$2.61152 \times 10^{- 10}$
0.8	2.2255409284924680	2.2255409261876826	$2.30479 \times 10^{- 9}$	$1.03561 \times 10^{- 9}$
0.9	2.4596031111569500	2.4596031026621000	$8.49485 \times 10^{- 9}$	$3.45375 \times 10^{- 9}$

FPS: fractional power series.

On the contrary, by applying the ADM,¹³ we have the following iteration

v_{0} (x, t) = x

and

v_{k + 1} (x, t) = J_{t}^{β} [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] v_{k} (x, t)

Using the property of equation (8), it follows that

v_{1} (x, t) = x \frac{t^{β}}{Γ (β + 1)}

v_{2} (x, t) = J_{t}^{β} [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] (x \frac{t^{β}}{Γ (β + 1)}) = x \frac{t^{2 β}}{Γ (2 β + 1)}

v_{3} (x, t) = J_{t}^{β} [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] (x \frac{t^{2 β}}{Γ (2 β + 1)}) = x \frac{t^{3 β}}{Γ (3 β + 1)}

Continuing this process, the nth approximate solution is $v_{n} (x, t) = x (t^{n β} / (Γ (n β + 1)))$

According to the ADM,¹³ the $k th$ ADM solution of IVPs (17) and (18) is given by

\begin{matrix} v_{k} (x, t) = x (1 + \frac{t^{β}}{Γ (β + 1)} + \frac{t^{2 β}}{Γ (2 β + 1)} + \frac{t^{3 β}}{Γ (3 β + 1)} \\ + \dots + \frac{t^{k β}}{Γ (k β + 1)}) = x \sum_{n = 0}^{k} \frac{t^{n β}}{Γ (n β + 1)} \end{matrix}

Hence, the solution is $v (x, t) = lim_{k \to \infty} v_{k} (x, t) = x \sum_{n = 0}^{\infty} (t^{n β} / (Γ (n β + 1))) = x E_{β} (t^{β})$

Obviously, the RPSM produced an identical analytical solution of the ADM solution for this example. Anyhow, to see the effect of the fractional derivative to Fokker–Planck equation, the tabulated and graphical results for the approximate solutions at different values of fractional order $β$ using the RPSM and (ADM)¹³ with $k = 20$ are summarized and listed in Table 2 and Figure 2. It is obvious from the current results that the RPS algorithm is of good agreements with earlier literature works, and the RPS solution of TF-FPE approaches to the solution of the classical case as soon as $β$ approaches to 1.

Table 2.

The FPS solution for Example 4.1 at $x = 1, k = 20$ with different values of $β$ .

t	$β = 1$	$β = 0.9$	$β = 0.8$	$β = 0.7$	$β = 0.6$	$β = 0.5$
0.1	1.051709180	1.140848723	1.189156774	1.255618802	1.349366562	1.486763398
0.2	1.221402758	1.280533157	1.357406208	1.459494954	1.599237102	1.799017244
0.3	1.349858808	1.430496420	1.533460260	1.668128885	1.850146792	2.107699204
0.4	1.491824698	1.593744892	1.722608099	1.889773347	2.114198816	2.430043141
0.5	1.648721271	1.772501193	1.928079927	2.128921777	2.397497008	2.774285980
0.6	1.822118800	1.968845658	2.152581489	2.389055235	2.704509965	3.146213032
0.7	2.013752707	2.184905988	2.398717509	2.673370205	3.039193101	3.550802681
0.8	2.225540928	2.422942168	2.669164047	2.985061731	3.405412950	3.992835823
0.9	2.459603111	2.685397789	2.966760222	3.327465540	3.807152834	4.477184772
1.0	2.718281828	2.974939075	3.294569235	3.704146145	4.248635002	5.008979960

FPS: fractional power series.

Figure 2.

Approximate RPS solution for different values of fractional order $β$ : (a) $x = 0.5, t \in [0, 10]$ and (b) $x \in [0, 2], t = 2$ .

Example 4.2

Consider the following TF-FPE

\begin{matrix} D_{t}^{β} v (x, t) = [- \frac{\partial}{\partial x} \frac{x}{6} + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{12}] v (x, t), \\ x \in [0, 2], t \geq 0, 0 < β \leq 1 \end{matrix}

(27)

with the initial conditions

v (x, 0) = x^{2}

(28)

The exact solution of IVPs (27) and (28) for standard case at $β = 1$ is given by $v (x, t) = x^{2} e^{\frac{t}{2}}$

In view of the RPS technique, by starting with $v_{0} (x) = x^{2}$ as the initial approximation, the $k th$ -residual function of IVPs (27) and (28) can be written as

\begin{matrix} {Res}_{v}^{k} (x, t) = D_{t}^{β} v_{k} (x, t) \\ + \frac{\partial}{\partial x} (\frac{x}{6} v_{k} (x, t)) - \frac{\partial^{2}}{\partial x^{2}} (\frac{x^{2}}{12} v_{k} (x, t)) \end{matrix}

(29)

where $v_{k} (x, t)$ is the $k th$ -truncated MFPS given by equation (13). Thus, for $k = 1$ , the first residual function is given by

\begin{matrix} {Res}_{v}^{1} (x, t) = v_{1} (x) \\ + \frac{1}{12} x (- 6 x - (2 {v'}_{1} (x) + x {v ″}_{1} (x)) \frac{t^{β}}{Γ (β + 1)}) \end{matrix}

(30)

Based on the result of equation (16), it yields that $v_{1} (x) = x^{2} / 2$ . Therefore, the first FPS approximate solution is

v_{1} (x, t) = x^{2} + \frac{x^{2}}{2} \frac{t^{β}}{Γ (β + 1)}

(31)

To determine the second coefficient, let $k = 2$ in the $k th$ -truncated MFPS (13), and substitute $v_{2} (x, t)$ into the ${Res}_{v}^{2} (x, t)$ of equation (30) such that

\begin{matrix} {Res}_{v}^{2} (x, t) = \frac{1}{12} (12 v_{1} (x) - 6 x^{2}) \\ + (12 v_{2} (x) - x (2 {v'}_{1} (x) + x {v ″}_{1} (x))) \frac{t^{β}}{Γ (β + 1)} \\ - x (2 {v'}_{2} (x) + x {v ″}_{2} (x)) \frac{t^{2 β}}{Γ (2 β + 1)} \end{matrix}

(32)

By considering the fact of equation (16) and solving $D_{t}^{β} {Res}_{v}^{2} (x, 0) = 0$ for $v_{2} (x)$ , then it can obtain that $v_{2} (x) = x^{2} / 4$ . Hence, the second FPS approximate solution can be written as

v_{2} (x, t) = x^{2} + \frac{x^{2}}{2} \frac{t^{β}}{Γ (β + 1)} + \frac{x^{2}}{4} \frac{t^{2 β}}{Γ (2 β + 1)}

(33)

For the third unknown coefficient, substitute $v_{3} (x, t)$ into ${Res}_{v}^{3} (x, t)$ of equation (30) to get

\begin{matrix} {Res}_{v}^{3} (x, t) = \frac{1}{12} (12 v_{1} (x) - 6 x^{2}) \\ + (12 v_{2} (x) - x (2 {v'}_{1} (x) + x {v ″}_{1} (x))) \frac{t^{β}}{Γ (β + 1)} \\ + (12 v_{3} (x) - x (2 {v'}_{2} (x) + x {v ″}_{2} (x))) \\ \frac{t^{2 β}}{Γ (2 β + 1)} + x (2 {v'}_{3} (x) + x {v ″}_{3} (x)) \frac{t^{3 β}}{Γ (3 β + 1)} \end{matrix}

(34)

Now, compute $D_{t}^{2 β} {Res}_{v}^{3} (x, t)$ for equation (34) and use $D_{t}^{2 β} {Res}_{v}^{3} (x, 0) = 0$ to get that $v_{3} (x) = x^{2} / 8$ . Therefore, the third FPS approximate solution is given by

\begin{matrix} v_{3} (x, t) = x^{2} + \frac{x^{2}}{2} \frac{t^{β}}{Γ (β + 1)} \\ + \frac{x^{2}}{4} \frac{t^{2 β}}{Γ (2 β + 1)} + \frac{x^{2}}{8} \frac{t^{3 β}}{Γ (3 β + 1)} \end{matrix}

(35)

Using the same process for $k \geq 4$ , the $k th$ unknown coefficient $v_{k} (x)$ can be obtained. Consequently, the solution $v (x, t)$ of IVPs (27) and (28) can be expressed in the form of an infinite series given by

\begin{matrix} v (x, t) = x^{2} \begin{matrix} (1 + \frac{t^{β}}{2 Γ (β + 1)} + \frac{t^{2 β}}{2^{2} Γ (2 β + 1)} + \frac{t^{3 β}}{2^{3} Γ (3 β + 1)} + \dots + \frac{t^{k β}}{2^{k} Γ (k β + 1)} + \dots) \end{matrix} \end{matrix} = x^{2} \sum_{n = 0}^{\infty} \frac{t^{n β}}{2^{n} Γ (n β + 1)}

(36)

For $β = 1$ , the result in equation (36) can be reduced to $v (x, t) = x^{2} e^{\frac{t}{2}}$ , which represents the exact solution for classical form of TF-FPEs.

In view of the previous discussion, the geometric behavior of the 10th FPS approximate solution of IVPs (27) and (28) has been constructed and presented in Figure 3 by drawing the 3D space graphs at different values of $β$ for each $x \in [0, 2]$ and $t \in [0, 10]$ . Furthermore, Figure 3 represents a comparison between the 10th FPS approximate solution for $β = {0.95, 0.75, 0.5}$ and the solution when $β = 1$ . Regarding these figures, it is worth mentioning that running the RPS algorithm via Mathematica 10.0 to get such graphs took $0.75, 79.66, 79.84$ , and 93.72 times in seconds for $β \in {1.0, 0.95, 0.75, 0.5}$ , respectively. It should be noted that the FPS approximate solutions are in good agreement with each other.

Figure 3.

Geometric behavior of 10th FPS of Example 4.2 at $x \in [0, 2]$ and $t \in [0, 10]$ : (a) the FPS approximate solution at $β = 1$ , (b) the FPS approximate solution at $β = 0.95$ , (c) the FPS approximate solution at $β = 0.75$ , and (d) the FPS approximate solution at $β = 0.5$ .

To illustrate the efficiency and accuracy of the fractional residual power series (FRPS) algorithm, some numerical results at fixed value of $x = 1$ and some selected grid points t with step size 0.1 on $[0, 1]$ and $k = 10$ are given in Table 3.

Table 3.

Numerical results of the 10th FPS approximate solution for Example 4.2 at $β = 1$ .

t	Exact solution	Approximation solution	Absolute error	Relative error
0.1	1.0512710963760241	1.0512710963760241	0.00000	0.00000
0.2	1.1051709180756477	1.1051709180756475	$\begin{matrix} 2.22045 \end{matrix} \times 10^{- 16}$	$\begin{matrix} 2.00914 \end{matrix} \times 10^{- 16}$
0.3	1.1618342427282830	1.1618342427282830	0.00000	0.00000
0.4	1.2214027581601699	1.2214027581601692	$6.66134 \times 10^{- 16}$	$5.45384 \times 10^{- 16}$
0.5	1.2840254166877414	1.2840254166877354	$5.99520 \times 10^{- 15}$	$4.66907 \times 10^{- 15}$
0.6	1.3498588075760032	1.3498588075759577	$\begin{matrix} 4.55191 \times 10^{- 14} \end{matrix}$	$\begin{matrix} 3.37214 \times 10^{- 14} \end{matrix}$
0.7	1.4190675485932573	1.4190675485930082	$2.49134 \times 10^{- 13}$	$1.75562 \times 10^{- 13}$
0.8	1.4918246976412703	1.4918246976401834	$1.08691 \times 10^{- 12}$	$7.28576 \times 10^{- 13}$
0.9	1.5683121854901690	1.5683121854861810	$3.98792 \times 10^{- 12}$	$2.54281 \times 10^{- 12}$

FPS: fractional power series.

On the contrary, by applying the ADM,¹³ we have the following iteration

\begin{matrix} v_{0} (x, t) = x^{2} and v_{k + 1} (x, t) \\ = J_{t}^{β} [- \frac{\partial}{\partial x} \frac{x}{6} + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{12}] v_{k} (x, t) \end{matrix}

Using the property of equation (8), we get

\begin{matrix} v_{1} (x, t) = x^{2} \frac{t^{β}}{2 Γ (β + 1)} \\ v_{2} (x, t) = J_{t}^{β} [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] \\ (x^{2} \frac{t^{β}}{2 Γ (β + 1)}) = x^{2} \frac{t^{2 β}}{4 Γ (2 β + 1)} \\ v_{3} (x, t) = J_{t}^{β} [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] \\ (x^{2} \frac{t^{2 β}}{4 Γ (2 β + 1)}) = x^{2} \frac{t^{3 β}}{8 Γ (3 β + 1)} \end{matrix}

and so on; therefore, in this manner, the rest of terms can be obtained. So, the ADM solution of IVPs (27) and (28) is

\begin{matrix} v (x, t) = lim_{k \to \infty} v_{k} (x, t) \end{matrix} = x^{2} * \lim_{k \to \infty} \sum_{n = 0}^{k} \frac{t^{n β}}{2^{n} Γ (n β + 1)}

The previous result is exactly in agreement with the result obtained by the RPSM. To see the effect of the fractional derivative to Fokker–Planck equation, the tabulated and graphical results for the approximate solutions at different values of fractional order $β$ by using the RPSM and ADM¹³ with $k = 10$ are given in Table 4 and Figure 4. It is obvious from these results that the RPS solution of TF-FPE approaches to the solution of the classical case as soon as $β$ approaches to 1.

Table 4.

The FPS approximate solutions for Example 4.2 at $x = 1$ and $k = 10$ .

t	$β = 1$	$β = 0.9$	$β = 0.8$	$β = 0.7$	$β = 0.6$	$β = 0.5$
0.1	1.051271096	1.067873047	1.089647419	1.118280169	1.156149183	1.206731021
0.2	1.105170918	1.130764608	1.162383287	1.201599744	1.250543786	1.312165429
0.3	1.161834243	1.194217917	1.232886661	1.279239821	1.335096142	1.402828108
0.4	1.221402758	1.259271171	1.303403550	1.354996692	1.415525830	1.486763385
0.5	1.284025417	1.326444560	1.374926420	1.430442617	1.494100595	1.567059193
0.6	1.349858808	1.396085151	1.448036274	1.506438647	1.572012045	1.645293775
0.7	1.419067549	1.468464873	1.523133711	1.583537184	1.649990828	1.722406977
0.8	1.491824698	1.543818988	1.600527385	1.662132979	1.728533408	1.799016640
0.9	1.568312185	1.622364270	1.680474811	1.742531667	1.808003639	1.875558788
1.0	1.648721271	1.704308722	1.763203672	1.824985031	1.888684476	1.952358374

FPS: fractional power series

Figure 4.

Approximate RPS solution for different values of fractional order $β$ : (a) $x = 0.5, t \in [0, 10]$ and (b) $x \in [0, 2], t = 2$ .

In addition, the VIM¹³ can be used to solve the IVPs (27) and (28) with the following iterations

\begin{array}{l} v_{0} (x, t) = x^{2} and v_{k + 1} (x, t) = v_{k} (x, t) - \int_{0}^{t} D_{τ}^{β} v_{k} (x, τ) \\ + [- \frac{\partial}{\partial x} x + \frac{\partial^{2}}{\partial x^{2}} \frac{x^{2}}{2}] v_{k} (x, τ) d τ \end{array}

Using this iteration with the help of property of equation (7), it follows that

\begin{matrix} v_{0} (x, t) = x^{2} \\ v_{1} (x, t) = x^{2} (1 + \frac{t}{2}) \\ v_{2} (x, t) = x (1 + t + \frac{t^{2}}{8} - \frac{t^{2 - β}}{2 Γ (3 - β)}) \end{matrix}

Continuing this process, the $k th$ -VIM solution can be given via Mathematica 10.0 for arbitrary k. Anyhow, Table 5 shows a comparison between approximate solutions using the RPSM, ADM,¹³ and VIM¹³ for different values of $β, x,$ and t. As it is evident, the approximate solutions obtained by the RPSM and ADM are identical and are in a good agreement with the VIM approximate solutions for all values of $β$ . It is to be noted that only the fourth term of each method is used in evaluating the results in Table 5.

Table 5.

A comparison of the RPS, ADM,¹³ and VIM¹³ of Example 4.2 for different values of $β$ .

t	x	$β = 0.5$			$β = 0.75$
t	x	RPS	ADM	VIM	RPS	ADM	VIM
0.2	0.25	0.0820103	0.0820103	0.07320989	0.0738078	0.0738078	0.07162351
	0.5	0.328041	0.328041	0.29283956	0.2952311	0.2952311	0.28649408
	0.75	0.738093	0.738093	0.65888901	0.664270	0.664270	0.64461167
	1	1.312165	1.312165	1.17135823	1.1809243	1.1809243	1.14597630
0.4	0.25	0.0929227	0.0929227	0.08280292	0.0830114	0.0830114	0.07997608
	0.5	0.371691	0.371691	0.33121168	0.3320456	0.3320456	0.31990430
	0.75	0.836304	0.836304	0.74522628	0.747103	0.747103	0.71978468
	1	1.486763	1.486763	1.32484671	1.3281823	1.3281823	1.27961720
0.6	0.25	0.1028309	0.1028309	0.09188702	0.0922740	0.0922740	0.88247584
	0.5	0.411323	0.411323	0.36754806	0.369096	0.369096	0.35299034
	0.75	0.925478	0.925478	0.82698314	0.830466	0.830466	0.79422825
	1	1.645294	1.645294	1.47019225	1.476385	1.476385	1.41196134
0.8	0.25	0.1124385	0.1124385	0.10067912	0.10191917	0.10191917	0.09663197
	0.5	0.449754	0.449754	0.40271647	0.407679	0.407679	0.38652789
	0.75	1.011947	1.011947	0.90611205	0.917277	0.917277	0.86968775
	1	1.799017	1.799017	1.61086586	1.630715	1.630715	1.54611160

RPS: residual power series; ADM: Adomain decomposition method; VIM: variational iteration method.

Concluding remarks

Developing analytical and numerical solutions for fractional mathematical models of physical and chemical phenomena are very essential in science. In this work, an analytic-approximate method, so-called RPS, has been employed effectively to solve a class of Fokker–Planck PDEs of fractional order with fitted initial conditions. The RPS algorithm has been applied directly to obtain the solution in rapidly convergent MFPS without being linearized, discretized, or exposed to perturbation. Graphs and numerical results show that the proposed method is complete reliability and performance with great potential for use in many scientific applications. The present results show that the RPS technique is a simple and quite powerful tool in finding the approximate solutions for different kinds of fractional PDEs. A comparison between the RPSM and those available in the literature are carried out through numerical examples. High agreements of numerical results are clear and remarkable.

Footnotes

Handling Editor: James Baldwin

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge Universiti Kebangsaan Malaysia for providing facilities for our research under the grant no. GP-K007788 and GP-K006926.

ORCID iDs

Mohammed Al-Smadi

Rokiah Rozita Ahmad

References

Podlubny

Fractional differential equations. San Diego, CA: Academic Press, 1999.

Millar

Ross

An introduction to the fractional calculus and fractional differential equations. New York: Wiley, 1993.

Caputo

Linear models of dissipation whose Q is almost frequency independent: part II. Geophys J Int 1967; 13: 529–539.

Al-Smadi

Simplified iterative reproducing kernel method for handling time-fractional BVPs with error estimation. Ain Shams Eng J 2018; 9: 2517–2525.

Al-Smadi

Solving fractional system of partial differential equations with parameters derivative by combining the GDTM and RDTM. Nonlin Stud 2019; 26: 587–601.

Al-Smadi

Abu Arqub

Computational algorithm for solving fredholm time-fractional partial integrodifferential equations of dirichlet functions type with error estimates. Appl Math Comput 2019; 342: 280–294.

Moaddy

Freihat

Al-Smadi

, et al. Numerical investigation for handling fractional-order Rabinovich–Fabrikant model using the multistep approach. Soft Comput 2018; 22: 773–782.

Al-Smadi

Freihat

Khalil

, et al. Numerical multistep approach for solving fractional partial differential equations. Int J Comput Method 2017; 14: 1750029.

Firoozjaee

Yousefi

Jafari

A numerical approach to Fokker-Planck equation with space- and time-fractional and non fractional derivatives. MATCH Commun Math Comput Chem 2015; 74: 449–464.

10.

Momani

Abu Arqub

Freihat

, et al. Analytical approximations for Fokker-Planck equations of fractional order in multistep schemes. Appl Comput Math 2016; 15: 319–330.

11.

Deng

Numerical algorithm for the time fractional Fokker-Planck equation. J Comput Phys 2007; 227: 1510–1522.

12.

Schneider

Wyss

Fractional diffusion and wave equations. J Math Phys 1989; 30: 134–144.

13.

Odibat

Momani

Numerical solution of Fokker-Planck equation with space-and time-fractional derivatives. Phys Lett A 2007; 369: 349–358.

14.

Abu Arqub

. Series solution of fuzzy differential equations under strongly generalized differentiability. J Adv Res Appl Math 2013; 5: 31–52.

15.

Freihet

Hasan

Al-Smadi

, et al. Construction of fractional power series solutions to fractional stiff system using residual functions algorithm. Adv Diff Equ 2019; 2019: 95.

16.

Hasan

Al-Smadi

Freihet

, et al. Two computational approaches for solving a fractional obstacle system in Hilbert space. Adv Diff Equat 2019; 2019: 55.

17.

Shammari

Al-Smadi

Abu Arqub

, et al. Adaptation of residual power series method to solve Fredholm fuzzy integro-differential equations. AIP Conf Proceed 2019; 2111: 020002.

18.

Komashynska

Al-Smadi

Abu Arqub

, et al. An efficient analytical method for solving singular initial value problems of nonlinear systems. Appl Math Inform Sci 2016; 10: 647–656.

19.

El-Ajou

Abu Arqub

Al-Smadi

A general form of the generalized Taylor’s formula with some applications. Appl Math Comput 2015; 256: 851–859.

20.

Alshammari

Al-Smadi

Hashim

, et al. Applications of fractional power series approach in solving fractional Volterra integro-differential equations. AIP Conf Proceed 2019; 2111: 020003.

21.

Alaroud

Al-Smadi

Ahmad

, et al. Computational optimization of residual power series algorithm for certain classes of fuzzy fractional differential equations. Int J Diff Equ 2018; 2018: 8686502.

22.

Hasan

Al-Smadi

Freihet

, et al. Application of power series method for solving obstacle problem of fractional order. In: Proceedings of the 2019 IEEE Jordan international joint conference on electrical engineering and information technology (JEEIT), Amman, Jordan, 9–11 April 2019, pp.513–518. New York: IEEE.

23.

Alaroud

Al-Smadi

Ahmad

, et al. Numerical computation of fractional Fredholm integro-differential equation of order 2β arising in natural sciences. J Phys Conf Ser 2019; 1212: 012022.

24.

Alaroud

Al-Smadi

Ahmad

, et al. An analytical numerical method for solving fuzzy fractional Volterra integro-differential equations. Symmetry 2019; 11: 205.

25.

Al-Smadi

Abu Arqub

Shawagfeh

, et al. Numerical investigations for systems of second-order periodic boundary value problems using reproducing kernel method. Appl Math Comput 2016; 291: 137–148.

26.

Zhang

Kumar

, et al. Residual power series method for time fractional Schrödinger equations. J Nonlinear Sci 2016; 9: 5821–5829.

27.

Kumar

Yan

SP.

Residual power series method for fractional diffusion equations. Fund Inform 2017; 151: 213–230.

28.

Hasan

Al-Zoubi

Freihet

, et al. Solution of fractional SIR epidemic model using residual power series method. Appl Math Inform Sci 2019; 13: 153–161.

29.

Kumar

Residual power series method for fractional Burger types equations. Nonlin Eng 2016; 5: 235–244.

30.

Abu Hammad

Alaroud

Abu Arqub

, et al. Solving fuzzy fractional IVPs of order 2β by residual power series algorithm. In: Proceedings of the 2019 IEEE Jordan international joint conference on electrical engineering and information technology (JEEIT), Amman, Jordan, 9–11 April 2019, pp.52–57. New York: IEEE.

Toward computational algorithm for time-fractional Fokker–Planck models

Abstract

Keywords

Introduction

Preliminaries and notations

Definition 2.1

Definition 2.2

Theorem 2.1

Definition 2.3

Theorem 2.2

Description of the RPS algorithm

Algorithm 3.1

Numerical experiments

Example 4.1

Example 4.2

Concluding remarks

Footnotes

Declaration of conflicting interests

Funding

ORCID iDs

References