Solution of Higher Order Nonlinear Time-Fractional Reaction Diffusion Equation

Abstract

The approximate analytical solution of fractional order, nonlinear, reaction differential equations, namely the nonlinear diffusion equations, with a given initial condition, is obtained by using the homotopy analysis method. As a demonstration of a good mathematical model, the present article gives graphical presentations of the effect of the reaction terms on the solution profile for various anomalous exponents of particular cases, to predict damping of the field variable. Numerical computations of the convergence control parameter, used to evaluate the convergence of approximate series solution through minimizing error, are also presented graphically for these cases.

1. Introduction

Nonlinear diffusion equations, an important class of parabolic equations, have come from a variety of diffusion phenomena which appear widely in nature. These are suggested as mathematical models of physical problems in many fields, such as filtration, phase transition, electromagnetism, acoustics, electrochemistry, cosmology, biochemistry and dynamics of biological groups. In many cases, the equations possess degeneracy or singularity, which makes the study more involved and challenging. Ultimately, this degeneracy or singularity has enriched the theory of partial differential equations. Science and engineering problems are generally nonlinear, therefore it is important to generate new efficient methods to solve such nonlinear problems. With the help of computerized symbolic computations, many researchers have implemented various methods to establish the solutions to different nonlinear differential equations e.g., the Exp-function method, the Jacobi elliptic function expansion method, the first integral method, the (G’/G)-expansion method, the direct algebraic method, the Cole-Hopf transformation method [1,2,3,4,5,6,7,8] and others. Fractional differential equations are considered as the general form of the differential equations that involved derivatives of any real or complex order. Over the years, fractional calculus has attracted the interest of engineers and scientists immensely.

Fractional differential equations are considered as the general form of differential equations, as they are involved with the derivatives of any real or complex order. In recent years, fractional calculus has attracted the interest of engineers and scientists immensely. The history and the comparative treatment of fractional order have been given by Oldham and Spanier [9], Miller and Ross [10] and Podlubny [11]. The phenomena, which occurs in engineering physics and other branches of science, can be described very successfully. Fractional partial differential equations have comprehensive applications in real world problems compared to integer order ones. Fractional partial differential equations are found to be an effective tool to describe certain physical phenomena, such as diffusion processes and visco-elasticity theories. In this article, we consider a type of time fractional partial differential equation that can be obtained from a standard diffusion equation by replacing the first order time derivative with the Caputo fractional order $\alpha$, $0<\alpha <1$, with higher order nonlinear diffusive and reaction terms. Fractional differential equations have been applied in the modeling of anomalous diffusive and sub-diffusive systems, description of fractional random walk, unification of diffusion. Due to the importance of differential equations of fractional order, many authors are working to find the exact or numerical solutions of the equations. Various methods for obtaining numerical solution have been applied, such as the artificial parameter method, the fractional sub-equation method, the $\delta$-expansion method, the homotopy perturbation method, the generalized differential transform method, the finite difference method, the power series method, and the Exp-function method etc. In 1992, the Chinese mathematician Liao [12] first introduced homotopy to propose an analytic method for the strongly nonlinear problem of the homotopy analysis method (HAM). Thereafter, the HAM has been improved step by step and has been widely applied for solving the nonlinear problems. The advantages of using the HAM as the solution for the nonlinear problems, are that the method is always valid no matter whether there are small physical parameters or not; it provides a simple way to guarantee the convergence of approximation series and it has the flexibility to choose the equation type of linear sub-problems and the base functions of solutions. The main advantage of the HAM is that it always provides a simple way to adjust and control the convergence radius of solution series. Another benefit is that, by studying the HAM in the wide context, we are actually able to obtain a more general iteration method which may exhibit better convergence properties [13,14,15,16,17,18].

The method is a useful analytical approach to get the series solution of fractional order linear and nonlinear partial differential equations (PDEs). This method has certain advantages over some existing numerical methods. Due to discretization used in numerical methods, loss of accuracy during the rounding of errors is possible and also computation requires much time. However, in HAM there is no need for a discretization of variables and it requires less time for computation. Many researchers [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39] have used the method for handling linear and nonlinear problems in fractional order systems.

