- Email: [email protected]
Jacobi–Gauss–Lobatto collocation method for solving nonlinear reaction–diffusion equations subject to Dirichlet boundary conditions
Accepted Manuscript
Jacobi-Gauss-Lobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions A.H. Bhrawy, E.H. Doha, M.A. Abdelkawy, R.A. Van Gorder PII: DOI: Reference:
S0307-904X(15)00495-3 10.1016/j.apm.2015.09.009 APM 10693
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
19 April 2013 4 February 2015 4 September 2015
Please cite this article as: A.H. Bhrawy, E.H. Doha, M.A. Abdelkawy, R.A. Van Gorder, Jacobi-GaussLobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.09.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Jacobi-Gauss-Lobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions A.H. Bhrawya,b , E.H. Dohac , M.A. Abdelkawyb , and R.A. Van Gorder
d
a
AN US
CR IP T
Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia (E-mail: [email protected]) b Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt (E-mail:[email protected]) c Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt (E-mail: [email protected]) d Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford OX2 6GG United Kingdom (E-mail: [email protected])
Introduction
PT
1
ED
M
Abstract: This paper extends the application of the spectral Jacobi-Gauss-Lobatto collocation (J-GL-C) method based on Gauss-Lobatto nodes to obtain semi-analytical solutions of nonlinear time-dependent reaction-diffusion equations (RDEs) subject to Dirichlet boundary conditions. This approach has the advantage of allowing us to obtain the solution in terms of the Jacobi parameters α and β, which therefore means that the method holds a number of collocation methods as a special case. In addition, the problem is reduced to the solution of systems of ordinary differential equations (SODEs) in the time variable, which may then be solved by any standard numerical technique. We consider five applications of the general method to concrete examples. In each of the examples considered, the numerical results show that the proposed method is of high accuracy and is efficient for solving nonlinear time-dependent RDEs. Keywords: Nonlinear reaction-diffusion equations; Collocation method; Jacobi-Gauss-Lobatto quadrature.
AC
CE
RDEs have a wide range of applications in many fields of science, and have recently attracted considerable attention, partly due to the interesting features and rich variety of properties of their solutions (see Ref. [1] and the references therein). The processes of diffusion and reaction play an essential role in the dynamics of many systems, e.g., in plasma [2], or semiconductor physics [3]. Spectral methods (see, e.g., [4]-[5] and the references therein) are techniques used in applied mathematics and scientific computing to numerically solve both linear and nonlinear differential equations. The Galerkin, tau and collocation methods are three well-known types of spectral methods. We mention that the spectral collocation method is very useful in providing highly accurate solutions to nonlinear differential equations [5]-[8], also it has become increasingly popular for solving fractional differential equations [8]-[10]. Guo and Yan [7] introduced a new Legendre-Gauss collocation method for solving nonlinear second-order ordinary differential equations. Saadatmandi and Dehghan [11] developed the Sinc-collocation approach for solving multi-point boundary value problem. In this approach, the computation of solution of such problems is reduced to solving algebraic equations. Bhrawy and Alofi [6] proposed the spectral shifted Jacobi-Gauss collocation method to find an accurate solution of the Lane-Emden type equation. Moreover, Doha et al. [12]
1
ACCEPTED MANUSCRIPT
Development of the J-GL-C method
CE
2
PT
ED
M
AN US
CR IP T
developed the shifted Jacobi-Gauss collocation method for solving nonlinear high-order multi-point boundary value problems. For time-dependent partial differential equations (PDEs), spectral methods have been studied for several decades. Glenn et al. [13] investigated the spectral methods for numerically solving time-dependent class of parabolic partial differential equations subject to periodic boundary conditions. Tal-Ezer [14, 15] introduced spectral methods by using the polynomial approximation of the evolution operator in Chebyshev least-squares sense for time-dependent parabolic and hyperbolic equations respectively. Recently, Bhrawy and Al-shomrani [16] proposed a new spectral method for treating the initial-boundary value Korteweg-de Vries equation. The scheme consists of a Jacobi dual-Petrov Galerkin-Jacobi collocation method in space combined with a Crank-Nicholson-leapfrog method in time such that at each time step only a sparse banded linear algebraic system needs to be solved. Solutions given in terms of the Jacobi parameters α and β are possible to obtain, and the use of Jacobi polynomials for solving differential equations has gained increasing popularity in recent years (see, [17]-[20]). To the best of our knowledge, there are no results on the J-GL-C method for solving nonlinear RDEs. The main aim of this paper is to investigate the J-GL-C method for the solution to nonlinear time-dependent RDEs subject to Dirichlet boundary conditions. Indeed, it is very useful to carry out a systematic study on J-GL-C method with general parameters α, β. The nonlinear time-dependent RDE is collocated only for the space variable at (N − 1) points, for suitable collocation points. We use the (N − 1) nodes of the Jacobi-Gauss-Lobatto interpolation which depends upon the two general parameters α, β. These equations together with the two-point boundary conditions constitute system of (N − 1) ordinary differential equations (ODEs) in time. This system can be solved by one of the standard numerical methods. In addition, this algorithm is applied to numerically solve two-dimensional RDEs. Finally, the accuracy of the proposed method is demonstrated by test problems. The paper has two main parts. In Section 2, some properties of Jacobi polynomials are presented. With this mathematical background, we are then able to develop a way of constructing the Gauss-Lobatto collocation technique for a general form of nonlinear time-dependent RDEs, using the Jacobi polynomials. Then, in Section 3, the proposed method is applied to five specific and relevant RDEs. Concluding remarks are made in Section 4.
AC
We start by outlining some properties of Jacobi polynomials. Then we discuss the general class of nonlinear reaction diffusion equation we shall consider. We finally develop the collocation method for this general class of reaction diffusion equation.
2
ACCEPTED MANUSCRIPT
2.1
Basic properties of Jacobi polynomials (α,β)
The standard Jacobi polynomial of degree k, Pk (x) (k = 0, 1, · · · ), with the parameters α > −1, β > −1 [21, 22, 23], satisfies the following relations: (α,β)
Pk
(α,β)
(−x) = (−1)k Pk
(x) ,
CR IP T
(−1)k Γ(k + β + 1) , k!Γ(β + 1) Γ(k + α + 1) (α,β) Pk (1) = . k!Γ(α + 1) (α,β)
Pk
(−1) =
(2.1)
(u, v)w(α,β) =
Z1
AN US
Let w(α,β) (x) = (1 − x)α (1 + x)β . Then we define the weighted space L2w(α,β) as usual, equipped with the following inner product and norm u(x) v(x)w(α,β) (x)dx,
−1
1
kukw(α,β) = (u, u)w2 (α,β) .
(2.2)
The set of Jacobi polynomials forms a complete L2w(α,β) -orthogonal system, and (α,β)
kw(α,β) = hk =
2α+β+1 Γ(k + α + 1)Γ(k + β + 1) . (2k + α + β + 1)Γ(k + 1)Γ(k + α + β + 1)
(2.3)
M
kPk
Let SN (−1, 1) be the set of polynomials of degree at most N on (−1, 1), and due to the property of the standard Jacobi-Gauss quadrature, it follows that for any φ ∈ S2N +1 (−1, 1),
ED
Z1
w(α,β) (x)φ(x)dx =
(α,β)
(α,β)
$N,j φ(xN,j ),
(2.4)
j=0
−1
(α,β)
N X
(α,β)
CE
PT
where xN,j (0 ≤ j ≤ N ) and $N,j (0 ≤ j ≤ N ) are the nodes and the corresponding Christoffel numbers of the Jacobi-Gauss-quadrature formula on the interval (−1, 1), respectively. Now, we introduce the following discrete inner product and norm (u, v)w(α,β) =
N X
(α,β)
(α,β)
(α,β)
u(xN,j ) v(xN,j ) $N,j ,
j=0
1
kukw(α,β) = (u, u)w2 (α,β) .
(2.5)
AC
For α = β one recovers the ultraspherical polynomials (symmetric Jacobi polynomials) and for α = β = ∓ 12 or α = β = 0 one has the Chebyshev of the first and second kinds or Legendre polynomials, respectively. For the nonsymmetric Jacobi polynomials, the two important special cases α = −β = ± 12 give the Chebyshev polynomials of the third and fourth kinds are also recovered.
2.2
Nonlinear reaction-diffusion model
The simple RDE concerning the concentration u of a single substance in one spatial dimension, ut = auxx + C(u), 3
(2.6)
ACCEPTED MANUSCRIPT
CR IP T
is also referred to as the Kolmogorov-Petrovsky-Piskounov (KPP) equation [24]. If the reaction term C(u) vanishes, then the equation represents a pure diffusion process. The corresponding equation is Fick’s second law. Making the choice C(u) = u(1−u2 ) yields the Newell-Whitehead-Segel equation, which describes the Rayleigh-Benard convection [25, 26]. Choosing C(u) = u(1 − u)(u − α) and 0 < α < 1, we obtain the more general Zeldovich equation that arises in combustion theory [27, 28]. Such equations with cubic reaction functions are also known as Nagumo or FitzhughNagumo models; these models appear in a variety of applications, and can include more general terms which allow for density dependence of the diffusion process [29, 30, 31, 32, 33]. In this paper, we consider nonlinear RDEs with a reaction term which can be written in the form ut = (A(u) ux )x + B(u) ux + C(u),
(2.7)
AN US
where u = u(x, t) is an unknown function, and A(u), B(u), C(u) are arbitrary smooth functions. The indices x and t denote differentiation with respect to these variables. Eq. (2.7) is a generalized form of a great number of well-known nonlinear second-order evolution equations describing various processes in biology [34]-[36]. Indeed, the term A(u) permits a density dependent manner of diffusion. Burgers’ equation is a particular case of Eq. (2.7) ut = uxx + λ1 u ux .
(2.8)
ED
M
Burgers simplified the Navier-Stokes equations by dropping the pressure term to obtain his one dimensional equation. This equation has many applications in applied mathematics, such as modeling of gas dynamics [37, 38], modeling of fluid dynamics, turbulence, boundary layer behavior, shock wave formation and traffic flow [39]. The second particular case of Eq. (2.7) is the well-known Fisher equation ut = uxx + λ2 u − λ3 u2 ,
(2.9)
AC
CE
PT
which firstly was introduced by Fisher [40] in order to describe the propagation of a mutant gene. Fisher equation has a wide application in a large number of the fields. It describes the spreading of biological populations [40], chemical kinetics [41], logistic population growth [42], flame propagation [43], population in one-dimensional habitual [44], neutron population in a nuclear reaction [45], neurophysiology [46], autocatalytic chemical reactions [47], branching Brownian motion processes [48], and nuclear reactor theory [49]. We also notice that the Murray equation [34, 35] represents a particular case of Eq. (2.7) ut = uxx + λ1 u ux + λ2 u − λ3 u2 ,
(2.10)
which can be considered as a generalization of the Fisher and Burgers equations. Due to the varied applications of nonlinear RDEs in science and engineering, solution methods for such equations are of great importance. We shall next develop the J-GL-C method for solving such equations under very general assumptions.
2.3
Jacobi spectral collocation method
Since collocation methods give good approximations to differential equations and integral equations in physical space, it is useful to implement and to adapt the methods to various problems, including 4
ACCEPTED MANUSCRIPT
variable coefficient and nonlinear PDEs (see, for instance [6, 8, 50]). In this section, without loss of generality, we derive the J-GL-C method to numerically solve the following general class of nonlinear time-dependent RDEs:
on the domain D = {x : −1 < x < 1} u(−1, t) = g1 (t),
(2.11)
u(1, t) = g2 (t)
(2.12)
AN US
subject to the Dirichlet boundary conditions
CR IP T
ut = A(u) uxx + B(u) ux + C(u)
(2.13)
and the initial condition u(x, 0) = f1 (x),
x ∈ D.
Here we assume A, B and C are sufficiently smooth functions of u. Now, assume that u(x, t) =
N X
(α,β)
aj (t)Pj
(x),
(2.14)
M
j=0
ED
and if we make use of the orthogonality condition of Jacobi polynomials (2.5), then we obtain form (2.14) that aj (t) =
N 1 X (α,β) (α,β) (α,β) (α,β) Pj (xN,i )$N,i u(xN,i , t), hj
(2.15)
i=0
PT
and accordingly, Eq. (2.14) takes the form
CE
u(x, t) =
N N X 1 X (α,β) (α,β) (α,β) (α,β) (α,β) Pj (xN,i )Pj (x)$N,i u(xN,i , t) , hj
(2.16)
N X N X 1 (α,β) (α,β) (α,β) (α,β) (α,β) Pj (xN,i )Pj (x)$N,i u(xN,i , t). hj
(2.17)
j=0
i=0
AC
or equivalently
u(x, t) =
i=0
j=0
The spatial partial derivatives in Eq. (2.11) can be computed at the collocation points to give (α,β)
ux (xN,n , t) =
N X N X 1 (α,β) (α,β) (α,β) (α,β) 0 (α,β) (α,β) Pj (xN,i )(Pj (xN,n )) $N,i u(xN,i , t) hj i=0
=
N X i=0
j=0
(α,β)
Ani u(xN,i , t),
n = 0, 1, · · · , N, 5
(2.18)
ACCEPTED MANUSCRIPT
N X N X 1 (α,β) (α,β) (α,β) (α,β) 00 (α,β) (α,β) Pj (xN,i )(Pj (xN,n )) $N,i u(xN,i , t) hj i=0
=
N X i=0
where
j=0
(α,β)
N X 1 (α,β) (α,β) (α,β) (α,β) 0 (α,β) P (xN,i )(Pj (xN,n )) $N,i , hj j
(2.20)
N X 1 (α,β) (α,β) (α,β) (α,β) 00 (α,β) = Pj (xN,i )(Pj (xN,n )) $N,i . hj
(2.21)
Ani =
j=0
Bni
(2.19)
Bni u(xN,i , t), n = 0, 1, · · · , N,
CR IP T
(α,β)
uxx (xN,n , t) =
j=0
N X
u˙ n (t) = A(un (t))
Bni ui (t) + B(un (t))
i=0
AN US
Making use of Eqs. (2.18)-(2.21), we rewrite Eq. (2.11) in the form: N X
Ani ui (t) + C(un (t)),
i=0
where
(α,β)
ui (t) = u(xN,i , t),
n = 1, · · · , N − 1,
(2.22)
i = 1, · · · , N − 1.
i=1
n = 1, · · · , N − 1,
CE
PT
where
N −1 X
fn (t)) + B(un (t)) ( Bni ui (t) + d
ED
N −1 X
u˙ n (t) = A(un (t)) (
M
Equation (2.22) with the two-point boundary conditions (2.12), constitutes a system of (N − 1) ODEs in the expansion coefficients aj (t), namely Ani ui (t) + dn (t)) + C(un (t)),
i=1
(2.23)
dn (t) = An0 g1 (t) + AnN g2 (t), fn (t) = Bn0 g1 (t) + BnN g2 (t). d
AC
This means that the problem (2.11)-(2.13) is transformed into the following system of ordinary differential equations: N −1 X
u˙ n (t) = A(un (t)) (
i=1
(α,β)
un (0) = f1 (xN,n ),
N −1 X
fn (t)) + B(un (t)) ( Bni ui (t) + d
n = 1, · · · , N − 1,
Ani ui (t) + dn (t)) + C(un (t)),
i=1
(2.24)
which can be written in the matrix form ˙ u(t) = F(t, u(t)) u(0) = f1 , where u(t) ˙ = [u˙ 1 (t), u˙ 2 (t), . . . , u˙ N −1 (t)]T , 6
(2.25)
ACCEPTED MANUSCRIPT
(α,β)
(α,β)
(α,β)
f1 = [f1 (xN,1 ), f1 (xN,2 ), . . . , f1 (xN,N −1 )]T , and F(t, u(t)) = [F1 (t, u(t)), F2 (t, u(t)), . . . , FN −1 (t, u(t))]T , (α,β)
Fn (t, u(t)) = f (t, xN,n , un (t),
N −1 X
Ani ui (t) + dn (t),
i=1
N −1 X
fn (t)). Bni ui (t) + d
CR IP T
where
i=1
The System of ODEs (2.25) can then be solved by using standard numerical solvers.
2.4
Two-dimensional RDEs
Now, we extend the above algorithm to solve the following two-dimensional RDEs: (x, y, t) ∈ Ω1 × Ω2 × Ω3 ,
AN US
ut = a1 uxx +b1 uyy + C(u), where Ω1 = [−1, 1],
Ω2 = [−1, 1]
subject to the initial condition
Ω3 = [0, T ],
(x, y) ∈ Ω1 × Ω2 ,
M
u(x, y, 0) = g1 (x, y),
and
and boundary conditions
u(1, y, t) = g3 (y, t),
u(x, −1, t) = g4 (x, t),
u(x, 1, t) = g5 (x, t)
ED
u(−1, y, t) = g2 (y, t),
(2.26)
(y, t) ∈ Ω2 × Ω3 ,
(x, t) ∈ Ω1 × Ω3 .
(2.27)
(2.28)
PT
The J-GL-C method will be extended to transform the previous PDE into SODEs. To this end, let us expand the dependent variable in the form
CE
u(x, y, t) =
N X M X
(α1 ,β1 )
ai,j (t)Pj
(α2 ,β2 )
(x)Pj
(y).
(2.29)
i=0 j=0
Employing Eqs. (2.14) and (2.5), the coefficients ai,j (t) can be approximated by N M 1 X X (α2 ,β2 ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) (α2 ,β2 ) Pj (yM,k )$M,k Pj (xN,l )$N,l u(xN,l1 1 , yM,k , t). hi hj
AC
ai,j (t) =
l=0 k=0
(2.30)
The expression (2.29) can be rewritten as u(x, y, t) =
N X M X
n,m τl,k ul,k (t),
l=0 k=0
where
(α ,β1 )
u(xN,l1
(α ,β2 )
2 , yM,k
7
, t) = ul,k (t),
(2.31)
ACCEPTED MANUSCRIPT
n,m τl,k =
N X M X i=0 j=0
(α2 ,β2 )
Pj
(α ,β2 )
2 (yM,k
(α ,β2 )
2 )$M,k
(α1 ,β1 )
Pi
(α ,β1 )
(xN,l1
(α ,β1 )
)$N,l1
hi hj
(α1 ,β1 )
Pi
(α2 ,β2 )
(x)Pj
(y).
ux (x, y, t) =
N X M X
CR IP T
In what follows, the approximation of the first-order spatial partial derivative with respect to x for u(x, y, t) can be computed at the Jacobi-Gauss-Lobatto interpolation nodes to give
n,m γl,k ul,k (t),
l=0 k=0
where n,m = γl,k
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X hi hj
0
(α2 ,β2 )
(x)) Pj
(y) ul,k (t),
AN US
l=0 k=0
(α1 ,β1 )
(Pi
(2.32)
(2.33)
also, the first-order spatial partial derivative with respect to y is given by uy (x, y, t) =
N X M X
n,m δl,k ul,k (t),
(2.34)
l=0 k=0
n,m = δl,k
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
M
where
hi hj
ED
l=0 k=0
(α1 ,β1 )
Pi
(α2 ,β2 )
(x)(Pj
0
(y)) ul,k (t). (2.35)
PT
In the same manner, the second spatial partial derivatives with respect to x and y are
n,m ζl,k =
n,m ηl,k =
(2.36)
n,m ηl,k ul,k (t),
(2.37)
l=0 k=0
uyy (x, y, t) =
N X M X l=0 k=0
(Pi
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
Pi
l=0 k=0
and
n,m ζl,k ul,k (t),
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
AC
where
CE
uxx (x, y, t) =
N X M X
l=0 k=0
hi hj
hi hj
8
(α1 ,β1 )
00
(α2 ,β2 )
(x)) Pj
(y) ul,k (t), (2.38)
(α1 ,β1 )
(α2 ,β2 )
(x)(Pj
00
(y)) ul,k (t). (2.39)
ACCEPTED MANUSCRIPT
In the proposed J-GL-C method the residual of (2.26) is set to be zero at (N − 1) × (M − 1) of Jacobi-Gauss-Lobatto interpolation points. Moreover, the values of u0,k (t), uN,k (t), ul,0 (t) and ul,N (t) can be given by
(α ,β1 )
ul,0 (t) = g4 (xN,l1
(α ,β2 )
, t),
2 uN,k (t) = g3 (yM,k
, t),
ul,N (t) = g5 (xN,l1
(α ,β1 )
, t),
k = 0, · · · , M,
, t),
l = 0, · · · , N.
CR IP T
(α ,β2 )
2 u0,k (t) = g2 (yM,k
(2.40)
Therefore, adapting (2.31)-(2.40), enable one to write (2.26)-(2.28) as (N − 1) × (M − 1) SODEs u˙ n,m (t) =a1
N X M X
n,m ζl,k ul,k (t) + b1
subject to the initial data (α ,β1 )
m = 1, · · · , M − 1,
AN US
n = 1, · · · , N − 1, 1 g1n,m = g1 (xN,n
n,m ηl,k ul,k (t) + C(un,m (t)),
l=0 k=0
l=0 k=0
un,m (0) = g1n,m ,
N X M X
(α ,β2 )
2 , yM,k
(α ,β )
2 2 )$M,m ) n = 1, · · · , N − 1,
u1,1 (0)
... u1,M −1 (0) .. u2,1 (0) . u2,M −1 (0) .. . ··· ··· .. . ··· ··· .. . ··· ··· .. uN −2,1 (0) . uN −2,M −1 (0) uN −1,1 (0) . . . uN −1,M −1 (0)
AC
CE
PT
ED
M
Finally, we can arrange (2.41)-(2.42) to their matrix formulation, Runge-Kutta scheme, H1,1 (t, u1 , · · · , uN −1 ) u˙ 1,1 (t) ... u˙ 1,M −1 (t) . .. u˙ 2,1 (t) u˙ 2,M −1 (t) H2,1 (t, u1 , · · · , uN −1 ) .. . ··· · · · · · · . . = . ··· ··· ··· . .. ··· ··· ··· .. u˙ N −2,1 (t) . u˙ N −2,M −1 (t) HN −2,1 (t, u1 , · · · , uN −1 ) HN −1,1 (t, u1 , · · · , uN −1 ) u˙ N −1,1 (t) . . . u˙ N −1,M −1 (t)
where
Hn,m (t, u1 , · · · , uN −1 ) =a1
N X M X
g 2,1 1 ··· = ··· ··· N −2,1 g1 g1N −1,1
n,m ζl,k ul,k (t) + b1
l=0 k=0
n = 1, · · · , N − 1,
g11,1
N X M X
m = 1, · · · , M − 1. (2.42)
which can be solved by implicit
... H1,M −1 (t, u1 , · · · , uN −1 ) .. . H2,M −1 (t, u1 , · · · , uN −1 ) .. . ··· .. . ··· .. . ··· .. . HN −2,M −1 (t, u1 , · · · , uN −1 ) . . . HN −1,M −1 (t, u1 , · · · , uN −1 ) (2.43) 1,M −1 ... g1 .. . g12,M −1 .. . ··· .. (2.44) , . ··· .. . ··· .. . g1N −2,M −1 . . . g1N −1,M −1
n,m ηl,k ul,k (t) + C(un,m (t)),
l=0 k=0
m = 1, · · · , M − 1. 9
(2.41)
(2.45)
,
ACCEPTED MANUSCRIPT
3
Numerical results
3.1
Murray equation
CR IP T
To illustrate the effectiveness of the proposed method, four test examples are considered in this section. Comparison of the results obtained by various choices of Jacobi parameters α and β reveal that the present method is very accurate and efficient for all choices of α and β.
Consider the nonlinear time-dependent one dimensional Murray equation ut = uxx + λ1 uux + λ2 u − λ3 u2 ,
(x, t) ∈ D × [0, 1],
AN US
subject to the boundary conditions λ2 λ2 2 2 1 + tanh (2 λ1 λ3 + (λ1 + 4 λ3 )t) , u(1, t) = 2λ3 8 λ23 λ2 λ2 2 2 u(−1, t) = 1 + tanh (−2 λ λ + (λ + 4 λ )t) , 1 3 1 3 2λ3 8 λ23 and the initial condition
λ2 (2 λ1 λ3 x) , 1 + tanh 8 λ23
x ∈ D.
(3.2)
(3.3)
M
λ2 u(x, 0) = 2λ3
(3.1)
ED
If we apply the generalized tanh method [51], then we find that the exact solution of Eq. (3.1) is λ2 λ2 2 2 u(x, t) = 1 + tanh (2 λ1 λ3 x + (λ1 + 4 λ3 )t) . (3.4) 2λ3 8 λ23
PT
The difference between the measured or inferred value of approximate solution and its actual value (absolute error) is given by E(x, t) = |u(x, t) − u e(x, t)|,
(3.5)
ME = max{E(x, t) : ∀(x, t) ∈ D × [0, 1]},
(3.6)
AC
CE
where u(x, t) and u e(x, t) are the exact solution and the approximate solution at the point (x, t), respectively. Moreover, the maximum absolute error is given by while the rate of convergence may be given by RE =
max{E(x, t) : ∀(x, t) ∈ D × [0, T ]}|N =Ni . max{E(x, t) : ∀(x, t) ∈ D × [0, T ]}|N =Ni+1
(3.7)
Maximum absolute errors and rates of convergence of (3.1) subject to (3.2) and (3.3) are introduced in Table 1 using the J-GL-C method for the two special values α = β = 0 and α = β = − 21 , respectively. We also fix λ1 = λ2 = λ3 = 1 with different choices of N . Meanwhile, the absolute errors of problem (3.1) are presented in Table 2 for α = β = − 12 , λ1 = λ2 = λ3 = 1 and N = 16 with different values of (x, t). 10
ACCEPTED MANUSCRIPT
Table 1 of convergence RE of solutions for equation (3.1) ME RE 4.00 × 10−5 8.00 × 10−8 2.00 × 10−3 2.00 × 10−8 2.50 × 10−1 1.00 × 10−8 5.00 × 10−1 8.00 × 10−6 8.00 × 10−8 1.00 × 10−2 6.00 × 10−8 7.50 × 10−1 1.50 × 10−8 2.50 × 10−1
CR IP T
Maximum absolute errors ME and rate N α β 4 8 0 0 12 16 4 8 − 12 − 21 12 16
A nonlinear RDE
M
3.2
AN US
In Figure 1, we display, in (a) the absolute error (b) the approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 8. Moreover in (c) the absolute error, (d) the approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 12 are plotted. From Figs. 1(c, d), we see that the curves of the approximate and exact solutions coincide for different values of t. This shows that the obtained numerical results are accurate and compare favorably with the exact solution.
Consider the nonlinear time-dependent RDE (a type of Nagumo equation) as given below
ED
ut = λ1 uxx + u(λ2 + λ3 u − λ4 u2 ),
(x, t) ∈ D × [0, 1],
CE
PT
together with the boundary conditions h i √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (1 + λ√3 2λλ1 t) 1 4 4 u(1, t) = , 2λ4 h i √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (−1 + λ√3 2λλ1 t) 1 4 4 u(−1, t) = , 2λ4
AC
and the initial condition
u(x, 0) =
λ3 +
h i p 4λ λ +λ2 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (x) 1 4
2λ4
,
x ∈ D.
(3.8)
(3.9)
(3.10)
Applying the generalized tanh method [51], we find that the exact solution of Eq. (3.8) is given by i h √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (x + λ√3 2λλ1 t) 1 4 4 (3.11) u(x, t) = . 2λ4 Maximum absolute errors and rates of convergence of (3.8) subject to (3.9) and (3.10) are presented in Table 3 using J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, and λ1 = λ2 = λ3 = λ4 = 1 with different choice of N . 11
ACCEPTED MANUSCRIPT
0
CR IP T
0.2
AN US
0
ED
0.1
M
Absolute t α
CE
PT
x - 0.5 -0.4 - 0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 - 0.5 -0.4 - 0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Table 2 errors of solutions for equation (3.1) β E α β E 5.05 × 10−9 3.18 × 10−9 4.93 × 10−9 2.72 × 10−9 −9 4.66 × 10 2.14 × 10−9 −9 4.28 × 10 1.53 × 10−9 −9 3.82 × 10 8.82 × 10−10 0 3.26 × 10−9 − 21 − 12 1.80 × 10−10 2.63 × 10−9 5.44 × 10−10 −9 1.95 × 10 1.25 × 10−9 1.21 × 10−9 1.94 × 10−9 −10 4.40 × 10 2.61 × 10−9 −10 3.18 × 10 3.18 × 10−9 7.31 × 10−10 3.32 × 10−9 −10 5.61 × 10 3.32 × 10−9 −10 3.52 × 10 3.18 × 10−9 −10 1.96 × 10 3.07 × 10−9 7.71 × 10−11 2.93 × 10−9 1 1 −11 0 3.60 × 10 − 2 − 2 2.68 × 10−9 −10 1.11 × 10 2.40 × 10−9 1.35 × 10−10 2.16 × 10−9 −10 1.65 × 10 1.87 × 10−9 −10 2.19 × 10 1.46 × 10−9 2.07 × 10−10 1.10 × 10−9
AC
Table 3 Maximum absolute errors and rate of convergence N α β ME 4 6.00 × 10−3 8 0 0 3.00 × 10−5 12 1.00 × 10−7 16 4.00 × 10−9 4 4.00 × 10−3 1 1 8 − 2 − 2 1.00 × 10−5 12 4.00 × 10−8 16 1.00 × 10−8
12
of the solution to equation (3.8) RE 5.00 × 10−3 3.33 × 10−3 4.00 × 10−2 2.50 × 10−3 4.00 × 10−3 2.50 × 10−1
ACCEPTED MANUSCRIPT
3.3
errors and rate of p ME 1.50 × 10−4 3 1.00 × 10−8 8.00 × 10−13 6.00 × 10−5 3 1.50 × 10−11 1.50 × 10−13
Table 4 convergence for the solution of RE p ME 3.00 × 10−3 6.67 × 10−5 5 6.00 × 10−5 8.00 × 10−5 4.00 × 10−7 2.00 × 10−2 −6 2.50 × 10 5 3.00 × 10−6 1.00 × 10−2 1.50 × 10−8
equation (3.12) RE 2.00 × 10−2 6.67 × 10−3
CR IP T
Maximum absolute N α β 4 8 0 0 12 4 8 − 21 − 12 12
1.50 × 10−4 5.00 × 10−3
Reaction-diffusion equation with power-law density dependent diffusion
AN US
In the following we analyse the nonlinear time-dependent one dimensional RDE which includes density dependent diffusion ut = (λ1 up−1 ux )x − λ2 u + λ3 up , The boundary conditions are given by
and subjected to the initial condition 2bp q(p + 1)
1 p−1
ED
u(x, 0) =
2bp q(p + 1)
1 p−1
λ3 > 0,
cos
p−1 2
(x, t) ∈ D × [0, 1] .
r
q ap
2 p−1
(3.12)
,
(3.13)
x ∈ D.
(3.14)
M
u(−1, t) = u(1, t) =
p > 1,
2 r p−1 p−1 q cos x , 2 ap
The exact solution [1] of Eq. (3.12) has the form
PT
u(x, t) =
2bp q(p + 1)
1 p−1
2 r p−1 p−1 q cos x . 2 ap
(3.15)
AC
CE
Maximum absolute errors and rates of convergence of (3.12) subject to (3.13) and (3.14) are depicted in Table 4 using the J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, and λ1 = λ2 = λ3 = 1 with various choice of N and p. In Figure 2, we present the absolute error of problem (3.12) where λ1 = λ2 = λ3 = 1, p = 5 α = β = −0.5 and N = 4, 8 and 12, respectively.
3.4
Diffusion- areaction equation with Fischer-Kolmogorov reaction term and density dependent diffusion
Finally, we investigate the nonlinear time-dependent one dimensional RDE presented below, which includes a Fischer-Kolmogorov reaction term and density dependent diffusion: ut = (u2 )xx + u(1 − u),
13
(x, t) ∈ D × [0, 1],
(3.16)
ACCEPTED MANUSCRIPT
for the solution of equation (3.16) RE 6.67 × 10−6 3.00 × 10−1 5.00 × 10−4 3.00 × 10−1
with to the boundary conditions
and the initial condition
2 , 1 + tanh −1+t 4 u(x, 0) =
u(1, t) =
2 , 1 + tanh 1+t 4
AN US
u(−1, t) =
CR IP T
Table 5 and rate of convergence α β ME 3.00 × 10−3 0 0 2.00 × 10−8 6.00 × 10−9 4.00 × 10−5 1 1 − 2 − 2 2.00 × 10−8 6.00 × 10−9
Maximum absolute errors N 4 8 12 4 8 12
2 , 1 + tanh x4
x ∈ D.
(3.17)
(3.18)
We notice that the exact solution [44] of Eq. (3.16) is
2 . 1 + tanh x+t 4
(3.19)
M
u(x, t) =
Two-dimensional nonlinear RDE
CE
3.5
PT
ED
Maximum absolute errors and rates of convergence of (3.16) subject to (3.17) and (3.18) can be seen in Table 5, where we have used the J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, with different values of N . On the other hand, the absolute errors of problem (3.16) are exhibited in Table 6 for α = β = − 12 and N = 12 with different values of (x, t). In Figure 3, we plot both the approximate and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.16) where N = 12. We remark that the curves of the approximate and exact solutions agree strongly for different values of t which are listed in their captions. This shows that the obtained numerical results are accurate and compare favorably with the exact solution.
AC
In order to confirm the high accuracy of our technique for the two-dimensional problem, we consider the 2+1 nonlinear RDE 1 ∂t u(x, y, t) = ∇2 u(x, y, t) + (u(x, y, t))2 (1 − u(x, y, t)), (x, y, t) ∈ [−1, 1] × [−1, 1] × [0, 20], 2 (3.20)
with the following initial-boundary conditions 1
u(−1, y, t) = 1+e u(x, −1, t) =
√1 (y−1− √1 t) 2 2
1 1+e
√1 (x−1− √1 t) 2 2
1+e 1+e
1+e
1 √ (x+y) 2
14
.
1 √ (y+1− √1 t) 2 2
1
, u(x, 1, t) = 1
u(x, y, 0) =
1
, u(1, y, t) =
1 √1 (x+1− √ t) 2 2
, ,
(3.21)
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Table 6 Absolute errors for the solution of equation (3.16) x t α β E α β E - 0.5 3.19 × 10−9 3.31 × 10−9 -0.4 3.50 × 10−9 3.66 × 10−9 −9 - 0.3 3.73 × 10 3.90 × 10−9 −9 -0.2 3.88 × 10 4.06 × 10−9 -0.1 3.96 × 10−9 4.14 × 10−9 1 1 −9 0 0.1 0 0 3.96 × 10 − 2 − 2 4.15 × 10−9 −9 0.1 3.90 × 10 4.09 × 10−9 −9 0.2 3.76 × 10 3.95 × 10−9 0.3 3.56 × 10−9 3.73 × 10−9 −10 0.4 3.29 × 10 3.46 × 10−9 −10 0.5 2.96 × 10 3.11 × 10−9 - 0.5 1.50 × 10−9 1.90 × 10−9 −9 -0.4 1.41 × 10 1.79 × 10−9 −9 - 0.3 1.23 × 10 1.57 × 10−9 -0.2 1.01 × 10−9 1.32 × 10−9 −10 -0.1 7.92 × 10 1.06 × 10−9 1 1 −10 0 0.2 0 0 5.96 × 10 − 2 − 2 8.35 × 10−10 0.1 4.32 × 10−10 6.46 × 10−10 −10 0.2 3.11 × 10 5.06 × 10−10 −10 0.3 2.44 × 10 4.21 × 10−10 −10 0.4 2.29 × 10 3.93 × 10−10 0.5 2.48 × 10−10 4.05 × 10−10
15
ACCEPTED MANUSCRIPT
The exact solution of this problem is 1
u(x, y, t) =
.
(3.22)
AN US
1+e
1 √ (x+y− √1 t) 2 2
ME 1.44 × 10−3 2.51 × 10−7 1.11 × 10−8 314 × 10−4 5.14 × 10−8 2.35 × 10−9
CR IP T
N =M =K 4 8 12 4 8 12
Table 7 Maximum absolute errors of problem (3.20) (α1 , β1 ) (α2 , β2 ) ME (α1 , β1 ) (α2 , β2 ) ( 21 , 21 ) (− 12 , 12 ) 1.14 × 10−3 (0, 12 ) (− 12 , 0) −7 3.49 × 10 2.89 × 10−9 1 1 (1, 0) (− 2 , − 2 ) 1.38 × 10−3 (− 12 , − 12 ) (− 12 , − 12 ) 4.19 × 10−7 3.67 × 10−9
Table 7 lists the maximum absolute errors of u(x, y, t) of problem (3.20) with various choices of α1 , α2 , β1 , β2 , , N, M, and K. The numerical results presented in this table show that the results are very accurate for small value of N, M, and K. Figure 4 demonstrates that the absolute errors E(x, y, t), at t = 1 and E(x, y, t), at y = t = 12 are very small for even the small number of grid points taken for all the Jacobi parameters.
Conclusions
M
4
AC
CE
PT
ED
We have investigated the J-GL-C method for the solution of nonlinear time-dependent RDEs subject to Dirichlet boundary conditions. The nonlinear time-dependent RDE has been collocated for the space variable at (N − 1) points, for suitable collocation points. We use the (N − 1) nodes of the Jacobi-Gauss-Lobatto interpolation which depends upon two general parameters α, β > −1. These equations, together with the two-point boundary conditions, yield a system of (N − 1) ODEs in the time variable. The resulting system can then be solved by one of the standard numerical methods. The solution method for a very general class of RDEs was considered in Section 2. This general class of equations holds a number of relevant physical and biological models as special cases. In Section 3, we have considered four specific nonlinear RDEs, in order to demonstrate the method. For each of the examples given, we were able to obtain accurate solutions with relatively few nodes used. Indeed, we often obtain very accurate results with fewer than ten nodes. As such, the method is quite efficient, requiring minimal computation to achieve a desired accuracy. Another facet of the solution method is that it permits solutions which appear to converge rapidly, as measured by the parameter RE , the rate of convergence. We find, in some cases, that the error decreases by as much as a factor of 10−3 due to an increase of only a single node. In the present paper, we have considered examples with power-law or polynomial type nonlinearity. However, the J-GL-C method is rather general, and one may consider other kinds of nonlinearity. In some applications [52], exponential nonlinearity occurs in reaction diffusion models, and such cases could also be addressed by the present method. Additionally, one sometimes encounters coupled RDEs. If the terms coupling the equations is strongly nonlinear, the solution of such equations can be rather difficult. Hence, for future work, it may be interesting to consider coupled RDEs with various forms of nonlinearity. 16
ACCEPTED MANUSCRIPT
CR IP T
The results presented here demonstrate that the J-GL-C method is both accurate and efficient for the computation of numerical solutions to nonlinear RDEs, even those with density-dependent diffusion terms and strongly nonlinear reaction functions. Therefore, the method warrants consideration from those involved in the numerical simulation of models arising in science or engineering which take the form of nonlinear RDEs.
References
[1] A.H. Khater, W. Malfliet, D.K. Callebaut, E.S. Kamel, The tanh method, a simple transformation and exact analytical solutions for nonlinear reaction-diffusion equations, Chaos Soliton. Fract., 14 (2002) 513-522.
AN US
[2] K. Mayawala, D. G. Vlachos, J. S. Edwards, Spatial modeling of dimerization reaction dynamics in the plasma membrane: Monte Carlo vs. continuum differential equations, Biophysical Chemistry, 121(3) (2006) 194-208. [3] P. A Markowich, P. Szmolyan, A system of convection-diffusion equations with small diffusion coefficient arising in semiconductor physics, J. Diff. Eqs., 81(2) (1989) 234-254. [4] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods: Fundamentals in Single Domains. Springer-Verlag, New York (2006).
M
[5] S. Nguyen and C. Delcarte, A spectral collocation method to solve Helmholtz problems with boundary conditions involving mixed tangential and normal derivatives, J. Comput. Phys., 200 (2004) 34-49.
ED
[6] A.H. Bhrawy and A.S. Alofi, A Jacobi-Gauss collocation method for solving nonlinear LaneEmden type equations, Commun. Nonlinear. Sci. Numer. Simulat., 17 (2012) 62-70.
PT
[7] B.-y. Guo and J.-p. Yan, Legendre-Gauss collocation method for initial value problems of second order ordinary differential equations, Appl. Numer. Math., 59 (2009) 1386-1408.
CE
[8] A.H. Bhrawy and Alghamdi, A shifted Jacobi-Gauss-Lobatto collocation method for solving nonlinear fractional Langevin equation involving two fractional orders in different intervals, Bound. Val. Prob., 2012 (2012) 62.
AC
[9] L. Zhu and Q. Fan, Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 2333-2341
[10] E.H. Doha, A.H. Bhrawy and S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equation, Appl. Math. Model., 36 (2012) 4931-4943.
[11] A. Saadatmandi and M. Dehghan, The use of Sinc-collocation method for solving multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 593-601. [12] E.H. Doha, A.H. Bhrawy and R.M. Hafez, On shifted Jacobi spectral method for high-order multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 3802-3810. 17
ACCEPTED MANUSCRIPT
[13] I. Glenn, S. Brian, W. Rodney, Spectral methods in time for a class of parabolic partial differential equations, J. Comput. Phys., 102 (1992) 88-97. [14] H. Tal-Ezer, Spectral methods in time for hyperbolic problems, SIAM J. Numer. Anal., 23 (1986) 11-26.
CR IP T
[15] H. Tal-Ezer, Spectral methods in time for parabolic problems, SIAM J. Numer. Anal., 26 (1989) 1-11. [16] A.H. Bhrawy, M. Al-shomrani, A Jacobi dual-Petrov Galerkin-Jacobi collocation method for solving Korteweg-de Vries equations, Abstract and Applied Analysis, 2012, (2012) Article ID 418943, 16
AN US
[17] E.H. Doha, A.H. Bhrawy, W.M. Abd-Elhameed, Jacobi spectral Galerkin method for elliptic Neumann problems, Numerical Algorithms, 50 (2009), 67-91 [18] E.H. Doha and A.H. Bhrawy, A Jacobi spectral Galerkin method for the integrated forms of fourth-order elliptic differential equations, Numerical Methods for Partial Differential Equations, 25 (2009) 712-739 [19] M. El-Kady, Jacobi discrete approximation for solving optimal control problems, J. Korean Math. Soc. 49 (2012) 99-112.
M
[20] E.H. Doha, A.H. Bhrawy and R.M. Hafez, A Jacobi-Jacobi dual-Petrov-Galerkin method for third- and fifth-order differential equations, Math. Comput. Model., 53 (2011) 1820-1832.
ED
[21] R. Beals, R. Wong, Special Functions, Cambridge University Press (2010). [22] G. Szeg¨ o, Orthogonal Polynomials, Colloquium Publications, XXIII, American Mathematical Society, ISBN 978-0-8218-1023-1, MR 0372517G. (1939).
PT
[23] Y. Luke, The special functions and their approximations, Academic Press, New York, 1969.
CE
[24] A. Kolmogorov, I. Petrovskii and N. Piscounov, A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem, In V. M. Tikhomirov, editor, Selected Works of A. N. Kolmogorov I, 248-270. Kluwer 1991. Translated by V. M. Volosov from Bull. Moscow Univ., Math. Mech., 1 (1937) 1-25.
AC
[25] A. C. Newell and J. A. Whitehead, Finite bandwidth, finite amplitude convection. J. Fluid Mech., 38 (1969) 279-303.
[26] L. A. Segel, Distant sidewalls cause slow amplitude modulation of cellular convection, J. Fluid Mech., 38 (1969) 203-224. [27] Y. B. Zeldovich and D. A. Frank-Kamenetsky, A theory of thermal propagation of flame, Acta Physicochim, 9 (1938) 341-350.
[28] B. H. Gilding and R. Kersner, Travelling Waves in Nonlinear Diffusion Convection Reaction, Nonlinear Differential Equations and their Applications, 60. Birkh¨auser Verlag, Basel (2004) 209 pp. 18
ACCEPTED MANUSCRIPT
[29] R. A. Van Gorder and K. Vajravelu, Analytical and numerical solutions of the density dependent diffusion Nagumo equation, Physics Letters A 372 (2008) 5152-5158.
CR IP T
[30] R. A. Van Gorder and K. Vajravelu, A variational formulation of the Nagumo reactiondiffusion equation and the Nagumo telegraph equation, Nonlinear Analysis Series B: Real World Applications 11 (2010) 2957-2962. [31] R. A. Van Gorder and K. Vajravelu, Analytical and Numerical Solutions for a Density Dependent Nagumo Telegraph Equation, Nonlinear Analysis Series B: Real World Applications 11 (2010) 3923-3929. [32] R. A. Van Gorder, Gaussian waves in the Fitzhugh-Nagumo equation demonstrate one role of the auxiliary function H(x) in the homotopy analysis method, Communications in Nonlinear Science and Numerical Simulation 17 (2012) 1233-1240.
AN US
[33] R. A. Van Gorder, Travelling waves for a density dependent diffusion Nagumo equation over the real line, Communications in Theoretical Physics 58 (2012) 5-11. [34] J.D. Murray, Nonlinear Differential Equation Models in Biology, Clarendon Press, Oxford, 1977. [35] J.D. Murray, Mathematical Biology, Springer, Berlin, 1989.
M
[36] P.C. Fife, Mathematical Aspects of Reacting and Diffusing Systems, Springer, Berlin, 1979.
ED
[37] A. Veksler and Y. Zarmi, Wave interactions and the analysis of the perturbed Burgers equation, Physica D, 211 (2005) 57-73. [38] A. Veksler, Y. Zarmi, Freedom in the expansion and obstacles to integrability in multiplesoliton solutions of the perturbed KdV equation, Physica D, 217 (2006) 77-87.
PT
[39] Z.-Y. Ma, X.-F. Wu and J.-M. Zhu, Multisoliton excitations for the Kadomtsev-Petviashvili equation and the coupled Burgers equation, Chaos Soliton. Fract., 31 (2007) 648-657.
CE
[40] R. A. Fisher, The wave of advance of advantageous genes. Ann. Eugen., 7 (1937) 353-369. [41] A. J. Khattak, A computational meshless method for the generalized Burgers’-Huxley equation, Appl. Math. Modell., 33 (2009) 3718-3729.
AC
[42] N. F. Britton, Reaction-diffusion equations and their applications to biology. London: Academic Press; 1986. [43] D. A. Frank, Diffusion and heat exchange in chemical kinetics. New Jersey: Princeton University Press; 1955. [44] A-M. Wazwaz, The extended tanh method for abundant solitary wave solutions of nonlinear wave equations, Appl. Math. Comput., 187 (2007) 1131-1142. [45] W. Malfliet, Solitary wave solutions of nonlinear wave equations, Am. J. Phys., 60 (7) (1992) 650-654.
19
ACCEPTED MANUSCRIPT
[46] Y. Tan, H. Xu, S-J. Liao, Explicit series solution of travelling waves with a front of Fisher equation, Chaos Soliton. Fract., 31 (2007) 462-472. [47] H. N.A. Ismail, K. Raslan, A. A. Abd Rabboh, Adomian decomposition method for Burgers’Huxley and Burgers’-Fisher equations, Appl. Math. Comput., 159 (2004) 291-301.
CR IP T
[48] A. J. Khattak, A computational meshless method for the generalized Burgers’-Huxley equation, Appl. Math. Modell., 33 (2009) 3718-3729. [49] J. Canosa, Diffusion in nonlinear multiplicate media, J. Math. Phys., 31 (1969) 1862-1869. [50] A.H. Bhrawy, An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system, Applied Mathematics and Computation, 247, (2014) 30-46.
AN US
[51] E. Fan and Y. C. Hon, Generalized tanh method extended to special types of nonlinear equations, Z. Naturforsch 57a (2002) 692-700.
AC
CE
PT
ED
M
[52] R. A. Van Gorder and K. Vajravelu, Nonlinear dispersion of a pollutant ejected into a channel flow, Central European Journal of Physics 9 (2011) 1182-1194.
20
ACCEPTED MANUSCRIPT
0.8
uHx,0L
0.7
u Hx,0L
CR IP T
u amd u
0.9
uHx,0.5L
0.6
u Hx,0.5L
0.5
uHx,0.9L
0.4
u Hx,0.9L
0.3 -1.0
0.0
0.5
1.0
AN US
-0.5
x
b
PT
0.8
uHx,0L
0.7
u Hx,0L uHx,0.5L
0.6
u Hx,0.5L
0.5
uHx,0.9L 0.4
CE AC
0.9
u amd u
ED
M
a
u Hx,0.9L
0.3 -1.0
-0.5
0.0
0.5
1.0
x
c
d
Figure 1. We plot the (a) absolute error (b) approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 8. Then, we plot the (c) absolute error (d) approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 12.
21
AN US
CR IP T
ACCEPTED MANUSCRIPT
b
AC
CE
PT
ED
M
a
c
Figure 2. The absolute error of problem (3.12) where λ1 = λ2 = λ3 = 1, p = 5 α = β = −0.5 and (a) N = 4, (b) N = 8 and (c) N = 12, respectively.
22
ACCEPTED MANUSCRIPT
2.4 2.3 uHx,0L
CR IP T
2.2
u Hx,0L
u amd u
2.1
uHx,0.5L
2.0 1.9
u Hx,0.5L
1.8
uHx,0.9L
1.7 1.6 -1.0
AN US
u Hx,0.9L
0.0
-0.5
0.5
1.0
x
ED
M
Figure 3. The approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.16) where N = 12, respectively.
EHx,0.5,0L
1.0
CE
1. ´ 10-9
6. ´ 10-12
PT
2. ´ 10-9 EHx,y,1L 1.5 ´ 10-9
8. ´ 10-12
4. ´ 10-12
0.5
5. ´ 10-10
2. ´ 10-12
0 -1.0
0.0
y
AC
-0.5
0.0 x
-0.5
0
0.5
-1.0 1.0
-0.5
-1.0
0.0
0.5
x
(b) E(x, 21 , 12 )
(a)E(x,y,1)
Figure 4. The absolute error of problem (3.20) with α1 = α2 = β1 = β2 = − 12 and N = M = K = 12.
23
1.0
Jacobi-Gauss-Lobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions A.H. Bhrawy, E.H. Doha, M.A. Abdelkawy, R.A. Van Gorder PII: DOI: Reference:
S0307-904X(15)00495-3 10.1016/j.apm.2015.09.009 APM 10693
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
19 April 2013 4 February 2015 4 September 2015
Please cite this article as: A.H. Bhrawy, E.H. Doha, M.A. Abdelkawy, R.A. Van Gorder, Jacobi-GaussLobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions, Applied Mathematical Modelling (2015), doi: 10.1016/j.apm.2015.09.009
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Jacobi-Gauss-Lobatto collocation method for solving nonlinear reaction-diffusion equations subject to Dirichlet boundary conditions A.H. Bhrawya,b , E.H. Dohac , M.A. Abdelkawyb , and R.A. Van Gorder
d
a
AN US
CR IP T
Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia (E-mail: [email protected]) b Department of Mathematics, Faculty of Science, Beni-Suef University, Beni-Suef, Egypt (E-mail:[email protected]) c Department of Mathematics, Faculty of Science, Cairo University, Giza, Egypt (E-mail: [email protected]) d Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford OX2 6GG United Kingdom (E-mail: [email protected])
Introduction
PT
1
ED
M
Abstract: This paper extends the application of the spectral Jacobi-Gauss-Lobatto collocation (J-GL-C) method based on Gauss-Lobatto nodes to obtain semi-analytical solutions of nonlinear time-dependent reaction-diffusion equations (RDEs) subject to Dirichlet boundary conditions. This approach has the advantage of allowing us to obtain the solution in terms of the Jacobi parameters α and β, which therefore means that the method holds a number of collocation methods as a special case. In addition, the problem is reduced to the solution of systems of ordinary differential equations (SODEs) in the time variable, which may then be solved by any standard numerical technique. We consider five applications of the general method to concrete examples. In each of the examples considered, the numerical results show that the proposed method is of high accuracy and is efficient for solving nonlinear time-dependent RDEs. Keywords: Nonlinear reaction-diffusion equations; Collocation method; Jacobi-Gauss-Lobatto quadrature.
AC
CE
RDEs have a wide range of applications in many fields of science, and have recently attracted considerable attention, partly due to the interesting features and rich variety of properties of their solutions (see Ref. [1] and the references therein). The processes of diffusion and reaction play an essential role in the dynamics of many systems, e.g., in plasma [2], or semiconductor physics [3]. Spectral methods (see, e.g., [4]-[5] and the references therein) are techniques used in applied mathematics and scientific computing to numerically solve both linear and nonlinear differential equations. The Galerkin, tau and collocation methods are three well-known types of spectral methods. We mention that the spectral collocation method is very useful in providing highly accurate solutions to nonlinear differential equations [5]-[8], also it has become increasingly popular for solving fractional differential equations [8]-[10]. Guo and Yan [7] introduced a new Legendre-Gauss collocation method for solving nonlinear second-order ordinary differential equations. Saadatmandi and Dehghan [11] developed the Sinc-collocation approach for solving multi-point boundary value problem. In this approach, the computation of solution of such problems is reduced to solving algebraic equations. Bhrawy and Alofi [6] proposed the spectral shifted Jacobi-Gauss collocation method to find an accurate solution of the Lane-Emden type equation. Moreover, Doha et al. [12]
1
ACCEPTED MANUSCRIPT
Development of the J-GL-C method
CE
2
PT
ED
M
AN US
CR IP T
developed the shifted Jacobi-Gauss collocation method for solving nonlinear high-order multi-point boundary value problems. For time-dependent partial differential equations (PDEs), spectral methods have been studied for several decades. Glenn et al. [13] investigated the spectral methods for numerically solving time-dependent class of parabolic partial differential equations subject to periodic boundary conditions. Tal-Ezer [14, 15] introduced spectral methods by using the polynomial approximation of the evolution operator in Chebyshev least-squares sense for time-dependent parabolic and hyperbolic equations respectively. Recently, Bhrawy and Al-shomrani [16] proposed a new spectral method for treating the initial-boundary value Korteweg-de Vries equation. The scheme consists of a Jacobi dual-Petrov Galerkin-Jacobi collocation method in space combined with a Crank-Nicholson-leapfrog method in time such that at each time step only a sparse banded linear algebraic system needs to be solved. Solutions given in terms of the Jacobi parameters α and β are possible to obtain, and the use of Jacobi polynomials for solving differential equations has gained increasing popularity in recent years (see, [17]-[20]). To the best of our knowledge, there are no results on the J-GL-C method for solving nonlinear RDEs. The main aim of this paper is to investigate the J-GL-C method for the solution to nonlinear time-dependent RDEs subject to Dirichlet boundary conditions. Indeed, it is very useful to carry out a systematic study on J-GL-C method with general parameters α, β. The nonlinear time-dependent RDE is collocated only for the space variable at (N − 1) points, for suitable collocation points. We use the (N − 1) nodes of the Jacobi-Gauss-Lobatto interpolation which depends upon the two general parameters α, β. These equations together with the two-point boundary conditions constitute system of (N − 1) ordinary differential equations (ODEs) in time. This system can be solved by one of the standard numerical methods. In addition, this algorithm is applied to numerically solve two-dimensional RDEs. Finally, the accuracy of the proposed method is demonstrated by test problems. The paper has two main parts. In Section 2, some properties of Jacobi polynomials are presented. With this mathematical background, we are then able to develop a way of constructing the Gauss-Lobatto collocation technique for a general form of nonlinear time-dependent RDEs, using the Jacobi polynomials. Then, in Section 3, the proposed method is applied to five specific and relevant RDEs. Concluding remarks are made in Section 4.
AC
We start by outlining some properties of Jacobi polynomials. Then we discuss the general class of nonlinear reaction diffusion equation we shall consider. We finally develop the collocation method for this general class of reaction diffusion equation.
2
ACCEPTED MANUSCRIPT
2.1
Basic properties of Jacobi polynomials (α,β)
The standard Jacobi polynomial of degree k, Pk (x) (k = 0, 1, · · · ), with the parameters α > −1, β > −1 [21, 22, 23], satisfies the following relations: (α,β)
Pk
(α,β)
(−x) = (−1)k Pk
(x) ,
CR IP T
(−1)k Γ(k + β + 1) , k!Γ(β + 1) Γ(k + α + 1) (α,β) Pk (1) = . k!Γ(α + 1) (α,β)
Pk
(−1) =
(2.1)
(u, v)w(α,β) =
Z1
AN US
Let w(α,β) (x) = (1 − x)α (1 + x)β . Then we define the weighted space L2w(α,β) as usual, equipped with the following inner product and norm u(x) v(x)w(α,β) (x)dx,
−1
1
kukw(α,β) = (u, u)w2 (α,β) .
(2.2)
The set of Jacobi polynomials forms a complete L2w(α,β) -orthogonal system, and (α,β)
kw(α,β) = hk =
2α+β+1 Γ(k + α + 1)Γ(k + β + 1) . (2k + α + β + 1)Γ(k + 1)Γ(k + α + β + 1)
(2.3)
M
kPk
Let SN (−1, 1) be the set of polynomials of degree at most N on (−1, 1), and due to the property of the standard Jacobi-Gauss quadrature, it follows that for any φ ∈ S2N +1 (−1, 1),
ED
Z1
w(α,β) (x)φ(x)dx =
(α,β)
(α,β)
$N,j φ(xN,j ),
(2.4)
j=0
−1
(α,β)
N X
(α,β)
CE
PT
where xN,j (0 ≤ j ≤ N ) and $N,j (0 ≤ j ≤ N ) are the nodes and the corresponding Christoffel numbers of the Jacobi-Gauss-quadrature formula on the interval (−1, 1), respectively. Now, we introduce the following discrete inner product and norm (u, v)w(α,β) =
N X
(α,β)
(α,β)
(α,β)
u(xN,j ) v(xN,j ) $N,j ,
j=0
1
kukw(α,β) = (u, u)w2 (α,β) .
(2.5)
AC
For α = β one recovers the ultraspherical polynomials (symmetric Jacobi polynomials) and for α = β = ∓ 12 or α = β = 0 one has the Chebyshev of the first and second kinds or Legendre polynomials, respectively. For the nonsymmetric Jacobi polynomials, the two important special cases α = −β = ± 12 give the Chebyshev polynomials of the third and fourth kinds are also recovered.
2.2
Nonlinear reaction-diffusion model
The simple RDE concerning the concentration u of a single substance in one spatial dimension, ut = auxx + C(u), 3
(2.6)
ACCEPTED MANUSCRIPT
CR IP T
is also referred to as the Kolmogorov-Petrovsky-Piskounov (KPP) equation [24]. If the reaction term C(u) vanishes, then the equation represents a pure diffusion process. The corresponding equation is Fick’s second law. Making the choice C(u) = u(1−u2 ) yields the Newell-Whitehead-Segel equation, which describes the Rayleigh-Benard convection [25, 26]. Choosing C(u) = u(1 − u)(u − α) and 0 < α < 1, we obtain the more general Zeldovich equation that arises in combustion theory [27, 28]. Such equations with cubic reaction functions are also known as Nagumo or FitzhughNagumo models; these models appear in a variety of applications, and can include more general terms which allow for density dependence of the diffusion process [29, 30, 31, 32, 33]. In this paper, we consider nonlinear RDEs with a reaction term which can be written in the form ut = (A(u) ux )x + B(u) ux + C(u),
(2.7)
AN US
where u = u(x, t) is an unknown function, and A(u), B(u), C(u) are arbitrary smooth functions. The indices x and t denote differentiation with respect to these variables. Eq. (2.7) is a generalized form of a great number of well-known nonlinear second-order evolution equations describing various processes in biology [34]-[36]. Indeed, the term A(u) permits a density dependent manner of diffusion. Burgers’ equation is a particular case of Eq. (2.7) ut = uxx + λ1 u ux .
(2.8)
ED
M
Burgers simplified the Navier-Stokes equations by dropping the pressure term to obtain his one dimensional equation. This equation has many applications in applied mathematics, such as modeling of gas dynamics [37, 38], modeling of fluid dynamics, turbulence, boundary layer behavior, shock wave formation and traffic flow [39]. The second particular case of Eq. (2.7) is the well-known Fisher equation ut = uxx + λ2 u − λ3 u2 ,
(2.9)
AC
CE
PT
which firstly was introduced by Fisher [40] in order to describe the propagation of a mutant gene. Fisher equation has a wide application in a large number of the fields. It describes the spreading of biological populations [40], chemical kinetics [41], logistic population growth [42], flame propagation [43], population in one-dimensional habitual [44], neutron population in a nuclear reaction [45], neurophysiology [46], autocatalytic chemical reactions [47], branching Brownian motion processes [48], and nuclear reactor theory [49]. We also notice that the Murray equation [34, 35] represents a particular case of Eq. (2.7) ut = uxx + λ1 u ux + λ2 u − λ3 u2 ,
(2.10)
which can be considered as a generalization of the Fisher and Burgers equations. Due to the varied applications of nonlinear RDEs in science and engineering, solution methods for such equations are of great importance. We shall next develop the J-GL-C method for solving such equations under very general assumptions.
2.3
Jacobi spectral collocation method
Since collocation methods give good approximations to differential equations and integral equations in physical space, it is useful to implement and to adapt the methods to various problems, including 4
ACCEPTED MANUSCRIPT
variable coefficient and nonlinear PDEs (see, for instance [6, 8, 50]). In this section, without loss of generality, we derive the J-GL-C method to numerically solve the following general class of nonlinear time-dependent RDEs:
on the domain D = {x : −1 < x < 1} u(−1, t) = g1 (t),
(2.11)
u(1, t) = g2 (t)
(2.12)
AN US
subject to the Dirichlet boundary conditions
CR IP T
ut = A(u) uxx + B(u) ux + C(u)
(2.13)
and the initial condition u(x, 0) = f1 (x),
x ∈ D.
Here we assume A, B and C are sufficiently smooth functions of u. Now, assume that u(x, t) =
N X
(α,β)
aj (t)Pj
(x),
(2.14)
M
j=0
ED
and if we make use of the orthogonality condition of Jacobi polynomials (2.5), then we obtain form (2.14) that aj (t) =
N 1 X (α,β) (α,β) (α,β) (α,β) Pj (xN,i )$N,i u(xN,i , t), hj
(2.15)
i=0
PT
and accordingly, Eq. (2.14) takes the form
CE
u(x, t) =
N N X 1 X (α,β) (α,β) (α,β) (α,β) (α,β) Pj (xN,i )Pj (x)$N,i u(xN,i , t) , hj
(2.16)
N X N X 1 (α,β) (α,β) (α,β) (α,β) (α,β) Pj (xN,i )Pj (x)$N,i u(xN,i , t). hj
(2.17)
j=0
i=0
AC
or equivalently
u(x, t) =
i=0
j=0
The spatial partial derivatives in Eq. (2.11) can be computed at the collocation points to give (α,β)
ux (xN,n , t) =
N X N X 1 (α,β) (α,β) (α,β) (α,β) 0 (α,β) (α,β) Pj (xN,i )(Pj (xN,n )) $N,i u(xN,i , t) hj i=0
=
N X i=0
j=0
(α,β)
Ani u(xN,i , t),
n = 0, 1, · · · , N, 5
(2.18)
ACCEPTED MANUSCRIPT
N X N X 1 (α,β) (α,β) (α,β) (α,β) 00 (α,β) (α,β) Pj (xN,i )(Pj (xN,n )) $N,i u(xN,i , t) hj i=0
=
N X i=0
where
j=0
(α,β)
N X 1 (α,β) (α,β) (α,β) (α,β) 0 (α,β) P (xN,i )(Pj (xN,n )) $N,i , hj j
(2.20)
N X 1 (α,β) (α,β) (α,β) (α,β) 00 (α,β) = Pj (xN,i )(Pj (xN,n )) $N,i . hj
(2.21)
Ani =
j=0
Bni
(2.19)
Bni u(xN,i , t), n = 0, 1, · · · , N,
CR IP T
(α,β)
uxx (xN,n , t) =
j=0
N X
u˙ n (t) = A(un (t))
Bni ui (t) + B(un (t))
i=0
AN US
Making use of Eqs. (2.18)-(2.21), we rewrite Eq. (2.11) in the form: N X
Ani ui (t) + C(un (t)),
i=0
where
(α,β)
ui (t) = u(xN,i , t),
n = 1, · · · , N − 1,
(2.22)
i = 1, · · · , N − 1.
i=1
n = 1, · · · , N − 1,
CE
PT
where
N −1 X
fn (t)) + B(un (t)) ( Bni ui (t) + d
ED
N −1 X
u˙ n (t) = A(un (t)) (
M
Equation (2.22) with the two-point boundary conditions (2.12), constitutes a system of (N − 1) ODEs in the expansion coefficients aj (t), namely Ani ui (t) + dn (t)) + C(un (t)),
i=1
(2.23)
dn (t) = An0 g1 (t) + AnN g2 (t), fn (t) = Bn0 g1 (t) + BnN g2 (t). d
AC
This means that the problem (2.11)-(2.13) is transformed into the following system of ordinary differential equations: N −1 X
u˙ n (t) = A(un (t)) (
i=1
(α,β)
un (0) = f1 (xN,n ),
N −1 X
fn (t)) + B(un (t)) ( Bni ui (t) + d
n = 1, · · · , N − 1,
Ani ui (t) + dn (t)) + C(un (t)),
i=1
(2.24)
which can be written in the matrix form ˙ u(t) = F(t, u(t)) u(0) = f1 , where u(t) ˙ = [u˙ 1 (t), u˙ 2 (t), . . . , u˙ N −1 (t)]T , 6
(2.25)
ACCEPTED MANUSCRIPT
(α,β)
(α,β)
(α,β)
f1 = [f1 (xN,1 ), f1 (xN,2 ), . . . , f1 (xN,N −1 )]T , and F(t, u(t)) = [F1 (t, u(t)), F2 (t, u(t)), . . . , FN −1 (t, u(t))]T , (α,β)
Fn (t, u(t)) = f (t, xN,n , un (t),
N −1 X
Ani ui (t) + dn (t),
i=1
N −1 X
fn (t)). Bni ui (t) + d
CR IP T
where
i=1
The System of ODEs (2.25) can then be solved by using standard numerical solvers.
2.4
Two-dimensional RDEs
Now, we extend the above algorithm to solve the following two-dimensional RDEs: (x, y, t) ∈ Ω1 × Ω2 × Ω3 ,
AN US
ut = a1 uxx +b1 uyy + C(u), where Ω1 = [−1, 1],
Ω2 = [−1, 1]
subject to the initial condition
Ω3 = [0, T ],
(x, y) ∈ Ω1 × Ω2 ,
M
u(x, y, 0) = g1 (x, y),
and
and boundary conditions
u(1, y, t) = g3 (y, t),
u(x, −1, t) = g4 (x, t),
u(x, 1, t) = g5 (x, t)
ED
u(−1, y, t) = g2 (y, t),
(2.26)
(y, t) ∈ Ω2 × Ω3 ,
(x, t) ∈ Ω1 × Ω3 .
(2.27)
(2.28)
PT
The J-GL-C method will be extended to transform the previous PDE into SODEs. To this end, let us expand the dependent variable in the form
CE
u(x, y, t) =
N X M X
(α1 ,β1 )
ai,j (t)Pj
(α2 ,β2 )
(x)Pj
(y).
(2.29)
i=0 j=0
Employing Eqs. (2.14) and (2.5), the coefficients ai,j (t) can be approximated by N M 1 X X (α2 ,β2 ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) (α2 ,β2 ) Pj (yM,k )$M,k Pj (xN,l )$N,l u(xN,l1 1 , yM,k , t). hi hj
AC
ai,j (t) =
l=0 k=0
(2.30)
The expression (2.29) can be rewritten as u(x, y, t) =
N X M X
n,m τl,k ul,k (t),
l=0 k=0
where
(α ,β1 )
u(xN,l1
(α ,β2 )
2 , yM,k
7
, t) = ul,k (t),
(2.31)
ACCEPTED MANUSCRIPT
n,m τl,k =
N X M X i=0 j=0
(α2 ,β2 )
Pj
(α ,β2 )
2 (yM,k
(α ,β2 )
2 )$M,k
(α1 ,β1 )
Pi
(α ,β1 )
(xN,l1
(α ,β1 )
)$N,l1
hi hj
(α1 ,β1 )
Pi
(α2 ,β2 )
(x)Pj
(y).
ux (x, y, t) =
N X M X
CR IP T
In what follows, the approximation of the first-order spatial partial derivative with respect to x for u(x, y, t) can be computed at the Jacobi-Gauss-Lobatto interpolation nodes to give
n,m γl,k ul,k (t),
l=0 k=0
where n,m = γl,k
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X hi hj
0
(α2 ,β2 )
(x)) Pj
(y) ul,k (t),
AN US
l=0 k=0
(α1 ,β1 )
(Pi
(2.32)
(2.33)
also, the first-order spatial partial derivative with respect to y is given by uy (x, y, t) =
N X M X
n,m δl,k ul,k (t),
(2.34)
l=0 k=0
n,m = δl,k
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
M
where
hi hj
ED
l=0 k=0
(α1 ,β1 )
Pi
(α2 ,β2 )
(x)(Pj
0
(y)) ul,k (t). (2.35)
PT
In the same manner, the second spatial partial derivatives with respect to x and y are
n,m ζl,k =
n,m ηl,k =
(2.36)
n,m ηl,k ul,k (t),
(2.37)
l=0 k=0
uyy (x, y, t) =
N X M X l=0 k=0
(Pi
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
Pi
l=0 k=0
and
n,m ζl,k ul,k (t),
(α ,β ) (α2 ,β2 ) (α2 ,β2 ) (α1 ,β1 ) (α1 ,β1 ) (α ,β ) N X M Pj 2 2 (yM,k )$M,k Pi (xN,l )$N,l1 1 X
AC
where
CE
uxx (x, y, t) =
N X M X
l=0 k=0
hi hj
hi hj
8
(α1 ,β1 )
00
(α2 ,β2 )
(x)) Pj
(y) ul,k (t), (2.38)
(α1 ,β1 )
(α2 ,β2 )
(x)(Pj
00
(y)) ul,k (t). (2.39)
ACCEPTED MANUSCRIPT
In the proposed J-GL-C method the residual of (2.26) is set to be zero at (N − 1) × (M − 1) of Jacobi-Gauss-Lobatto interpolation points. Moreover, the values of u0,k (t), uN,k (t), ul,0 (t) and ul,N (t) can be given by
(α ,β1 )
ul,0 (t) = g4 (xN,l1
(α ,β2 )
, t),
2 uN,k (t) = g3 (yM,k
, t),
ul,N (t) = g5 (xN,l1
(α ,β1 )
, t),
k = 0, · · · , M,
, t),
l = 0, · · · , N.
CR IP T
(α ,β2 )
2 u0,k (t) = g2 (yM,k
(2.40)
Therefore, adapting (2.31)-(2.40), enable one to write (2.26)-(2.28) as (N − 1) × (M − 1) SODEs u˙ n,m (t) =a1
N X M X
n,m ζl,k ul,k (t) + b1
subject to the initial data (α ,β1 )
m = 1, · · · , M − 1,
AN US
n = 1, · · · , N − 1, 1 g1n,m = g1 (xN,n
n,m ηl,k ul,k (t) + C(un,m (t)),
l=0 k=0
l=0 k=0
un,m (0) = g1n,m ,
N X M X
(α ,β2 )
2 , yM,k
(α ,β )
2 2 )$M,m ) n = 1, · · · , N − 1,
u1,1 (0)
... u1,M −1 (0) .. u2,1 (0) . u2,M −1 (0) .. . ··· ··· .. . ··· ··· .. . ··· ··· .. uN −2,1 (0) . uN −2,M −1 (0) uN −1,1 (0) . . . uN −1,M −1 (0)
AC
CE
PT
ED
M
Finally, we can arrange (2.41)-(2.42) to their matrix formulation, Runge-Kutta scheme, H1,1 (t, u1 , · · · , uN −1 ) u˙ 1,1 (t) ... u˙ 1,M −1 (t) . .. u˙ 2,1 (t) u˙ 2,M −1 (t) H2,1 (t, u1 , · · · , uN −1 ) .. . ··· · · · · · · . . = . ··· ··· ··· . .. ··· ··· ··· .. u˙ N −2,1 (t) . u˙ N −2,M −1 (t) HN −2,1 (t, u1 , · · · , uN −1 ) HN −1,1 (t, u1 , · · · , uN −1 ) u˙ N −1,1 (t) . . . u˙ N −1,M −1 (t)
where
Hn,m (t, u1 , · · · , uN −1 ) =a1
N X M X
g 2,1 1 ··· = ··· ··· N −2,1 g1 g1N −1,1
n,m ζl,k ul,k (t) + b1
l=0 k=0
n = 1, · · · , N − 1,
g11,1
N X M X
m = 1, · · · , M − 1. (2.42)
which can be solved by implicit
... H1,M −1 (t, u1 , · · · , uN −1 ) .. . H2,M −1 (t, u1 , · · · , uN −1 ) .. . ··· .. . ··· .. . ··· .. . HN −2,M −1 (t, u1 , · · · , uN −1 ) . . . HN −1,M −1 (t, u1 , · · · , uN −1 ) (2.43) 1,M −1 ... g1 .. . g12,M −1 .. . ··· .. (2.44) , . ··· .. . ··· .. . g1N −2,M −1 . . . g1N −1,M −1
n,m ηl,k ul,k (t) + C(un,m (t)),
l=0 k=0
m = 1, · · · , M − 1. 9
(2.41)
(2.45)
,
ACCEPTED MANUSCRIPT
3
Numerical results
3.1
Murray equation
CR IP T
To illustrate the effectiveness of the proposed method, four test examples are considered in this section. Comparison of the results obtained by various choices of Jacobi parameters α and β reveal that the present method is very accurate and efficient for all choices of α and β.
Consider the nonlinear time-dependent one dimensional Murray equation ut = uxx + λ1 uux + λ2 u − λ3 u2 ,
(x, t) ∈ D × [0, 1],
AN US
subject to the boundary conditions λ2 λ2 2 2 1 + tanh (2 λ1 λ3 + (λ1 + 4 λ3 )t) , u(1, t) = 2λ3 8 λ23 λ2 λ2 2 2 u(−1, t) = 1 + tanh (−2 λ λ + (λ + 4 λ )t) , 1 3 1 3 2λ3 8 λ23 and the initial condition
λ2 (2 λ1 λ3 x) , 1 + tanh 8 λ23
x ∈ D.
(3.2)
(3.3)
M
λ2 u(x, 0) = 2λ3
(3.1)
ED
If we apply the generalized tanh method [51], then we find that the exact solution of Eq. (3.1) is λ2 λ2 2 2 u(x, t) = 1 + tanh (2 λ1 λ3 x + (λ1 + 4 λ3 )t) . (3.4) 2λ3 8 λ23
PT
The difference between the measured or inferred value of approximate solution and its actual value (absolute error) is given by E(x, t) = |u(x, t) − u e(x, t)|,
(3.5)
ME = max{E(x, t) : ∀(x, t) ∈ D × [0, 1]},
(3.6)
AC
CE
where u(x, t) and u e(x, t) are the exact solution and the approximate solution at the point (x, t), respectively. Moreover, the maximum absolute error is given by while the rate of convergence may be given by RE =
max{E(x, t) : ∀(x, t) ∈ D × [0, T ]}|N =Ni . max{E(x, t) : ∀(x, t) ∈ D × [0, T ]}|N =Ni+1
(3.7)
Maximum absolute errors and rates of convergence of (3.1) subject to (3.2) and (3.3) are introduced in Table 1 using the J-GL-C method for the two special values α = β = 0 and α = β = − 21 , respectively. We also fix λ1 = λ2 = λ3 = 1 with different choices of N . Meanwhile, the absolute errors of problem (3.1) are presented in Table 2 for α = β = − 12 , λ1 = λ2 = λ3 = 1 and N = 16 with different values of (x, t). 10
ACCEPTED MANUSCRIPT
Table 1 of convergence RE of solutions for equation (3.1) ME RE 4.00 × 10−5 8.00 × 10−8 2.00 × 10−3 2.00 × 10−8 2.50 × 10−1 1.00 × 10−8 5.00 × 10−1 8.00 × 10−6 8.00 × 10−8 1.00 × 10−2 6.00 × 10−8 7.50 × 10−1 1.50 × 10−8 2.50 × 10−1
CR IP T
Maximum absolute errors ME and rate N α β 4 8 0 0 12 16 4 8 − 12 − 21 12 16
A nonlinear RDE
M
3.2
AN US
In Figure 1, we display, in (a) the absolute error (b) the approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 8. Moreover in (c) the absolute error, (d) the approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 12 are plotted. From Figs. 1(c, d), we see that the curves of the approximate and exact solutions coincide for different values of t. This shows that the obtained numerical results are accurate and compare favorably with the exact solution.
Consider the nonlinear time-dependent RDE (a type of Nagumo equation) as given below
ED
ut = λ1 uxx + u(λ2 + λ3 u − λ4 u2 ),
(x, t) ∈ D × [0, 1],
CE
PT
together with the boundary conditions h i √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (1 + λ√3 2λλ1 t) 1 4 4 u(1, t) = , 2λ4 h i √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (−1 + λ√3 2λλ1 t) 1 4 4 u(−1, t) = , 2λ4
AC
and the initial condition
u(x, 0) =
λ3 +
h i p 4λ λ +λ2 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (x) 1 4
2λ4
,
x ∈ D.
(3.8)
(3.9)
(3.10)
Applying the generalized tanh method [51], we find that the exact solution of Eq. (3.8) is given by i h √ p 4λ λ +λ2 λ3 + 4λ2 λ4 + λ23 tanh 2√2 2λ4 λ 3 (x + λ√3 2λλ1 t) 1 4 4 (3.11) u(x, t) = . 2λ4 Maximum absolute errors and rates of convergence of (3.8) subject to (3.9) and (3.10) are presented in Table 3 using J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, and λ1 = λ2 = λ3 = λ4 = 1 with different choice of N . 11
ACCEPTED MANUSCRIPT
0
CR IP T
0.2
AN US
0
ED
0.1
M
Absolute t α
CE
PT
x - 0.5 -0.4 - 0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 - 0.5 -0.4 - 0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Table 2 errors of solutions for equation (3.1) β E α β E 5.05 × 10−9 3.18 × 10−9 4.93 × 10−9 2.72 × 10−9 −9 4.66 × 10 2.14 × 10−9 −9 4.28 × 10 1.53 × 10−9 −9 3.82 × 10 8.82 × 10−10 0 3.26 × 10−9 − 21 − 12 1.80 × 10−10 2.63 × 10−9 5.44 × 10−10 −9 1.95 × 10 1.25 × 10−9 1.21 × 10−9 1.94 × 10−9 −10 4.40 × 10 2.61 × 10−9 −10 3.18 × 10 3.18 × 10−9 7.31 × 10−10 3.32 × 10−9 −10 5.61 × 10 3.32 × 10−9 −10 3.52 × 10 3.18 × 10−9 −10 1.96 × 10 3.07 × 10−9 7.71 × 10−11 2.93 × 10−9 1 1 −11 0 3.60 × 10 − 2 − 2 2.68 × 10−9 −10 1.11 × 10 2.40 × 10−9 1.35 × 10−10 2.16 × 10−9 −10 1.65 × 10 1.87 × 10−9 −10 2.19 × 10 1.46 × 10−9 2.07 × 10−10 1.10 × 10−9
AC
Table 3 Maximum absolute errors and rate of convergence N α β ME 4 6.00 × 10−3 8 0 0 3.00 × 10−5 12 1.00 × 10−7 16 4.00 × 10−9 4 4.00 × 10−3 1 1 8 − 2 − 2 1.00 × 10−5 12 4.00 × 10−8 16 1.00 × 10−8
12
of the solution to equation (3.8) RE 5.00 × 10−3 3.33 × 10−3 4.00 × 10−2 2.50 × 10−3 4.00 × 10−3 2.50 × 10−1
ACCEPTED MANUSCRIPT
3.3
errors and rate of p ME 1.50 × 10−4 3 1.00 × 10−8 8.00 × 10−13 6.00 × 10−5 3 1.50 × 10−11 1.50 × 10−13
Table 4 convergence for the solution of RE p ME 3.00 × 10−3 6.67 × 10−5 5 6.00 × 10−5 8.00 × 10−5 4.00 × 10−7 2.00 × 10−2 −6 2.50 × 10 5 3.00 × 10−6 1.00 × 10−2 1.50 × 10−8
equation (3.12) RE 2.00 × 10−2 6.67 × 10−3
CR IP T
Maximum absolute N α β 4 8 0 0 12 4 8 − 21 − 12 12
1.50 × 10−4 5.00 × 10−3
Reaction-diffusion equation with power-law density dependent diffusion
AN US
In the following we analyse the nonlinear time-dependent one dimensional RDE which includes density dependent diffusion ut = (λ1 up−1 ux )x − λ2 u + λ3 up , The boundary conditions are given by
and subjected to the initial condition 2bp q(p + 1)
1 p−1
ED
u(x, 0) =
2bp q(p + 1)
1 p−1
λ3 > 0,
cos
p−1 2
(x, t) ∈ D × [0, 1] .
r
q ap
2 p−1
(3.12)
,
(3.13)
x ∈ D.
(3.14)
M
u(−1, t) = u(1, t) =
p > 1,
2 r p−1 p−1 q cos x , 2 ap
The exact solution [1] of Eq. (3.12) has the form
PT
u(x, t) =
2bp q(p + 1)
1 p−1
2 r p−1 p−1 q cos x . 2 ap
(3.15)
AC
CE
Maximum absolute errors and rates of convergence of (3.12) subject to (3.13) and (3.14) are depicted in Table 4 using the J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, and λ1 = λ2 = λ3 = 1 with various choice of N and p. In Figure 2, we present the absolute error of problem (3.12) where λ1 = λ2 = λ3 = 1, p = 5 α = β = −0.5 and N = 4, 8 and 12, respectively.
3.4
Diffusion- areaction equation with Fischer-Kolmogorov reaction term and density dependent diffusion
Finally, we investigate the nonlinear time-dependent one dimensional RDE presented below, which includes a Fischer-Kolmogorov reaction term and density dependent diffusion: ut = (u2 )xx + u(1 − u),
13
(x, t) ∈ D × [0, 1],
(3.16)
ACCEPTED MANUSCRIPT
for the solution of equation (3.16) RE 6.67 × 10−6 3.00 × 10−1 5.00 × 10−4 3.00 × 10−1
with to the boundary conditions
and the initial condition
2 , 1 + tanh −1+t 4 u(x, 0) =
u(1, t) =
2 , 1 + tanh 1+t 4
AN US
u(−1, t) =
CR IP T
Table 5 and rate of convergence α β ME 3.00 × 10−3 0 0 2.00 × 10−8 6.00 × 10−9 4.00 × 10−5 1 1 − 2 − 2 2.00 × 10−8 6.00 × 10−9
Maximum absolute errors N 4 8 12 4 8 12
2 , 1 + tanh x4
x ∈ D.
(3.17)
(3.18)
We notice that the exact solution [44] of Eq. (3.16) is
2 . 1 + tanh x+t 4
(3.19)
M
u(x, t) =
Two-dimensional nonlinear RDE
CE
3.5
PT
ED
Maximum absolute errors and rates of convergence of (3.16) subject to (3.17) and (3.18) can be seen in Table 5, where we have used the J-GL-C method for the two special values α = β = 0 and α = β = − 12 , respectively, with different values of N . On the other hand, the absolute errors of problem (3.16) are exhibited in Table 6 for α = β = − 12 and N = 12 with different values of (x, t). In Figure 3, we plot both the approximate and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.16) where N = 12. We remark that the curves of the approximate and exact solutions agree strongly for different values of t which are listed in their captions. This shows that the obtained numerical results are accurate and compare favorably with the exact solution.
AC
In order to confirm the high accuracy of our technique for the two-dimensional problem, we consider the 2+1 nonlinear RDE 1 ∂t u(x, y, t) = ∇2 u(x, y, t) + (u(x, y, t))2 (1 − u(x, y, t)), (x, y, t) ∈ [−1, 1] × [−1, 1] × [0, 20], 2 (3.20)
with the following initial-boundary conditions 1
u(−1, y, t) = 1+e u(x, −1, t) =
√1 (y−1− √1 t) 2 2
1 1+e
√1 (x−1− √1 t) 2 2
1+e 1+e
1+e
1 √ (x+y) 2
14
.
1 √ (y+1− √1 t) 2 2
1
, u(x, 1, t) = 1
u(x, y, 0) =
1
, u(1, y, t) =
1 √1 (x+1− √ t) 2 2
, ,
(3.21)
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Table 6 Absolute errors for the solution of equation (3.16) x t α β E α β E - 0.5 3.19 × 10−9 3.31 × 10−9 -0.4 3.50 × 10−9 3.66 × 10−9 −9 - 0.3 3.73 × 10 3.90 × 10−9 −9 -0.2 3.88 × 10 4.06 × 10−9 -0.1 3.96 × 10−9 4.14 × 10−9 1 1 −9 0 0.1 0 0 3.96 × 10 − 2 − 2 4.15 × 10−9 −9 0.1 3.90 × 10 4.09 × 10−9 −9 0.2 3.76 × 10 3.95 × 10−9 0.3 3.56 × 10−9 3.73 × 10−9 −10 0.4 3.29 × 10 3.46 × 10−9 −10 0.5 2.96 × 10 3.11 × 10−9 - 0.5 1.50 × 10−9 1.90 × 10−9 −9 -0.4 1.41 × 10 1.79 × 10−9 −9 - 0.3 1.23 × 10 1.57 × 10−9 -0.2 1.01 × 10−9 1.32 × 10−9 −10 -0.1 7.92 × 10 1.06 × 10−9 1 1 −10 0 0.2 0 0 5.96 × 10 − 2 − 2 8.35 × 10−10 0.1 4.32 × 10−10 6.46 × 10−10 −10 0.2 3.11 × 10 5.06 × 10−10 −10 0.3 2.44 × 10 4.21 × 10−10 −10 0.4 2.29 × 10 3.93 × 10−10 0.5 2.48 × 10−10 4.05 × 10−10
15
ACCEPTED MANUSCRIPT
The exact solution of this problem is 1
u(x, y, t) =
.
(3.22)
AN US
1+e
1 √ (x+y− √1 t) 2 2
ME 1.44 × 10−3 2.51 × 10−7 1.11 × 10−8 314 × 10−4 5.14 × 10−8 2.35 × 10−9
CR IP T
N =M =K 4 8 12 4 8 12
Table 7 Maximum absolute errors of problem (3.20) (α1 , β1 ) (α2 , β2 ) ME (α1 , β1 ) (α2 , β2 ) ( 21 , 21 ) (− 12 , 12 ) 1.14 × 10−3 (0, 12 ) (− 12 , 0) −7 3.49 × 10 2.89 × 10−9 1 1 (1, 0) (− 2 , − 2 ) 1.38 × 10−3 (− 12 , − 12 ) (− 12 , − 12 ) 4.19 × 10−7 3.67 × 10−9
Table 7 lists the maximum absolute errors of u(x, y, t) of problem (3.20) with various choices of α1 , α2 , β1 , β2 , , N, M, and K. The numerical results presented in this table show that the results are very accurate for small value of N, M, and K. Figure 4 demonstrates that the absolute errors E(x, y, t), at t = 1 and E(x, y, t), at y = t = 12 are very small for even the small number of grid points taken for all the Jacobi parameters.
Conclusions
M
4
AC
CE
PT
ED
We have investigated the J-GL-C method for the solution of nonlinear time-dependent RDEs subject to Dirichlet boundary conditions. The nonlinear time-dependent RDE has been collocated for the space variable at (N − 1) points, for suitable collocation points. We use the (N − 1) nodes of the Jacobi-Gauss-Lobatto interpolation which depends upon two general parameters α, β > −1. These equations, together with the two-point boundary conditions, yield a system of (N − 1) ODEs in the time variable. The resulting system can then be solved by one of the standard numerical methods. The solution method for a very general class of RDEs was considered in Section 2. This general class of equations holds a number of relevant physical and biological models as special cases. In Section 3, we have considered four specific nonlinear RDEs, in order to demonstrate the method. For each of the examples given, we were able to obtain accurate solutions with relatively few nodes used. Indeed, we often obtain very accurate results with fewer than ten nodes. As such, the method is quite efficient, requiring minimal computation to achieve a desired accuracy. Another facet of the solution method is that it permits solutions which appear to converge rapidly, as measured by the parameter RE , the rate of convergence. We find, in some cases, that the error decreases by as much as a factor of 10−3 due to an increase of only a single node. In the present paper, we have considered examples with power-law or polynomial type nonlinearity. However, the J-GL-C method is rather general, and one may consider other kinds of nonlinearity. In some applications [52], exponential nonlinearity occurs in reaction diffusion models, and such cases could also be addressed by the present method. Additionally, one sometimes encounters coupled RDEs. If the terms coupling the equations is strongly nonlinear, the solution of such equations can be rather difficult. Hence, for future work, it may be interesting to consider coupled RDEs with various forms of nonlinearity. 16
ACCEPTED MANUSCRIPT
CR IP T
The results presented here demonstrate that the J-GL-C method is both accurate and efficient for the computation of numerical solutions to nonlinear RDEs, even those with density-dependent diffusion terms and strongly nonlinear reaction functions. Therefore, the method warrants consideration from those involved in the numerical simulation of models arising in science or engineering which take the form of nonlinear RDEs.
References
[1] A.H. Khater, W. Malfliet, D.K. Callebaut, E.S. Kamel, The tanh method, a simple transformation and exact analytical solutions for nonlinear reaction-diffusion equations, Chaos Soliton. Fract., 14 (2002) 513-522.
AN US
[2] K. Mayawala, D. G. Vlachos, J. S. Edwards, Spatial modeling of dimerization reaction dynamics in the plasma membrane: Monte Carlo vs. continuum differential equations, Biophysical Chemistry, 121(3) (2006) 194-208. [3] P. A Markowich, P. Szmolyan, A system of convection-diffusion equations with small diffusion coefficient arising in semiconductor physics, J. Diff. Eqs., 81(2) (1989) 234-254. [4] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods: Fundamentals in Single Domains. Springer-Verlag, New York (2006).
M
[5] S. Nguyen and C. Delcarte, A spectral collocation method to solve Helmholtz problems with boundary conditions involving mixed tangential and normal derivatives, J. Comput. Phys., 200 (2004) 34-49.
ED
[6] A.H. Bhrawy and A.S. Alofi, A Jacobi-Gauss collocation method for solving nonlinear LaneEmden type equations, Commun. Nonlinear. Sci. Numer. Simulat., 17 (2012) 62-70.
PT
[7] B.-y. Guo and J.-p. Yan, Legendre-Gauss collocation method for initial value problems of second order ordinary differential equations, Appl. Numer. Math., 59 (2009) 1386-1408.
CE
[8] A.H. Bhrawy and Alghamdi, A shifted Jacobi-Gauss-Lobatto collocation method for solving nonlinear fractional Langevin equation involving two fractional orders in different intervals, Bound. Val. Prob., 2012 (2012) 62.
AC
[9] L. Zhu and Q. Fan, Solving fractional nonlinear Fredholm integro-differential equations by the second kind Chebyshev wavelet, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 2333-2341
[10] E.H. Doha, A.H. Bhrawy and S.S. Ezz-Eldien, A new Jacobi operational matrix: an application for solving fractional differential equation, Appl. Math. Model., 36 (2012) 4931-4943.
[11] A. Saadatmandi and M. Dehghan, The use of Sinc-collocation method for solving multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 593-601. [12] E.H. Doha, A.H. Bhrawy and R.M. Hafez, On shifted Jacobi spectral method for high-order multi-point boundary value problems, Commun. Nonlinear Sci. Numer. Simulat., 17 (2012) 3802-3810. 17
ACCEPTED MANUSCRIPT
[13] I. Glenn, S. Brian, W. Rodney, Spectral methods in time for a class of parabolic partial differential equations, J. Comput. Phys., 102 (1992) 88-97. [14] H. Tal-Ezer, Spectral methods in time for hyperbolic problems, SIAM J. Numer. Anal., 23 (1986) 11-26.
CR IP T
[15] H. Tal-Ezer, Spectral methods in time for parabolic problems, SIAM J. Numer. Anal., 26 (1989) 1-11. [16] A.H. Bhrawy, M. Al-shomrani, A Jacobi dual-Petrov Galerkin-Jacobi collocation method for solving Korteweg-de Vries equations, Abstract and Applied Analysis, 2012, (2012) Article ID 418943, 16
AN US
[17] E.H. Doha, A.H. Bhrawy, W.M. Abd-Elhameed, Jacobi spectral Galerkin method for elliptic Neumann problems, Numerical Algorithms, 50 (2009), 67-91 [18] E.H. Doha and A.H. Bhrawy, A Jacobi spectral Galerkin method for the integrated forms of fourth-order elliptic differential equations, Numerical Methods for Partial Differential Equations, 25 (2009) 712-739 [19] M. El-Kady, Jacobi discrete approximation for solving optimal control problems, J. Korean Math. Soc. 49 (2012) 99-112.
M
[20] E.H. Doha, A.H. Bhrawy and R.M. Hafez, A Jacobi-Jacobi dual-Petrov-Galerkin method for third- and fifth-order differential equations, Math. Comput. Model., 53 (2011) 1820-1832.
ED
[21] R. Beals, R. Wong, Special Functions, Cambridge University Press (2010). [22] G. Szeg¨ o, Orthogonal Polynomials, Colloquium Publications, XXIII, American Mathematical Society, ISBN 978-0-8218-1023-1, MR 0372517G. (1939).
PT
[23] Y. Luke, The special functions and their approximations, Academic Press, New York, 1969.
CE
[24] A. Kolmogorov, I. Petrovskii and N. Piscounov, A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem, In V. M. Tikhomirov, editor, Selected Works of A. N. Kolmogorov I, 248-270. Kluwer 1991. Translated by V. M. Volosov from Bull. Moscow Univ., Math. Mech., 1 (1937) 1-25.
AC
[25] A. C. Newell and J. A. Whitehead, Finite bandwidth, finite amplitude convection. J. Fluid Mech., 38 (1969) 279-303.
[26] L. A. Segel, Distant sidewalls cause slow amplitude modulation of cellular convection, J. Fluid Mech., 38 (1969) 203-224. [27] Y. B. Zeldovich and D. A. Frank-Kamenetsky, A theory of thermal propagation of flame, Acta Physicochim, 9 (1938) 341-350.
[28] B. H. Gilding and R. Kersner, Travelling Waves in Nonlinear Diffusion Convection Reaction, Nonlinear Differential Equations and their Applications, 60. Birkh¨auser Verlag, Basel (2004) 209 pp. 18
ACCEPTED MANUSCRIPT
[29] R. A. Van Gorder and K. Vajravelu, Analytical and numerical solutions of the density dependent diffusion Nagumo equation, Physics Letters A 372 (2008) 5152-5158.
CR IP T
[30] R. A. Van Gorder and K. Vajravelu, A variational formulation of the Nagumo reactiondiffusion equation and the Nagumo telegraph equation, Nonlinear Analysis Series B: Real World Applications 11 (2010) 2957-2962. [31] R. A. Van Gorder and K. Vajravelu, Analytical and Numerical Solutions for a Density Dependent Nagumo Telegraph Equation, Nonlinear Analysis Series B: Real World Applications 11 (2010) 3923-3929. [32] R. A. Van Gorder, Gaussian waves in the Fitzhugh-Nagumo equation demonstrate one role of the auxiliary function H(x) in the homotopy analysis method, Communications in Nonlinear Science and Numerical Simulation 17 (2012) 1233-1240.
AN US
[33] R. A. Van Gorder, Travelling waves for a density dependent diffusion Nagumo equation over the real line, Communications in Theoretical Physics 58 (2012) 5-11. [34] J.D. Murray, Nonlinear Differential Equation Models in Biology, Clarendon Press, Oxford, 1977. [35] J.D. Murray, Mathematical Biology, Springer, Berlin, 1989.
M
[36] P.C. Fife, Mathematical Aspects of Reacting and Diffusing Systems, Springer, Berlin, 1979.
ED
[37] A. Veksler and Y. Zarmi, Wave interactions and the analysis of the perturbed Burgers equation, Physica D, 211 (2005) 57-73. [38] A. Veksler, Y. Zarmi, Freedom in the expansion and obstacles to integrability in multiplesoliton solutions of the perturbed KdV equation, Physica D, 217 (2006) 77-87.
PT
[39] Z.-Y. Ma, X.-F. Wu and J.-M. Zhu, Multisoliton excitations for the Kadomtsev-Petviashvili equation and the coupled Burgers equation, Chaos Soliton. Fract., 31 (2007) 648-657.
CE
[40] R. A. Fisher, The wave of advance of advantageous genes. Ann. Eugen., 7 (1937) 353-369. [41] A. J. Khattak, A computational meshless method for the generalized Burgers’-Huxley equation, Appl. Math. Modell., 33 (2009) 3718-3729.
AC
[42] N. F. Britton, Reaction-diffusion equations and their applications to biology. London: Academic Press; 1986. [43] D. A. Frank, Diffusion and heat exchange in chemical kinetics. New Jersey: Princeton University Press; 1955. [44] A-M. Wazwaz, The extended tanh method for abundant solitary wave solutions of nonlinear wave equations, Appl. Math. Comput., 187 (2007) 1131-1142. [45] W. Malfliet, Solitary wave solutions of nonlinear wave equations, Am. J. Phys., 60 (7) (1992) 650-654.
19
ACCEPTED MANUSCRIPT
[46] Y. Tan, H. Xu, S-J. Liao, Explicit series solution of travelling waves with a front of Fisher equation, Chaos Soliton. Fract., 31 (2007) 462-472. [47] H. N.A. Ismail, K. Raslan, A. A. Abd Rabboh, Adomian decomposition method for Burgers’Huxley and Burgers’-Fisher equations, Appl. Math. Comput., 159 (2004) 291-301.
CR IP T
[48] A. J. Khattak, A computational meshless method for the generalized Burgers’-Huxley equation, Appl. Math. Modell., 33 (2009) 3718-3729. [49] J. Canosa, Diffusion in nonlinear multiplicate media, J. Math. Phys., 31 (1969) 1862-1869. [50] A.H. Bhrawy, An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system, Applied Mathematics and Computation, 247, (2014) 30-46.
AN US
[51] E. Fan and Y. C. Hon, Generalized tanh method extended to special types of nonlinear equations, Z. Naturforsch 57a (2002) 692-700.
AC
CE
PT
ED
M
[52] R. A. Van Gorder and K. Vajravelu, Nonlinear dispersion of a pollutant ejected into a channel flow, Central European Journal of Physics 9 (2011) 1182-1194.
20
ACCEPTED MANUSCRIPT
0.8
uHx,0L
0.7
u Hx,0L
CR IP T
u amd u
0.9
uHx,0.5L
0.6
u Hx,0.5L
0.5
uHx,0.9L
0.4
u Hx,0.9L
0.3 -1.0
0.0
0.5
1.0
AN US
-0.5
x
b
PT
0.8
uHx,0L
0.7
u Hx,0L uHx,0.5L
0.6
u Hx,0.5L
0.5
uHx,0.9L 0.4
CE AC
0.9
u amd u
ED
M
a
u Hx,0.9L
0.3 -1.0
-0.5
0.0
0.5
1.0
x
c
d
Figure 1. We plot the (a) absolute error (b) approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 8. Then, we plot the (c) absolute error (d) approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.1) where λ1 = λ2 = λ3 = 1, α = β = 0 and N = 12.
21
AN US
CR IP T
ACCEPTED MANUSCRIPT
b
AC
CE
PT
ED
M
a
c
Figure 2. The absolute error of problem (3.12) where λ1 = λ2 = λ3 = 1, p = 5 α = β = −0.5 and (a) N = 4, (b) N = 8 and (c) N = 12, respectively.
22
ACCEPTED MANUSCRIPT
2.4 2.3 uHx,0L
CR IP T
2.2
u Hx,0L
u amd u
2.1
uHx,0.5L
2.0 1.9
u Hx,0.5L
1.8
uHx,0.9L
1.7 1.6 -1.0
AN US
u Hx,0.9L
0.0
-0.5
0.5
1.0
x
ED
M
Figure 3. The approximate solution and the exact solution for different values of t (0, 0.5 and 0.9) of problem (3.16) where N = 12, respectively.
EHx,0.5,0L
1.0
CE
1. ´ 10-9
6. ´ 10-12
PT
2. ´ 10-9 EHx,y,1L 1.5 ´ 10-9
8. ´ 10-12
4. ´ 10-12
0.5
5. ´ 10-10
2. ´ 10-12
0 -1.0
0.0
y
AC
-0.5
0.0 x
-0.5
0
0.5
-1.0 1.0
-0.5
-1.0
0.0
0.5
x
(b) E(x, 21 , 12 )
(a)E(x,y,1)
Figure 4. The absolute error of problem (3.20) with α1 = α2 = β1 = β2 = − 12 and N = M = K = 12.
23
1.0