Solving Fractional Partial Differential Equations with Variable Coefficients by the Reconstruction of Variational Iteration Method

1 Introduction

Fractional differential equations are generalised form of integer order ones, which are achieved by replacing the integer order derivatives by fractional ones. Compared with integer order differential equations, fractional differential equations are more capable for simulating natural physical process and dynamic system [1, 2]. Fractional calculus has become the focus of interest for many researchers in different fields of science and engineering. A lot of research has shown the advantages of using the fractional calculus in the modelling and control of many dynamical systems [3, 4]. Many physics and engineering systems can be elegantly modeled with the help of the fractional derivative, such as dielectric polarisation [5], electrolyte–electrolyte polarisation [6] and viscoelastic systems [7]. Owing to the increasing applications, there has been important interest in developing numerical methods for the solution of fractional differential equations. These methods include Adomian decomposition method (ADM) [8], generalised differential transform method [9], variational iteration method (VIM) [10], finite difference method [11], homotopy analysis method [12] and wavelet method [13]. In this article, our study focuses are on the fractional partial differential equations (FPDEs). FPDEs are generalisations of classical partial differential equations, and they can be classified into three principal kinds: space-time fractional differential equation, space-fractional differential equation and time-fractional differential equation. One of the simplest examples of the former is fractional order diffusion equations, which are generalisations of classical diffusion equations. In the present article, we consider the following FPDE with variable coefficients:

(1)∂αu(x, t)∂tα + Lu(x, t) = q(x, t), a < x < b, 0 < t ≤ τ,0 < α < 1, (1)

where L is the linear operator

(2)L = a0(x) + a1(x)∂β∂xβ + a2(x)∂∂x + a3(x)∂2∂x2, 0 ≤ β ≤ α, (2)

subject to the initial condition

(3)u(x, 0) = f(x), a < x < b, 0 < t ≤ τ, (3)

where ∂^α/∂t^α and ∂^β/∂x^β are fractional derivatives of the Caputo sense, a_i(x) for i=0, …, 3 and q(x, t) are known continuous functions and u(x, t) is the unknown function. The analytic results on existence and uniqueness of solutions to fractional differential equations have been investigated by many authors, see, say for example [14]. Uddin and Haq [15] used radial basis functions for the solution of (1) with constant coefficients. Also, Garrappa and Popolizio [16] used the matrix functions for solving this problem. We construct the solution of problem (1–3) using the reconstruction of variational iteration method (RVIM). The RVIM was first suggested by Hesameddini and Latifizadeh for solving ordinary differential equations (ODE) [17]. In this work, we extend this method to solve FPDEs. It does not only simplify the problem but also speeds up the computation. Our method contains an iterative formula that can provide rapidly convergent successive approximations of the exact solution if such a closed-form solution exists. The rest of this article is as follows. In Section 2, some preliminaries are presented that will be used in later sections. In Section 3, the method is discussed. In Section 4, study on the convergence of this method is illustrated. Section 5 is devoted to numerical experiments, and the results are compared with the exact solutions and some previous works. Section 6 is the conclusion.

2 Preliminaries

For the reader’s convenience, we present some necessary definitions that are used further in this article.

Definition 1 A real function f(t), t>0, is said to be in the space C_μ, μ∈R if there exists a real number p(>μ), such that f(t)=t^pv(t), where v(t)∈C[0, ∞), and it is said to be in the space Cμm if f^m∈C_μ, m∈N.

Definition 2 The Riemann–Liouville fractional integral operator of order α ≥ 0, of a function f(t)∈C_μ, μ ≥ −1, is defined as

Jαf(t) = 1Γ(α)∫0t(t − τ)α − 1f(τ)dτ, α > 0, t > 0,

such that J⁰f(t)=f(t). Note that Γ is the gamma function [14].

Definition 3 For a continuous function f(t), the Caputo fractional partial derivative of order α, where n − 1<α<n, is defined as [18]

∂αf(t)∂tα = Jn − α(∂nf(τ)∂tn) = 1Γ(n − α)∫0t(t − τ)n − α − 1∂nu(x, τ)∂tndτ.

Definition 4 The Laplace transform of the Caputo fractional derivative ∂^αu(x, t)/∂t^α is given by

ℓ{∂αu(x, t)∂tα;s} = sαU(x, s) − ∑k = 0n − 1s(α − k − 1)∂ku(x, t)∂tk(0),

where U(x, s)=ℓ{u(x, t); s}, n − 1 ≤ α<n.

3 Implementation of the Method

To illustrate the basic idea of our technique, we rewrite the general problem (1–3) as follows:

(4)∂αu(x, t)∂tα = g(x, t, a1(x)∂βu(x, t)∂xβ, a2(x)∂u(x, t)∂x,a3(x)∂2u(x, t)∂x2), 0 ≤ β ≤ α, (4)

where

g(x, t, a1(x)∂βu(x, t)∂xβ, a2(x)∂u(x, t)∂x, a3(x)∂2u(x, t)∂x2) = q(x, t) − Lu(x, t),

with the initial conditions

(5)u(x, 0) = f(x), a < x < b, 0 < t ≤ τ. (5)

By taking the Laplace transform on both sides of (4), in the usual way and using the zero artificial initial conditions, the following result is obtained:

(6)sαℓ{u(x, s)} = ℓ{g(x, s, a1(x)∂βu(x, s)∂xβ, a2(x)∂u(x, s)∂x,a3(x)∂2u(x, s)∂x2)}. (6)

Therefore, we can conclude that

(7)ℓ{u(x, s)} = 1sαℓ{g(x, s, a1(x)∂βu(x, s)∂xβ, a2(x)∂u(x, s)∂x,a3(x)∂2u(x, s)∂x2)}. (7)

Suppose that 1/s^α=H(s), then by using the convolution theorem, one obtains

(8)ℓ{u(x, s)} = ℓ{h ∗ g}, (8)

where ℓ⁻¹{H(s)}=h(t). Taking the inverse Laplace transform to both sides of (7), the result is as follows:

(9)u(x, t) = ∫0th(t − τ)g(x, τ, a1(x)∂βu(x, τ)∂xβ, a2(x)∂u(x, τ)∂x,a3(x)∂2u(x, τ)∂x2)dτ. (9)

Now, we must impose the actual initial conditions to obtain the solution of (4). Thus, the following iteration formulation is resulted:

(10)uk + 1(x, t) = u0(x, t) + ∫0th(t − τ)g(x, τ, a1(x)∂βuk(x, τ)∂xβ,a2(x)∂uk(x, τ)∂x, a3(x)∂2uk(x, τ)∂x2)dτ, (10)

where the value of u₀(x, t) is obtained by the actual initial conditions (5) and satisfied in

(11)∂αu0(x, t)∂tα = 0. (11)

The RVIM assumes a series solution for u(x, t) by an infinite sum of components u(x, t)=lim_k→∞u_k(x, t).

Now, define the operator φ[u] as

(12)φ[u] = ∫0th(t − τ)g(x, τ, a1(x)∂βuk(x, τ)∂xβ, a2(x)∂uk(x, τ)∂x,a3(x)∂2uk(x, τ)∂x2)dτ, (12)

and also the components v_k, k=0, 1, … as the following:

(13)v0 = u0, v1 = φ[v0],v2 = φ[v0 + v1] − v1,⋮vk + 1 = φ[∑i = 0kvi] − (∑i = 0kvi). (13)

Using the relations of (10), (12) and (13), result in

(14)u0 = v0,u1 = u0 + φ[u0] = v0 + v1,⋮uk + 1 = u0 + φ[uk] = v0 + φ[∑i = 0kvi] = v0 + ∑i = 1k + 1vi = ∑i = 0k + 1vi. (14)

Then, consequently, we have u(x, t) = limk→∞uk(x, t) = ∑i = 0∞vi(x, t). Therefore, as a result, the solution of problem (4) can be derived in the series form

(15)u(x, t) = ∑i = 0∞vi(x, t). (15)

Corollary 1If 0<α<1, then

(16)∂α∂tαφ[u(x,t)] = g(x, t, a1(x)∂β∂xβu(x, t), a2(x)∂∂xu(x, t),a3(x)∂2∂x2u(x, t)), (16)

where the operatorφ[u] is defined in (12).

Proof. For clarifying the corollary, it is sufficient to take the Laplace transform of both sides of (16), then we have

(17)sαℓ{φ[u(x, s)]} − sα − 1φ[u(x, t)]|t = 0 = sαℓ{φ[u(x,s)]}=ℓ{g(x, s, a1(x)∂βu(x, s)∂xβ, a2(x)∂u(x, s)∂x,a3(x)∂2u(x, s)∂x2)}. (17)

Suppose that 1/s^α=H(s) and ℓ⁻¹{H(s)}=h(t), then by taking the inverse Laplace transform from both sides of (17), one obtains

(18)φ[u(x, t)] = ∫0th(t − τ)g(x, τ, a1(x)∂βu(x, τ)∂xβ,a2(x)∂u(x, τ)∂x, a3(x)∂2u(x, τ)∂x2)dτ. (18)

Therefore, according to the definition of the operator φ[u], the proof is complete.