In this article, an effort has been made to solve a nonlinear fractional order reaction diffusion equation using the HAM. Here, diffusive term is taken as a cubic order polynomial and the reaction term is taken as cubic order. The approximate analytical solution of the equation is found using the powerful technique HAM. The salient feature of the article is the graphical presentation of the probability density function u(x,t) with and without the presence of highly nonlinear reaction term for different Brownian motions and also for the standard motion.

The paper is organized as follows. In Section 2, some necessary definitions of the fractional calculus are provided. In Section 3, the HAM is introduced. In Section 4, we have successfully applied the HAM on time nonlinear fractional order diffusion equation with higher order reaction term. Section 5 consists of the discussion of numerical calculations of the obtained result, followed by a conclusion consisting of the overall work and a future study given in Section 6.

2. Definitions

(a)

A real function $f\left(x\right)$, $x>0$ is said to be in the space $C_{\lambda }$, $\lambda \in \Re$ if $\exists$ a real number $q\left(>\lambda \right)$, such that $f\left(x\right)=x^{q}g\left(x\right)$, where $g\left(x\right)\in C[0,\infty )$. Clearly $C_{\lambda }\subset C_{\mu }$ if $\mu \le \lambda$.

(b)

A function $f\left(x\right)$, $x>0$ is said to be in the space ${C_{\lambda }}^m$, m $\in N\cup \left\{0\right\}$ if $f^{\left(m\right)}\in C_{\lambda }.$

(c)

The (left sided) Riemann–Liouville fractional integral of order $\delta >0$ of a function $f$ $\in C_{\lambda }$, $\lambda \ge -1$ is defined as $J^{\delta }f\left(t\right)=\frac{1}{\mathsf{\Gamma }\left(\delta \right)}{\displaystyle \underset{0}{\overset{t}{\int }}\frac{f\left(\tau \right)}{{\left(t-\tau \right)}^{1-\delta }}}d\tau$, $\delta >0$, $t>0$ where $J^{0}f\left(t\right)=f\left(t\right)$.

(d)

The (left sided) Riemann–Liouville fractional derivative of $f$, $f\in C_{-1}^m$, $m\in N\cup \left\{0\right\}$ of order $\delta >0$, is defined by $D^{\delta }f\left(t\right)=J^{m-\delta }\frac{d^m}{d{t}^m}f\left(t\right)$, $m-1<\delta \le m$, $m\in N$.

(e)

The (left sided) Caputo fractional derivative of $f$, $f\in C_{-1}^m$, $m\in N\cup \left\{0\right\}$ is defined as

