Exact Solution of Fractional Black-Scholes European Option Pricing Equations

Maryeme Ouafoudi, Fei Gao

Mathematics Department, Wuhan University of Technology, Wuhan, China.

DOI: 10.4236/am.2018.91006                

Abstract

We introduce two algorithms in order to find the exact solution of the nonlinear Time-fractional Partial differential equation, in this research work. Those algorithms are proposed in the following structure: The Modified Homotopy Perturbation Method (MHPM), The Homotopy Perturbation and Sumudu Transform Method. The results achieved using the both methods are the same. However, we calculate the approached theoretical solution of the Black-Scholes model in the form of a convergent power series with a regularly calculated element. Finally, we propose a descriptive example to demonstrate the efficiency and the simplicity of the methods.

Keywords

Homotopy Perturbation Method, Modified Homotopy Perturbation Method, Sumudu Transform, Black-Scholes Equations

1. Introduction

In recent years, fractional calculus has been increasingly used for numerous applications in many scientific and technical fields such as medical sciences, biological research, as well as various chemical, biochemical and physical fields. Fractional calculus can be, for instance, employed to solve a lot of problems within the biomedical research field. Such an important application is studying membrane biophysics and polymer viscoelasticity [1] .

Numerous methods and approaches have been presented in recent years. Some of these methods are analytical such as Fourier transform method [2] , the fractional Green function method [3] , the popular Laplace method [4] [5] , the iteration method [6] , the Mellin transform method and the method of the orthogonal polynomial [4] .

Numerical methods and approaches are also popular and used to obtain approximate solutions of FPDEs. Examples of such a numerical method for solving FPDEs are the Homotopy Perturbation Method (HPM) [7] [8] [9] [10] [11] , the Differential Transform Method (DTM) [12] , the Variational Iteration Method (VIM) [13] , the New Iterative Method (NIM) [3] [14] , the Homotopy Analysis Method (HAM) [1] [6] [15] and the Adomian Decomposition Method (ADM) [16] .

Among these numerical methods, the VIM and the ADM are the most popular ones that are used to solve differential and integral equations of integer and fractional order. The HPM is a universal approach which can be used to solve both fractional and ordinary differential equations FODs as well as fractional partial differential equations FPDEs. The HPM method was originally proposed by He [8] [17] . The HPM is a coupling of homotopy in the topology and the perturbation method. The method is used to solve various types of equations such us the Hellmholtz equation, the fifth order KdV [6] , the Kleein-Gorden Equation [4] , the Fokker-Plank equation [9] [18] , the nonlinear Kolmogorov- Petrovskii-Piskunov Equation [3] as well as other types of equations as proposed and used in [10] [19] .

In this paper, we present two approaches to derive the exact solution of various types of FPDEs. Those methods are: The Modified Homotopy Perturbation Method (MHPM) and the Homotopy Perturbation and Sumudu Transform Method (HPSTM). The MHPM is a fast approach that’s based on designing and utilizing a proper initial approximation which satisfies the initial condition of the HPM. However, the HPSTM is a combination of the homotopy perturbation method and Sumudu transform. The paper is structured in six sections. In Section 2, we begin with an introduction to some necessary definitions of fractional calculus theory. In Section 3, we describe the basic idea of the HPM. In Section 4, we describe the MHPM. In Section 5, we describe the HPSTM. In Section 6, we present two examples to show the efficiency of using the MHPM and HPSTM to solve the fractional Black-Scholes equation. Finally, relevant conclusions are drawn in Section 7.

2. Basic Definitions of Fractional Calculus

In this section, we present the basic definitions and properties of the fractional calculus theory, which are used further in this paper.

Definition 2.1

A real function $f\left(t\right),t>0$ is said to be in the space $C_{\sigma },\sigma \in ℝ$ , if there exists a real number $p>\sigma$ such that $f\left(t\right)=t^p{f}_1\left(t\right)$ where $f_1\left(t\right)\in C\left[0,\infty \right]$ and it is said to be in the space $C_{\sigma }^m$ if $f^m\in C_{\sigma },m\in \mathbb{N}$ .

Definition 2.2

The left sided Riemann-Louiville fractional integral of order $\alpha \ge 0$ , of a function $f\in C_{\sigma },\sigma \ge -1$ is defined as:

$J_t^{\alpha }f\left(t\right)=\frac{1}{\text{Γ}\left(\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\alpha \right)}^{\alpha -1}f\left(\zeta \right)\text{d}\zeta }$ (1)

where $\alpha >0$ , $t>0$ , and $\Gamma \left(m\right)$ is the Gamma function. Properties of the operator $J_t^{\alpha }$ for $f\in C_{\mu }$ , $\mu \ge -1$ , $\alpha ,\beta \ge 0$ , $\gamma \ge -1$ , can be founded in [20] , and are defined as follow

1) $J_t^{0}f\left(t\right)=f(\; t\; )$

2) $J_t^{\alpha }f\left(t\right)J_t^{\beta }f\left(t\right)=J_t^{\alpha +\beta }f(\; t\; )$

3) $J_t^{\alpha }f\left(t\right)J_t^{\beta }f\left(t\right)=J_t^{\beta }f\left(t\right)J_t^{\alpha }f(\; t\; )$

4) $J_t^{\beta }{t}^{\gamma }=\frac{\text{Γ}\left(\gamma +1\right)}{\text{Γ}\left(\alpha +\gamma +1\right)}{t}^{\alpha +\gamma +1}$

Definition 2.3

Let $f\in C_{\mu }^m$ , $n\in \mathbb{N}\cup \left\{0\right\}$ . The left-Sided Caputo fractional of $f$ in the Caputo sense is defined by [4] as follows

$D_t^{\alpha }f\left(t\right)=f\left(x\right)=\{\begin{array}{l}\frac{1}{\text{Γ}\left(n-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\alpha \right)}^{n-1-\alpha }{f}^n\left(\zeta \right)\text{d}\zeta },n-1<\alpha \le n\\ D_t^{n}f\left(t\right),\alpha =n\end{array}$ (2)

According to the [20] , Equations (1) and (2) becomes, for $n-1<\alpha \le n$

