Numerical Methods for Fractional Order Singular Partial Differential Equations with Variable Coefficients

Abstract

We implement relatively analytical methods, the homotopy perturbation method and the variational iteration method, for solving singular fractional partial differential equations of fractional order. The process of the methods which produce solutions in terms of convergent series is explained. The fractional derivatives are described in Caputo sense. Some examples are given to show the accurate and easily implemented of these methods even with the presence of singularities.

1. Introduction

In the last decades, fractional calculus found many applications in various fields of engineering and physical sciences such as physics, chemistry, biology, economy, viscoelasticity, electrochemistry, electromagnetic, relaxation processes, diffusion, control, porous media, and many more; see, for example, [1–6].

The numerical solution of singular differential equations of integer order has been a hot topic in numerical and computational mathematics for a long time [7, 8]. Singular partial differential equations of fractional order, as generalizations of classical singular partial differential equations of integer order, are increasingly used to model problems in physics and engineering. Consequently, considerable attention has been given to the solution of singular partial differential equations of fractional order.

In most of these equations analytical solutions are either quite difficult or impossible to achieve, so approximations and numerical techniques must be used. Several methods have been used to solve singular differential equations such as variational iteration method [8], homotopy perturbation method [7], and homotopy analysis method [9].

The homotopy perturbation method [10–12] and variational iteration method [13–16] are relatively new approach to provide an analytical approximation to linear and nonlinear problems, and they are particularly valuable as tools for scientists and applied mathematicians, because it provides immediate and visible symbolic terms of analytic solutions, as well as numerical approximate solutions to both linear and nonlinear differential equations.

Recently, the application of these methods extended for fractional differential equations. He [15] was the first to apply the variational iteration method to fractional differential equations.

The objective of the present paper is to extend the application of homotopy perturbation method to provide approximate solutions and to make comparison with that obtained by the variational iteration method for singular partial differential equations of fractional order: 𝜕^𝛼𝑢𝜕𝑡^𝛼(𝑥)𝜕+𝜇⁴𝑢𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢𝜕𝑧⁴=0,𝑎<𝑥,𝑦,𝑧<𝑏,𝑡>0, (1) where 𝜇(𝑥), 𝜆(𝑦), ℎ(𝑧) subject to initial conditions 𝑢𝑥,𝑦,𝑧,0=𝑓₀𝑥,𝑦,𝑧,𝜕𝑢𝜕𝑡𝑥,𝑦,𝑧,0=𝑓₁𝑥,𝑦,𝑧 (2) and boundary conditions 𝑢𝑎,𝑦,𝑧,𝑡=𝑔₀,𝑢𝑦,𝑧,𝑡𝑏,𝑦,𝑧,𝑡=𝑔₁,𝑢()𝑦,𝑧,𝑡𝑥,𝑎,𝑧,𝑡=𝑔₂(),𝑢()𝑥,𝑧,𝑡𝑥,𝑏,𝑧,𝑡=𝑔₃(),𝑢𝑥,𝑧,𝑡𝑥,𝑦,𝑎,𝑡=𝑔₄(),𝑢𝑥,𝑧,𝑡𝑥,𝑦,𝑏,𝑡=𝑔₅(),𝜕𝑥,𝑧,𝑡₂𝑢𝜕𝑥²𝑎,𝑦,𝑧,𝑡=𝑘₀,𝜕𝑦,𝑧,𝑡₂𝑢𝜕𝑥²𝑏,𝑦,𝑧,𝑡=𝑘₀,𝜕𝑦,𝑧,𝑡₂𝑢𝜕𝑦²(𝑥,𝑎,𝑧,𝑡)=𝑘₂(),𝜕𝑥,𝑧,𝑡₂𝑢𝜕𝑦²(𝑥,𝑎,𝑧,𝑡)=𝑘₃(),𝜕𝑥,𝑧,𝑡₂𝑢𝜕𝑧²𝑥,𝑦,𝑎,𝑡=𝑘₄,𝜕𝑥,𝑦,𝑡₂𝑢𝜕𝑧²𝑥,𝑦,𝑏,𝑡=𝑘₅𝑦,𝑧,𝑡, (3) where 𝜕^𝛼(⋅)/𝜕𝑡^𝛼 is the fractional derivative in the Caputo sense, 1<𝛼≤2, and 𝑔_𝑖 and 𝑘_𝑖, 𝑖=0,…,5 are continuous. The paper is organized as follows. In Section 2, some basic definitions and properties of fractional calculus theory are given. In Section 3, the basic idea of HPM is given. In Section 4, the basic idea of variational iteration method is given. In Sections 5 and 6 the analysis of HPM and VIM for singular fractional partial differential equations, respectively, is given. Some examples are given in Section 7. Concluding remarks are listed in Section 8.

