# A Novel Effective Approach for Solving Fractional Nonlinear PDEs

Authors: Hossein Aminikhah; Nasrin Malekzadeh; Hadi Rezazadeh
Journal: International Scholarly Research Notices
 (2014)

---

## Abstract

The present work introduces an effective modification of homotopy perturbation method for the solution of nonlinear time-fractional biological population model and a system of three nonlinear time-fractional partial differential equations. In this approach, the solution is considered a series expansion that converges to the nonlinear problem. The new approximate analytical procedure depends only on two iteratives. The analytical approximations to the solution are reliable and confirm the ability of the new homotopy perturbation method as an easy device for computing the solution of nonlinear equations.


---

## Body

## 1. Introduction

In recent years, fractional differential equations have played an important role in different research areas such as mechanics, electricity, biology, economics, notably control theory, and signal and image processing [1–5]. Many phenomena can be described very successfully by models using mathematical tools from fractional calculus, that is, the theory of derivatives and integrals of fractional order. The most important advantage of using fractional differential equations in these and other applications is their nonlocal property. It is well known that the integer order differential operator is a local operator, but the fractional-order differential operator is nonlocal. This means that the next state of a system depends not only upon its current state but also upon all of its historical states. In fact, this is the main reason why differential operators of fractional order provide an excellent instrument for description of memory and hereditary properties of various mathematical, physical, and engineering processes. The reader is asked to refer to [6–9] in order to know more details about the fractional DEs, including their history and kinds, their existence and uniqueness of solutions, and their applications and methods of solutions.During the last few years, the numerical methods and exact solution methods have been proposed to solve fractional differential equations, for example, the Adomian decomposition method [10], the homotopy perturbation method [11, 12], the variational iteration method [13, 14], the differential transform method [15, 16], the G
′
/
G method [17, 18], the first integral method [19], and the exp-function method [20]. Aminikhah and Hemmatnezhad [21] and Aminikhah and Biazar [22] proposed a new homotopy perturbation method to obtain the exact and numerical solutions of ordinary differential equation.In this paper, we will implement the new homotopy perturbation method to obtain the approximate solution for the following time-fractional derivative nonlinear partial differential equation.(i)
The fractional-order biological population equation as a nonlinear model is as follows [23, 24]:
(1)
D
t
α
u
=
(
u
2
)
x
x
+
(
u
2
)
y
y
+
f
(
u
)
,
t
>
0
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiii
x
,
y
∈
R
,
0
<
α
≤
1
f
(
u
)
=
h
u
a
(
1
-
r
u
b
)
,
u
(
x
,
y
,
0
)
=
g
(
x
,
y
)
,where u denotes population density, f represents the population supply due to births and deaths, h, a, r, and b are real numbers, g is given initial condition, and D
t
α denotes the differential operator in the sense of Caputo.(ii)
The system of three nonlinear time-fractional partial equations which have been widely discussed in the literature [25, 26] is as follows:
(2)
D
t
α
u
=
-
v
x
w
y
+
v
y
w
x
-
u
,
D
t
α
v
=
-
w
x
u
y
-
w
y
u
x
+
v
,
D
t
α
w
=
-
u
x
v
y
-
u
y
v
x
+
w
,
iiiiiiiiiiiiii
t
>
0
,
0
<
α
≤
1
,with initial conditions
(3)
u
(
x
,
y
,
0
)
=
h
1
(
x
,
y
)
,
v
(
x
,
y
,
0
)
=
h
2
(
x
,
y
)
,
w
(
x
,
y
,
0
)
=
h
3
(
x
,
y
)
,
where D
t
α denotes the differential operator in the sense of Caputo.The rest of this paper is organized as follows. In Section2, the properties of fractional derivative are described. In Section 3, basic ideas of new homotopy perturbation method for solving fractional partial differential equations are presented. In Section 4, the application of new homotopy perturbation method for solving the nonlinear fractional biological population equation and the system of three nonlinear time-fractional partial equations are presented. Finally, we give a brief conclusion in the last section.


## 2. Background on Fractional Derivatives