$J_t^{\alpha }f\left(t\right)=\frac{1}{\text{Γ}\left(\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\zeta \right)}^{\alpha -1}f\left(\zeta \right)\text{d}\zeta }$ (3)

$D_t^{\alpha }f\left(t\right)=\frac{1}{\text{Γ}\left(n-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\alpha \right)}^{n-1-\alpha }{f}^n\left(\zeta \right)\text{d}\zeta }$ (4)

Definition 2.4

The single parameter and the two parameters variants of the Mitting-Leffler function are denoted by $E_{\alpha }\left(t\right)$ and $E_{\alpha ,\beta }\left(t\right)$ , respectively, which are relevant for their connection with fractional calculus, and are defined as:

$E_{\alpha }\left(t\right)={\displaystyle {\sum }_{j=0}^{\infty }\frac{t^j}{\Gamma \left(\alpha j+1\right)}},\alpha >0,t\in ℂ$ (5)

$E_{\alpha ,\beta }\left(t\right)={\displaystyle {\sum }_{j=0}^{\infty }\frac{t^j}{\Gamma \left(\alpha j+\beta \right)}},\alpha ,\beta >0,t\in ℂ$ (6)

the k-th derivatives are

$E_{\alpha }\left(t\right)=\frac{{\text{d}}^k}{\text{d}{t}^k}{E}_{\alpha ,\beta }\left(t\right)={\displaystyle {\sum }_{j=0}^{\infty }\frac{\left(k+j\right)!t^j}{j!\Gamma \left(\alpha j+\alpha k+1\right)}},k=0,1,2,\cdots ,$ (7)

$E_{\alpha ,\beta }\left(t\right)=\frac{{\text{d}}^k}{\text{d}{t}^k}{E}_{\alpha ,\beta }\left(t\right)={\displaystyle {\sum }_{j=0}^{\infty }\frac{\left(k+j\right)!t^j}{\Gamma \left(\alpha j+\alpha k+\beta \right)}},k=0,1,2,\cdots ,$ (8)

some special cases of the Mitting-Leffler function are as follows

1) $E_1\left(t\right)={\text{e}}^t$

2) $E_{\alpha ,1}\left(t\right)=E_{\alpha }(\; t\; )$

3) $\frac{{\text{d}}^k}{\text{d}{t}^k}\left[t^{\beta -1}{E}_{\alpha ,\beta }\left(a{t}^{\alpha }\right)\right]=t^{\beta -\alpha -1}{E}_{\alpha ,\beta -k}\left(a{t}^{\alpha }\right)$

Other properties of the Mitting-Leffler functions can be found in [20] . These functions, are generalizations of the exponential function, because, most linear differential equations of fractional order have solutions that are expressed in terms of these functions.

Definition 2.5

Sumudu transform over the following set of functions

$A=\left\{f\left(t\right)|\exists M>0,{\tau }_1,{\tau }_2>0,\left|f\left(t\right)\right|>M{\text{e}}^{\frac{\left|t\right|}{{\tau }_j}}\text{if}t\in {\left(-1\right)}^j\times \left[0,\infty \right]\right\}$ (9)

is defined by

$\mathbb{S}\left[f\left(t\right)\right]=G\left(u\right){\displaystyle {\int }_0^{\infty }f\left(ut\right){\text{e}}^{-t}\text{d}t}$ (10)

where $u\in \left({\tau }_1,{\tau }_2\right)$ . Some special properties of the Sumudu transform are as follows

1) $\mathbb{S}\left[1\right]=1$

2) $\mathbb{S}\left[\frac{t^n}{\text{Γ}\left(n+1\right)}\right]=u^n,n>0$

3) $\mathbb{S}\left[{\text{e}}^{at}\right]=\frac{1}{1-au}$

4) $\mathbb{S}\left[\alpha f\left(x\right)+\beta g\left(x\right)\right]=\alpha \mathbb{S}\left[f\left(x\right)\right]+\beta \mathbb{S}\left[g\left(x\right)\right]$

Other properties of the Sumudu transform can be found in [21] .

Definition 2.6

$G\left(u\right)$ is the Sumudu transform of $f\left(t\right)$ . And according to [22] we have:

1) $\frac{G\left(1/s\right)}{s}$ is a meromorphic function, with singularities having $Re\left(s\right)<\gamma$ .

2) There exists a circular region $\Gamma$ with radius R and positive constants, M and k, with

$\left|\frac{G\left(1/s\right)}{s}\right|<M{R}^{-k}$

then the function $f\left(t\right)$ is given by

${\mathbb{S}}^{-1}\left[G\left(s\right)\right]=\frac{1}{2\text{π}i}{\displaystyle {\int }_{\gamma -i\infty }^{\gamma +i\infty }{\text{e}}^{st}G\left(\frac{1}{s}\right)\frac{\text{d}s}{s}}={{\displaystyle \sum }}^{\text{}}\text{residuse}\left[{\text{e}}^{st}\frac{G\left(\frac{1}{s}\right)}{s}\right]$ (11)

Definition 2.6

The Sumudu transform $\mathbb{S}$ , of the Caputo functional integral is defined as [22]

$\mathbb{S}\left[D_t^{\alpha }f\left(t\right)\right]=\frac{G\left(u\right)}{u^{\alpha }}-\underset{k=0}{\overset{n-1}{{\displaystyle \sum }}}\frac{f^{\left(k\right)}\left(0\right)}{u^{\alpha -k}}$ (12)

then it can be easily understood that

$\mathbb{S}\left[D_t^{\alpha }f\left(t\right)\right]=\frac{\mathbb{S}\left[f\left(x,t\right)\right]}{u^{\alpha }}-\underset{k=0}{\overset{n-1}{{\displaystyle \sum }}}\frac{f^{\left(k\right)}\left(0\right)}{u^{\alpha -k}},n-1\le \alpha <n$ (13)

Theorem 2.1

Consider the following n-term liner fractional differential equation [7] theorem

$\left(a_n{D}_t^{{\beta }_n}+a_{n-1}{D}_t^{{\beta }_{n-1}}+\cdots +a_0{D}_t^{{\beta }_0}\right){\omega }_t=f\left(t\right)$ (14)

with the constant initial condition

${\omega }^{j_i}=C_{i{j}_i},i=0,1,\cdots ,n,j=1,2,\cdots ,l_i$