2. Preliminaries

In this section, we give some basic definitions and properties of fractional calculus theory which is used in this paper.

Definition 1. A real function 𝑓(𝑥), 𝑥>0 is said to be in space 𝐶𝜇, 𝜇∈𝑅 if there exists a real number 𝑝>𝜇, such that 𝑓(𝑥)=𝑥^𝑝𝑓₁(𝑥) where 𝑓₁(𝑥)∈𝐶(0,∞), and it is said to be in the space 𝐶^𝑛_𝜇 if 𝑓^𝑛∈𝑅_𝜇, 𝑛∈𝑁.

Definition 2. The Riemann-Liouville fractional integral operator of order 𝛼≥0 of a function 𝑓∈𝐶𝜇, 𝜇≥−1 is defined as 𝐽^𝛼𝑓(𝑥)=1(Γ(𝛼)^𝑥₀(𝑥−𝑡)^𝛼−1𝑓(𝑡)𝑑𝑡,𝛼>0,𝑡>0; (4) in particular 𝐽⁰𝑓(𝑥)=𝑓(𝑥).

For 𝛽≥0 and 𝛾≥−1, some properties of the operator 𝐽^𝛼(i)𝐽^𝛼𝐽^𝛽𝑓(𝑥)=𝐽^𝛼+𝛽𝑓(𝑥),(ii)𝐽^𝛼𝐽^𝛽𝑓(𝑥)=𝐽^𝛽𝐽^𝛼𝑓(𝑥),(iii)𝐽^𝛼𝑥^𝛾=(Γ(𝛾+1)/Γ(𝛼+𝛾+1))𝑥^𝛼+𝛾.

Definition 3. The Caputo fractional derivative of 𝑓∈𝐶^𝑚₋₁, 𝑚∈𝑁 is defined as 𝐷^𝛼𝑓(𝑥)=1Γ(𝑚−𝛼)^𝑥₀(𝑥−𝑡)^{𝑚−𝛼−1}𝑓^𝑚(𝑡)𝑑𝑡,𝑚−1<𝛼≤𝑚. (5)

Lemma 4. if 𝑚−1<𝛼≤𝑚, 𝑚∈𝑁, 𝑓∈𝐶^𝑚_𝜇, 𝜇>−1 then the following two properties hold: (i)𝐷^𝛼[𝐽^𝛼𝑓(𝑥)]=𝑓(𝑥), (ii)𝐽^𝛼[𝐷^𝛼𝑓(𝑥)]=𝑓(𝑥)−∑^𝑚−1_𝑘=1𝑓^𝑘(0)(𝑥^𝑘/𝑘!).

3. Homotopy Perturbation Method

To illustrate the basic idea of this method, we consider the following nonlinear differential equation: 𝐴(𝑢)−𝑓(𝑟)=0,𝑟∈Ω, (6) with boundary conditions 𝐵𝑢,𝜕𝑢𝜕𝑛=0,𝑟∈Γ, (7) where 𝐴 is a general differential operator, 𝐵 is a boundary operator, 𝑓(𝑟) is a known analytic function, and Γ is the boundary of the domain Ω.

In general, the operator 𝐴 can be divided into two parts 𝐿 and 𝑁, where 𝐿 is linear, while 𝑁 is nonlinear. Equation (6) therefore can be rewritten as follows: 𝐿(𝑢)+𝑁(𝑢)−𝑓(𝑟)=0. (8) By the homotopy technique [10, 11], we construct a homotopy 𝑣(𝑟,𝑝)∶Ω×[0,1]→𝑅 which satisfies 𝐻𝑣,𝑝=1−𝑝𝐿(𝑣)−𝐿𝑢₀𝐴(𝑣)(𝑟[]+𝑝−𝑓)=0𝑝∈0,1,𝑟∈Ω (9) or 𝐻𝑣,𝑝=𝐿(𝑣)−𝐿𝑢₀+𝑝𝐿𝑢₀+𝑝𝑁(𝑣)−𝑓(𝑟)=0, (10) where 𝑝∈[0,1] is an embedding parameter and 𝑢₀ is an initial approximation of (6) which satisfies the boundary conditions.

From (9) and (10) we have 𝐻(𝑣,0)=𝐿(𝑣)−𝐿𝑢₀𝐻()(𝑣)(𝑟)=0,𝑣,1=𝐴−𝑓=0. (11) The change in the process of 𝑝 from zero to unity is just that of 𝑣(𝑟,𝑝) from 𝑢₀(𝑟) to 𝑢(𝑟). In topology this is called deformation and 𝐿(𝑣)−𝐿(𝑢₀) and 𝐴(𝑣)−𝑓(𝑟) are called homotopic.

Now, assume that the solution of (9) and (10) can be expressed as 𝑣=𝑣₀+𝑝𝑣₁+𝑝²𝑣₂+⋯. (12) Setting 𝑝=1 results in the approximate solution of (6). Consider 𝑢=lim𝑣_𝑝→1=𝑣₀+𝑣₁+𝑣₂+⋯. (13)

4. The Variational Iteration Method

To illustrate the basic concepts of VIM, we consider the following differential equation: 𝐿(𝑢)+𝑁(𝑢)=𝑔(𝑥), (14) where 𝐿 is a linear operator, 𝑁 is a nonlinear operator, and 𝑔(𝑥) is a nonhomogeneous term.

According to VIM, one constructs a correction functional as follows: 𝑦_𝑛+1=𝑦_𝑛+^𝑥₀𝜆𝐿𝑦_𝑛(𝑠)−𝑁𝑦_𝑛(𝑠)𝑑𝑠, (15) where 𝜆 is a general Lagrange multiplier, and 𝑦_𝑛 denotes restricted variation; that is, 𝛿𝑦_𝑛=0.

5. Analysis of Homotopy Perturbation Method

To illustrate the basic concepts of HPM for (1) with (2) and (3), we use the view of He [10, 11], where the following homotopy was constructed for (1): 𝜕1−𝑝^𝛼𝑢𝜕𝑡^𝛼𝜕+𝑝^𝛼𝑢𝜕𝑡^𝛼(𝑥)𝜕+𝜇⁴𝑢𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢𝜕𝑧⁴=0, (16) or 𝜕^𝛼𝑢𝜕𝑡^𝛼𝜇(𝑥)𝜕+𝑝⁴𝑢𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢𝜕𝑧⁴=0, (17) where 𝑝∈[0,1] is an embedding parameter. If 𝑝=0, (17) becomes linear fractional differential equation 𝜕^𝛼𝑢𝜕𝑡^𝛼=0, (18) and when 𝑝=1, (17) turns out to be the original equation. In view of basic assumption of HPM, solution of (1) can be expressed as a power series in 𝑝𝑢(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑝₁𝑢₁(𝑥,𝑡)+𝑝₂𝑢₂(𝑥,𝑡)+⋯. (19) When 𝑝=1, we get the approximate solution of (1) 𝑢(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+𝑢₂(𝑥,𝑡)+⋯. (20) The convergence of series (20) has been proved in [12].

Substituting (19) into (17) and then equating the terms with identical power of 𝑝, we obtain the following series of linear equations: 𝑝⁰∶𝜕^𝛼𝑢₀𝜕𝑡^𝛼𝑝=0₁∶𝜕^𝛼𝑢₁𝜕𝑡^𝛼(𝑥)𝜕=𝜇⁴𝑢₀𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢₀𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢₀𝜕𝑧⁴𝑝²∶𝜕^𝛼𝑢₂𝜕𝑡^𝛼(𝑥)𝜕=𝜇⁴𝑢₁𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢₁𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢₁𝜕𝑧⁴𝑝³∶𝜕^𝛼𝑢₃𝜕𝑡^𝛼(𝑥)𝜕=𝜇⁴𝑢₂𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢₂𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢₂𝜕𝑧⁴⋮⋮ (21) with the initial conditions (2) and boundary conditions (3).

The initial approximation can be chosen in the following manner. 𝑢₀=¹_𝑗=0𝛾_𝑗𝑡^𝑗𝑗!=𝛾₀+𝛾₁𝑡,where𝛾₀=𝑓₀𝑥,𝑦,𝑧,𝛾₁=𝑓₁𝑥,𝑦,𝑧. (22) The equations in the system (21) can be solved by applying the operator 𝐽^𝛼_∗ and by simple computation; we approximate the series solution of HPM by the following N-term truncated series: 𝑢(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+⋯+𝑢_𝑛(𝑥,𝑡), (23) which is the approximate solution of (1) with (2) and (3).

6. Analysis of Variational Iteration Method

To solve the fractional singular partial differential equations by using the variational iteration method, with initial and boundary conditions (2) and (3), we construct the following correction functional: 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘(𝑥,𝑡)+𝐽^𝛼_𝑡𝜆𝜕^𝛼𝑢_𝑘(𝑥,𝑡)𝜕𝑡^𝛼+𝐑𝑢_𝑘(𝑥,𝑡), (24) or 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘()+1𝑥,𝑡Γ(𝛼)𝛿×^𝑡₀(𝑡−𝑠)^𝛼−1𝜆(𝑠𝜕)^𝛼𝑢_𝑘(𝑥,𝑠)𝜕𝑠^𝛼+𝐑𝑢_𝑘(𝑥,𝑠)𝑑𝑠, (25) where 𝐑𝑢_𝑘()(𝑥)𝜕𝑥,𝑠=𝜇⁴𝑢_𝑘(𝑥,𝑠)𝜕𝑥⁴𝑦𝜕+𝜆⁴𝑢_𝑘(𝑥,𝑠)𝜕𝑦⁴(𝑧)𝜕+ℎ⁴𝑢_𝑘(𝑥,𝑠)𝜕𝑧⁴, (26) and 𝐽^𝛼_𝑡 is the Riemann-Liouville fractional integral operator of order 𝛼, with respect to variable 𝑡, 𝜆 is a general Lagrange multiplier which can be identified optimally by variational theory [17], and 𝑢_𝑘(𝑥,𝑡) are considered as restricted variation; that is, 𝛿𝑢_𝑘(𝑥,𝑡)=0.

Making the above correction functional stationary, the following condition can be obtained: 𝛿𝑢_𝑘+1(𝑥,𝑡)=𝛿𝑢_𝑘()+1𝑥,𝑡Γ(𝛼)𝛿×^𝑡₀(𝑡−𝑠)^𝛼−1𝜆(𝑠)𝜕^𝛼𝑢_𝑘(𝑥,𝑠)𝜕𝑠^𝛼+𝐑𝑢_𝑘(𝑥,𝑠)𝑑𝑠, (27) which yields to Lagrange multiplier 𝜆(𝑠)=𝑠−𝑡. (28) We obtain the following iteration formula by substitution of (28) in (25): 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘()+1𝑥,𝑡Γ()×𝛼−1_𝑡0(𝑡−𝑠)^𝛼−2(𝜕𝑡−𝑠)^𝛼𝑢_𝑘(𝑥,𝑠)𝜕𝑠^𝛼+𝑅𝑢_𝑘(𝑥,𝑠)𝑑𝑠. (29) That is, 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘()−()𝑥,𝑡𝛼−1Γ(𝛼)×^𝑡₀(𝑡−𝑠)^𝛼−1𝜕^𝛼𝑢_𝑘(𝑥,𝑠)𝜕𝑠^𝛼+𝑅𝑢_𝑘(𝑥,𝑠)𝑑𝑠. (30) This yields the following iteration formula: 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢_𝑘(𝑥,𝑡)𝜕𝑡^𝛼+𝑅𝑢_𝑘(𝑥,𝑡). (31) The initial approximation 𝑢₀ can be chosen by the following manner which satisfies initial conditions: 𝑢₀=¹_𝑗=0𝛾_𝑗𝑡^𝑗𝑗!=𝛾₀+𝛾₁𝑡where𝛾₀=𝑓₀𝑥,𝑦,𝑧,𝛾₁=𝑓₁𝑥,𝑦,𝑧. (32) We can obtain the following first-order approximation by substitution of (32) in (31): 𝑢₁(𝑥,𝑡)=𝑢₀(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡𝜕^𝛼𝑢₀(𝑥,𝑡)𝜕𝑡^𝛼+𝑅𝑢₀(𝑥,𝑡). (33) Finally, by substituting the constant values of 𝛾₀ and 𝛾₁ in (33), we have the results as the approximate solutions of (1) with (2) and (3).

7. Applications

In this section we have applied homotopy perturbation method and variational iteration method to fractional singular partial differential equations with known exact solution.

Example 5. Consider the following fourth-order fractional singular partial differential equation: 𝜕^𝛼𝑢𝜕𝑡^𝛼+𝑥𝜕sin𝑥−1⁴𝑢𝜕𝑥⁴=0,0<𝑥<1,𝑡>0,1<𝛼≤2. (34) With initial conditions 𝑢(𝑥,0)=𝑥−sin𝑥,0<𝑥<1𝜕𝑢()()𝜕𝑡𝑥,0=−𝑥−sin𝑥,0<𝑥<1 (35) and with boundary conditions 𝑢(0,𝑡)=0,𝑢(1,𝑡)=𝑒^−𝑡()𝜕1−sin1,𝑡>0,₂𝑢𝜕𝑥²()𝜕0,𝑡=0,²𝑢𝜕𝑥²(1,𝑡)=𝑒^−𝑡sin1,𝑡>1. (36) The exact solution in special case 𝛼=2 is 𝑢(𝑥,𝑡)=(𝑥−sin𝑥)𝑒^−𝑡. (37) First, we construct 𝜕^𝛼𝑢𝜕𝑡^𝛼𝑥=−𝑝𝜕sin𝑥−1⁴𝑢𝜕𝑥⁴. (38) Substituting (19) in (38) and then equating the terms with same powers of 𝑝, we get the series 𝑝⁰∶𝜕^𝛼𝑢₀𝜕𝑡^𝛼=0,𝑢₀(𝑥,𝑡)=𝑥−sin𝑥,𝜕𝑢₀()(),𝑝𝜕𝑡𝑥,0=−𝑥−sin𝑥₁∶𝜕^𝛼𝑢₁𝜕𝑡^𝛼𝑥=−𝜕sin𝑥−1⁴𝑢₀𝜕𝑥⁴,𝑝²∶𝜕^𝛼𝑢₂𝜕𝑡^𝛼𝑥=−𝜕sin𝑥−1⁴𝑢₁𝜕𝑥⁴,⋮ (39) Now applying the operator 𝐽^𝛼_𝑡 to (39) and using initial and boundary conditions yield 𝑢₀()=()−()𝑢𝑥,𝑡𝑥−sin𝑥𝑥−sin𝑥𝑡,₁()=()𝑡𝑥,𝑡𝑥−sin𝑥^𝛼Γ()−()𝑡𝛼+1𝑥−sin𝑥^𝛼+1Γ(),𝑢𝛼+2₂()=()𝑡𝑥,𝑡𝑥−sin𝑥^2𝛼Γ()−()𝑡2𝛼+1𝑥−sin𝑥^2𝛼+1Γ(),⋮2𝛼+2 (40) so that the solution 𝑢(𝑥,𝑡) of the problem (34) is given by 𝑢(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+𝑢₂()𝑢()=()−()()×𝑡𝑥,𝑡+⋯𝑥,𝑡𝑥−sin𝑥𝑥−sin𝑥𝑡+𝑥−sin𝑥_𝛼Γ()−()𝑡𝛼+1𝑥−sin𝑥^𝛼+1Γ()+()𝑡𝛼+2𝑥−sin𝑥^2𝛼Γ()−()𝑡2𝛼+1𝑥−sin𝑥^2𝛼+1Γ(2𝛼+2)+⋯. (41) Now, we solve by variational iteration method. According to variational iteration method, the formula (31) for (34) can be expressed in the following form: 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢_𝑘(𝑥,𝑡)𝜕𝑡^𝛼+𝑥𝜕sin𝑥−1⁴𝑢𝜕𝑥⁴. (42) Suppose that an initial approximation has the following form which satisfies the initial conditions: 𝑢₀(𝑥,𝑡)=(𝑥−sin𝑥)−(𝑥−sin𝑥)𝑡. (43) Now by iteration formula (42), we obtain the first approximation: 𝑢₁(𝑥,𝑡)=𝑢₀(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢₀(𝑥,𝑡)𝜕𝑡^𝛼+𝑥𝜕sin𝑥−1⁴𝑢₀𝜕𝑥⁴=()−()()×𝑡𝑥−sin𝑥𝑥−sin𝑥𝑡+𝛼−1)(𝑥−sin𝑥_𝛼Γ()−()𝑡𝛼+1𝛼−1)(𝑥−sin𝑥^𝛼+1Γ(𝛼+2) (44) and second approximation: 𝑢₂(𝑥,𝑡)=𝑢₁(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢₁(𝑥,𝑡)𝜕𝑡^𝛼+𝑥𝜕sin𝑥−1⁴𝑢₁𝜕𝑥⁴=()−()𝑡−()𝑡𝑥−sin𝑥𝑥−sin𝑥𝛼−3)(𝛼−1)(𝑥−sin𝑥^𝛼Γ()+()𝑡𝛼+1𝛼−3)(𝛼−1)(𝑥−sin𝑥^𝛼+1Γ()+()𝛼+2𝛼−1²()𝑡𝑥−sin𝑥^2𝛼Γ()−()2𝛼+1𝛼−1²()𝑡𝑥−sin𝑥^2𝛼+1Γ().2𝛼+2 (45)
Table 1 shows the approximate solution for (34) obtained for different values of 𝛼 using the homotopy perturbation method and the variational iteration method. The value 𝛼=2 is the only case for which we know the exact solution 𝑢(𝑥,𝑡)=(𝑥−sin𝑥)𝑒^−𝑡 and our approximate solution using the variational iteration method is more accurate than approximate solution obtained using homotopy perturbation method. It is to be noted that only the third-order term series was used in evaluating the approximate solutions for Table 1.

Table 1 

Numerical values when 𝛼=1.5, 1.75, and 2.0 for (34).

Example 6. Consider the following singular two-dimensional partial differential equation of fractional order: 𝜕^𝛼𝑢𝜕𝑡^𝛼1+2𝑥²+𝑥⁴𝜕6!⁴𝑢𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢𝜕𝑦⁴,0<𝑥,𝑦<1,𝑡>0,1<𝛼≤2. (46) With initial conditions 𝑢𝑥,𝑦,0=0,0<𝑥<1𝜕𝑢𝑥𝜕𝑡𝑥,𝑦,0=2+⁶+𝑦6!6!⁶,0<𝑥<1 (47) and with boundary conditions 𝑢=()0.5,𝑦,𝑡2+0.5⁶+𝑦6!⁶𝑢=16!sin𝑡,1,𝑦,𝑡2++𝑦6!⁶𝜕6!sin𝑡,𝑡>0,₂𝑢𝜕𝑥²=()0.5,𝑦,𝑡0.5⁴𝜕6!sin𝑡,₂𝑢𝜕𝑥²=11,𝑦,𝑡𝜕6!sin𝑡,𝑡>1,₂𝑢𝜕𝑦²()=()𝑥,0.5,𝑡0.5⁴𝜕6!sin𝑡,₂𝑢𝜕𝑦²()=1𝑥,1,𝑡6!sin𝑡,𝑡>1, (48) and the exact solution in special case 𝛼=2 is 𝑢=𝑥𝑥,𝑦,𝑡2+⁶+𝑦6!6!⁶sin𝑡. (49) According to HPM, we construct the following homotopy: 𝜕^𝛼𝑢𝜕𝑡^𝛼21=𝑝𝑥²+𝑥⁴𝜕6!⁴𝑢𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢𝜕𝑦⁴. (50) Substituting (19) in (50) and then equating the terms with same powers of 𝑝, we get the series 𝑝⁰∶𝜕^𝛼𝑢₀𝜕𝑡^𝛼𝑢=0,₀(𝑥,𝑡)=0,𝜕𝑢₀()𝑥𝜕𝑡𝑥,0=2+⁶+𝑦6!⁶𝑝6!₁∶𝜕^𝛼𝑢₁𝜕𝑡^𝛼21=−𝑥²+𝑥⁴𝜕6!⁴𝑢₀𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢₀𝜕𝑦⁴,𝑝²∶𝜕^𝛼𝑢₂𝜕𝑡^𝛼21=−𝑥²+𝑥⁴𝜕6!⁴𝑢₁𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢₁𝜕𝑦⁴,⋮ (51) Now applying the operator 𝐽^𝛼_𝑡 to (51) and using initial and boundary conditions yield 𝑢₀()=𝑥𝑥,𝑡2+⁶+𝑦6!⁶𝑢6!𝑡,₁()𝑥𝑥,𝑡=−2+⁶+𝑦6!⁶𝑡6!Γ(^𝛼+1),𝑢𝛼+2₂()=𝑥𝑥,𝑡2+⁶+𝑦6!⁶𝑡6!Γ(^2𝛼+1),⋮2𝛼+2 (52) so that the solution 𝑢(𝑥,𝑡) of the problem (46) is given by 𝑢(𝑥,𝑡)=𝑢₀(𝑥,𝑡)+𝑢₁(𝑥,𝑡)+𝑢₂()𝑢()=𝑥𝑥,𝑡+⋯𝑥,𝑡2+⁶+𝑦6!⁶𝑥6!𝑡−2+⁶+𝑦6!⁶𝑡6!Γ(^𝛼+1)+𝑥𝛼+22+⁶+𝑦6!⁶𝑡6!Γ(^2𝛼+1)=𝑥2𝛼+2+⋯2+⁶+𝑦6!⁶𝑡6!𝑡−^𝛼+1Γ()+𝑡𝛼+2Γ(^2𝛼+1).2𝛼+2+⋯ (53) Now, we solve the problem (46) by variational iteration method. According to variational iteration method, the formula (31) for (46) can be expressed in the following form: 𝑢_𝑘+1(𝑥,𝑡)=𝑢_𝑘(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢_𝑘(𝑥,𝑡)𝜕𝑡^𝛼1+2𝑥²+𝑥⁴𝜕6!⁴𝑢𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢𝜕𝑦⁴. (54) Suppose that an initial approximation has the following form which satisfies the initial conditions 𝑢₀()=𝑥𝑥,𝑡2+⁶+𝑦6!6!⁶𝑡. (55) Now by iteration formula (42), we obtain the following approximations: 𝑢₁(𝑥,𝑡)=𝑢₀(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢₀(𝑥,𝑡)𝜕𝑡^𝛼1+2𝑥²+𝑥⁴𝜕6!⁴𝑢₀𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢₀𝜕𝑦⁴=𝑥2+⁶+𝑦6!⁶()𝑥6!𝑡−𝛼−12+⁶+𝑦6!⁶𝑡6!Γ(^𝛼+1).𝛼+2 (56) The second approximation takes the following form: 𝑢₂(𝑥,𝑡)=𝑢₁(𝑥,𝑡)−(𝛼−1)𝐽^𝛼_𝑡×𝜕^𝛼𝑢₁(𝑥,𝑡)𝜕𝑡^𝛼1+2𝑥²+𝑥⁴𝜕6!⁴𝑢₁𝜕𝑥⁴1+2𝑦²+𝑦⁴𝜕6!⁴𝑢₁𝜕𝑦⁴=𝑥2+⁶+𝑦6!⁶()𝑥6!𝑡−2𝛼−12+⁶+𝑦6!⁶×𝑡6!_𝛼+1Γ()+()𝛼+2𝛼−1²𝑥2+⁶+𝑦6!⁶𝑡6!Γ(^2𝛼+1)=𝑥2𝛼+12+⁶+𝑦6!⁶×()𝑡6!𝑡−2𝛼−1^𝛼+1Γ()+()𝛼+2𝛼−1²𝑡^2𝛼+1Γ().2𝛼+1 (57)
Table 2 shows the approximate solution for (46) obtained for different values of 𝛼 using the homotopy perturbation method and the variational iteration method. The value 𝛼=2 is the only case for which we know the exact solution 𝑢(𝑥,𝑡)=(2+𝑥⁶/6!+𝑦⁶/6!)sin𝑡, and our approximate solution using homotopy perturbation method is more accurate than approximate solution obtained using the variational iteration method. It is to be noted that only the third-order term series was used in evaluating the approximate solutions for Table 2.

Table 2 

Numerical values when 𝛼=1.5, 1.75, and 2.0 for (46).

8. Conclusion

The essential goal of this work has been to construct an approximate solution of singular partial differential equations of fractional order. The goal has been done by using the homotopy perturbation method and the variational iteration method. The methods provide the solutions in terms of convergent series with easily computable component even with the presence of singularities. Two examples were presented in Section 7 to show the accuracy of these methods; in Example 5 it seems that the approximate solution using the variational iteration method converges faster than approximate solution using homotopy perturbation method while the approximate solution in Example 6 using homotopy perturbation method converges faster than the approximate solution obtained using the variational iteration method to the exact solution. The fact that the proposed methods solve nonlinear problems without using Adomian’s polynomial can be considered as advantage of these techniques over the decomposition method.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors express their sincere thanks to the referees for the careful and noteworthy reading of the paper and for their very helpful suggestions that improved the paper substantially. The authors also gratefully acknowledge that this work is partially supported by the University Putra Malaysia Grant no. 5527179.