$$D_{•}^{\delta }f\left(t\right)=\{\begin{array}{ll}{J}^{m-\delta }{f}^{\left(m\right)}\left(t\right),& m-1<\delta <m,m\in N,\\ \frac{d^m}{d{t}^m}f\left(t\right),& \delta =m.\end{array}$$

3. Basic Idea of HAM

Perturbation techniques are widely applied to obtain analytic approximations of nonlinear equations. However, perturbation methods are essentially based on small physical parameters known as perturbation quantity, but unfortunately many nonlinear problems have no such kind of small physical parameters at all. In addition, neither perturbation techniques, nor the traditional non-perturbation techniques, can provide a way to guarantee the convergence of approximation series. Therefore, both perturbation techniques and the traditional non-perturbation methods mentioned above are, in essence, valid only for weakly nonlinear problems.

We consider the following differential equation

$$N\left[U\left(x,t\right)\right]=0,$$

(1)

where N is a nonlinear operator, r and t are independent variables, $u\left(x,t\right)$ is an unknown function. By means of generalizing the traditional HAM, Liao [14] constructed the so-called zero-order deformation equation as

$$\left(1-q\right)\pounds \left[U\left(x,t;q\right)-u_0\left(x,t\right)\right]=qhH\left(x,t\right)N\left[U\left(x,t;q\right)\right],$$

(2)

where £ is an auxiliary linear operator, $q\in [0,1]$ is the embedding parameter, $h$ and $H\left(x,t\right)$ are a nonzero auxiliary parameter and function, respectively. When $q=0,U\left(x,t;0\right)=u_0\left(x,t\right)$; and when $q=1$, $U\left(x,t;1\right)=u\left(x,t\right)$.

As $q$ increases from 0 to 1, the solution $U\left(x,t;q\right)$ varies from the initial guess $u_0\left(x,t\right)$ to the solution $u\left(x,t\right)$. Assuming that the auxiliary function $H\left(x,t\right)$ and auxiliary parameter h are properly chosen so that $U\left(x,t;q\right)$ can be expressed by the Taylor series as

$$U\left(x,t;q\right)=u_0\left(x,t\right)+{\displaystyle \sum _{m=1}^{\infty }{u}_m}\left(x,t\right)q^m$$

(3)

where $u_m\left(x,t\right)=\frac{1}{m!}{\left(\frac{{\partial }^{m}U\left(x,t;q\right)}{\partial q^m}\right)}_{q=0}$, $m\ge 1$.

The above series is convergent for q = 1. Then by Equation (3), we have

$$U\left(x,t;1\right)=u_0\left(x,t\right)+{\displaystyle \sum _{m=1}^{\infty }{u}_m\left(x,t\right)}.$$

(4)

Differentiating the zero-order deformation Equation (2) m-times with respect to the embedding parameter $q$ and setting $q=0$, and dividing by (m-1)!, we get the m-th order deformation equation as

$$\mathsf{\pounds }\left[u_m\left(x,t\right)-{\chi }_m{u}_{m-1}\left(x,t\right)\right]=hH\left(x,t\right){\mathrm{R}}_m$$

(5)

where the vectors ${\overrightarrow{u}}_m\left(x,t\right)$ are defined by ${\overrightarrow{u}}_m\left(x,t\right)=\left\{u_0\left(x,t\right),u_1\left(x,t\right),\mathrm{\dots },u_m\left(x,t\right)\right\}$, with the initial condition

$$u_m(x,0)=0,$$

(6)

where

$${\mathrm{R}}_m[\overrightarrow{u_{m-1}}(x,t)]=\frac{1}{\left(m-1\right)!}{\left\{\frac{{\partial }^{m-1}}{\partial q^{m-1}}{\rm N}\left[U\left(r,t;q\right)\right]\right\}}_{q=0},$$

(7)

and

$${\chi }_m=0,m\le 1\mathrm{and}{\chi }_m=1,m>1.$$

(8)

On solving the m-th order deformation Equation (5), we get

$$u_m(x,t)={\chi }_m{u}_{m-1}(x,t)+h{J}_t^{\alpha }[{\mathrm{R}}_m(\overrightarrow{u_{m-1}}(x,t))]+c,$$

where $c$ is the constant of integration determined by the initial condition as given by (6).

In this way the N-th order approximation of $u(x,t)$ is obtained by

$$u(x,t)=\underset{N\to \infty }{\lim }{\phi }_N(x,t),$$

(9)

where

$${\phi }_N(x,t)={\displaystyle \sum _{m=0}^{N-1}{u}_m}(x,t).$$

(10)

4. Solution of the Problem

Consider the following fractional diffusion-reaction equation

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{\partial }{\partial x}\left(\left(a+bu+c{u}^2+d{u}^3\right)\frac{\partial u}{\partial x}\right)+ku\left(1-u^3\right),0<\alpha \le 1$$

(11)

with

$$u\left(x,0\right)=x,x\in R,$$

(12)

where u(x,t) is a field variable, which is assumed to vanish for t < 0. The equation will represent a death process for the sink term k > 0 and a birth process for the source term k < 0.

Let the linear operator £ $=D_t^{\alpha }\left[\phi \left(x,t;q\right)\right]$ with the property £ $\left[c\right]=0$, where c is a constant.

We now define a nonlinear operator as

$$\begin{array}{l}N[\mathsf{\phi }\left(x,t;q\right)]=-\frac{{\partial }^{\alpha }\left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial t^{\alpha }}+a\frac{{\partial }^2\left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x^2}+b\left[\mathsf{\phi }\left(x,t;q\right)\right]\frac{{\partial }^2\left[\phi \left(x,t;q\right)\right]}{\partial x^2}+b{\left(\frac{\partial \left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x}\right)}^2\\ +c{\left[\mathsf{\phi }\left(x,t;q\right)\right]}^2\frac{{\partial }^2\left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x^2}+2c\left[\mathsf{\phi }\left(x,t;q\right)\right]{\left(\frac{\partial \left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x}\right)}^2+d{\left[\mathsf{\phi }\left(x,t;q\right)\right]}^3\frac{{\partial }^2\left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x^2}\\ +3d{\left[\mathsf{\phi }\left(x,t;q\right)\right]}^2{\left(\frac{\partial \left[\mathsf{\phi }\left(x,t;q\right)\right]}{\partial x}\right)}^2+k\left[\mathsf{\phi }\left(x,t;q\right)\right]-k{\left[\mathsf{\phi }\left(x,t;q\right)\right]}^4.\end{array}$$

Choosing H(x,t) = 1, the zero-th-order deformation Equation (5) becomes

$$\mathsf{\pounds }\left[u_m\left(x,t\right)-{\chi }_m{u}_{m-1}\left(x,t\right)\right]=\hslash {\mathrm{R}}_m\left(\overline{u_{m-1}}\left(x,t\right)\right),$$

(13)

which gives rise to

$$u_m\left(x,t\right)-{\chi }_m{u}_{m-1}\left(x,t\right)=\hslash J_t^{\alpha }\left[{\mathrm{R}}_m\left({\overrightarrow{u}}_{m-1}\left(x,t\right)\right)\right]$$

where

$$\begin{array}{l}{\mathrm{R}}_m\left(\overline{u_{m-1}}\left(x,t\right)\right)=\frac{1}{\left(m-1\right)!}\left\{\frac{{\partial }^{m-1}}{\partial q^{m-1}}{\rm N}\left[\phi \left(x,t;q\right)\right]\right\}\\ =-\frac{{\partial }^{\alpha }{u}_{m-1}}{\partial t^{\alpha }}+a\frac{{\partial }^2{u}_{m-1}}{\partial x^2}+b{\displaystyle \sum _{i=0}^{m-1}{u}_i\frac{{\partial }^2{u}_{m-1-i}}{\partial x^2}}+b{\displaystyle \sum _{i=0}^{m-1}\frac{\partial u_i}{\partial x}}\frac{\partial u_{m-1-i}}{\partial x}+c{\displaystyle \sum _{i=0}^{m-1}\left({\displaystyle \sum _{j=0}^i{u}_j{u}_{i-j}}\right)}\frac{{\partial }^2{u}_{m-1-i}}{\partial x^2}\\ +2c{\displaystyle \sum _{i=0}^{m-1}\left({\displaystyle \sum _{j=0}^i\frac{\partial u_j}{\partial x}}\frac{\partial u_{i-j}}{\partial x}\right)}{u}_{m-1-i}+d{\displaystyle \sum _{i=0}^{m-1}\left({\displaystyle \sum _{j=0}^i\left({\displaystyle \sum _{k=0}^j{u}_k}{u}_{j-k}\right)u_{i-j}}\right)}\frac{{\partial }^2{u}_{m-1-i}}{\partial x^2}\\ +3d{\displaystyle \sum _{i=0}^{m-1}\left({\displaystyle \sum _{j=0}^i\left({\displaystyle \sum _{k=0}^j{u}_k{u}_{j-k}}\right)\frac{\partial u_{i-j}}{\partial x}}\right)}\frac{\partial u_{m-1-i}}{\partial x}+k{u}_{m-1}-k{\displaystyle \sum _{i=0}^{m-1}\left({\displaystyle \sum _{j=0}^i{\displaystyle \sum _{k=0}^j}{u}_k{u}_{j-k}}\right)}{u}_{i-j}{u}_{m-1-i}.\end{array}$$

(14)

Taking $u_0\left(x,0\right)=x$, we get

$$\begin{array}{l}{u}_1\left(x,t\right)=h\left(b+2cx+3d{x}^2+kx-k{x}^4\right)\frac{t^{\mathsf{\alpha }}}{\mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)},\\ u_2\left(x,t\right)=h(1-h)\left(b+2cx+3d{x}^2+kx-k{x}^4\right)\frac{t^{\alpha }}{\mathsf{\Gamma }\left(1+\alpha \right)}\\ +h^2[\left(6ad+6bc+3bk\right)+\left(24bd+12c^2+8ck+k^2\right)x+\left(60cd+15kd-12ak\right)x^2\\ +\left(-24bk+60d^2\right)x^3+\left(-30ck-5k^2\right)x^4-54dk{x}^4+4k^2{x}^7]\frac{t^{2\alpha }}{\mathsf{\Gamma }\left(1+2\alpha \right)}.\end{array}$$