where

$a_i,C_{i{j}_i}\in ℝ,n_i-1<{\beta }_i\le n_i,n_i\in \mathbb{N}$ and ${\beta }_0<{\beta }_1<{\beta }_2<\cdots <{\beta }_{n-1}<n\le {\beta }_n<n+1$ .

Then we see that the analytical general solution of Equation (14) is

$\omega \left(t\right)={\displaystyle {\int }_0^t{G}_n\left(t-\zeta \right)f\left(\zeta \right)\text{d}\zeta }+\underset{i=0}{\overset{\infty }{{\displaystyle \sum }}}\underset{j_i=0}{\overset{l_i-1}{{\displaystyle \sum }}}{a}_i{C}_{i{j}_i}{G}_n^{B_i-j_i-1}\left(t\right)$ (15)

where $G_n$ is the Green function and it’s defined by

$\begin{array}{c}{G}_n\left(t\right)=\frac{1}{a_n}\underset{m=0}{\overset{\infty }{{\displaystyle \sum }}}\frac{{\left(-1\right)}^m}{m!}\times \underset{k_0,k_1,\cdots ,k_{n-2}\le 0,k_0+k_1+\cdots +k_{n-2}=m}{{\displaystyle \sum }}\left(m;k_0,k_1,\cdots ,k_{n-2}\right)\\ \times \underset{p=0}{\overset{n-2}{{\displaystyle \prod }}}{\left(\frac{a_p}{a_n}\right)}^{k_p}{t}^{\left({\beta }_n-{\beta }_{n-1}\right)m+{\beta }_n+{\displaystyle {\sum }_{j=0}^{n-2}\left({\beta }_{n-1}-{\beta }_j\right)k_j}-1}\\ \times E_{{\beta }_n-{\beta }_{n-1},{\beta }_n+{\displaystyle {\sum }_{j=0}^{n-2}\left({\beta }_{n-1}-{\beta }_j\right)k_j}-1}^m\left(-\frac{a_{n-1}}{a_n}{D}^{{\beta }_n-{\beta }_{n-1}}\right)\end{array}$ (16)

where $\left(m;k_0,k_1,\cdots ,k_{n-2}\right)=\frac{m!}{k_0!k_{01}!\cdots k_{n-2}!}$ .

And $E_{(.),(.)}^m$ is the m-th derivative of the Mitting-Leffler function.

In a special case of the latter theorem, the following relaxation-oscillation [7] is solved:

$D_t^{\alpha }\omega \left(t\right)+A\omega \left(t\right)=f\left(t\right),t>0,{\omega }^i\left(0\right)=b_i,i=1,2,\cdots ,n-1$ (17)

where $b_i$ are real constants, and $n-1<\alpha \le n$ . By utilizing theorem 2.1, we obtain the solution of Equation (14) as follows

$\omega \left(t\right)={\displaystyle {\int }_0^t{G}_2\left(t-\zeta \right)f\left(\zeta \right)}+\underset{j=0}{\overset{n-1}{{\displaystyle \sum }}}{b}_j{D}_t^{\alpha -j-1}{G}_2\left(t\right)$ (18)

where $G_2\left(t\right)=t^{\alpha }{E}_{\alpha ,\alpha }\left(-A{t}^{\alpha }\right)$ . It is easy to see that if $0<\alpha \le 1$ , then the solution of the Equation (14) is given as follows

$\omega \left(t\right)={\displaystyle {\int }_0^t{G}_2\left(t-\zeta \right)f\left(\zeta \right)\text{d}\zeta }+b_0{D}_t^{\alpha -1}G\left(t\right)$ (19)

which will be used in the coming examples, discussed in this work.

3. The Basic Idea of the HPM

In this section, we will briefly present the main idea of the HPM. At first, we will consider the following nonlinear differential equation

$A\left(u\right)-f\left(x\right)=0,x\in \Omega$ (20)

where $A,B,f\left(x\right)$ and Gamma are a general differential function operator, a boundary operator, a known analytical function and the boundary of the domain $\Omega$ , respectively. The operator A can be decomposed into a linear operator denoted by L, and nonlinear operator denoted by N. Hence, the equation can be written as follows

$L\left(u\right)+N\left(u\right)-f\left(x\right)=0,x\in \Omega$ (21)

As a result, we construct a homotopy $v\left(x,p\right):\Omega \times \left[0,1\right]\to ℝ$ that satisfies:

$\mathcal{H}\left(v,p\right)=\left(1-p\right)\left[L\left(u\right)+N\left(u_0\right)\right]+p\left[A\left(u\right)-f\left(x\right)\right]=0,0\le p<1$ (22)

which equivalent to

$\mathcal{H}\left(v,p\right)=L\left(v\right)-N\left(u_0\right)+pL\left(u_0\right)+p\left[N\left(v\right)+f\left(x\right)\right]=0,0\le p<1$ (23)

when the value of p is varied from zero to one, we can see easily that

$\mathcal{H}\left(v,p\right)=L\left(v\right)-N\left(u_0\right)=0$ (24)

$\mathcal{H}\left(v,p\right)=L\left(v\right)-N\left(u\right)-f\left(x\right)=A\left(u\right)-f\left(x\right)=0$ (25)

If the parameter p is assumed as small, then the solution of the Equation (23) can be expressed as a power series in p as follows

$v=v_0+p{v}_1+p^2{v}_2+p^3{v}_3+\cdots$ (26)

The best approximation for the solution of Equation (23) is

$\omega ={\lim }_{p\to 1}v={\displaystyle {\sum }_{i=0}^{\infty }v}=v_0+v_1+v_2+v_3+\cdots$ (27)

4. The Modified Homotopy Perturbation Method (MHPM)

Meanwhile, we are able to apply the HPM to solve the class of time-fractional partial differential equations defined by

$D_t^{\alpha }\left(\omega \left(x,t\right)\right)=L\left(\omega \left(x,t\right)\right)+N\left(\omega \left(x,t\right)\right)+f\left(\omega \left(x,t\right)\right)$ (28)

subject to the initial condition

$\omega \left(x,0\right)=h\left(x\right)$ (29)