4 Study on the Convergence of Method

In this section, according to the approach that is proposed in the pervious section, the sufficient conditions for convergence of the method and the error estimate are presented. The main results are proposed in the following theorems.

Theorem 1Let φ, defined in (12), be an operator from a Hilbert space H to H. The series solution u(x, t) = ∑i = 0∞vi(x, t), defined in (15), converges if there exists 0<θ<1, such that||vk + 1|| < θ||vk||, k∈N∪{0}.

Proof. Define a sequence {Sn}n = 0∞ as

s0 = v0, s1 = v0 + v1, …, sn = ∑i = 0nvi.

We show that {Sn}n = 0∞ is a Cauchy sequence in the Hilbert space H. For this purpose, consider

(19)||Sn + 1 − Sn|| = ||vn + 1|| ≤ θ||vn|| ≤ θ2||vn − 1|| ≤ … ≤ θn + 1||v0||. (19)

For every n, m∈N, n ≥ m, we have

(20)||Sn − Sm|| = ||(Sn − Sn − 1) + (Sn − 1 − Sn − 2) + … + (Sm + 1 − Sm)|| ≤||(Sn − Sn − 1)|| + ||(Sn − 1 − Sn − 2)|| + … + ||(Sm + 1 − Sm)|| ≤θn||v0|| + θn − 1||v0|| + … + θm + 1||v0|| =θm + 1(1 + θ + … + θn − m − 1)||v0|| = 1 − θn − m1 − θθm + 1||v0||. (20)

Since 0<θ<1, we get lim_m,n→||S_n − S_m||=0. Therefore, {Sn}n = 0∞ is a Cauchy sequence in the Hilbert space H, and it implies the convergence of series solution u(x, t) = ∑i = 0∞vi(x, t) defined in (15).

Theorem 2If the series solution given by (15) converges, then it is an exact solution of the FPDE with variable coefficients (4).

Proof. Suppose that the series solution (15) converges, then we have

limk→∞ vk = 0, ∑n = 0k(vn + 1 − vn) = vk + 1 − v0.

Also,

(21)∑n = 0∞(vn + 1 − vn) = limk→∞ vk + 1 − v0 = −v0. (21)

Applying the linear operator ∂^α/∂t^α to both sides of (21) and according to relations (11) and (13), we obtain

(22)∂α∂tα(∑n = 0∞(vn + 1 − vn)) = ∂α∂tα(−v0) = 0. (22)

On the other hand, definition (12) implies

(23)∂α∂tα∑n = 0∞(vn + 1 − vn) = ∑n = 0∞∂α∂tα[φ[∑i = 0nvi] − φ[∑i = 0n − 1vi] − vn]. (23)

For n>1, using Corollary 1, one obtains

(24)∂α∂tα[φ[∑i = 0nvi] − φ[∑i = 0n − 1vi] − vn] = ∂α∂tαφ[∑i = 0nvi] − ∂α∂tαφ[∑i = 0n − 1vi] − ∂α∂tαvn =g(x, t, a1(x)∂β∂xβ∑i = 0nvi, a2(x)∂∂x∑i = 0nvi, a3(x)∂2∂x2∑i = 0nvi) −g(x, t, a1(x)∂β∂xβ∑i = 0n − 1vi, a2(x)∂∂x∑i = 0n − 1vi, a3(x)∂2∂x2∑i = 0n − 1vi) − ∂α∂tαvn, (24)

also for n=0, we have

(25)∂α∂tα[φ[v0] − v0] = g(x, t, a1(x)∂β∂xβv0, a2(x)∂∂xv0,a3(x)∂2∂x2v0) − ∂α∂tαv0. (25)

Now substituting relations (24) and (25) in (23), result in

(26)∑n = 0m∂α∂tα[φ[∑i = 0nvi]−φ[∑i = 0n − 1vi]−vn]=[g(x, t, a1(x)∂β∂xβv0, a2(x)∂∂xv0, a3(x)∂2∂x2v0) − ∂α∂tαv0]+[g(x, t, a1(x)∂β∂xβv0 + v1, a2(x)∂∂xv0 + v1, a3(x)∂2∂x2v0 + v1)−g(x, t, a1(x)∂β∂xβv0, a2(x)∂∂xv0, a3(x)∂2∂x2v0) − ∂α∂tαv1]+[g(x, t, a1(x)∂β∂xβv0 + v1 + v2, a2(x)∂∂xv0 + v1 + v2,a3(x)∂2∂x2v0+v1+v2) − g(x, t, a1(x)∂β∂xβv0 + v1, a2(x)∂∂xv0 + v1, a3(x)∂2∂x2v0 + v1) − ∂α∂tαv2]⋮+[g(x, t, a1(x)∂β∂xβ∑i = 0mvi, a2(x)∂∂x∑i = 0mvi, a3(x)∂2∂x2∑i = 0mvi)−g(x, t, a1(x)∂β∂xβ∑i = 0mvi − 1, a2(x)∂∂x∑i = 0mvi − 1, a3(x)∂2∂x2∑i = 0mvi − 1) − ∂α∂tαvm]. (26)