Proceeding in a similar manner, we can calculate the other components $u_n$ ($n\ge 3$) and finally substituting those expressions in Equations (9) and (10), we get the required approximate series solution of u(x,t).

As given by Liao [12,40] at the m-th order of approximation, one can define the exact square residual error as

$$E_m={\displaystyle {\int }_{\mathsf{\Gamma }}{\left(N{\displaystyle \sum _{i=0}^m{u}_i}(t)\right)}^{2}dt}.$$

(15)

During numerical computation, the limits of the equation will be taken from 0 to 1.The optimal value will be obtained by minimizing the exact residual error defined by the above expression (15) corresponding to the nonlinear algebraic equation

$$\frac{d{E}_m}{dh}=0.$$

5. Numerical Results and Discussion

In this section, numerical results of the probability density function u(x,t) for different Brownian motion $\alpha =0.8,0.9$ and also for standard motion $\alpha =1$ are calculated for three specific cases: (I) $a=b=c=1,d=0$, (II) $a=b=c=0,d=1$ and (III) $a=b=c=d=1.$

The variations of u(x,t) vs. t at $x=1$ for the above three cases at $\alpha =0.8,0.9,1.0$ are depicted through Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6. During numerical computation, only three terms of the series solution are considered and the accuracy can be increased by introducing more terms in the solution.