where $0<\alpha \le 1$ . The constructed homotopy satisfies

$\begin{array}{c}\mathcal{H}\left(\omega \left(x,t\right)\right)=\left(1-p\right)\left[D_t^{\alpha }\left(\omega \left(x,t\right)\right)-D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)\right]\\ +p\left[D_t^{\alpha }\left(\omega \left(x,t\right)\right)-L\left(\omega \left(x,t\right)\right)-N\left(\omega \left(x,t\right)\right)-f\left(\omega \left(x,t\right)\right)\right]\end{array}$ (30)

By simplifying Equation (30) we get

$\begin{array}{c}{D}_t^{\alpha }\left(\omega \left(x,t\right)\right)=D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)+p[D_t^{\alpha }\left(\omega \left(x,t\right)\right)-L\left(\omega \left(x,t\right)\right)\\ -N\left(\omega \left(x,t\right)\right)-f\left(\omega \left(x,t\right)\right)]\end{array}$ (31)

where the embedding parameter p is considered to be small, and applied to the classical perturbation technique. The next step is to use the homotopy parameter p to expand the solution into the following form

$\omega \left(x,t\right)=\underset{k=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^k{\omega }_k\left(x,t\right)$ (32)

For the sake of clarifying the solution procedure of this method, we consider a general nonlinear time-fractional partial differential equation:

$D_t^{\alpha }\omega \left(x,t\right)=L\omega \left(x,t\right)+N\omega \left(x,t\right)+f\left(x,t\right)$ (33)

By substituting Equation (32) into Equation (31) and equating the terms with identical power of p, we can obtain a series of equations as the following

$\begin{array}{l}{p}^0:D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)=D_t^{\alpha }\left(h\left(x\right)\right)\\ p^1:D_t^{\alpha }\left({\omega }_1\left(x,t\right)\right)=D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)-L{\omega }_0\left(x,t\right)-N{\omega }_0\left(x,t\right)-f\left(x,t\right)\\ p^2:D_t^{\alpha }\left({\omega }_2\left(x,t\right)\right)=D_t^{\alpha }\left({\omega }_1\left(x,t\right)\right)-L{\omega }_1\left(x,t\right)-N{\omega }_1\left(x,t\right)-f\left(x,t\right)\\ ⋮\end{array}$ (34)

Applying the Riemann-Louiville fractional integral of order $\alpha \ge 0,J_t^{\alpha }$ in the both sides

$\begin{array}{l}{\omega }_0\left(x,t\right)=h\left(x\right)\\ {\omega }_1\left(x,t\right)=J_t^{\alpha }\left[D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)-L{\omega }_0\left(x,t\right)\left(x,t\right)-N{\omega }_0\left(x,t\right)-f\left(x,t\right)\right]\\ {\omega }_2\left(x,t\right)=J_t^{\alpha }\left[D_t^{\alpha }\left({\omega }_1\left(x,t\right)\right)-L{\omega }_1\left(x,t\right)\left(x,t\right)-N{\omega }_1\left(x,t\right)-f\left(x,t\right)\right]\\ ⋮\end{array}$ (35)

Substitute the results of the Equation (35) into the Equation (32), and applying the Perturbation technique [15] [22] we obtain an authentic n-th approximation of the exact solution, given by

$\omega \left(x,t\right)=\underset{i=0}{\overset{n}{{\displaystyle \sum }}}{\omega }_i\left(x,t\right)$ (36)

If there exist some terms $vn=0,n>0$ , then the exact solution can be written in the following form

$\omega \left(x,t\right)=\underset{i=0}{\overset{n-1}{{\displaystyle \sum }}}{\omega }_i\left(x,t\right)$ (37)

We Assume that the initial approximation of Equation (28) is given by

${\omega }_0\left(x,t\right)=\omega \left(x,0\right)c_1\left(t\right)+\omega \left(x\right)c_2\left(t\right)$ (38)

where $\omega \left(x,0\right)$ is the initial condition of the Equation (28), and $\omega \left(x\right)=\frac{\partial \omega }{\partial t}\left(x,0\right)$ .

The aim of this algorithm is to find the terms $c_1$ and $c_2$ . We assume that ${\omega }_n\equiv 0$ , $\forall n\ge 0$ , which means that the exact solution is given by

$\omega \left(x,t\right)={\omega }_0\left(x,t\right)$ (39)

As $\omega \left(x,t\right)$ satisfies also the initial condition, we get

$\omega \left(x,0\right)={\omega }_0\left(x,0\right)=\omega \left(x,0\right)c_1\left(0\right)+\omega \left(x\right)c_2\left(0\right)=h\left(x\right)$ (40)

and $c_1\left(0\right)=1,c_2\left(0\right)=0$ , on the other part we have,

$D_t^{\alpha }\left({\omega }_1\left(x,t\right)\right)=D_t^{\alpha }\left({\omega }_0\left(x,t\right)\right)-L{\omega }_0\left(x,t\right)-N{\omega }_0\left(x,t\right)-f\left(x,t\right)\equiv 0$ (41)

By substitution the Equation (38) into the Equation (28), we obtain

$\begin{array}{l}\omega \left(x,0\right)D_t^{\alpha }\left(c_1\left(t\right)\right)+\omega \left(x\right)D_t^{\alpha }\left(c_2\left(t\right)\right)\\ =L\left(\omega \left(x,0\right)c_1\left(t\right)+\omega \left(x\right)c_2\left(t\right)\right)\\ +N\left(\omega \left(x,0\right)c_1\left(t\right)+\omega \left(x\right)c_2\left(t\right)\right)+f\left(x,t\right)\end{array}$ (42)

Therein, the FPDEs are changed to a fractional order differential equation, which simplifies the problem. When $\alpha$ is an integer, the equation is transformed to an ordinary differential equation.

5. The Homotopy Perturbation Sumudu Transform Method (HPSTM)

We illustrate the basic idea of HPSTM by considering the general time-fractional nonlinear non-homogeneous partial differential equation with the initial condition of the general form defined by:

$D_t^{\alpha }\omega \left(x,t\right)=L\omega \left(x,t\right)+N\omega \left(x,t\right)+g\left(x,t\right)$ (43)

where $\alpha >0$ and $n-1\le \alpha <n$ .