The history of fractional calculus is more than three centuries old, yet only in the past 20 years has the field received much attention and interest. The reader may refer to [8, 9]. There are several definitions for fractional differential equations. These definitions include Grunwald-Letnikov, Riemann-Liouville, Caputo, Weyl, Marchaud, and Riesz fractional derivatives; among those Riemann-Liouville and Caputo fractional derivatives are the most popular. The differential equations defined in terms of Riemann-Liouville derivatives require fractional initial conditions whereas the differential equations defined in terms of Caputo derivatives require boundary conditions involving integer order derivatives, which have clear physical interpretations. For this reason, Caputo fractional derivatives are popular among scientists and engineers.Now, we give some basic definitions and properties of the fractional calculus theory used in this work.Definition 1.
The Riemann-Liouville fractional integral operatorJ
μ of order μ on the usual Lebesgue space L
1
[
a
,
b
] is given by
(4)
J
μ
f
(
x
)
=
1
Γ
(
μ
)
∫
0
x
(
x
-
t
)
μ
-
1
f
(
t
)
d
t
,
μ
>
0
,
J
0
f
(
x
)
=
f
(
x
)
.
It has the following properties:(i)
J
μ exists for any x
∈
[
a
,
b
],
(ii)
J
μ
J
β
=
J
μ
+
β,
(iii)
J
μ
J
β
=
J
β
J
μ,
(iv)
J
μ
J
β
f
(
x
)
=
J
β
J
μ
f
(
x
),
(v)
J
μ
(
x
-
a
)
γ
=
(
Γ
(
γ
+
1
)
/
Γ
(
μ
+
γ
+
1
)
)
(
x
-
a
)
μ
+
γ,wheref
∈
L
1
[
a
,
b
], μ
,
β
≥
0, and γ
>
-
1.Definition 2.
The Caputo definition of fractal derivative operator is given by(5)
D
μ
f
(
x
)
=
J
m
-
μ
D
n
f
(
x
)
=
1
Γ
(
m
-
μ
)
∫
0
t
(
x
-
τ
)
m
-
μ
-
1
f
m
(
τ
)
d
τ
,
where m
-
1
<
μ
≤
m,  m
∈
N, x
>
0, and f
∈
L
1
[
a
,
b
].Remark 3.
Forn to be the smallest integer that exceeds α, the Caputo time-fractional derivative operator of α
>
0 is defined as
(6)
D
t
α
u
(
x
,
y
,
t
)
=
∂
α
u
(
x
,
y
,
t
)
∂
t
α
=
{
1
Γ
(
n
-
α
)
∫
0
t
(
t
-
τ
)
n
-
α
-
1
∂
n
u
(
x
,
y
,
τ
)
∂
t
n
d
τ
,
n
-
1

<

α
<

n
∂
n
u
(
x
,
y
,
t
)
∂
t
n
,
α
=
n
∈
N
.


## 3. Analysis of New Homotopy Perturbation Method

In this section, we construct the solution of system of partial differential equations with time-fractional derivative by extending the idea of [21, 22]. Let us consider the system of nonlinear fractional differential equations:
(7)
D
t
α
u
i
(
x
,
t
)
=
A
i
(
u
i
)
+
f
i
(
t
,
x
)
,
x
,
t
∈
Ω
,
i
=
1,2
,
…
,
n
,
with the following initial conditions:
(8)
u
i
(
x
,
0
)
=
α
i
,
i
=
1,2
,
…
,
n
,
where A
i are the operators, f
i are known functions, and u
i are sought functions. Assume that operators A
i can be written as A
i
(
u
i
)
=
L
i
(
u
i
)
+
N
i
(
u
i
), where L
i are the linear operators and N
i are the nonlinear operators. Hence, (7) can be rewritten as follows:
(9)
D
t
α
u
i
(
x
,
t
)
=
L
i
(
u
i
)
+
N
i
(
u
i
)
+
f
i
(
x
,
t
)
.
For solving system (7) by NHPM [22] we construct the following homotopy:
(10)
H
(
U
i
;
p
)
=
(
1
-
p
)
(
D
t
α
U
i
(
x
,
t
)
-
u
i
,
0
)
+
p
(
D
t
α
u
i
(
x
,
t
)
-
L
i
(
U
i
)
-
N
i
(
U
i
)
-
f
i
(
t
,
x
)
)
=
0
,
where p
∈
[
0,1
]is an embedding or homotopy parameter, H
(
t
,
x
;
p
)
:
Ω
×
[
0,1
]
→
R, and u
i
,
0 are the initial approximation of solution of the problem in (9).Clearly, the homotopy equationsH
(
U
i
;
0
)
=
0 and H
i
(
U
i
;
1
)
=
0 are equivalent to the equations D
t
α
U
i
(
x
,
t
)
-
u
i
,
0
=
0 and D
t
α
U
i
(
x
,
t
)
-
L
i
(
U
i
)
+
N
i
(
U
i
)
+
f
i
(
t
,
x
)
=
0, respectively. Thus, a monotonous change of parameter p from zero to one corresponds to a continuous change of the trivial problem D
t
α
U
i
(
x
,
t
)
-
u
i
,
0
=
0 to the original problem. Next, we assume that the solution of equation H
(
U
i
,
p
) can be written as a power series in embedding parameter p, as follows:
(11)
U
i
=
U
i
,
0
+
p
U
i
,
1
,
i
=
1,2
,
…
,
n
.
Now, let us write (11) in the following form:
(12)
D
t
α
U
i
(
x
,
t
)
=
u
i
,
0
+
p
(
L
i
(
U
i
)
+
N
i
(
U
i
)
+
f
i
(
t
,
x
)
)
.
Applying the inverse operator J
t
α which is the Riemann-Liouville fractional integral of order α
≥
0, on both sides of (12), we have
(13)
U
i
(
x
,
t
)
=
U
i
(
x
,
0
)
+
J
t
α
u
i
,
0
+
p
J
t
α
(
L
i
(
U
i
)
+
N
i
(
U
i
)
+
f
i
(
t
,
x
)
)
.Suppose that the initial approximation of (9) has the form
(14)
u
i
,
0
(
x
,
t
)
=
∑
n
=
0
∞
a
i
,
n
(
x
)
p
n
(
t
)
,
i
=
1,2
,
…
,
n
,
where a
i
,
n
(
x
)
,
n
=
0,1
,
2
,
…, are unknown coefficients and p
n
(
t
)
,
n
=
0,1
,
2
,
…, are specific functions on the problem. By substituting (11) and (14) into (13), we get
(15)
U
i
,
0
+
p
U
i
,
1
=
U
i
(
x
,
0
)
+
J
t
α
(
∑
n
=
0
∞
a
i
,
n
(
x
)
p
n
(
t
)
)
+
p
J
t
α
(
L
i
(
U
i
,
0
+
p
U
i
,
1
)
+
N
i
(
U
i
,
0
+
p
U
i
,
1
)
+
f
i
(
t
,
x
)
)
.Equating the coefficients of like powers ofp, we get the following set of equations:
(16)
coefficient
of
p
0
:
U
i
,
0
(
x
,
t
)
=
U
i
(
x
,
0
)
+
J
t
α
(
∑
n
=
0
∞
a
i
,
n
(
x
)
p
n
(
t
)
)
,
coefficient
of
p
1
:
U
i
,
1
(
x
,
t
)
=
J
t
α
(
L
i
(
U
i
,
0
)
+
N
i
(
U
i
,
0
)
+
f
i
(
t
,
x
)
)
.Now, we solve these equations in such a way thatU
i
,
1
(
x
,
t
)
=
0. Therefore, the approximate solution may be obtained as
(17)
u
i
(
x
,
t
)
=
U
i
,
0
(
x
,
t
)
=
U
i
(
x
,
0
)
+
J
t
α
(
∑
n
=
0
∞
a
i
,
n
(
x
)
p
n
(
t
)
)
.


## 4. Applications

In this section, we present some examples with analytical solution to show the efficiency of methods described in the previous section for solving (1) and (2).Example 1.
Let us consider (1) with a
=
1, r
=
0, and g
(
x
,
y
)
=
x
y, then we have
(18)
D
t
α
u
(
x
,
y
,
t
)
=
∂
2
(
u
2
(
x
,
y
,
t
)
)
∂
x
2
+
∂
2
(
u
2
(
x
,
y
,
t
)
)
∂
y
2
+
h
u
(
x
,
y
,
t
)
,
u
(
x
,
y
,
0
)
=
x
y
.
The exact solution of (18), for the special case α
=
1, is u
(
x
,
y
,
t
)
=
x
y
e
h
t.
To obtain the solution of (18) by new homotopy perturbation method, we construct the following homotopy:
(19)
(
1
-
p
)
(
D
t
α
U
(
x
,
y
,
t
)
-
u
0
(
x
,
y
,
t
)
)
+
p
(
D
t
α
U
(
x
,
y
,
t
)
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
-
h
U
(
x
,
y
,
t
)
)
=
0
.
Applying the inverse operator J
t
α of D
t
α on both sides of (19), we obtain
(20)
U
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
-
p
J
t
α
(
u
0
(
x
,
y
,
t
)
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
iiiiiiiiiiiiiiiii
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
-
h
U
(
x
,
y
,
t
)
)
.
The solution of (21) has the following form:
(21)
U
(
x
,
y
,
t
)
=
U
0
(
x
,
y
,
t
)
+
p
U
1
(
x
,
y
,
t
)
.
Substituting (21) into (20) and equating the coefficients of like powers of p, we get the following set of equations:
(22)
U
0
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
,
U
1
(
x
,
y
,
t
)
=
J
t
α
(
-
u
0
(
x
,
y
,
t
)
+
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
+
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
+
h
U
(
x
,
y
,
t
)
)
.
Assumingu
0
(
x
,
y
,
t
)
=
∑
n
=
0
∞
a
n
(
x
,
y
)
p
n
(
t
), p
n
(
t
)
=
t
n
α, and U
(
x
,
y
,
0
)
=
u
(
x
,
y
,
0
) and solving (22) for U
1
(
x
,
y
,
t
) lead to the following result:
(23)
U
1
(
x
,
y
,
t
)
=
(
-
a
0
(
x
,
y
)
+
h
x
y
)
t
α
Γ
(
α
+
1
)
+
(
2
(
∂
a
0
(
x
,
y
)
/
∂
y
)
x
x
y
-
1
2
a
0
(
x
,
y
)
y
2
(
x
y
)
3
/
2
+
2
(
∂
a
0
(
x
,
y
)
/
∂
x
)
y
x
y
+
h
a
0
(
x
,
y
)
+
2
x
y
∂
2
a
0
(
x
,
y
)
∂
y
2
-
a
1
(
x
,
y
)
-
1
2
a
0
(
x
,
y
)
x
2
(
x
y
)
3
/
2
+
2
x
y
∂
2
a
0
(
x
,
y
)
∂
x
2
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
(
(
∂
a
1
(
x
,
y
)
/
∂
y
)
x
x
y
-
1
4
a
1
(
x
,
y
)
y
2
(
x
y
)
3
/
2
1
4
a
1
(
x
,
y
)
x
2
(
x
y
)
3
/
2
iiiiiiiiiiiiiii
+
2
(
∂
a
0
(
x
,
y
)
∂
y
)
2
+
2
a
0
(
x
,
y
)
∂
2
a
0
(
x
,
y
)
∂
y
2
iiiiiiiiiiiiiii
+
x
y
∂
2
a
1
(
x
,
y
)
∂
y
2
+
2
(
∂
a
0
(
x
,
y
)
∂
x
)
2
iiiiiiiiiiiiiii
+
1
2
h
a
1
(
x
,
y
)
+
2
a
0
(
x
,
y
)
∂
2
a
0
(
x
,
y
)
∂
x
2
iiiiiiiiiiiiiii
+
(
∂
a
1
(
x
,
y
)
/
∂
x
)
y
x
y
+
x
y
∂
2
a
1
(
x
,
y
)
∂
x
2
iiiiiiiiiiiiiii
-
a
2
(
x
,
y
)
-
1
4
a
1
(
x
,
y
)
x
2
(
x
y
)
3
/
2
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
(
a
0
(
x
,
y
)
∂
2
a
1
(
x
,
y
)
∂
y
2
iiiiiiiiiiiiiiii
+
2
3
x
y
∂
2
a
2
(
x
,
y
)
∂
y
2
⋯
-
a
3
(
x
,
y
)
)
×
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
.
Vanishing U
1
(
x
,
y
,
t
) lets the coefficients a
i
,
i
=
0,1
,
2
,
…, have the following values:
(24)
a
0
(
x
,
y
)
=
h
x
y
,
a
1
(
x
,
y
)
=
h
2
x
y
,
a
2
(
x
,
y
)
=
1
2
!
h
3
x
y
,
a
3
(
x
,
y
)
=
1
3
!
h
4
x
y
,
a
4
(
x
,
y
)
=
1
4
!
h
5
x
y
,
…
.
Therefore, we obtain the solutions of (18) as
(25)
u
(
x
,
y
,
t
)
=
x
y
+
h
x
y
t
α
Γ
(
α
+
1
)
+
h
2
x
y
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
1
2
!
h
3
x
y
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
1
3
!
h
4
x
y
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
=
x
y
(
1
+
∑
n
=
1
∞
Γ
(
(
n
-
1
)
α
+
1
)
(
h
t
α
)
n
(
n
-
1
)
!
Γ
(
n
α
+
1
)
)
.
If α
=
1, we have
(26)
u
(
x
,
y
,
t
)
=
x
y
(
1
+
h
t
+
h
2
t
2
2
!
+
h
3
t
3
3
!
+
⋯
)
=
x
y
e
h
t
,
which is an exact solution.Example 2.
Consider (1) with a
=
1
,
b
=
1, and g
(
x
,
y
)
=
e
h
r
/
8
(
x
+
y
), then we have
(27)
D
t
α
u
(
x
,
y
,
t
)
=
∂
2
(
u
2
(
x
,
y
,
t
)
)
∂
x
2
+
∂
2
(
u
2
(
x
,
y
,
t
)
)
∂
y
2
+
h
u
(
x
,
y
,
t
)
(
1
-
r
u
(
x
,
y
,
t
)
)
;
subject to the initial condition
(28)
u
(
x
,
y
,
0
)
=
e
h
r
/
8
(
x
+
y
)
.
The exact solution of (27), for the special case α
=
1, is
(29)
u
(
x
,
y
,
t
)
=
e
h
r
/
8
(
x
+
y
)
+
h
t
.
To obtain the solution of (27) by the new homotopy perturbation method, we construct the following homotopy:
(30)
(
1
-
p
)
(
D
t
α
U
(
x
,
y
,
t
)
-
u
0
(
x
,
y
,
t
)
)
+
p
(
D
t
α
U
(
x
,
y
,
t
)
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
iiiiiiiiiiiiiiiiii
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
iiiiiiiiiiiiiiiiii
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
h
U
(
x
,
y
,
t
)
(
1
-
r
U
(
x
,
y
,
t
)
)
)
=
0
.
Applying the inverse operator J
t
α of D
t
α on both sides of (30), we obtain
(31)
U
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
-
p
J
t
α
(
u
0
(
x
,
y
,
t
)
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
iiiiiiiiiiiiiiiii
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
-
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
iiiiiiiiiiiiiiiii
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
-
h
U
(
x
,
y
,
t
)
(
1
-
r
U
(
x
,
y
,
t
)
)
)
.
The solution of (31) has the following form:
(32)
U
(
x
,
y
,
t
)
=
U
0
(
x
,
y
,
t
)
+
p
U
1
(
x
,
y
,
t
)
.
Substituting (32) into (33) and equating the coefficients of like powers of p, we get the following set of equations:
(33)
U
0
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
,
U
1
(
x
,
y
,
t
)
=
J
t
α
(
-
u
0
(
x
,
y
,
t
)
+
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
x
2
iiiiiiiiiiii
+
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
iiiiiiiiiiii
∂
2
(
U
2
(
x
,
y
,
t
)
)
∂
y
2
+
h
U
(
x
,
y
,
t
)
(
1
-
r
U
(
x
,
y
,
t
)
)
)
.
Assuming u
0
(
x
,
y
,
t
)
=
∑
n
=
0
10
a
n
(
x
,
y
)
p
n
(
t
), p
n
(
t
)
=
t
n
α, U
(
x
,
y
,
0
)
=
u
(
x
,
y
,
0
), and solving (33) for U
1
(
x
,
y
,
t
) lead to the following result:
(34)
U
1
(
x
,
y
,
t
)
=
(
-
a
0
(
x
,
y
)
+
(
e
h
r
/
8
(
x
+
y
)
)
2
h
r
+
h
(
e
h
r
/
8
(
x
+
y
)
-
r
(
e
h
r
/
8
(
x
+
y
)
)
2
)
)
t
α
Γ
(
α
+
1
)
+
(
2
e
h
r
/
8
(
x
+
y
)
∂
2
a
0
(
x
,
y
)
∂
x
2
+
2
h
r
e
h
r
/
8
(
x
+
y
)
∂
a
0
(
x
,
y
)
∂
x
+
h
(
-
2
e
h
r
/
8
(
x
+
y
)
r
a
0
(
x
,
y
)
+
a
0
(
x
,
y
)
)
+
2
h
r
e
h
r
/
8
(
x
+
y
)
∂
a
0
(
x
,
y
)
∂
y
+
1
2
h
r
e
h
r
/
8
(
x
+
y
)
a
0
(
x
,
y
)
-
a
1
(
x
,
y
)
+
2
e
h
r
/
8
(
x
+
y
)
∂
2
a
0
(
x
,
y
)
∂
y
2
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
(
e
h
r
/
8
(
x
+
y
)
∂
2
a
1
(
x
,
y
)
∂
x
2
-
a
2
(
x
,
y
)
+
2
(
∂
a
0
(
x
,
y
)
∂
y
)
2
+
h
(
-
e
h
r
/
8
(
x
+
y
)
r
a
1
(
x
,
y
)
-
(
a
0
(
x
,
y
)
)
2
r
+
1
2
iiiiiiiiiiiiiiiiiiiii
-
(
a
0
(
x
,
y
)
)
2
r
+
1
2
a
1
(
x
,
y
)
)
+
1
4
h
r
e
h
r
/
8
(
x
+
y
)
a
1
(
x
,
y
)
+
2
a
0
(
x
,
y
)
∂
2
a
0
(
x
,
y
)
∂
y
2
+
e
h
r
/
8
(
x
+
y
)
∂
2
a
1
(
x
,
y
)
∂
y
2
+
2
(
∂
a
0
(
x
,
y
)
∂
x
)
2
+
2
a
0
(
x
,
y
)
∂
2
a
0
(
x
,
y
)
∂
x
2
+
1
2
2
h
r
e
h
r
/
8
(
x
+
y
)
∂
a
1
(
x
,
y
)
∂
y
+
1
2
2
h
r
e
h
r
/
8
(
x
+
y
)
∂
a
1
(
x
,
y
)
∂
x
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
(
1
3
2
h
r
e
h
r
/
8
(
x
+
y
)
∂
a
2
(
x
,
y
)
∂
y
-
a
3
(
x
,
y
)
+

·
·
·
+
2
3
e
h
r
/
8
(
x
+
y
)
∂
2
a
2
(
x
,
y
)
∂
y
2
)
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)

+

⋯
.
Vanishing U
1
(
x
,
y
,
t
) lets the coefficients a
i
,
i
=
0,1
,
2
,
…, have the following values:
(35)
a
0
(
x
,
y
)
=
h
e
h
r
/
8
(
x
+
y
)
,
a
1
(
x
,
y
)
=
h
2
e
h
r
/
8
(
x
+
y
)
,
a
2
(
x
,
y
)
=
1
2
!
h
3
e
h
r
/
8
(
x
+
y
)
,
a
3
(
x
,
y
)
=
1
3
!
h
4
e
h
r
/
8
(
x
+
y
)
,
a
4
(
x
,
y
)
=
1
4
!
h
5
e
h
r
/
8
(
x
+
y
)
,
…
.
Therefore, we obtain the solution of (27) as
(36)
u
(
x
,
y
,
t
)
=
e
h
r
/
8
(
x
+
y
)
+
h
e
h
r
/
8
(
x
+
y
)
t
α
Γ
(
α
+
1
)
+
h
2
e
h
r
/
8
(
x
+
y
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
1
2
!
h
3
e
h
r
/
8
(
x
+
y
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
1
3
!
h
4
e
h
r
/
8
(
x
+
y
)
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
-
·
·
·
=
e
h
r
/
8
(
x
+
y
)
(
1
+
∑
n
=
1
∞
Γ
(
(
n
-
1
)
α
+
1
)
(
h
t
α
)
n
(
n
-
1
)
!
Γ
(
n
α
+
1
)
)
.
If we put α
→
1 in (36) or solve (30) with α
=
1, we obtain the exact solution
(37)
u
(
x
,
y
,
t
)
=
e
h
r
/
8
(
x
+
y
)
(
1
+
h
t
+
h
2
t
2
2
!
+
h
3
t
3
3
!
+
⋯
)
=
e
h
r
/
8
(
x
+
y
)
+
h
t
.Example 3.
Let us consider (2) with h
1
(
x
,
y
)
=
e
(
x
+
y
), h
2
(
x
,
y
)
=
e
(
x
-
y
), and h
3
(
x
,
y
)
=
e
(
-
x
+
y
). Therefore, we have
(38)
D
t
α
u
(
x
,
y
,
t
)
=
-
∂
v
(
x
,
y
,
t
)
∂
x
∂
w
(
x
,
y
,
t
)
∂
y
+
∂
v
(
x
,
y
,
t
)
∂
y
∂
w
(
x
,
y
,
t
)
∂
x
-
u
(
x
,
y
,
t
)
,
D
t
α
v
(
x
,
y
,
t
)
=
-
∂
w
(
x
,
y
,
t
)
∂
x
∂
u
(
x
,
y
,
t
)
∂
y
-
∂
w
(
x
,
y
,
t
)
∂
y
∂
u
(
x
,
y
,
t
)
∂
x
+
v
(
x
,
y
,
t
)
,
D
t
α
w
(
x
,
y
,
t
)
=
-
∂
u
(
x
,
y
,
t
)
∂
x
∂
v
(
x
,
y
,
t
)
∂
y
-
∂
u
(
x
,
y
,
t
)
∂
y
∂
v
(
x
,
y
,
t
)
∂
x
+
w
(
x
,
y
,
t
)
subject to the initial conditions
(39)
u
(
x
,
y
,
0
)
=
e
x
+
y
,
v
(
x
,
y
,
0
)
=
e
x
-
y
,
w
(
x
,
y
,
0
)
=
e
-
x
+
y
.
The exact solution of (38), for the special case α
=
1, is
(40)
u
(
x
,
y
,
t
)
=
e
x
+
y
-
t
,
v
(
x
,
y
,
t
)
=
e
x
-
y
+
t
,
w
(
x
,
y
,
t
)
=
e
-
x
+
y
+
t
.
To obtain the solution of (38) by new homotopy perturbation method, we construct the following homotopy:
(41)
(
1
-
p
)
(
D
t
α
U
(
x
,
y
,
t
)
-
u
0
(
x
,
y
,
t
)
)
+
p
(
D
t
α
U
(
x
,
y
,
t
)
+
∂
V
(
x
,
y
,
t
)
∂
x
∂
W
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiii
-
∂
V
(
x
,
y
,
t
)
∂
y
∂
W
(
x
,
y
,
t
)
∂
x
+
U
(
x
,
y
,
t
)
)
=
0
,
(
1
-
p
)
(
D
t
α
V
(
x
,
y
,
t
)
-
v
0
(
x
,
y
,
t
)
)
+
p
(
D
t
α
V
(
x
,
y
,
t
)
+
∂
W
(
x
,
y
,
t
)
∂
x
∂
U
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiii
+
∂
W
(
x
,
y
,
t
)
∂
y
∂
U
(
x
,
y
,
t
)
∂
x
-
V
(
x
,
y
,
t
)
)
=
0
,
(
1
-
p
)
(
D
t
α
W
(
x
,
y
,
t
)
-
w
0
(
x
,
y
,
t
)
)
+
p
(
D
t
α
W
(
x
,
y
,
t
)
+
∂
U
(
x
,
y
,
t
)
∂
x
∂
V
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiii
+
∂
U
(
x
,
y
,
t
)
∂
y
∂
V
(
x
,
y
,
t
)
∂
x
-
W
(
x
,
y
,
t
)
)
=
0
.
Applying the inverse operator J
t
α of D
t
α on both sides of the above equation, we obtain
(42)
U
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
-
p
J
t
α
(
u
0
(
x
,
y
,
t
)
+
∂
V
(
x
,
y
,
t
)
∂
x
∂
W
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiiiii
+
∂
V
(
x
,
y
,
t
)
∂
x
∂
W
(
x
,
y
,
t
)
∂
y
-
∂
V
(
x
,
y
,
t
)
∂
y
∂
W
(
x
,
y
,
t
)
∂
x
+
U
(
x
,
y
,
t
)
)
,
V
(
x
,
y
,
t
)
=
V
(
x
,
y
,
0
)
+
J
t
α
v
0
(
x
,
y
,
t
)
-
p
J
t
α
(
v
0
(
x
,
y
,
t
)
+
∂
W
(
x
,
y
,
t
)
∂
x
∂
U
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiiiiiii
+
∂
W
(
x
,
y
,
t
)
∂
y
∂
U
(
x
,
y
,
t
)
∂
x
-
V
(
x
,
y
,
t
)
)
,
W
(
x
,
y
,
t
)
=
W
(
x
,
y
,
0
)
+
J
t
α
w
0
(
x
,
y
,
t
)
-
p
J
t
α
(
w
0
(
x
,
y
,
t
)
+
∂
U
(
x
,
y
,
t
)
∂
x
∂
V
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiiiiiii
+
∂
U
(
x
,
y
,
t
)
∂
y
∂
V
(
x
,
y
,
t
)
∂
x
-
W
(
x
,
y
,
t
)
)
.
The solution of (42) has the following form:
(43)
U
(
x
,
y
,
t
)
=
U
0
(
x
,
y
,
t
)
+
p
U
1
(
x
,
y
,
t
)
,
V
(
x
,
y
,
t
)
=
V
0
(
x
,
y
,
t
)
+
p
V
1
(
x
,
y
,
t
)
,
W
(
x
,
y
,
t
)
=
W
0
(
x
,
y
,
t
)
+
p
W
1
(
x
,
y
,
t
)
.
Substituting (43) into (42) and equating the coefficients of like powers of p, we get the following set of equations:
(44)
U
0
(
x
,
y
,
t
)
=
U
(
x
,
y
,
0
)
+
J
t
α
u
0
(
x
,
y
,
t
)
,
V
0
(
x
,
y
,
t
)
=
V
(
x
,
y
,
0
)
+
J
t
α
v
0
(
x
,
y
,
t
)
,
W
0
(
x
,
y
,
t
)
=
W
(
x
,
y
,
0
)
+
J
t
α
w
0
(
x
,
y
,
t
)
,
U
1
(
x
,
y
,
t
)
=
J
t
α
(
-
u
0
(
x
,
y
,
t
)
-
∂
V
(
x
,
y
,
t
)
∂
x
∂
W
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiii
+
∂
V
(
x
,
y
,
t
)
∂
y
∂
W
(
x
,
y
,
t
)
∂
x
-
U
(
x
,
y
,
t
)
)
,
V
1
(
x
,
y
,
t
)
=
J
t
α
(
-
v
0
(
x
,
y
,
t
)
-
∂
W
(
x
,
y
,
t
)
∂
x
∂
U
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiii
-
∂
W
(
x
,
y
,
t
)
∂
y
∂
U
(
x
,
y
,
t
)
∂
x
+
V
(
x
,
y
,
t
)
)
,
W
1
(
x
,
y
,
t
)
=
J
t
α
(
-
w
0
(
x
,
y
,
t
)
-
∂
U
(
x
,
y
,
t
)
∂
x
∂
V
(
x
,
y
,
t
)
∂
y
iiiiiiiiiiiiiiiii
∂
U
(
x
,
y
,
t
)
∂
x
-
∂
U
(
x
,
y
,
t
)
∂
y
∂
V
(
x
,
y
,
t
)
∂
x
+
W
(
x
,
y
,
t
)
)
.
Assuming(45)
u
0
(
x
,
y
,
t
)
=
∑
n
=
0
∞
a
n
(
x
,
y
)
p
n
(
t
)
,
v
0
(
x
,
y
,
t
)
=
∑
n
=
0
∞
b
n
(
x
,
y
)
p
n
(
t
)
,
w
0
(
x
,
y
,
t
)
=
∑
n
=
0
∞
c
n
(
x
,
y
)
p
n
(
t
)
,
p
n
(
t
)
=
t
n
α
,
U
(
x
,
y
,
0
)
=
u
(
x
,
y
,
0
)
,
V
(
x
,
y
,
0
)
=
v
(
x
,
y
,
0
)
,
W
(
x
,
y
,
0
)
=
w
(
x
,
y
,
0
)
and solving the above equation for U
1
(
x
,
y
,
t
) and V
1
(
x
,
y
,
t
) and W
1
(
x
,
y
,
t
) lead to the following result:
(46)
U
1
(
x
,
y
,
t
)
=
(
-
a
0
(
x
,
y
)
-
e
x
+
y
)
t
α
Γ
(
α
+
1
)
+
(
-
e
x
-
y
∂
c
0
(
x
,
y
)
∂
y
-
∂
b
0
(
x
,
y
)
∂
x
e
-
x
+
y
-
a
1
(
x
,
y
)
iiiiiiIiiiiiii
-
a
0
(
x
,
y
)
-
∂
b
0
(
x
,
y
)
∂
y
e
-
x
+
y
iiiiiiiiiiiiiii
-
e
x
-
y
∂
c
0
(
x
,
y
)
∂
x
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
(
-
1
2
e
x
-
y
∂
c
1
(
x
,
y
)
∂
y
-
1
2
∂
b
1
(
x
,
y
)
∂
y
e
-
x
+
y
iiiiiiiiiiiiiii
-
1
2
e
x
-
y
∂
c
1
(
x
,
y
)
∂
x
-
1
2
∂
b
1
(
x
,
y
)
∂
x
e
-
x
+
y
iiiiiiiiiiiiiii
-
1
2
a
1
(
x
,
y
)
-
∂
b
0
(
x
,
y
)
∂
x
∂
c
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
+
∂
b
0
(
x
,
y
)
∂
y
∂
c
0
(
x
,
y
)
∂
x
-
a
2
(
x
,
y
)
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
(
-
1
3
∂
b
2
(
x
,
y
)
∂
x
e
-
x
+
y
+
1
2
∂
b
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
×
∂
c
1
(
x
,
y
)
∂
x
+
⋯
-
a
3
(
x
,
y
)
)
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
,
V
1
(
x
,
y
,
t
)
=
(
-
b
0
(
x
,
y
)
+
e
x
-
y
)
t
α
Γ
(
α
+
1
)
+
(
-
e
-
x
+
y
∂
a
0
(
x
,
y
)
∂
x
-
∂
c
0
(
x
,
y
)
∂
y
e
x
+
y
iiiiiiiiiiiiiii
-
∂
c
0
(
x
,
y
)
∂
x
e
x
+
y
+
b
0
(
x
,
y
)
-
b
1
(
x
,
y
)
iiiiiiiiiiiiiii
+
∂
a
0
(
x
,
y
)
∂
y
e
-
x
+
y
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
(
1
2
e
-
x
+
y
∂
a
1
(
x
,
y
)
∂
y
-
1
2
∂
a
1
(
x
,
y
)
∂
x
e
-
x
+
y
iiiiiiiiiiiiiii
-
1
2
e
x
+
y
∂
c
1
(
x
,
y
)
∂
x
-
1
2
∂
c
1
(
x
,
y
)
∂
y
e
x
+
y
iiiiiiiiiiiiiii
-
1
2
b
1
(
x
,
y
)
-
∂
c
0
(
x
,
y
)
∂
x
∂
a
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
-
∂
c
0
(
x
,
y
)
∂
y
∂
a
0
(
x
,
y
)
∂
x
-
b
2
(
x
,
y
)
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
(
-
1
3
∂
a
2
(
x
,
y
)
∂
x
e
-
x
+
y
-
1
2
∂
c
1
(
x
,
y
)
∂
x
∂
a
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
1
2
∂
c
1
(
x
,
y
)
∂
x
+
⋯
-
b
3
(
x
,
y
)
)
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
,
W
1
(
x
,
y
,
t
)
=
(
-
c
0
(
x
,
y
)
+
e
-
x
+
y
)
t
α
Γ
(
α
+
1
)
+
(
e
x
-
y
∂
a
0
(
x
,
y
)
∂
x
-
∂
a
0
(
x
,
y
)
∂
y
e
x
-
y
iiiiiiiiiiiiiii
-
∂
b
0
(
x
,
y
)
∂
x
e
x
+
y
+
c
0
(
x
,
y
)
-
c
1
(
x
,
y
)
iiiiiiiiiiiiiii
-
∂
b
0
(
x
,
y
)
∂
y
e
x
+
y
)
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
(
1
2
e
x
-
y
∂
a
1
(
x
,
y
)
∂
x
-
1
2
∂
b
1
(
x
,
y
)
∂
x
e
x
+
y
iiiiiiiiiiiiiii
-
1
2
e
x
+
y
∂
b
1
(
x
,
y
)
∂
y
-
1
2
∂
a
1
(
x
,
y
)
∂
y
e
x
-
y
iiiiiiiiiiiiiii
+
1
2
c
1
(
x
,
y
)
-
∂
a
0
(
x
,
y
)
∂
x
∂
b
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
-
∂
a
0
(
x
,
y
)
∂
y
∂
b
0
(
x
,
y
)
∂
x
-
c
2
(
x
,
y
)
)
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
(
-
1
3
∂
a
2
(
x
,
y
)
∂
x
e
x
-
y
-
1
2
∂
b
1
(
x
,
y
)
∂
x
∂
a
0
(
x
,
y
)
∂
y
iiiiiiiiiiiiiii
-
1
3
∂
a
2
(
x
,
y
)
∂
x
+
⋯
-
c
3
(
x
,
y
)
)
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
.
Vanishing U
1
(
x
,
y
,
t
), V
1
(
x
,
y
,
t
), and W
1
(
x
,
y
,
t
) lets the coefficients a
i
,
b
i
,
c
i
,
i
=
0,1
,
2
,
…, have the following values:
(47)
a
0
(
x
,
y
)
=
-
e
x
+
y
,
a
1
(
x
,
y
)
=
e
x
+
y
,
a
2
(
x
,
y
)
=
-
1
2
!
e
x
+
y
,
a
3
(
x
,
y
)
=
1
3
!
e
x
+
y
,
a
4
(
x
,
y
)
=
-
1
4
!
e
x
+
y
,
…
,
b
0
(
x
,
y
)
=
e
x
-
y
,
b
1
(
x
,
y
)
=
e
x
-
y
,
b
2
(
x
,
y
)
=
1
2
!
e
x
-
y
,
b
3
(
x
,
y
)
=
1
3
!
e
x
-
y
,
b
4
(
x
,
y
)
=
1
4
!
e
x
-
y
,
…
c
0
(
x
,
y
)
=
e
-
x
+
y
,
c
1
(
x
,
y
)
=
e
-
x
+
y
,
c
2
(
x
,
y
)
=
1
2
!
e
-
x
+
y
,
c
3
(
x
,
y
)
=
1
3
!
e
-
x
+
y
,
c
4
(
x
,
y
)
=
1
4
!
e
-
x
+
y
,
…
.
Therefore, we obtain the solutions of (38) as
(48)
u
(
x
,
y
,
t
)
=
e
x
+
y
-
e
x
+
y
t
α
Γ
(
α
+
1
)
+
e
x
+
y
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
-
1
2
!
e
x
+
y
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
1
3
!
e
x
+
y
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
-
⋯
=
e
x
+
y
(
1
+
∑
n
=
1
∞
(
-
1
)
n
Γ
(
(
n
-
1
)
α
+
1
)
t
n
α
(
n
-
1
)
!
Γ
(
n
α
+
1
)
)
,
v
(
x
,
y
,
t
)
=
e
x
-
y
+
e
x
-
y
t
α
Γ
(
α
+
1
)
+
e
x
-
y
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
1
2
!
e
x
-
y
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
1
3
!
e
x
-
y
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
=
e
x
-
y
(
1
+
∑
n
=
1
∞
Γ
(
(
n
-
1
)
α
+
1
)
t
n
α
(
n
-
1
)
!
Γ
(
n
α
+
1
)
)
,
w
(
x
,
y
,
t
)
=
e
-
x
+
y
+
e
-
x
+
y
t
α
Γ
(
α
+
1
)
+
e
-
x
+
y
Γ
(
α
+
1
)
t
2
α
Γ
(
2
α
+
1
)
+
1
2
!
e
-
x
+
y
Γ
(
2
α
+
1
)
t
3
α
Γ
(
3
α
+
1
)
+
1
3
!
e
-
x
+
y
Γ
(
3
α
+
1
)
t
4
α
Γ
(
4
α
+
1
)
+
⋯
=
e
-
x
+
y
(
1
+
∑
n
=
1
∞
Γ
(
(
n
-
1
)
α
+
1
)
t
n
α
(
n
-
1
)
!
Γ
(
n
α
+
1
)
)
.
If we solve (38) with α
=
1, we obtain the exact solution
(49)
u
(
x
,
y
,
t
)
=
e
x
+
y
(
1
-
t
+
t
2
2
!
-
t
3
3
!
+
⋯
)
=
e
x
+
y
-
t
,
v
(
x
,
y
,
t
)
=
e
x
+
y
(
1
+
t
+
t
2
2
!
+
t
3
3
!
+
⋯
)
=
e
x
-
y
+
t
,
w
(
x
,
y
,
t
)
=
e
-
x
+
y
(
1
+
t
+
t
2
2
!
+
t
3
3
!
+
⋯
)
=
e
-
x
+
y
+
t
.


## 5. Conclusion

In the present work, we proposed a new homotopy perturbation method to solve nonlinear time-fractional biological population model and a system of three nonlinear time-fractional partial differential equations. The new method for solving fractional-order differential equations is based on two component procedures and polynomial initial condition. The new homotopy perturbation method is very simple in application and is less computational and more accurate in comparison with other mentioned methods. Unlike the Adomian decomposition method, the new homotopy perturbation method is free from the need to use Adomian polynomials. In this method, we do not need the Lagrange multiplier, correction functional, stationary conditions, and calculating integrals, which eliminate the complications that exist in the variational iteration method. The obtained results show that these approaches can solve the problem effectively and can be applied to many nonlinear differential equations. The computations associated with the examples discussed above were performed by MAPLE.



---

*Source: 647492-2014-10-29.xml*