Therefore,

(27)∑n = 0∞∂α∂tα[φ[∑i = 0nvi] − φ[∑i = 0n − 1vi] − vn]= g(x, t, a1(x)∂β∂xβ(∑i = 0∞vi),a2(x)∂∂x(∑i = 0∞vi), a3(x)∂2∂x2(∑i = 0∞vi)) − ∂α∂tα(∑i = 0∞vi). (27)

From (22), (23) and (27), we can observe that u(x, t) = ∑i = 0∞vi(x, t) is an exact solution of the problem (4).

Theorem 3If the series solutionu(x, t) = ∑i = 0∞vi(x, t)defined in (15) converges, then it is an exact solution of the FPDE with variable coefficients (4). Also, if the truncated series∑i = 0mvi(x, t)is used as an approximation to the solution u(x, t) of the problem (4), then the maximum error, E_m(t), is estimated as

Em(t) = ||u(x, t) − ∑i = 0mvi(x, t)|| ≤ 11 − θθm + 1||v0||.

Proof. Theorem 1 and inequality (20) imply that

(28)||Sn − Sm|| = 1 − θn − m1 − θθm + 1||v0||, n ≥ m. (28)

Now, if n→∞, we have S_n→u(x, t). So,

(29)||u(x, t) − ∑i = 0mvi(x, t)|| ≤ 1 − θn − m1 − θθm + 1||v0||. (29)

Also, since 0<θ<1, we have (1 − θ^(n–m)) ≤ 1. Therefore, the preceding inequality becomes ||u(x, t) − ∑i = 0mvi(x, t)|| ≤ 11 − θθm+1||v0||.

In summary, Theorems 1 and 2 indicate that the solution of the problem (4) is obtained by using the iteration formula (10). This approximate formula converges to its exact solution under the condition of existing 0<θ<1 such that ||vn + 1|| ≤ θ||vn||, ∀n∈N∪{0}. In other words, if we define, for every i∈N∪{0}, the parameters