With the initial condition

$D^m\left(\omega \right)\left(0\right)=\{\begin{array}{l}{f}_m\left(x\right),\left(m=0,1,\cdots ,n-1\right)\\ 0,m=n\equiv \left[\alpha \right]\end{array}$ (44)

where $\omega :=\omega \left(x,t\right)$ .

Here by applying Sumudu Transform on the both sides of the Equation (43), we get

$\mathbb{S}\left[D_t^{\alpha }\omega \left(x,t\right)\right]=\mathbb{S}\left[L\omega \left(x,t\right)\right]+\mathbb{S}\left[N\omega \left(x,t\right)\right]+\mathbb{S}\left[g\left(x,t\right)\right]$ (45)

which, upon using a property of the Sumudu transform, yields

$\mathbb{S}\left[\omega \right]=f\left(x\right)+u^{\alpha }\mathbb{S}\left[L\omega \right]+u^{\alpha }\mathbb{S}\left[N\omega \right]+u^{\alpha }\mathbb{S}\left[g\left(x,t\right)\right]$ (46)

where f is a function of x.

Taking the inverse of Sumudu on the both sides, we obtain

$\begin{array}{c}\omega \left(x,t\right)={\mathbb{S}}^{-1}\left[{\displaystyle {\sum }_{k=0}^{n-1}{u}^k{\omega }^{\left(k\right)}f\left(x,0\right)}\right]\\ +p{\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[L\omega \left(x,t\right)+N\omega \left(x,t\right)+g\left(x,t\right)\right]\right]\end{array}$ (47)

Now, we apply the HPM to the Equation (47) to get the solution

$\omega \left(x,t\right)=\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\omega }_n\left(x,t\right)$ (48)

where $p\in \left[0,1\right]$ is the embedding (or homotopy) parameter. ${\omega }_0$ is an initial approximation which satisfies the boundary conditions, and ${\omega }_n,n\in \mathbb{N}$ are the nth order approximations which are functions yet to be determined. It is noted that setting $p=1$ in the Equation (48) gives an approximate solution of the given nonlinear differential equation.

The nonlinear term can be decomposed as follows

$N\left(\omega \left(x,t\right)\right)=\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{H}_n\left({\omega }_0,{\omega }_1,\cdots ,{\omega }_n\right)$ (49)

where $H_n$ are the He’s polynomials given by,

$H_n\left({\omega }_0,{\omega }_1,\cdots ,{\omega }_n\right)=\frac{1}{n!}\frac{{\partial }^n}{\partial p^n}N\left({\displaystyle {\sum }_{k=0}^n{p}^k{\omega }_k\left(x,t\right)}\right)$ (50)

By substituting Equation (48) and Equation (49) into Equation (47), we obtain the following equation

$\begin{array}{c}{\displaystyle {\sum }_{n=0}^{\infty }{p}^n{\omega }_n\left(x,t\right)}={\mathbb{S}}^{-1}\left[{\displaystyle {\sum }_{k=0}^{n-1}{u}^k{f}_k\left(x\right)}\right]+p{\mathbb{S}}^{-1}[u^{\alpha }\mathbb{S}[L\left({\displaystyle {\sum }_{n=0}^{\infty }{p}^n{\omega }_n\left(x,t\right)}\right)\\ +{\displaystyle {\sum }_{n=0}^{\infty }{p}^n{H}_n}+g\left(x,t\right)]]\end{array}$ (51)

In the meantime, we have completed the coupling of the Sumudu transform and the HPM using He’s polynomials. By equating the terms of the same power of p in the Equation (51), we obtain the given approximation for ${\omega }_n,n\in \mathbb{N}$ .

$\begin{array}{l}{p}^0:{\omega }_0\left(x,t\right)={\mathbb{S}}^{-1}\left[\underset{k=0}{\overset{n-1}{{\displaystyle \sum }}}{u}^k{f}_k\left(x\right)\right]\\ p^1:{\omega }_1\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[L\left({\omega }_0\right)+H_0\left({\omega }_0\right)\right]\right]\\ p^2:{\omega }_2\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[L\left({\omega }_1\right)+H_1\left({\omega }_0,{\omega }_1\right)\right]\right]\\ ⋮\\ p^n:{\omega }_n\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[L\left({\omega }_{n-1}\right)+H_1\left({\omega }_0,{\omega }_1,{\omega }_2,\cdots ,{\omega }_{n-1}\right)\right]\right]\end{array}$ (52)

Finally, the approximate and analytical solution of the Equation (43) is given by truncating the following series

$\omega \left(x,t\right)=\underset{p\to 1}{\lim }\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\omega }_n\left(x,t\right)=\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{\omega }_n\left(x,t\right)$ (53)

6. Applications

In order to apply the MHPM and the HPSTM in Sections 4 and 5, respectively, consider the following fractional Black-Scholes option pricing equations as follows:

6.1. Example 1

$\frac{{\partial }^{\alpha }\psi }{\partial x^{\alpha }}=\frac{{\partial }^2\psi }{\partial x^2}+\left(k-1\right)\frac{\partial \psi }{\partial x}-k\psi$ (51)

$\psi \left(x,0\right)=\max \left({\text{e}}^x-1,0\right)$ (55)

6.1.1. Use the MHPM

By applying the MHPM, the initial approximation will be given by

${\psi }_0\left(x,t\right)=\psi \left(x,0\right)c_1\left(t\right)+\psi \left(x\right)c_2\left(t\right)$ (56)

where

$\frac{\partial \psi }{\partial x}\left(x,0\right)=\max \left({\text{e}}^x,0\right)$ (57)

therefore

${\psi }_0\left(x,t\right)=\max \left({\text{e}}^x-1,0\right)c_1\left(t\right)+\max \left({\text{e}}^x,0\right)c_2\left(t\right)$ (58)

Hence,