It is seen from the Figure 1 that if the diffusivity term is considered as a quadratic polynomial as $1+u+u^2$, then it is found that u(x,t) decreases with the increases in t and $\alpha$ due to absence, as well as presence, of reaction terms. The rate of decrease is less (Figure 2),when there will be reaction term ($k=1$) compared to when there will be no reaction term ($k=0$), which reveals that there is the possibility of a highly nonlinear diffusion process to be stabilized for both fractional, as well as standard case, using suitable reaction term.

The nature for the other two considered cases, described through Figure 3, Figure 4, Figure 5 and Figure 6, will be similar. It is also seen from the figures that a decrease of the values of u(x,t) is maximum for the first case and minimum for the third case. This implies that the effect of reaction term will be more if the diffusivity term has quadratic order nonlinearity rather than cubic order nonlinearity.

It is seen from Figure 7, Figure 8 and Figure 9 that the magnitudes of the errors are at a minimum at $h=-0.198372,-0.102184,-0.17583$ for case I, case II, and case III respectively during the presence of reaction term at the standard order case $\alpha =1$ taking three terms of the series solutions.

6. Conclusions

The present article has achieved two important goals. The first one is to demonstrate damping of the probability density function, through the use of proper reaction term on the nonlinear diffusion equation, with a diffusivity term consisting of cubic polynomial for different Brownian motion and also for standard motion. The second one is how to accelerate the convergence of approximate solution using the convergence control parameter through error analysis, which reveals the feasibility and the potential of the applied mathematical tool, HAM, during the solution of highly nonlinear partial differential equations in fractional order system. This makes room for future research work where we will study the effect of the advection term in our considered model.