Variational Iteration Method for Solving Fractional Integro-Differential Equations with Conformable Differointegration

Abstract

Multidimensional integro-differential equations are obtained when the unknown function of several independent variable and/or its derivatives appear under an integral sign. When the differentiation or integration operators or both are of fractional order, the integral equation in this case is called a multidimensional fractional integro-differential equation. Such equations are difficult to solve analytically; therefore, as the main objective of this paper, an approximate method—which is the variational iteration method—will be used to solve this type of equation with conformable fractional-order derivatives and integrals. First, we drive the iterative sequence of approximate solutions using the proposed method, and then, under certain conditions over the kernel of the integro-differential equation, prove its convergence to the exact solution. Two illustrative examples, linear and nonlinear, are given, and their approximated solutions are simulated using computer programs in order to verify from the reliability and applicability of the proposed method.

1. Introduction

Recently, fractional differential equations, fractional calculus and fractional differential integral equations became highly important in several branches of science and engineering because many mathematical models are used to formulate different phenomena, such as mechanics, physics, chemical kinetics, astronomy, biology, economics, potential theory and electrostatistics, which are modeled using integro-differential equations [1,2,3,4,5,6]. In 2016, Shi [7] introduced a formula of mild solutions for impulsive fractional evolution equations.

Many academic researchers continue to be interested in the use of fractional differential equations and/or integral equations, which are based on the development and applications of fractional calculus [8,9]. Nonlinear fractional integral equations and integro-differential equations are notoriously difficult to solve analytically. Moreover, accurate solutions for these equations are extremely rare. As a result, various authors have taken an interest in numerically solving these problems, particularly after the big revolution in computer application. Among the techniques used to solve integro-differential equations are the multi-step methods [10], Adomian decomposition method (ADM) [11], homotopy perturbation method (HPM) [12,13,14], homotopy analysis method (HAM) [15], variational iteration method (VIM) [16] and so on.

The VIM has been successfully applied to solve many problems in different fields of mathematics and its applications. For example, He was the first researcher to propose the use of VIM to solve linear and nonlinear differential and integral equations [17]. In 1998, He used VIM to solve the classical Blasius equation, ref. [18] and in 1999, he provided the approximate solutions for some well-known nonlinear problems [19]. In 2000, He used the VIM to solve autonomous ordinary differential equations. Moreover, in 2006, Soliman applied the VIM to solve the kdv-Burger’s and Lax’s seventh-order kdv equations. In the same year, VIM was applied by Abulwafa and Momani [20] to solve a nonlinear coagulation problem with mass loss and Odibat et al. employed the VIM to solve nonlinear differential equations of fractional order in 2006. In 2006, the VIM was also utilised by Bildiki et al. to solve a variety of problems, including nonlinear partial differential equations, Dehghan and Tateri’s Fokker–Plank equation and quadratic Riccati differential equations with constant coefficients. Wang [21] used the VIM to solve integro-differential equations in 2009, while Sweilam used the VIM to solve both linear and nonlinear boundary value problems of the fourth-order integro-differential equations.

In 2009, Wen-Hua Wang used the VIM to solve certain types of fractional integro-differential equations [21]. In 2011, Muhammet and Adin used the VIM to solve the problem of nonlinear fractional integro-differential equations [20].

If an exact solution exists, the VIM provides rapidly convergent consecutive approximations to the precise solution; otherwise, a few approximations might be employed for numerical results.

In this paper, we shall present the VIM used to solve integro–differential equations with conformable fractional order differointegration of the form:

$$T_x^{\alpha }u(x,y)=g(x,y)+I_x^{\beta }{I}_y^{\gamma }K(x,y,s,t,u(x,y))$$

(1)

where K is given continuous function, $0<\alpha \le 1,\beta ,\gamma >0,x,y\in [a,b]\times [c,d]$ here, $T^{\alpha }$ is understood as a conformable fractional derivative of order $\alpha$, while $I^{\beta }$ and $I^{\gamma }$ stands for conformable fractional-order integrals of order $\beta$ and $\gamma$, respectively.

2. Main Concepts of Factional Calculus

Among the most important definitions of fractional-order derivatives or integrals which will be used next in this paper is the conformable type, which is more simple than other definitions and more stable in comparison with the nonfractional (or integer order) derivatives and integrals.

Definition 1

(Conformable Fractional-Order Derivative [22]). Given a function $f:[a,\infty )\to \mathbb{R}$, then the left conformable fractional derivative of order α can be defined as:

$$\left(T_{\alpha }^a\left(f\right)\right)\left(x\right)=\underset{\epsilon \to 0}{lim}\frac{f\left(x+\epsilon {(x-a)}^{1-\alpha }\right)-f\left(x\right)}{\epsilon }$$