$\begin{array}{l}{D}_t^{\alpha }{\psi }_1=D_t^{\alpha }\left(\max \left({\text{e}}^x-1,0\right)c_1\left(t\right)+\max \left({\text{e}}^x,0\right)c_2\left(t\right)\right)\\ -\frac{{\partial }^2}{\partial x^2}\left(\max \left({\text{e}}^x-1,0\right)c_1\left(t\right)+\max \left({\text{e}}^x,0\right)c_2\left(t\right)\right)\\ +\left(k-1\right)\frac{\partial }{\partial x}\left(\max \left({\text{e}}^x-1,0\right)c_1\left(t\right)+\max \left({\text{e}}^x,0\right)c_2\left(t\right)\right)\\ -k\left(\max \left({\text{e}}^x-1,0\right)c_1\left(t\right)+\max \left({\text{e}}^x,0\right)c_2\left(t\right)\right)\equiv 0.\end{array}$ (59)

We obtain the fractional system

$\{\begin{array}{l}{D}_t^{\alpha }{c}_1\left(t\right)+k{c}_1\left(t\right)=0\\ c_1\left(0\right)=0\end{array}$ (60)

$\{\begin{array}{l}{D}_t^{\alpha }{c}_1\left(t\right)+D_t^{\alpha }{c}_2\left(t\right)=0\\ c_2\left(0\right)=0\end{array}$ (61)

Solving the Equation (60) and Equation (61) by applying the theorem 2.1, we obtain

$\{\begin{array}{l}{c}_1\left(t\right)=E_{\alpha }\left(-k{t}^{\alpha }\right)\\ c_2\left(t\right)=1-E_{\alpha }\left(-k{t}^{\alpha }\right)\end{array}$ (62)

and the exact solution is

$\psi \left(x,t\right)=\max \left({\text{e}}^x-1,0\right)E_{\alpha }\left(-k{t}^{\alpha }\right)+\max \left({\text{e}}^x,0\right)\left(1-E_{\alpha }\left(-k{t}^{\alpha }\right)\right)$ (63)

If we put $\alpha \to 1$ in Equation (63) or solving Equation (60) and Equation (61) with $\alpha =1$ , we obtain the exact solution

$\psi \left(x,t\right)=\max \left({\text{e}}^x-1,0\right){\text{e}}^{\left(-k{t}^{\alpha }\right)}+\max \left({\text{e}}^x,0\right)\left(1-{\text{e}}^{\left(-k{t}^{\alpha }\right)}\right)$ (64)

6.1.2. Use the HPSTM

Applying the Sumudu transform to the both sides of the Equation (54), we get

$\mathbb{S}\left[\psi \left(x,t\right)\right]=\max \left({\text{e}}^x-1,0\right)+u^{\alpha }\mathbb{S}\left[\frac{{\partial }^2\psi }{\partial x^2}+\left(k-1\right)\frac{\partial \psi }{\partial x}-k\psi \right]$ (65)

Operating the inverse Sumudu transform on both sides in the Equation (65), we have

$\psi \left(x,t\right)=\max \left({\text{e}}^x-1,0\right)+p{\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[\frac{{\partial }^2\psi }{\partial x^2}+\left(k-1\right)\frac{\partial \psi }{\partial x}-k\psi \right]\right]$ (66)

Now, applying the HPM

${\displaystyle {\sum }_{n=0}^{\infty }{p}^n{\psi }_n\left(x,t\right)}=\max \left({\text{e}}^x-1,0\right)+p\left({\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\displaystyle {\sum }_{n=0}^{\infty }{p}^n{\Psi }_n\left(x,t\right)}\right]\right]\right)$ (67)

where

${\Psi }_n\left(x,t\right)=\frac{{\partial }^2{\psi }_n}{\partial x^2}+\left(k-1\right)\frac{\partial {\psi }_n}{\partial x}-k{\psi }_n$ (68)

Equating the corresponding power where,
of p on both sides in Equation (67), we obtain

$p^0:{\psi }_0\left(x,t\right)=\max \left({\text{e}}^x-1,0\right)$

$\begin{array}{l}{p}^1:{\psi }_1\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\Psi }_0\left(x,t\right)\right]\right]\\ ={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[k\max \left({\text{e}}^x,0\right)-k\max \left({\text{e}}^x-1,0\right)\right]\right]\\ =k\frac{t^{\alpha }}{\Gamma \left(\alpha +1\right)}\max \left({\text{e}}^x,0\right)-k\frac{t^{\alpha }}{\Gamma \left(\alpha +1\right)}\max \left({\text{e}}^x-1,0\right)\end{array}$

$\begin{array}{l}{p}^2:{\psi }_2\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\text{Ψ}}_1\left(x,t\right)\right]\right]\\ =\frac{{\left(kt\right)}^{2\alpha }}{\text{Γ}\left(2\alpha +1\right)}\max \left({\text{e}}^x,0\right)-k\frac{{\left(-kt\right)}^{2\alpha }}{\text{Γ}\left(2\alpha +1\right)}\max \left({\text{e}}^x-1,0\right)\end{array}$

$⋮$

$\begin{array}{l}{p}^n:{\psi }_n\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\Psi }_{n-1}\left(x,t\right)\right]\right]\\ =\frac{{\left(kt\right)}^{n\alpha }}{\Gamma \left(n\alpha +1\right)}\max \left({\text{e}}^x,0\right)-k\frac{{\left(-kt\right)}^{n\alpha }}{\Gamma \left(n\alpha +1\right)}\max \left({\text{e}}^x-1,0\right)\end{array}$ (69)

So that the solution $\psi \left(x,t\right)$ of the problem is given by

$\begin{array}{c}\psi \left(x,t\right)=\underset{p\to 1}{\lim }\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\psi }_n\left(x,t\right)\\ =\max \left({\text{e}}^x-1,0\right)E_{\alpha }\left(-k{t}^{\alpha }\right)+\max \left({\text{e}}^x,0\right)\left(1-E_{\alpha }\left(-k{t}^{\alpha }\right)\right)\end{array}$ (70)

where $E_{\alpha }\left(x\right)$ is the Mittag-Leffler function in one parameter. For special case, $\alpha =1$ , we get

$\psi \left(x,t\right)=\max \left({\text{e}}^x-1,0\right){\text{e}}^{\left(-k{t}^{\alpha }\right)}+\max \left({\text{e}}^x,0\right)\left(1-{\text{e}}^{\left(-k{t}^{\alpha }\right)}\right)$ (71)

6.2. Example 2