(30)λn = {||vi + 1||||vi||, ||vi|| ≠ 0,0,||vi|| = 0, (30)

then the series solution ∑i = 0∞vi(x, t) of the problem (4) converges to its exact solution, u(x, t), when 0 ≤ λn < 1, ∀n∈N∪{0}. Moreover, as stated in Theorem 3, the maximum absolute truncation error is estimated:

||u(x, t) − ∑i = 0mvi(x, t)|| ≤ 11 − λλm + 1||v0||,

where λ=max{λ_i, i=0, 1, …, m}.

5.1 Numerical Examples

In this section, we apply the proposed approach of RVIM to solve some FPDEs with variable coefficients. We also examine the conditions for convergent RVIM solution of these problems to their exact solutions.

Example 1Consider the following fractional heat-like equation

(31)∂αu(x, t)∂tα = x22∂2u(x, t)∂x2, 0 < α ≤ 1, 0 < t ≤ 1, 0 < x < 1, (31)

subject to the initial condition

(32)u(x, 0) = x2. (32)

The exact solution is u(x, t) = ∑k = 0∞(x2tkα/Γ(kα + 1)). So, for α=1, the exact solution of the problem is u(x, t)=x²e^t.

To obtain the RVIM solution, we apply the Laplace transform on both sides of (31)

(33)ℓ{u(x, s)} = 1sαℓ{x22∂2u(x, s)∂x2}. (33)

Applying the inverse Laplace transform to both sides of (33), result in

(34)u(x, t) = 1Γ(α + 1)∫0t(t − τ)α − 1(x22∂2u(x, τ)∂x2)dτ. (34)

Considering the initial condition (32), the following iterative relation is obtained

(35)un + 1(x, t) = u0(x, t) + 1Γ(α + 1)∫0t(t − τ)α − 1(x22∂2un(x, τ)∂x2)dτ, (35)

where u₀(x, t)=x² and u_n(x, t) indicate the nth approximation of u(x, t). So, φ[u] is defined as

(36)φ[u] = 1Γ(α + 1)∫0t(t − τ)α − 1(x22∂2u(x, τ)∂x2)dτ. (36)

According to (36) and (13), after some simplification and substitution, the following relations are resulted:

v0(x,t) = x2,v1(x, t) = x2tαΓ(α + 1),v2(x, t) = x2t2αΓ(2α + 1),⋮vn + 1(x, t) = x2t(n + 1)αΓ((n + 1)α + 1),

therefore, the solution of (31) will be obtained as

(37)u(x, t) = ∑i = 0∞vi(x, t). (37)

For the approximation purpose, we approximate the solution u(x, t) by the nth-order truncated series u(x, t)=∑i = 0nvi(x, t). Using (37), the closed-form solution of (31) is as follows:

(38)u(x, t) = limn→∞∑k = 0nvk(x, t) = limn→∞∑k = 0nx2tkαΓ(kα + 1) = x2Eα(tα) (38)

that leads to its exact solution, where E_α(z) is the one parameter Mittag–Leffler function defined by [19]

Eα(z) = ∑k = 0∞zkαΓ(kα + 1), α > 0, z∈C.

By using this method, the exact solution of the problem is obtained. To illustrate the effectiveness of our method, we compare it to other numerical methods. For α=0.75, 0.9 and α=1, a comparison between our method and the other methods such as VIM, ADM and Sinc–Legendre collocation method [20] is demonstrated in Tables 1–3. Also, Figure 1 shows the comparison between the absolute error function |u(x, t)|−u_n(x, t)|, which is obtained by the presented method with n=10, α=1 and Sinc–Legendre collocation method. This figure shows that our method can provide more accurate solution than the Sinc–Legendre collocation method [20]. Also, Table 3 shows that our method is clearly reliable if it is compared with the other mentioned methods. In addition, the λ_i for this problem will be obtained as

Figure 1:

Graphs of the resulting absolute error, left (RVIM) and right (Sinc–Legendre), for α=1, in Example 1.

(39)λi = ||vi+1||||vi|| ≤ Γ(iα + 1)Γ((i + 1)α + 1) < 1. (39)

Figure 2 shows that the relation (39) is correct for i ≥ 10. This confirms that the reconstruction of variational approach for problem (31) converges to its exact solution.

Figure 2:

Plot of the function Γ(iα + 1)/Γ((i + 1)α + 1) with i=10.

Example 2Consider the following FPDE of order 0<α<1

(40)∂αu(x, t)∂tα + x∂u(x, t)∂x+∂2u(x, t)∂x2 = 2tα + 2x2 + 2, 0 < x, t ≤ 1, (40)

subject to the initial condition

(41)u(x, 0) = x2. (41)

The exact solution of this problem is u(x, t) = x² + 2(Γ(α + 1)/Γ(2α+1))t^2α [21, 22].

According to RVIM, the recursive relation is given by

(42){un + 1(x, t) = u0(x, t)+1Γ(α + 1)∫0t(t − τ)α − 1(−x∂un(x, τ)∂x − ∂2un(x, τ)∂x2 + 2τα + 2x2 + 2)dτ,u0(x, t) = x2. (42)

Thus, φ[u] is defined as

(43)φ[u] = 1Γ(α + 1)∫0t(t − τ)α − 1(−x∂un(x, τ)∂x − ∂2un(x, τ)∂x2+ 2τα + 2x2 + 2)dτ. (43)

Using (13) and (43), the components v_i can be obtained as follows:

v0(x, t) = x2, v1(x, t)=2Γ(α + 1)t2αΓ(2α + 1), v2(x, t) = 0.

Therefore, the exact solutions is

(44)u(x, t) = ∑i = 0∞vi(x, t) = x2 + 2Γ(α + 1)Γ(2α + 1)t2α. (44)

Furthermore, λ_i=||v_{i+ 1}||/||v_i||=0<1 for i ≥ 1. This confirms that the reconstruction of variational approach for FPDE (40) converges to its exact solution.

Figure 3 shows that the solution of this problem for α=0.6, which is obtained by the presented method, is the same as its exact solution.

Figure 3:

Plot of the RVIM solution for α=0.6 in Example 2.

Example 3 Finally, consider the following FPDE

(45)∂13u(x, t)∂x13 + ∂12u(x, t)∂t12 = 2x53Γ(83) + 2t32Γ(52), 0 ≤ t ≤ 1, 0 ≤ x ≤ 1, (45)

subject to the initial condition u(x, 0)=x².

The exact solution of this problem is u(x, t)=x² + t² [23].

According to the RVIM, the recursive relation is given by

(46){un + 1(x, t) = u0(x, t)  + 1Γ(12)∫0t(t − τ)12 − 1(2x53Γ(83) + 2τ32Γ(52) − ∂13un(x, τ)∂x13)dτ, u0(x, t) = x2. (46)