for all $x>0,\alpha \in (0,1]$. When $a=0$, we write $T_{\alpha }$. If $\left(T_{\alpha }^a\left(f\right)\right)\left(x\right)$ exist on the interval $(a,b)$, then define:

$$\left(T_{\alpha }^a\left(f\right)\right)\left(a\right)=\underset{x\to a^{+}}{lim}{f}^{\left(\alpha \right)}\left(x\right)$$

The right conformable fractional derivative of order $\alpha \in (0,1]$ terminating at b of f, is defined by:

$${}_{\alpha }^{b}T\left(f\right)\left(x\right)=lim\frac{f\left(x+\epsilon {(b-x)}^{1-\alpha }\right)-f\left(x\right)}{\epsilon }$$

If ${}_{\alpha }^b\left(T\left(f\right)\left(x\right)\right)$ exist on the interval $(a,b)$, then define:

$${}_{\alpha }^b\left(T\left(f\right)\right)\left(b\right)=\underset{x\to b^{-}}{lim}{}^b{T}_{\alpha }\left(f\right)\left(x\right)$$

Definition 2

(Conformable Fractional-Order Integral [23]). Given a continuous function $f:[a,\infty )\to R$, then the left conformable fractional integral of order α of f is:

$$\left(I_{\alpha }^a\left(f\right)\right)\left(x\right)={\int }_a^{x}f\left(s\right)d\alpha (s,a)={\int }_a^x\frac{f\left(s\right)}{{(s-a)}^{1-\alpha }}ds$$

(2)

where the integral is considered as the usual Riemann improper integral and $a\ge 0$. On the other hand, in the right case, we have:

$${}_{\alpha }^b{I}_x\left(f\right)\left(x\right)={\int }_x^{b}f\left(s\right)d\alpha (b,s)={\int }_x^b\frac{f\left(s\right)}{{(b-s)}^{1-\alpha }}ds$$

(3)

Among some different properties concerning fractional integrals and derivative which are very useful in applications are the following, where $T^{\alpha }$ and $I^{\alpha }$ refer to fractional-order conformable and integral, respectively:

• $T^{\alpha }(c_{1}f+b{c}_2)=c_1{T}^{\alpha }\left(f\right)+c_2{T}^{\alpha }\left(g\right),\mathrm{for}\mathrm{all}{c}_1,c_2\in \mathbb{R}.$
• $T^{\alpha }(c_{1}f+c_{2}g)=c_1{T}^{\alpha }f+c_2{T}^{\alpha }g.$
• $T^{\alpha }(c_{1}f+b{c}_2)=c_1{T}^{\alpha }\left(f\right)+c_2{T}^{\alpha }\left(g\right),\mathrm{for}\mathrm{all}{c}_1,c_2\in \mathbb{R}.$
• $T^{\alpha }\left(x^p\right)=p{x}^{p-\alpha },\mathrm{for}\mathrm{all}p\in \mathbb{R}.$
• $T^{\alpha }\left(\lambda \right)=0,\mathrm{for}\mathrm{all}\mathrm{constant}\mathrm{functions}f\left(x\right)=\lambda .$
• $T^{\alpha }\left(fg\right)=f{T}^{\alpha }\left(g\right)+g{T}^{\alpha }\left(f\right)$
• $T^{\alpha }\left(\frac{f}{g}\right)=\frac{g{T}^{\alpha }\left(f\right)-f{T}^{\alpha }\left(g\right)}{g^2},g\ne 0.$
• If, in addition, f is differentiable, then $T^{\alpha }\left(f\right)\left(x\right)=x^{1-\alpha }\frac{df}{dt}\left(x\right).$
• $T_{\alpha x}^a{I}_{\alpha }^a\left(f\right)\left(x\right)=f\left(x\right)$, for $x>a$, where f is any continuous function in the domain of $I_{\alpha }^a$ and ${}_{\alpha }^b{T}_{\alpha }^b{I}_x\left(f\right)\left(x\right)=f\left(x\right)$, for $x<b$, where f is any continuous function in the domain of ${}_{\alpha }^{b}I$.

3. Variational Iteration Method

The essential aspect of the VIM, as previously stated in the literature, is that the solution of a mathematical problem under the linearization assumption is utilized as a starting approximation or trial function for the next successive approximate solution to the problem under certain conditions [2].

Consider the following general nonlinear equation in operator form to demonstrate the VIM’s essential concept [3]:

$$Au\left(x\right)=g\left(x\right)$$

(4)

and suppose that Equation (1) may be decomposed as:

$$L\left(u\right(x\left)\right)+N\left(u\right(x\left)\right)=g\left(x\right),x\in [a,b]$$