$\frac{{\partial }^{\alpha }\psi }{\partial x^{\alpha }}+0.08{\left(2+\sin x\right)}^{2}x\frac{{\partial }^2\psi }{\partial x^2}+0.06x\frac{\partial \psi }{\partial x}-0.06\psi =0,0<\alpha \le 1$ (72)

Subject to the initial condition,

$\psi \left(x,0\right)=\max \left(x-25{\text{e}}^x,0\right)$ (73)

6.2.1. Use the MHPM

By applying the MHPM, the initial approximation will be given by

${\psi }_0\left(x,t\right)=\psi \left(x,0\right)c_1\left(t\right)+\psi \left(x\right)c_2\left(t\right)$ (74)

where,

$\frac{\partial \psi }{\partial x}\left(x,0\right)=1$ (75)

Therefore,

${\psi }_0\left(x,t\right)=\max \left(x-25{\text{e}}^x,0\right)c_1\left(t\right)+c_2\left(t\right)$ (76)

Thus,

$\begin{array}{l}\max \left(x-25{\text{e}}^x,0\right)+D_t^{\alpha }{c}_1\left(t\right)+D_t^{\alpha }{c}_2\left(t\right)+0.06c_1\left(t\right)\\ -0.06\max \left(x-25{\text{e}}^x,0\right)c_1\left(t\right)-0.06c_2\left(t\right)=0\end{array}$ (77)

We obtain the final system

$\{\begin{array}{l}{D}_t^{\alpha }{c}_1\left(t\right)-0.06c_1\left(t\right)=0\\ c_1\left(0\right)=0\end{array}$ (78)

$\{\begin{array}{l}{D}_t^{\alpha }{c}_2\left(t\right)+0.06x{c}_1\left(t\right)-0.06c_2\left(t\right)=0\\ c_2\left(0\right)=0\end{array}$ (79)

Solving the Equation (49) and Equation (50), we get

$\{\begin{array}{l}{c}_1\left(t\right)=E_{\alpha }\left(0.06t^{\alpha }\right)\\ c_2\left(t\right)=x\left(1-E_{\alpha }\left(0.06t^{\alpha }\right)\right)\end{array}$ (80)

Hence,

$\psi \left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right)E_{\alpha }\left(0.06t^{\alpha }\right)+x\left(1-E_{\alpha }\left(0.06t^{\alpha }\right)\right)$ (81)

This is the exact solution of the given option pricing Equation (72). The solution at the special case $\alpha =1$ , is given as follows

$\psi \left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right){\text{e}}^{0.06t}+x\left(1-{\text{e}}^{0.06t}\right)$ (82)

6.2.2. Use the HPSTM

Applying the Sumudu transform on both sides of the Equation subject to initial condition given in Equation , we get

$\begin{array}{c}\mathbb{S}\left[\psi \left(x,t\right)\right]=\max \left(x-25{\text{e}}^{-0.06},0\right)\\ -u^{\alpha }\mathbb{S}\left[0.08{\left(2+\sin x\right)}^{2}x\frac{{\partial }^2\psi }{\partial x^2}+0.06x\frac{\partial \psi }{\partial x}-0.06\psi \right]\end{array}$ (83)

Operating the inverse Sumudu transform on both sides in the Equation (83), we have

$\begin{array}{c}\psi \left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right)\\ -p{\mathbb{S}}^{-1}\left\{u^{\alpha }\mathbb{S}\left[0.08{\left(2+\sin x\right)}^{2}x\frac{{\partial }^2\psi }{\partial x^2}+0.06x\frac{\partial \psi }{\partial x}-0.06\psi \right]\right\}\end{array}$ (84)

Now, applying the homotopy perturbation method we have

$\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\psi }_n\left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right)-p{\mathbb{S}}^{-1}\left\{u^{\alpha }\mathbb{S}\left[\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\Psi }_n\left(x,t\right)\right]\right\}$ (85)

where,

${\Psi }_n\left(x,t\right)=0.08{\left(2+\sin x\right)}^{2}x\frac{{\partial }^2{\psi }_n}{\partial x^2}+0.06x\frac{\partial {\psi }_n}{\partial x}-0.06{\psi }_n$ (86)

Equating the corresponding power of p on both sides in Equation (85), we obtain

$p^0:{\psi }_0\left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right)$

$\begin{array}{l}{p}^1:{\psi }_1\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\Psi }_0\left(x,t\right)\right]\right]\\ ={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[0.06x-0.06\max \left(x-25{\text{e}}^{-0.06},0\right)\right]\right]\\ =x\frac{-0.06t^{\alpha }}{\Gamma \left(\alpha +1\right)}+\frac{0.06t^{\alpha }}{\Gamma \left(\alpha +1\right)}\max \left(x-25{\text{e}}^{-0.06},0\right)\end{array}$

$\begin{array}{l}{p}^2:{\psi }_2\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\Psi }_1\left(x,t\right)\right]\right]\\ =x\frac{{\left(-0.06t\right)}^{2\alpha }}{\text{Γ}\left(2\alpha +1\right)}+\frac{{\left(0.06t\right)}^{2\alpha }}{\text{Γ}\left(2\alpha +1\right)}\frac{{\left(0.06t\right)}^{2\alpha }}{\text{Γ}\left(2\alpha +1\right)}\max \left(x-25{\text{e}}^{-0.06},0\right)\end{array}$

$⋮$

$\begin{array}{l}{p}^n:{\psi }_n\left(x,t\right)={\mathbb{S}}^{-1}\left[u^{\alpha }\mathbb{S}\left[{\Psi }_{n-1}\left(x,t\right)\right]\right]\\ =x\frac{{\left(-0.06t\right)}^{n\alpha }}{\text{Γ}\left(n\alpha +1\right)}+\frac{{\left(0.06t\right)}^{n\alpha }}{\text{Γ}\left(n\alpha +1\right)}\frac{{\left(0.06t\right)}^{n\alpha }}{\text{Γ}\left(n\alpha +1\right)}\max \left(x-25{\text{e}}^{-0.06},0\right)\end{array}$ (87)

So that the solution $\psi \left(x,t\right)$ of the problem is given by

$\begin{array}{c}\psi \left(x,t\right)=\underset{p\to 1}{\lim }\underset{n=0}{\overset{\infty }{{\displaystyle \sum }}}{p}^n{\psi }_n\left(x,t\right)\\ =\max \left(x-25{\text{e}}^{-0.06},0\right)E_{\alpha }\left(0.06t^{\alpha }\right)+x\left(1-E_{\alpha }\left(0.06t^{\alpha }\right)\right)\end{array}$ (88)

This is the exact solution of the given option pricing equation Equation (72). The exact solution at the special case $\alpha =1$ is

$\psi \left(x,t\right)=\max \left(x-25{\text{e}}^{-0.06},0\right){\text{e}}^{0.06t}+x\left(1-{\text{e}}^{0.06t}\right)$ (89)

7. Conclusion

It is widely known that various scientific models in the life sciences end up with a PDEs or FPDEs. Solving these equations is very important to understand many phenomena in life sciences. In this work, two methods are applied: MHPM and HPSTM, in order to find the solution of the Fractional Black-Scholes model. The both methods are based on the HPM. However, they are effectively employed for getting the solution. To conclude, MPHM and HPSTM are exceptionally effective and productive methods to discover the approached solution, also the numerical solution for the Time-Fractional Partial Differential Equations.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Hilfer R. (2000) Applications of Fractional Calculus in Physics: World Scientific. https://doi.org/10.1142/3779
[2]	Magin, R. and Ovadia, M. (2006) Modeling the Cardiac Tissue Electrode Interface Using Fractional Calculus. IFAC Proceedings, 39, 302-307. https://doi.org/10.3182/20060719-3-PT-4902.00056
[3]	Gepreel, K.A. (2011) The Homotopy Perturbation Method Applied to the Nonlinear Fractional Kolmogorov-Petrovskii-Piskunov Equations. Applied Mathematics Letters, 24, 1428-1434. https://doi.org/10.1016/j.aml.2011.03.025
[4]	Odibat, Z. and Momani, S. (2007) A Reliable Treatment of Homotopy Perturbation Method for Klein-Gordon Equations. Physics Letters A, 365, 351-357. https://doi.org/10.1016/j.physleta.2007.01.064
[5]	Podlubny, I. (1998) Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications. Academic Press.
[6]	Rafei, M. and Ganji, D. (2006) Explicit Solutions of Helmholtz Equation and Fifth-Order KdV Equation Using Homotopy Perturbation Method. International Journal of Nonlinear Sciences and Numerical Simulation, 7, 321-328. https://doi.org/10.1515/IJNSNS.2006.7.3.321
[7]	Cheng, J.-F. and Chu, Y.-M. (2011) Solution to the Linear Fractional Differential Equation Using Adomian Decomposition Method. Mathematical Problems in Engineering, 2011, Article ID: 587068. https://doi.org/10.1155/2011/587068
[8]	He, J.-H. (1998) Approximate Analytical Solution for Seepage Flow with Fractional Derivatives in Porous Media. Computer Methods in Applied Mechanics and Engineering, 167, 57-68. https://doi.org/10.1016/S0045-7825(98)00108-X
[9]	He, J.-H. (1999) Homotopy Perturbation Technique. Computer Methods in Applied Mechanics and Engineering, 178, 257-262. https://doi.org/10.1016/S0045-7825(99)00018-3
[10]	He, J.-H. (2003) Homotopy Perturbation Method: A New Nonlinear Analytical Technique. Applied Mathematics and Computation, 135, 73-79. https://doi.org/10.1016/S0096-3003(01)00312-5
[11]	He, J.-H. (2005) Limit Cycle and Bifurcation of Nonlinear Problems. Chaos, Solitons & Fractals, 26, 827-833. https://doi.org/10.1016/j.chaos.2005.03.007
[12]	He, J.H. (2006) New Interpretation of Homotopy Perturbation Method. International Journal of Modern Physics B, 20, 25661-25668. https://doi.org/10.1142/S0217979206034819
[13]	Miller, K.S. and Ross, B. (1993) An Introduction to the Fractional Calculus and Fractional Differential Equations.
[14]	Bhalekar, S. and Daftardar-Gejji, V. (2008) New Iterative Method: Application to Partial Differential Equations. Applied Mathematics and Computation, 203, 778-783. https://doi.org/10.1016/j.amc.2008.05.071
[15]	He, J.-H. (2004) Solution of Nonlinear Equations by an Ancient Chinese Algorithm. Applied Mathematics and Computation, 151, 293-297. https://doi.org/10.1016/S0096-3003(03)00348-5
[16]	Daftardar-Gejji, V. and Bhalekar, S. (2010) Solving Fractional Boundary Value Problems with Dirichlet Boundary Conditions using a New Iterative Method. Computers & Mathematics with Applications, 59, 1801-1809. https://doi.org/10.1016/j.camwa.2009.08.018
[17]	Adomian, G. (1994) Solving Frontier Problems of Physics: The Decomposition Method. Kluwer, Boston.
[18]	Y?ld?r?m, A. (2010) Application of the Homotopy Perturbation Method for the Fokker-Planck Equation. International Journal for Numerical Methods in Biomedical Engineering, 26, 1144-1154. https://doi.org/10.1002/cnm.1200
[19]	Y?ld?r?m, A. (2010) Analytical Approach to Fokker-Planck Equation with Space-and Time-Fractional Derivatives by Means of the Homotopy Perturbation Method. Journal of King Saud University—Science, 22, 257-264. https://doi.org/10.1016/j.jksus.2010.05.008
[20]	Arafa, A., Rida, S. and Mohamed, H. (2011) Homotopy Analysis Method for Solving Biological Population Model. Communications in Theoretical Physics, 56, 797-800. https://doi.org/10.1088/0253-6102/56/5/01
[21]	Kilbas, A.A., Saigo, M. and Saxena, R. (2004) Generalized Mittag-Leffler Function and Generalized Fractional Calculus Operators. Integral Transforms and Special Functions, 15, 31-49. https://doi.org/10.1080/10652460310001600717
[22]	Belgacem, F.B.M. and Karaballi, A.A. (2006) Sumudu Transform Fundamental Properties Investigations and Applications. International Journal of Stochastic Analysis, 2006, Article ID: 91083. https://doi.org/10.1155/JAMSA/2006/91083