(5)

where A is any operator that may be decomposed into linear and nonlinear operators L and N, respectively, and $g\left(x\right)$ is any function that is referred to as the nonhomogeneous term. Equation (5) may be solved iteratively using the VIM by using the correction functional defined by:

$$u_{n+1}\left(x\right)=u_n\left(x\right)+{\int }_{x_a}^a\lambda (x,s)\{L(u_n\left(s\right)+N\left({\tilde{u}}_n\left(s\right)\right)-g\left(s\right)\}ds,n=0,1,\dots$$

(6)

where $\lambda$ is the general Lagrange multiplier that may ideally be discovered using variational theory the $nth$ approximation of the subscript n denotes the solution u, and ${\tilde{u}}_n$ is considered a restricted variation, i.e., $\delta {\tilde{u}}_n=0$ where the $\delta$ is the first variation [3].

4. Applications of the Vim for Multidimensional Integro-Differential Equations of Fractional Order

Consider the fractional integro-differential Equation (1), which may be written as:

$$T_x^{\alpha }u(x,y)-I_x^{\beta }{I}_y^{\gamma }K(x,y,s,t,u(s,t))-g(x,y)=0$$

(7)

Multiplying Equation (7) by a general Lagrange multiplier $\lambda$, yields to:

$$\lambda (s,t)\{T_x^{\alpha }u(x,y)-I_x^{\beta }{I}_y^{\gamma }K(x,y,s,t,u(s,t))-g(x,y)\}=0$$

(8)

Now, take $I_x^{\alpha }$ to both sides of Equation (8), which give:

$$I_x^{\alpha }[\lambda (s,t)\{T_x^{\alpha }u(x,y)+I_x^{\beta }{I}_y^{\alpha }K(x,y,s,t,u(s,t))\}-g(x,y)]=0$$

(9)

Then, the correction functional with respect to x will be read as follows:

$$u_{n+1}(x,y)=u_n(x,y)+I_x^{\alpha }\left[\lambda (w,y)\{T_w^{\alpha }u(w,y)-I_w^{\beta }{I}_y^{\gamma }K(w,y,s,t,u(s,t))-g(w,y)\}\right]=0$$

(10)

and the problem now is to evaluate $\lambda$. The problem of evaluating $\lambda$ is difficult, since Equation (10) consists of functional derivatives and integrals, so to avoid this difficulty, approximate $I_x^{\alpha }$ and $T_w^a$ for $0<\alpha ⩽1$ by the first integral and derivative. Hence:

$$u_{n+1}(x,y)=u_n(x,y)+{\int }_a^x\left[\lambda \left(x,y\left)\right(\frac{\partial u_n}{\partial w}(w,y)-g(w,y)-I_w^{\beta }{I}_y^{\gamma }k\left(w,y,s,t,\widetilde{u_n}(s,t)\right)\right)\right]dw$$

(11)

Now, taking the first variation of Equation (11) with respect to $u_n$, give:

$$\begin{array}{cc}\hfill \delta u_{n+1}(x,y)& =\delta u_n(x,y)+\delta {\int }_{x_0}^x\left[\lambda \left(\frac{\partial u_n}{\partial s}-g-I_x^{\beta }{I}_y^{y}k\left(s,t,u_n(s,t)\right)\right)\right]dw\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & =\delta u_n(x,y)+\delta {\int }_{x_a}^x\lambda (w,y)\left(\frac{\partial u_n(w,y)}{\partial w}\right)dw\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & =\delta u_n(x,y)+{\int }_a^x\lambda (w,y)\delta \left(\frac{\partial u_n(w,y)}{\partial w}\right)dw\hfill \end{array}$$

(14)

Using integration by parts,

$$\begin{array}{cc}\hfill \delta u_{n+1}(x,y)& =\delta u_n(x,y)+{\left.\delta u_n(w,y)\right|}_{w=x}-\int \frac{\partial \lambda (w,y)}{\partial w}\delta u_{n}dw\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & ={\left.(1+\lambda )\delta u_n(w,y)\right|}_{w=x}-\int \frac{\partial \lambda }{\partial w}\delta u_n(w,y)ds\hfill \end{array}$$

(16)

$$\frac{\partial \lambda }{\partial w}{(w,y)|}_{w=x}=0$$

(17)

$${1+\lambda (w,y)|}_{w=x}=0$$

(18)

Solving this equation will give $\lambda (w,y)=-1$, and substituting in Equation (10), we get:

$$u_{n+1}(x,y)=u_n(x,y)-I_x^{\alpha }\left[T_w^{\alpha }{u}_n(w,y)-g(w,y)-I_w^{\beta }{I}_y^{\gamma }K\right(w,y,s,t,u_n(s,t)\left)\right]$$

(19)

Theorem 1.

Let $u,u_n\in C_x^m\left(\left[a,b\right]\times \left[c,d\right]\right)$, which is a Banach space with a m-th-order continuous partial derivative with respect to x, be the exact and approximate solutions of the integro-differential equation of fractional order (1). If $E_n\left(x,y\right)=u_n\left(x,y\right)-u(x,y)$ and the kernel K satisfy Lipschitz with respect to u, with constant L satisfying

$$L<{\left[\frac{\theta \left(2\theta +1\right)\left(3\theta +1\right)\cdots \left(n\theta +1\right)\gamma \left(2\gamma +1\right)\left(3\gamma +1\right)\cdots \left(n\gamma +1\right)}{{(b-a)}^{n\theta +1}{(d-c)}^{n\gamma +1}}\right]}^{1/n}$$

then, the sequences solutions of approximation $\left\{u_n\right\}$ converge to the exact solution u.

Proof. 

Consider the integro-differential equation of fractional order:

$$T_x^{\alpha }u\left(x,y\right)=g\left(x,y\right)+I_x^{\beta }{I}_y^{\gamma }K(x,y,s,t,u\left(s,t\right)),whereu(0,y)=u_0$$

The approximate solution using the VIM is given by:

$$u_{n+1}\left(x,y\right)=u_n(x,y)-I_x^{\alpha }\{T_w^{\alpha }{u}_n\left(w,y\right)-g\left(w,y\right)-I_w^{\beta }{I}_y^{\gamma }k\left(x,y,s,t,u\left(s,t\right)\right\}$$

(20)

Since u is the exact solution of the integro-differential of the fractional order,

$$u\left(x,y\right)=u(x,y)-I_x^{\alpha }\{T_w^{\alpha }u\left(w,y\right)-g\left(x,y\right)-I_w^{\beta }{I}_y^{\gamma }k\left(x,y,s,t,u\left(s,t\right)\right\}$$

(21)

Then, subtracting Equation (21) from Equation (20), we get:

$$\begin{array}{cc}\hfill E_{n+1}(x,y)=E_n(x,y)-I_x^{\alpha }& T_w^{\alpha }{u}_n(w,y)-g(x,y)+g(x,y)\hfill \\ \hfill -I_w^{\beta }{I}_y^{\gamma }k{u}_n(x,y,s,t,u_n(s,t)-& T_w^{\beta }u(w,y)+I_w^{\beta }{I}_y^{\gamma }K(x,y,s,t,u(s,t)+g(x,y))\hfill \end{array}$$

(22)

$$=E_n\left(x,y\right)-I_x^{\alpha }\{T_w^{\beta }{E}_n\left(w,y\right)-I_w^{\beta }{I}_y^{\gamma }\left\{K\right(w,y,s,t,u_n\left(s,t\right)-K(x,y,s,t,u_n\left(s,t\right)\}$$

From property (3), $I_x^{\alpha }{I}_w^{\alpha }{E}_n\left(x,y\right)=E_n\left(x,y\right)-E_n\left(0,y\right)$ and since $E_n\left(0,y\right)=0$, then

$$E_{n+1}\left(x,y\right)=E_n\left(x,y\right)-E_n\left(x,y\right)+I_x^{\alpha }{I}_x^{\beta }{I}_y^{\gamma }\left[K\left(x,y,s,t,u_n\left(s,t\right)\right)-K\left(x,y,s,t,u_n\left(s,t\right)\right)\right]$$

If $\theta =\alpha +\beta$, then:

$$E_{n+1}\left(x,y\right)=I_x^{\theta }{I}_y^{\gamma }\left[K\left(x,y,s,t,u_n\left(s,t\right)\right)-K\left(x,y,s,t,u\left(s,t\right)\right)\right]$$

(23)

Now, taking the supermum norm for both sides of Equation (23)

$$∥E_{n+1}(x,y)∥\le I_x^{\theta }{I}_y^{\gamma }∥K(x,y,s,t,u_n\left(s,t\right)-K(x,y,s,t,u\left(s,t\right)∥$$

where the K (kernel function) satisfies the Lipschitz condition with constant L, then:

$$∥E_{n+1}(x,y)∥\le L{I}_x^{\theta }{I}_y^{\gamma }∥u_n-u∥$$

$$=L{I}_x^{\theta }{I}_y^{\gamma }\parallel E_n(x,y)\parallel$$

(24)

Using the conformable definition of integrals in Equation (24) implies:

$$∥E_{n+1}(x,y)∥\le L{\int }_a^x{(s-a)}^{\theta -1}{\int }_c^y{(t-c)}^{\gamma -1}∥E_n(s,t)∥dsdt$$

$$=L{\int }_c^x{\int }_a^y{(s-a)}^{\theta -1}{(t-c)}^{\gamma -1}\parallel E_n(x,y)\parallel dsdt$$

(25)

Now, applying mathematical induction over the last inequality:

If $n=0$

$$∥E_1(x,y)∥\le L{\int }_a^x{\int }_c^y{(s-a)}^{\theta -1}{(t-c)}^{\gamma -1}∥E_0(x,y)∥dsdt$$

$$=L{\left.\frac{{(s-a)}^{\theta }}{\theta }\right|}_a^x{\left.\frac{{(t-c)}^{\gamma }}{\gamma }\right|}_c^y∥E_0(x,y)∥$$

$$\le L\frac{{(x-a)}^{\theta }}{\theta }\frac{{(y-c)}^{\gamma }}{\gamma }∥E_0(x,y)∥$$

If $n=1$, then:

$$∥E_2(x,y)∥\le L{\int }_a^x{\int }_c^y{(s-a)}^{\theta -1}{(t-c)}^{\gamma -1}∥E_1(x,y)∥dsdt$$

$$\le L{\int }_a^x{\int }_c^y{\left(s-a\right)}^{\theta -1}{\left(t-c\right)}^{\gamma -1}L\frac{{(s-a)}^{\theta }}{\theta }\frac{{(t-c)}^{\gamma }}{\gamma }∥E_0(x,y)∥dsdt$$

$$=\frac{L^2}{\theta \gamma }{\int }_a^x{\int }_c^y{(s-a)}^{2\theta }{(t-c)}^{2\gamma -1}∥E_0(x,y)∥dsdt$$

$$=\frac{L^2}{\theta \gamma }{\left.\frac{{(s-a)}^{2\theta +1}}{2\theta +1}\right|}_a^x{\left.\frac{{(t-c)}^{2\gamma +1}}{2\gamma +1}\right|}_c^y∥E_0(x,y)∥$$

$$=\frac{L^2}{\theta (2\theta +1)\gamma (2\gamma +1)}{(x-a)}^{2\theta +1}{(y-c)}^{2\gamma +1}∥E_0(x,y)∥$$

If $n=2$

$$∥E_3(x,y)∥\le L{\int }_a^x{\int }_c^y{(s-a)}^{\theta -1}{(t-c)}^{\gamma -1}∥E_2(x,y)∥dsdt$$

$$\begin{array}{c}\hfill \le L{\int }_a^x{\int }_c^y{\left(s-a\right)}^{\theta -1}{\left(t-c\right)}^{\gamma -1}\frac{L^2}{w(2\theta +1)\gamma (2\gamma +1)}{(s-a)}^{2\theta +1}{(t-c)}^{2\gamma +1}∥E_0(x,y)∥dsdt\end{array}$$

$$=\frac{L^3}{\theta (2\theta +1)\gamma (2\gamma +1)}{\int }_a^x{\int }_c^y{\left(s-a\right)}^{3\theta }{\left(t-c\right)}^{3\gamma }∥E_0(x,y)∥dsdt$$

$$=\frac{L^3}{\theta (2\theta +1)\gamma (2\gamma +1)}{\left.\frac{{(s-a)}^{3\theta +1}}{3\theta +1}\right|}_a^x{\left.\frac{{(t-c)}^{3\gamma +1}}{3\gamma +1}\right|}_c^y∥E_0(x,y)∥$$

$$=\frac{L^3}{\theta (2\theta +1)(3\theta +1)\gamma (2\gamma +1)(3\gamma +1)}{(x-a)}^{3\theta +1}{(y-c)}^{3\gamma +1}∥E_0(x,y)∥$$

Hence, by induction, we have:

$$\begin{array}{c}\hfill E_{n+1}(x,y)∥\le \\ & \frac{L^n}{\theta \left(2\theta +1\right)\left(3\theta +1\right)\cdots \left(n\theta +1\right)\gamma \left(2\gamma +1\right)\left(3\gamma +1\right)\cdots \left(n\gamma +1\right)}\hfill & \\ & {(x-a)}^{n\theta +1}{(y-c)}^{n\gamma +1}∥E_0(x,y)∥\hfill \end{array}$$

and upon taking the supremum values of x and y, over $[a,b]\times [c.d]$, getting:

$$\begin{array}{c}\hfill E_{n+1}(x,y)∥\le \\ & \frac{L^n}{\theta \left(2\theta +1\right)\left(3\theta +1\right)\cdots \left(n\theta +1\right)\gamma \left(2\gamma +1\right)\left(3\gamma +1\right)\cdots \left(n\gamma +1\right)}\hfill & \\ & {(b-a)}^{n\theta +1}{(d-c)}^{n\gamma +1}∥E_0(x,y)∥\hfill \end{array}$$

Since

$$\frac{L^n{(b-a)}^{n\theta +1}{(d-c)}^{n\gamma +1}}{\theta \left(2\theta +1\right)\left(3\theta +1\right)\cdots \left(n\theta +1\right)\gamma \left(2\gamma +1\right)\left(3\gamma +1\right)\cdots \left(n\gamma +1\right)}<1$$

Because

$$L<{\left[\frac{\theta \left(2\theta +1\right)\left(3\theta +1\right)\cdots \left(n\theta +1\right)\gamma \left(2\gamma +1\right)\left(3\gamma +1\right)\cdots \left(n\gamma +1\right)}{{(b-a)}^{n\theta +1}{(d-c)}^{n\gamma +1}}\right]}^{1/n}$$

Hence as $asn\to \infty$, we have $∥E_n(x,y)∥\to 0$, i.e., $u_n\left(x,y\right)\to u\left(x,y\right)asn\to \infty$. □

5. Illustrative Examples

Two examples of using the VIM to solve linear and nonlinear integro-differential equations with conformable fractional order differo-integration are presented in this section.

Example 1.

Consider the linear integral equation of fractional order:

$$\begin{array}{c}\hfill T_x^{\alpha }u\left(x,y\right)=g\left(x,y\right)+I_x^{\beta }{I}_y^{\gamma }\left[\left(x-y\right)u(x,y)\right]\end{array}$$

where $g\left(x,y\right)=T_x^{\alpha }u\left(x,y\right)-I_x^{\beta }{I}_y^{\gamma }\left[\left(x-y\right)u(x,y)\right]$.

The exact solution is given for comparison purpose by $u_e(x,y)=x^{3}y$

Hence, starting with the initial guess solution:

$$u_0\left(x,y\right)=g(x,y)=3y{x}^{3-\alpha }-\frac{x^{\beta +4}}{\beta +4}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}.\frac{y^{\gamma +2}}{\gamma +2}$$

then, to find $u_1(x,y)$, if $n=1$, getting:

$$u_1\left(x,y\right)=u_0\left(x,y\right)-I_x^{\alpha }[T_x^{\alpha }{u}_0\left(x,y\right)-g\left(x,y\right)-I_x^{\beta }{I}_y^{\gamma }\left[\left(x-y\right)u(x,y)\right]$$

$$\begin{array}{cc}\hfill & =3y{x}^{2-\alpha }-\frac{x^{\beta +4}}{\beta +4}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}\hfill \\ \hfill & .\frac{y^{\gamma +2}}{\gamma +2}-I_x^{\alpha }[\left(9-3\alpha \right)y^{3-2\alpha }-x^{\beta +4-\alpha }.\frac{y^{\gamma +1}}{\gamma +1}+x^{\beta +3-\alpha }.\hfill \\ \hfill & \frac{y^{\gamma +2}}{\gamma +2}-3y{x}^{2-\alpha }-\frac{x^{\beta +4}}{\beta +4}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}.\frac{y^{\gamma +2}}{\gamma +2}-I_x^{\beta }{I}_y^{\gamma }\hfill \\ \hfill & \left[3y{x}^{2-\alpha }+\frac{x^{\beta +5}}{\beta +4}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +4}}{\beta +3}.\frac{y^{\gamma +2}}{\gamma +2}+3y^2{x}^{3-\alpha }+\frac{x^{\beta +4}}{\beta +4}.\frac{y^{\gamma +2}}{\gamma +1}-\frac{x^{\beta +3}}{\beta +3}.\frac{y^{\gamma +3}}{\gamma +2}\right]\hfill \end{array}$$

(26)

and carrying out recursively fractional-order integrals of the order γ represent to y, of order β with respect to x and of order of α represent to y, we get the final form of $u_1\left(x,y\right)$, which is as follows: when $\alpha =0.8,\beta =0.5and\gamma =0.75$ substitution in Equation (26).

We get the following result of $u_1\left(x,y\right)$ approximate up to six decimals

$$\begin{array}{c}\hfill u_1\left(x,y\right)\cong 0.024161885x^{4.3}{y}^{2.75}+\\ \hfill 1.0x^{3.0}y+0.1029601x^{4.5}{y}^{1.75}& -0.0012449424x^{6.8}{y}^{2.5}-0.023959269x^{5.3}{y}^{1.75}\hfill \\ \hfill -0.11544012x^{3.5}{y}^{2.75}+0.0022746821x^{5.8}{y}^{3.5}& +2.1684043e^{-19}{x}^{2.2}y-0.0012025012x^{4.8}{y}^{4.5}\hfill \end{array}$$

Similarly to the calculations $u_1\left(x,y\right)$, we may find new approximation solution up to three iterations, which are found to be

$$\begin{array}{cc}\hfill u_2\left(x,y\right)=9.09495e-10x^{4.3}{y}^{2.75}& -0.000199751x^{5.6}{y}^{4.5}+1.0x^{3.0}y\hfill \\ \hfill -0.00000507162x^{9.1}{y}^{3.25}+& 2.96859e-9x^{4.5}{y}^{1.75}+0.00100941x^{6.8}{y}^{2.5}+\hfill \\ \hfill 5.67525e-10x^{5.3}{y}^{1.75}& +0.00000595113x^{6.1}{y}^{6.25}-0.0000148071x^{7.1}{y}^{5.25}+4.5693e-\hfill \\ \hfill 9x^{3.5}{y}^{2.75}-0.000185443x^{7.6}{y}^{2.5}& -0.00215173x^{5.8}{y}^{3.5}+0.000359167x^{6.6}{y}^{3.5}+\hfill \\ \hfill 0.0000140055x^{8.1}{y}^{4.25}+& 9.31323e-10x^{2.2}y+\hfill \\ \hfill 0.00133611x^{4.8}{y}^{4.5}\end{array}$$

$$\begin{array}{cc}\hfill u_3\left(x,y\right)=0.0241619x^{4.3}{y}^{2.75}+& 1.0x^{3.0}y+2.00089e-11x^{6.8}{y}^{2.5}-0.0239593\hfill \\ \hfill x^{5.3}{y}^{1.75}+9.41681e-12x^{3.5}{y}^{2.75}+6.36646e-& 12x^{5.8}{y}^{3.5}-\left(-0.386473x^{2.3}{y}^{0.75}+0.879121x^{1.3}{y}^{1.75}\right).\hfill \end{array}$$

$$\begin{array}{cc}\hfill (0.00000525264x^{8.1}{y}^{3.25}+& 4.0829e-11x^{5.8}{y}^{2.5}-7.91503e-7x^{8.9}{y}^{3.25}\hfill \\ \hfill +0.0379687x^{4.3}{y}^{1.75}& +0.00000787187x^{6.1}{y}^{5.25}-0.0000113188x^{7.1}{y}^{4.25}+\hfill \\ \hfill 5.946e-12x^{6.6}{y}^{2.5}& +6.78554e-11x^{4.8}{y}^{3.5}-9.03963e-7x^{6.9}{y}^{5.25}+\hfill \\ \hfill 0.00000150668x^{7.9}{y}^{4.25}-& 3.86567e-8x^{8.4}{y}^6+9.63085e-\hfill \\ \hfill 12x^{5.6}{y}^{3.5}+3.465e-8x^{9.4}{y}^5& -1.26994e-8x^{10.4}{y}^4+\hfill \\ \hfill 1.3112e-20x^{3.5}{y}^{1.75}+& 1.74071e-8x^{7.4}{y}^7)+9.31323e-10x^{2.2}y+4.77485e-\hfill & \\ \hfill 11x^{4.8}{y}^{4.5}\end{array}$$

Figure 1 shows the comparison between the exact and the approximated solution for different values of y.

Example 2.

Let the non-linear equation is:

$$T_x^{\alpha }u\left(x,y\right)=g\left(x,y\right)+I_x^{\beta }{I}_y^{\gamma }\left[\left(xy\right)e^{u(x,y}\right]$$

When the exact solution $u_e(x,y)=x^2{y}^2$.

For the simplicity of calculations and in order to use the properties of conformable and integrations $e^u$ by the Tayler series

$$e^u=1+u+\frac{u^2}{2!}+\cdots$$

After tow terms, we have $e^u\cong 1+u$.

By calculus of variation, the initial condition to

$$u_0\left(x,y\right)=g(x,y)=2y^2{x}^{2-\alpha }-\frac{x^{\beta +1}}{\beta +1}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}.\frac{y^{\gamma +3}}{\gamma +3}$$

Now, using the variation iteration method to find the next approximation solution as follows:

$$\begin{array}{c}\hfill u_1\left(x,y\right)=u_0\left(x,y\right)-I_x^{\alpha }[T_x^{\alpha }{u}_0\left(x,y\right)-g\left(x,y\right)-I_x^{\beta }{I}_y^{\gamma }\left[xy(1+u_0(x,y)\right]\end{array}$$

$$\begin{array}{cc}\hfill =2y^2{x}^{2-\alpha }-\frac{x^{\beta +1}}{\beta +1}& .\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}\hfill \\ \hfill .\frac{y^{\gamma +3}}{\gamma +3}-I_y^{\alpha }[T_x^{\alpha }2y^2{x}^{2-\alpha }& -\frac{x^{\beta +1}}{\beta +1}.\frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}.\hfill \\ \hfill \frac{y^{\gamma +3}}{\gamma +3}-[2y^2{x}^{2-\alpha }-\frac{x^{\beta +1}}{\beta +1}.& \frac{y^{\gamma +1}}{\gamma +1}+\frac{x^{\beta +3}}{\beta +3}.\hfill \\ \hfill \frac{y^{\gamma +3}}{\gamma +3}]& -I_x^{\beta }{I}_y^{\gamma }\left[xy+2y^3{x}^{1-\alpha }-\frac{x^{\beta +2}}{\beta +1}.\frac{y^{\gamma +2}}{\gamma +1}+\frac{x^{\beta +4}}{\beta +3}.\frac{y^{\gamma +4}}{\gamma +3}\right]\hfill \end{array}$$

(27)

Carrying out recursively order integrals of order γ to represent y, of order β to represent x and of order α to represent y, we get the final from of $u_1\left(x,y\right)$, which is as follows.

When $\alpha =0.8,\beta =0.5and\gamma =0.75$ substitution in Equation (27)

We get the following result of $u_1\left(x,y\right)$, which is approximate up to six decimals

$$\begin{array}{cc}\hfill u_1\left(x,y\right)=1.81826e-8x^{5.5}{y}^{5.75}& -0.000477683x^{5.8}{y}^{5.5}\hfill \\ \hfill -0.0177187x^{4.3}{y}^{3.75}+1.0x^{2.0}{y}^2& -0.00250957x^{6.3}{y}^{5.75}+\hfill \\ \hfill 9.31323e^{-10}{x}^{1.2}{y}^2+0.0564374x^{3.5}& y^{3.75}-0.00954771x^{3.8}{y}^{3.5}\hfill \\ \hfill -0.0000386088x^{7.8}{y}^{7.5}\end{array}$$

If n = 1, we get

$$\begin{array}{cc}\hfill u_2\left(x,y\right)=7.27596e-12x^{5.5}{y}^{5.75}+& 0.000353838x^{5.8}{y}^{5.5}\hfill \\ \hfill -0.00000111428x^{8.1}{y}^{7.25}-0.0000562515x^{6.1}{y}^{5.25}& -4.65079e-8x^{4.3}{y}^{3.75}\hfill \\ \hfill -0.00000498821x^{8.6}{y}^{7.5}& +1.0x^{2.0}{y}^2-0.00250957\hfill \\ \hfill x^{6.3}{y}^{5.75}-0.0000841583x^{6.6}{y}^{5.5}& -2.42514e^{-15}{x}^{1.2}{y}^2\hfill \\ \hfill +5.82077e^{-11}{x}^{3.5}{y}^{3.75}& -4.94072e^{-8}{x}^{9.6}{y}^{9.25}\hfill \end{array}$$

If n = 2, we get

$$\begin{array}{cc}\hfill u_3\left(x,y\right)=-4.89841e^{-9}{x}^{10.9}{y}^{9.25}& -7.20333e^{-12}{x}^{5.5}{y}^{5.75}+2.27374e^{-13}\hfill \\ \hfill x^{5.8}{y}^{5.5}+8.25388e^{-7}{x}^{8.1}{y}^{7.25}& -1.61021e^{-7}{x}^{8.9}{y}^{7.25}-7.27596e^{-12}{x}^{4.3}\hfill \\ \hfill y^{3.75}-0.00000498821x^{8.6}{y}^{7.5}+& 1.0x^{2.0}{y}^2-0.00250957x^{6.3}{y}^{5.75}\hfill \\ \hfill -3.03032e^{-10}{x}^{6.6}{y}^{5.5}& -1.25876e^{-7}{x}^{8.4}{y}^7-1.7528e^{-21}{x}^{1.2}{y}^2-\hfill \\ \hfill 1.24008e^{-9}{x}^{10.4}{y}^9-3.40038e^{-11}& x^{11.9}{y}^{11}+1.77679e^{-14}{x}^{7.8}{y}^{7.5}\hfill \end{array}$$

Figure 2 compares the exact and the approximated solutions for different values of y.

6. Conclusions

The present study shows that VIM is a very accurate method that gives the exact solution in a few steps. In some cases, however, it requires more calculations, which will add some difficulties to the problem under consideration. This work may be improved in future by including integro-differential equations with kernels including fractional-order derivatives of the unknown function, in addition to considering fractional-order derivatives greater than 1.