- Email: [email protected]
Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation
Ain Shams Engineering Journal xxx (2016) xxx–xxx
Contents lists available at ScienceDirect
Ain Shams Engineering Journal journal homepage: www.sciencedirect.com
Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation H.S. Shekarabi, J. Rashidinia ⇑ Department of Mathematics, Karaj Branch, Islamic Azad University, Karaj, Iran School of Mathematics, Iran University of Science and Technology, Narmak, Tehran, Iran
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 4 June 2016 Revised 22 September 2016 Accepted 30 October 2016 Available online xxxx
In this work, the numerical approximation of Convection-Reaction-Diffusion equation is investigated using the method based on tension spline function and finite difference approximation. For nonlinear term, nonstandard finite difference method by nonlocal approximation is utilized. We describe the mathematical formulation procedure in detail and also analyze the stability and convergence of the method. Numerical results are provided to justify the good performance of the proposed scheme. Ó 2016 Ain Shams University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
MSC: 65M06 12 99 Keywords: Convection-Reaction-Diffusion equation Tension spline Nonstandard finite difference Implicit scheme Convergence analysis
1. Introduction We consider the Convection-Reaction-Diffusion equation:
ut þ aðx; tÞux ¼ bðx; tÞuxx þ gðuÞ;
ðx; tÞ 2 ½d; e ½t 0 ; T;
ð1Þ
with initial conditions
uðx; t 0 Þ ¼ /ðxÞ;
d 6 x 6 e;
ð2Þ
and Dirichlet boundary conditions
uðd; tÞ ¼ p0 ðtÞ;
uðe; tÞ ¼ p1 ðtÞ;
t0 6 t 6 T
ð3Þ
where /ðxÞ, p0 ðtÞ and p1 ðtÞ are known functions while the function uðx; tÞ is unknown. bðx; tÞ and aðx; tÞ are sufficiently smooth and positive functions. The Convection-Reaction-Diffusion equation implies three processes: Peer review under responsibility of Ain Shams University. ⇑ Corresponding author. E-mail addresses: [email protected] (H.S. Shekarabi), [email protected]. ir (J. Rashidinia).
Convection, Diffusion and Reaction. The term Convection is because of movement of molecules from a region (fluid) to another, while the term Diffusion means the spread of particles through random motion from a region with high concentration to a low concentration region. Finally, Reaction is due to decay, adsorption and chemical reactions of substance with another components, or accounting for biodegradation of the contaminants. aðx; tÞ denotes the velocity of the medium in the x direction and bðx; tÞ is the Diffusion coefficient [3]. Convection-Reaction-Diffusion equation (CRD) reveals itself in the mathematical modeling of many important fields of applied science such as physics, mechanics, biology, and hydrology [4]. For example Eq. (1) has been used to describe the following: transport phenomena in porous media [6], theory of combustion and detonation [7], population dynamics and mathematical biology [8,10,11,13], chemical reaction-diffusion [14,15], heat transfer in a draining film [16], finance [17], fluid flow [18,19] as well as transport chemistry in the atmosphere [20]. In hydrology, this type of equation models the transport of adsorbing contaminants and microbe-nutrient systems in groundwater [21]. The solutions of nonlinear Convection-Reaction-Diffusion equation are important in many applications aforesaid. The difficulties
http://dx.doi.org/10.1016/j.asej.2016.10.005 2090-4479/Ó 2016 Ain Shams University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
2
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
in obtaining analytical solutions made many researchers to use numerical methods in order to approximate the solution of these types of equations. There are a lot of studies on the numerical solutions of nonlinear (CRD) equations. Finite difference method, is one of the simplest and the oldest ones to deal with (CRD) equation. In [22,23] essentially non-oscillatory (ENO) and weighted ENO (WENO) schemes have been developed for solutions of partial differential equations. Application and comparison of the above mentioned schemes have been used for solution of ConvectionDiffusion in [24]. Standard finite difference in some problems produces numerical instabilities, like oscillations and numerical dispersion that appear in the solution approximation [3]. Nonstandard finite difference method for overcoming numerical instabilities is used in [25]. Meanwhile, new Eulerian-Lagrangian numerical method that combined the idea of the exact timestepping scheme with the finite difference method has been proposed in [26]. Recently in [27], a new finite difference method based on a special grid has been provided, the approach contains discretization in time, by using Crank-Nicolson and BackwardEuler finite difference schemes, while for the spatial discretization we used the method in [28]. The streamline-upwind PetrovGalerkin (SUPG) method is proposed in [29]. In new research, it is shown that SUPG-SC (streamline-upwind Petrov-Galerkin and shock capturing) technique is capable of reducing the unphysical oscillation in cross-wind direction [30–32,48]. Furthermore, several finite element methods have been developed for solving (CRD) equation [45–47]. In [33] some non-linear Convection dominated problems are solved with Discontinuous Galerkin method with the shock-capturing technique. Adaptive Discontinuous Galerkin method has been extended in [34]. Residual-free bubble method (RFB) which is based on enriching the finite elements space, is described in [35,36]. In [37,38] researchers developed the schemes that preserve some properties of exact solution like positivity. Another approach is the method of upper and a lower solution coupled with its associated monotone iterative techniques [39]. Compact exponential finite difference scheme for solving the time fractional (CRD) equation with variable coefficients is considered in [49]. In [12] a new monotone iterative technique known as modified accelerated monotone iterative method based on improving the rate of convergence, has been suggested. In [40,41] nonpolynomial cubic tension spline scheme has been used for solution of partial differential equation. In this paper we will focus on developing a new method using tension spline function and nonstandard finite difference scheme. We will demonstrate how to construct an implicit scheme for a (CRD) equation. This scheme is unconditionally stable. The arising system is tri-diagonal. Motivation behind the proposed scheme is to present a solution which can be obtained efficiently by the well-known algorithms. The rest of the paper is organized as follows: In Section 2, we describe the formulation of the method and achieve the system of equations. Section 3, is devoted to dealing with stability and convergence of the method. For this purpose, we use energy method. In Section 4, the method has been tested for some problems to demonstrate the effectiveness and accuracy of the proposed scheme. The paper ends with some conclusions in Section 5.
0 such that h ¼ ed , and tj ¼ t 0 þ jk, j ¼ 0; 1; . . . ; m such that j ¼ Tt , in m n which m and n are integers. For each segment ½xl ; xlþ1 , ðl ¼ 0; 1; . . . ; n 1Þ, and point t j , ðj ¼ 0; . . . ; mÞ we define the non-polynomial spline in tension as
j
j
sðx; tj Þ s j ðxÞ ¼ f l þ bl ðx xl Þ þ clj sinh xðx xl Þ þ j
j dl
ð4Þ
cosh xðx xl Þ;
j
j
where the f l , bl , clj , dl are unknown coefficients and x > 0 is a free parameter. The derivatives of non-polynomial tension spline function s j ðxÞ are given by 00
j
j
sj ðxÞ ¼ bl þ xclj cosh xðx xl Þ þ xdl sinh xðx xl Þ; l ¼ 1; 2; . . . ; n ð5Þ 00
j
sj ðxÞ ¼ x2 ðclj sinh xðx xl Þ þ dl cosh xðx xl ÞÞ; l ¼ 1; 2; . . . ; n ð6Þ Adapting the following notation:
s j ðxl Þ ¼ ulj ;
j s j ðxlþ1 Þ ¼ ulþ1 ;
0
0
sj ðxl Þ ¼ mlj ;
j sj ðxlþ1 Þ ¼ mlþ1 ;
00
00
sj ðxl Þ ¼ Mlj ;
j sj ðxlþ1 Þ ¼ Mlþ1 :
ð7Þ
The notation ulj be a grid function which is used in the discrete approximation value of uðlh; jkÞ, l ¼ 0; 1; . . . ; n and j ¼ 0; 1; . . . ; m. Now using (7) and after some algebraic manipulation, we can determine the unknown coefficients in (4) as M
j
j
f l ¼ ulj x2l ; M lj ; 2x2
j
dl ¼
j
bl ¼
u j u j lþ1
h
l
þ
M j M j l
lþ1
xh:
;
clj ¼
2M j ðeh þeh ÞM j lþ1
2x2 ðeh eh Þ
l
;
h ¼ xh:
Noticing the continuity of the first derivative at ðxl ; ul Þ, that is 0
0
jþ j sj ðxj l Þ ¼ s ðxl Þ, we obtain the following equations for l ¼ 1; . . . ; n
j j ulþ1 2ulj þ ul1
h
aþb h
j j ¼ aM lþ1 þ 2bM lj þ aMl1 ;
j j j ðulþ1 ul1 Þ ¼ amlj þ 2bmlj þ aml1 ;
where
a¼
2
1
1
h2
2h ; eh eh
b¼
ð8Þ
ð9Þ
h hðe þ eh Þ 1 : eh eh h2 1
When x ! 0, that h ! 0, then ða; bÞ ¼ 16 ; 13 , and the relations defined by (8) reduced to ordinary cubic spline. In the other words, 0
mlj ¼ sj ðxl Þ ¼ uxjl ¼
j j ulþ1 ul1 ah j j ½Mlþ1 M l1 2 2h
0
j ¼ sj ðxlþ1 Þ ¼ uxjlþ1 ¼ mlþ1
0
j ¼ sj ðxl1 Þ ¼ uxjl1 ¼ ml1
ð10Þ
j ulþ1 ulj j þ aMlj þ h½bM lþ1 h
ð11Þ
j ulj ul1 j þ aMlj h½bM l1 h
ð12Þ
2. Description of the method
Considering finite difference approximation
2.1. Non-polynomial tension spline
xj ¼ u lþ1
In order to develop the numerical method for approximation solution of (1)–(3), the domain of equation,½d; e ½t 0 ; T is divided into an n m mesh. Using the grid points xl ¼ d þ lh, l ¼ 0; 1; . . . ; n
xj ¼ u l
j j 3ulþ1 4ulj þ ul1 2 ¼ uxjlþ1 þ Oðh Þ 2h
j j ulþ1 ul1 2 ¼ uxjl þ Oðh Þ 2h
ð13Þ
ð14Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
3
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx j 3ulþ1 þ 4ulj ulþ1 2 ¼ uxjl1 þ Oðh Þ 2h
xj ¼ u l1
ujþ1 l
tj ¼ u l
2k
uj1 l
ð15Þ
2
¼ utjl þ Oðk Þ
ð16Þ
the spatial derivatives by non-polynomial tension spline would be 2 j xx ¼ s00 ðxl ; t j Þ ¼ M lj þ Oðh Þ u l
ð17Þ
We obtain ðj 1Þth and ðj þ 1Þth time level similar to (8) as follows: j1 uj1 þ uj1 lþ1 2ul l1
h ujþ1 lþ1
2
j1 ¼ aM j1 þ aMj1 lþ1 þ 2bM l l1
ð18Þ
þ
ujþ1 l1
jþ1 ¼ aM jþ1 þ aMjþ1 lþ1 þ 2bM l l1
ð19Þ
By using the relations of (10), (16) and(17), we can obtain a new approximation for the solutions of Eqs. (1)–(3) as follows:
jþ1
ul
j1
ul
2k
l
þ alj mlj g lj
¼
ð20Þ
gMj1 þ ð1 2gÞMlj þ gMjþ1 l l ; j
where 0 6 g 6 1 and bl ¼ bðlh; jkÞ, alj ¼ aðlh; jkÞ and g lj ¼ gðulj Þ. We multiply equations (18) and (19) by g and equation (8) by ð1 2gÞ and add together. We have
1h 2
h
1 h
2
½gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l
! ! jþ1 jþ1 ulþ1 uj1 ul uj1 2b lþ1 l ¼ þ ðA CÞ þ ðB FÞ þ j j 2k 2k bl blþ1 ! ! j jþ1 j j ulþ1 ul1 ul1 uj1 2al b a l1 þ þ þ ðD EÞ þ ðB FÞ j j 2k 2h bl bl1 ! ! j j a j a ulþ1 ulj a j a ulj ul1 þ lþ1 þ l1 j j h h b b
a
lþ1
2ujþ1 l 2 h
1 j b
To simplify Eq. (23), we use (13)–(16); therefore, we achieve
gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l
j þ ðA CÞalþ1
2h j j j 3u a j j l1 þ 4ul ulþ1 þ ðD EÞal1 þ ðA CÞ g lþ1 j 2h blþ1 2b a j j þ ðB FÞ g þ ðD EÞ g l1 ; l j j bl bl1
A¼
D¼
j j þ að1 2gÞM lþ1 þ 2ð1 2gÞbMlj þ að1 2gÞMl1 j1 þ agM j1 þ agMj1 lþ1 þ 2gbM l l1 ;
ð21Þ
We know that
j alþ1 abh j
ðblþ1 Þ
2
alj abh j j bl1 bl
;
;
B¼
j alþ1 a2 h
E¼
j
j
j al1 abh j
ðblþ1 Þ
j j j þ alþ1 mlþ1 g lþ1 jþ1 j1 u u l l þ alj mlj g lj þ 2bj 2k b l jþ1 j1 ul1 ul1 j j j þ a m g þ aj l1 l1 l1 ; 2k
jþ1
j1
ulþ1 ulþ1 2k
X¼h
Z¼h
2
h
Þd2x ulj
½g ¼
a
ð
j
blþ1 þ
þ ð1 2g
ujþ1 lþ1
j alþ1 hb j
2k
uj1 lþ1
þ
j Þ þ alþ1
h
ulj
d1 ¼
j g lþ1
j j ðutjlþ1 þ alþ1 uxjlþ1 g lþ1 Þþ
j alþ1 ha j
ðutjl þ alj uxjl g lj Þ
j
:
Y ¼h
2
ð25Þ
2b j
bl
! þ ðB FÞ ;
! þ ðD EÞ :
a
a
h
j j al1 uxjl1 g l1 Þ
j al1 ha j
bl
d4 ¼
! j h 2al b þ ðB FÞ ; 2 bj l j al1 ha
d2 ¼
j alþ1 ha j
blþ1
;
;
j
bl1
j alþ1 h ðA CÞ; 2
d5 ¼
j al1 h ðD EÞ; 2
and 2
a jh j j ðutjlþ1 þ alþ1 uxjlþ1 g lþ1 Þ þ l j ðutjl1 2bl1 jþ1 j j ul1 uj1 u u j j l1 l1 ð Þ þ al1 l g l1
j 2k bl1 j al1 hb j ðutl1 þ j bl1
j
bl1
d3 ¼
alj h j 2blþ1
a
a
2
þg
j ulþ1
þ ðA CÞ ;
j
blþ1
d2x uj1 l
bl blþ1 jþ1 j j j1 u u 2b ul ul l1 þ j ð Þ þ alj lþ1 g lj 2k 2h bl
j
bl1 bl
If we suppose
where l ¼ 1ð1Þn 1. Now we replace (10)–(12) in (22), and get
d2x ujþ1 l
;
j al1 a2 h
F¼
!
a
ð22Þ
bl1
1
j
j Zg l1 ;
2
j b lþ1
j
bl blþ1
jþ1 j j j j1 k1 ujþ1 þ k3 ujþ1 lþ1 þ k2 ul l1 ¼ k4 ulþ1 þ k5 ul þ k6 ul1 þ k7 ulþ1
where
gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l a
;
alj abh
The simplest form of Eq. (24) is
i
¼
2
C¼
j j þ k8 uj1 þ k9 uj1 l l1 Xg lþ1 Yg l
Using Eq. (20), and then eliminating M’s in above equations, lead us to
h
;
bl blþ1
j j d2x ulj ¼ ul1 2ulj þ ulþ1
1 h2
ð24Þ
where
i
jþ1 ¼ agMjþ1 þ agMjþ1 lþ1 þ 2gbM l l1
l1
j j 3ulþ1 4ulj þ ul1
j þ al1 uxjl1
j g l1 Þ
k1 ¼ g h2k ð
a
bj
þ ðA CÞÞ;
lþ1 2
k2 ¼ 2g h2k ð2bj þ ðB FÞÞ; bl
! 2 h a þ ðD EÞ ; k3 ¼ g 2k b j l1
ðutjl þ alj uxjl g lj Þ
then
ð23Þ
k4 ¼ ð1 2gÞ d1 d2 3d4 þ d5 ;
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
4
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
k5 ¼ 2ð1 2gÞ þ d2 d3 þ 4d4 4d5 ;
j j u j ulj 2ulj þ ul1 ujþ1 ulj ju l þ gðulj Þ þ alj lþ1 ¼ bl lþ1 2 k h h
k6 ¼ ð1 2gÞ þ d1 þ d3 d4 3d5 ;
Non-standard finite difference scheme was defined as follows by using two of Mickens rules:
and
k7 ¼ g
2 h a þ ðA CÞ ; 2k blþ1
k8 ¼ 2g
2 h 2b þ ðB FÞ ; 2k bl
k9 ¼ g
Definition 1. A finite difference scheme (27) is called nonstandard if at least one of the following conditions is met: (a) The classical denominator k or h of the discrete derivatives ðut Þlj or ðux Þlj or ðuxx Þlj is replaced by a non-negative function mðkÞ or mðhÞ where
2 h a þ ðD EÞ ; 2k bl1
2
s¼
h ; 2bk
mðzÞ ¼ z þ oðz2 Þ as 0 < z ! 0 mðzÞ is an increasing function of z:
k k¼ : h
mðzÞ may depend on the parameter appearing in the
It is obvious that, if aðx; tÞ and bðx; tÞ are constant functions, we have
A ¼ C;
B ¼ F;
d2 ¼ d3 ;
differential equation:
D ¼ E;
For example the denominator of simple equation ut ¼ ax ah 1
given by Mickens rules is mðkÞ ¼ e
d4 ¼ d5 ¼ 0; k1 ¼ k3 ;
X ¼ Z;
k7 ¼ k9 ;
jþ1 j j j j1 T lj ¼ k1 U jþ1 þ k3 U jþ1 lþ1 þ k2 U l l1 k4 U lþ1 k5 U l k6 U l1 k7 U lþ1
(b) In the expression nonlinear terms that appear in gðuÞ are approximated in a non-local way, by a suitable function of several points of the mesh. The nonlocal terms in the right-hand side of gðuÞ in (1) can be approximated non- locally in many different ways. For example
j j j k8 U j1 k9 U j1 l l1 þ Xg lþ1 þ Yg l þ Zg l1 ;
j
j j ðu2 Þl ’ ul1 ulþ1 ;
By using Taylor series expansion at the grid point ðxl ; tj Þ and
¼
such
aðxlj Þ
j Þ ¼ aðxlj Þ h aðxl1
approximation
þ OðhÞ
aðxlj Þ
j bðxl1 Þ
and
¼
bðxlj Þ
j
@bðx Þ h @xl
j @aðx Þ l
@x
j ulj ulþ1 ;
þ
(
¼
2
h
ð2a þ 2bÞ (
þ
2
h
1 j bl
alj
alj
bl
bl
þ j
j
j
ð2a þ 2bÞ þ
1 j bl
j
ul1 ulþ1 j u l
j
2
;
j j ul1 ulj ulþ1 ;
2
;
ð28Þ
j 1 ðulþ1 2
j þ ul1 Þðulj Þ
ð29Þ
2
This paper focuses on the use of nonlocal approximation strategy (second item in Definition 1) for the construction of nonstandard finite difference. Therefore we choose mðkÞ ¼ k and mðhÞ ¼ h.
@ulj @t
2 j @ ul d2 d3 2 þ h g þ ð1 2gÞ þ þ d4 d5 þ g 1 2 2 @x2 ( !) 3 @ 3 ulj k 2a 2b 2 þ j þ h j 3 @t 3 2bl k 2bl k 4 h 1 1 þ g þ ð1 2gÞ þ d2 d3 þ 2d4 2d5 þ g 12a 2 2 12 @ 4 ulj þ ; @x4
ð26Þ
Consistency condition dictates that 2a þ 2b ¼ 1 (where 2a þ 2b ¼ 1, the error tends to zero as h, k tend to zero simultane2
2
j
ul1 þulþ1
j þ ulj þ ul1 Þulj
2
j
!)
@ulj @x !)
j 1 ðulþ1 3
ulj
j ðu3 Þl ’ ðulj Þ ulþ1 ;
þ ¼ bðxlj Þþ
OðhÞ bðxlj Þ and substituting coefficient defined above, after simplification, we obtain
T lj
.
a
gðulj Þ,
In order to determine the local truncation error in (25), we suppose the solution of the Eqs. (1)–(3) defined by U then, we have
using
ð27Þ
2 2
4
ously). We can obtain the scheme of Oðk þ k h þ h Þ by Simplification of the above equation due to the consistency condition and 1 5 choosing a ¼ 12 and b ¼ 12 correspondingly. 2.2. Nonstandard finite difference Nonstandard finite difference schemes were introduced by Mickens in [42] as powerful numerical methods that preserve significant properties of exact solution of the involved differential equation. It is known that the finite difference scheme for Eq. (1) can be defined as
3. Stability of the method In this section, by following [40], we will investigate the stability and convergence of scheme (25) for numerical solutions of the Convection-Reaction-Diffusion equation using the initial and boundary conditions (2) and (3). For this purpose, we suppose that: first,aðx; tÞ ¼ a and bðx; tÞ ¼ b are constants and second: we introduce the difference operators and some lemmas.
dþt ulj ¼ ujþ1 ulj ; l
d2t ulj ¼ ujþ1 uj1 l l ; dt ulj ¼ ulj uj1 l ;
d2t ulj ¼ ujþ1 2ulj þ uj1 l l ;
d2t ulj ¼ ðdþt þ dt Þulj ; j ulj ; dþx ulj ¼ ulþ1
utjl ¼
1 þ j d u; k t l
d2t ulj ¼ ðdþt dt Þulj ;
j dx ulj ¼ ulj ul1 ;
uxjl ¼
1 þ j d u h x l
If fv i j0 6 i 6 ng and fwi j0 6 i 6 ng are two grid functions on Xh , we can define the inner product as follows: n X hv ; wi ¼ h v i wi ; i¼0
kv k ¼ hv ; v i2 ¼ 1
h
n X
v 2i
!12 :
i¼0
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
5
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
and we have
hdþx
v
; dþx wi
¼h
where n1 X
dþx
v
2
þ i dx wi
k1 ¼
i¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdþx v k ¼ hdþx v ; dþx v i ¼
h
n1 X
ðdþx v l Þ
2
!12
k2 ¼
2
To prove the stability and convergence of the method the following lemmas are required (see Cui[9] and Ben-Yu [5]). Lemma 3.1. For any two mesh function fv i j0 6 i 6 ng and fwi j0 6 i 6 ng, there holds
k01
¼
ah þ 2b þ 10h 20bg þ 2
k02
¼
2
2
8abh þ ð5h þ 8gbÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ah þ 2b þ 10h 20bg
2
2
2
2
2
2
16ab þ ðah þ 2b þ 10h 20bgÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 16ab þ ðah þ 2b þ 10h 20bgÞ : 4ab ð33Þ
W 4;þ1 denotes 4-periodic Sobolev spaces [44].
n1 n1 X X ðd2x v i Þwi ¼ ðdþx v i Þðdþx wi Þ w1 dþx v 0 þ wn dx v n : i¼1
Proof. We rewrite (25) in the following form:
gd2x d2t u j þ ð1 2gÞd2x u j þ sd2t u j 2 þaksd2x u j þ hb ð1 þ ad2x Þg j ¼ 0
It is easy to see that if w0 ¼ wn ¼ 0, then n1 n1 X X ðd2x v i Þwi ¼ ðdþx v i Þðdþx wi Þ: i¼1
5h 8gb þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 8abh þ ð5h þ 8gbÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
i¼0
i¼1
5h 8gb
ð34Þ
Multiplying (34) by ðhd2t ulj Þ and summing up for l from 0 to n yields to
i¼1
this lemma is the discrete Green formula.
hgd2x d2t u j ; d2t u j i þ ð1 2gÞhd2x u j ; d2t u j i þ shd2t u j ; d2t u j i Lemma 3.2. For any mesh function fv i j0 6 i 6 ng and defined on the interval ½L0 ; L1 with v 0 ¼ v n ¼ 0, the inequality holds:
1 þ 2 1 2 kd v k 6 kv k2 6 2 ðb aÞ kdþx v k2 : 4 x 8h
ð35Þ
2
þakshd2x u j ; d2t u j i þ hb hð1 þ ad2x Þg j ; d2t u j i ¼ 0
Using difference operators and Lemma 3.1, the first term of (35) becomes
hgd2x d2t u j ; d2t u j i ¼ ghdþx d2t u j ; dþx d2t u j i 6 gðkdþx d2t u j k2 Þ
ð36Þ
Similarly, for the second term of (35), we obtain Lemma 3.2 is analogous to Poincare inequality.
hd2x u j ; d2t u j i ¼ hdþx u j ; dþx ðujþ1 uj1 Þi
Lemma 3.3. Let xðsÞ and mðsÞ be nonnegative mesh functions. If C P 0 and mðsÞ is nondecreasing and
ð38Þ
Fourth term of (35), by the use of inequalities 2pq 6 p2 þ q2 becomes
i¼0
hd2x u j ; d2t u j i ¼ hdþx u j ; d2t u j i þ hdx u j ; d2t u j i
s
2
this lemma is the discrete Gronwall lemma. We verify the stability analysis of the scheme (25) by using the energy method described in the next theorem. Theorem 3.1. We assume that the exact solution of (1)–(3) satisfy U 2 W 4;þ1 ð0; þ1; L1 ðL0 ; L1 ÞÞ \ L1 ð0; þ1 : W 4;þ1 ðL0 ; L1 Þ. Let u be the solution of the difference scheme (25) and there exist constants ci > 0, i ¼ 1; 2 such that
jgðuÞj 6 c1 ð1 þ jujÞ @gðuÞ @u 6 c2
ð30Þ
Then for sufficiently small k there exits a constant C > 0 (independent of h and k) such that mþ1 kumþ1 k þ kum k6C t k þ kux
ð31Þ
provided that
k1 < k < k2
if g > 2bþah 4b
k01 < k < k02
if g < 2bþah 4b
2
2
2
12 ðkdþx u j k þ kd2t u j k þ kdx u j k þ kd2t u j k Þ
xðsÞ 6 mðsÞeCsk :
(
and for the third term of (35), we have
hd2t u j ; d2t u j i ¼ kd2t u j k2
j1 X xðsÞ 6 mðsÞ þ Ck xðiÞ;
then for all
ð37Þ
¼ ðhdþx u j ; dþx ujþ1 i hdþx u j ; dþx uj1 iÞ
ð32Þ
ð39Þ
For fifth term of (35) the inequalities ðp þ q þ u þ v Þ2 6 4ðp2 þ q2 þ u2 þ v 2 Þ and 2pq 6 p2 þ q2 imply that
1 2 2 hð1 þ ad2x Þg j ; d2t u j i kgðu j Þk þ kd2t u j k2 þ a2 kdþx gðu j Þk 2 1 2 þ kdþx d2t u j k2 kc1 ð1þ j u j jÞk 2 1 þ kd2t u j k2 þ a2 kc2 j dþx u j j k2 2 1 1 þ kdþx d2t u j k2 kd2t u j k2 2 2 1 þ j 2 þ kdx d2t u k þ 4c21 ð1 þ ku j k2 Þ 2 þ 4a2 ðc22 kdþx u j k2
ð40Þ
Substituting (36)–(40) into (35), we get
gkdþx d2t u j k2 þ ð1 2gÞðhdþx u j ; dþx ujþ1 i hdþx u j ; dþx uj1 iÞ þ skd2t u j k2 1 þ aksð ðkdþx u j k2 þ kdx u j k2 Þ þ kd2t u j k2 Þ 2 2 i h h1 1
kd2t u j k2 þ kdþx d2t u j k2 þ 4c21 ð1 þ ku j k2 Þ þ 4a2 ðc22 kdþx u j k2 Þ 2 b 2 ð41Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
6
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
Applying the sum j from 1 to m, we can obtain
g
m X
h ðkdþx d2t u j k2 Þ
m X
2
2b
j¼1
þ ðs þ aks
It can be proved that 2
kdþx d2t u j k2 þ ð1 2
kdþx dþt um k ¼ hdþx ðumþ1 um Þ; dþx ðumþ1 um Þi
Þðhdþx um ; dþx umþ1 iÞ
g
2
2 m m m X h X aks X kd2t u j k2 þ kdþx u j k2 þ kdx u j k2 Þ Þ ð 2b j¼1 2 j¼1 j¼1
ð42Þ Note that
P
kþ1 k i ui1 ui
¼
P
k kþ1 i ui uiþ1 ,
therefore
m h X þ ðs þ aks Þ kd2t u j k2 2b j¼1 2
þ
m m X aks X kdþx ujþ1 k2 kdþx uj1 k2 Þ ð1 2gÞhdþx u1 ; dþx u0 i ð 2 j¼1 j¼1
2 m m i X X h h 2 þ ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð43Þ
by simplifying Eq. (43), then we have m h X ðkdþx d2t u j k2 Þ þ ð1 2gÞðhdþx um ; dþx umþ1 iÞ Þ 2b j¼1 2
ðg
m h X aks kd2t u j k2 þ Þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2b j¼1 2 2
þ ðs þ aks
aks ðkdþx u0 k2 þ kdþx u1 k2 Þ 2 2 m m i X X h h 2 þ ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð1 2gÞhdþx u1 ; dþx u0 i þ
ð44Þ
ð45Þ
and similarly
2hdþx u1 ; dþx u0 i ¼ kdþx u1 k2 þ kdþx u0 k2 kdþx dþt u0 k2
ð46Þ
it follows that m h X ð1 2gÞ aks ðkdþx d2t u j k2 Þ þ ½ þ Þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2b j¼1 2 2
ðg
m ð1 2gÞ þ þ m 2 h X kdx dt u k þ ðs þ aks Þ kd2t u j k2 2 2b j¼1 2
ð1 2gÞ aks ð1 2gÞ þ þ 0 2 þ kdx dt u k ðkdþx u1 k2 þ kdþx u0 k2 Þ 2 2 2 2 m m X X þ 2 h þ ½4c21 ðm þ ku j k2 Þ þ 4a2 c22 kdx u j k ð47Þ b j¼1 j¼1
½
Now, we use Lemma 3.2 m 2h h X ð4g ðkd2t u j k2 Þ þ s þ aks Þ b 2b j¼1 2
2
ð1 2gÞ aks ð1 2gÞ þ þ m 2 þ kdx dt u k ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2 2 ð1 2gÞ aks ð1 2gÞ þ þ 0 2 þ kdx dt u k
½ ðkdþx u1 k2 þ kdþx u0 k2 Þ 2 2 2 2h m m i X X þ 2 h þ ku j k2 Þ þ 4a2 c22 kdx u j k ð48Þ 4c21 ðm þ b j¼1 j¼1 þ½
ð1 2gÞ þ þ m 2 ð1 2gÞ aks kdx dt u k þ ½ þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2 2 ð1 2gÞ aks ½ð1 2gÞ aks þ ð þ Þkdþt um k2 2 2 ð1 2gÞ aks þ þ ð ð50Þ Þðkdþx umþ1 k2 þ kdþx um k2 Þ 4 4
ð1 2gÞ aks þ f1 þ 2g þ aks ð Þgkdþt um k2 2 2 2 m X 5h þ f4g ðkd2t u j k2 þ s þ aksg 2b j¼1 ð1 2gÞ aks þ Þðkdþx umþ1 k2 þ kdþx um k2 Þ 4 4 ð1 2gÞ þ 2ðkdþx u1 k2 þ kdþx u0 k2 Þ
½ 2 aks 2 m m X X h ku j k2 Þ þ 4a2 c22 kdþx u j k2 þ 2ð1 2gÞkdþt u0 k2 þ ½4c21 ðm þ b j¼1 j¼1 þ ð
aks dþx u1 2 dþ u0 dþ u 0 2 1 ð1 2gÞ ¼k f 2ð þ Þðk k þ k x k2 Þ þ 2ð1 2gÞk t k2 2 2 h h k k m m X X 1 2 2 þ þ 2 ½4c21 ðm þ ku j k Þ þ 4a2 c22 kdx u j k g ð51Þ bk j¼1 j¼1 Regarding the arbitrary sufficiently small of , if 1 þ 2g þ aks > 0, i.e. g > 2bþah , it requires 4b 2
4g
With the following identity
2hdþx umþ1 ; dþx um i ¼ kdþx umþ1 k2 þ kdþx um k2 kdþx dþt um k2
Where is a sufficiently small coefficient, the equation may be rewritten as
m h X ðkdþx d2t u j k2 Þ þ ð1 2gÞðhdþx um ; dþx umþ1 iÞ Þ 2b j¼1 2
ðg
ð49Þ
Noticing to (49) and Lemma 3.2, for the second and the third terms on the left-hand side of (48), we have
2 m m i X X h h 2 ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð1 2gÞhdþx u1 ; dþx u0 i þ
2
2ðkdþx umþ1 k þ kdþx um k Þ
j¼1
2
5h h þ ak > 0 2b 2bk
ð52Þ
On the other hand, we know that
kd2t u j k2 ¼ kdþt u j þ dt u j k2 6 2ðkdþt u j k2 þ kdt u j k2 Þ
ð53Þ
We need k to satisfy the following condition (for j ¼ m): 2
2
5h h 1 þ 2g þ aks þ 2 4g þ ak 2b 2bk
!
>0
ð54Þ
Then by the use of Eqs. (52) and (54), we proved the stability condition (32). 0 1 þ 1 0 We need to show that 1h dþ and 1k dþ are all x u , h dx u t u bounded. It means that, all terms on the right-hand side of (51) have upper bound. 2 n1 2 n1 X X uðxlþ1 ; t 0 Þ uðxl ; t0 Þ 2 u0 u0l 1 k dþx u0 k ¼ ð lþ1 ð Þ h Þ h h h h l¼0 l¼0 Z b 2 n1 Z xlþ1 1X @/ @/ dxÞ
j dx C ð j ¼ h l¼0 @x xl a @x
ð55Þ
Note that from (25), u1 satisfies 2
h 1 ah þ u þ ðd þ dt Þu0l þ s0l ; 2bk l 2b t s0l ¼ k7 ½/ððl þ 1ÞhÞkðb/xx ððl þ 1ÞhÞ a/x ððl þ 1ÞhÞ þ gð/ðl þ 1ÞhÞÞ þ /ððl 1ÞhÞ kðb/xx ððl 1ÞhÞ a/x ððl 1ÞhÞ þ gð/ðl 1ÞhÞÞ þ k8 ð/ððlÞhÞ kðb/xx ððlÞhÞ a/x ððlÞhÞ þ gð/ðlÞhÞÞÞ
k1 d2x dþt u0l ¼ ð1 þ 2g k1 Þd2x u0l þ
2
h ðaðgðu0lþ1 Þ þ gðu0l1 ÞÞ þ bgðu0l ÞÞ b l ¼ 1ð1Þn 1
ð56Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
7
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx T
Simplifying Eq. (56), we have
ð
2bk1 2 þ 0 2b Þð1 þ 2g k1 Þd2x u0l Þdx dt ul þ dþt u0l ¼ ð ah ah h u1l u0l þ u1 þ s0l l ak
Taking inner product with
0 dþ t u
j where T j ¼ ðT 1j ; T 2j ; . . . ; T n1 Þ is the local truncation error, and we know that
2
ð57Þ
on both sides of (57), we get
2bk1 2 þ 0 þ 0 2b Þð1 þ 2g k1 Þhd2x u0 ; dþt u0 i kdþt u0 k2 ð Þhdx dt u ; dt u i ¼ ð ah ah h h hdþt u0 ; dþt u0 i þ ð1 Þhu0 ; dþt u0 i þ hdþt u0 ; s0 i ð58Þ ak ak By the use of Lemma 3.1, we can show that
2bk1 kdþt u0 k2 þ kdþx dþt u0 k 6 kCkdþt u0 k ah
ð60Þ
Therefore we have
kdþt u0 k þ 4
gd2x d2t E j þ ð1 2gÞd2x E j þ sd2t E j þ 2 2 aksd2x E j þ hb ð1 þ ad2x Þg j ¼ hb T j Multiplying (67) by ðhd2t E j Þ and taking the inner product, we get to
hgd2x d2t E j ; d2t E j i þ ð1 2gÞhd2x E j ; d2t E j i þ shd2t E j ; d2t E j i þ akshd2x E j ; d2t E j i þ
h hð1 þ ad2x Þ½gðU j Þ gðu j Þ; b
2
h j T ; d2t E j i b
ð68Þ
We know the inequalities ðp þ q þ uÞ2 3ðp2 þ q2 þ u2 Þ and pq 6 p2 þ 41 q2 ð8 > 0Þ, by applying these inequalities, imply that
hgd2x d2t E j ; d2t E j i þ ð1 2gÞhd2x E j ; d2t E j i þ shd2t E j ; d2t E j i þ akshd2x E j ; d2t E j i
1 þ 0 1 þ þ 0 d u þ d d u 6C k t k x t Provided that
ð66Þ
Now we subtract (65) from (34), and let E ¼ U u, then we have
d2t E j i ¼ h
g 1 Provided that 2bk > 0, we get k2 < 2b ah ak . Dividing the above inequality þ 0 by kdt u k, we obtain
2bk1 ah
4
2
ð59Þ
kdþt u0 k þ kdþx dþt u0 k 6 kC
2 2
jT j jl1 6 Cðk þ k h þ h Þ
ð61Þ
< 0, by the use of Lemma 3.2, we achieve
bk1 4bk1 kdþt u0 k 6 kC 1 þ kdþt u0 k 6 kC ah ah
ð62Þ
2
1 1 þ 0 1 þ þ 0 know that 1h dþ x u ¼ h dx u þ h dt dx u . It can also be proved that
ð63Þ
Regarding Eq. (63), we proved that all terms on the right-hand side of (51) are bounded. Hence from Lemma 3.2 in both cases we obtain þ mþ1 m mþ1 that kdþ k 6 kC, i.e.,kum k 6 C, under the t k þ kkux t u k þ kkdx u P þ j mþ1 0 ¼u þ m gives condition (32). The relation u j¼1 dt u , kumþ1 k 6 ku0 k þ mkC. So, the proof is completed. h Convergence of the spline difference scheme (25) by using the energy method is the issue of the next theorem.
ðg þ 12Þðkdþx d2t E j k2 Þ þ ak2s ðkdþx E j k2 þ kdx E j k2 Þ
difference scheme (25) converges to U in the discrete L2 -norm and we have
ð72Þ
By the use of Theorem 3.1 and summing up for j from 1 to m and using Lemma 3.3, we obtain (64) and the proof is completed. h 4. Numerical experiment In this section we present the numerical result of the new method on several problems. We can obtain u1 by using Taylor l about u0l and using (1): series expansion of u1 l 2
u1 ¼ u0l kðut Þ0l þ Oðk Þ l 2
¼ u0l kfbl ðuxx Þ0l a0l ðux Þ0l þ gðu0l Þg þ Oðk Þ u1 l By using initial values, we have
n o 0 u1 ¼ /ðlhÞ k bl /xx ðlhÞ a0l /x ðlhÞ þ gð/ðlhÞÞ l Example 1. We consider the following equation: [1]
mþ1 kðU uÞmþ1 k þ kðU uÞm k t k þ kðU uÞx 2 2
ð70Þ
ð71Þ
þð1 2gÞðhdþx E j ; dþx Ejþ1 i hdþx E j ; dþx Ej1 iÞ ! 2 1 h kd2t u j k2 þ s þ aks þ 2 b 2 h b a 1 6 2 3 þ C 22 kE j k2 þ kT j k2 4 b 2 h
0
Theorem 3.2. Suppose that the solution of (1)–(3) is sufficiently regular and defined by U and the assumptions of Theorem 3.1 are also valid. For sufficiently small k, the approximation solution of the spline
2
ð69Þ
Similar to Theorem 3.1, we get the following inequalities
1 Therefore we must have 1 þ 4bk 0, take k, h sufficiently small proah
2 a h ah 0 vided g > 2bk 4b , and then we have 1k dþ 6 C hence t u 1 þ þ 0 d d u 6 4 1 dþ u0 6 C. k x t k t 1 þ þ 0 0 In all cases we obtain 1k dþ þ k dx dt u 6 C. t u 1 is bounded. We For the last step, we must prove that 1h dþ xu
1 þ 1 1 þ 0 d u 6 d u þ k 1 dþ dþ u0 6 C h x h x h t x
2
h h hð1 þ ad2x Þ½gðU j Þ gðu j Þ; d2t E j i ¼ h T j ; d2t E j i b b * + h2 h2 1 kd2t E j k2 þ kT j k2 T j ; d2t E j 6 b 4 b þ
ut þ aux ¼ buxx ;
4
Cðk þ k h þ h Þ
0 6 t 6 T;
06x61
with a ¼ 3:5, b ¼ 0:022 and the following initial condition: Proof. By substituting U by u in (34), we obtain
/ðxÞ ¼ expðcxÞ
gd2x d2t U j þ ð1 2gÞd2x U j
The exact solution of above equation is taken as 2
þ sd2t U j þ aksd2x U j þ
2
h h ð1 þ ad2x Þg j ¼ T j b b
ð65Þ
uðx; tÞ ¼ expðcx þ ntÞ where c ¼ 0:02854797991928 and n ¼ 0:0999.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
8
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
Table 1 Maximum absolute error for Example 1. x
Ismail et al. [1]
0.10 0.50 0.90
t¼1
t¼2
t¼5
2.56(10) 8.93(10) 1.33(9)
2.38(10) 1.38(9) 2.83(9)
5.65(10) 1.91(9) 3.97(9)
Douglas method
0.10 0.50 0.90
t¼1
t¼2
t¼5
2.56(7) 8.37(7) 1.33(6)
2.37(7) 1.38(6) 2.82(6)
5.63(7) 1.90(6) 3.95(6)
Present method
0.10 0.50 0.90
t¼1
t¼2
t¼5
1.32898(11) 6.72048(11) 1.99724(9)
1.20256(11) 6.08158(11) 7.01701(11)
8.91154(12) 4.50668(11) 4.72124(10)
The boundary conditions can be obtained from the exact solution. We solve Example 1 with h ¼ 0:1 and k ¼ 0:001. The result are computed in different time levels. The maximum absolute errors in the solution are reported in Table 1. The computed results by our method are compared with[1]. The table shows our method is more accurate. Example 2. We consider the following equation:
ut þ aux ¼ buxx ;
t P 0;
06x61
Example 3. We consider the following model for the recolonization by oaks of Europe after the last glaciation; therefore, uðx; tÞ is the local density of oaks at time [38]
ut þ aux ¼ buxx þ ru;
t P 0;
0 6 x 6 10
If the population size at time 0 is M and is concentrated at the origin, the exact solution of this equation is
M ðx atÞ2 uðx; tÞ ¼ pffiffiffiffiffiffiffiffi exp rt 4bt 2 pbt
!
The boundary conditions can be obtained from the exact solution. In this example, we take a ¼ b ¼ 1, r ¼ 0:1 with h ¼ 0:1 and k ¼ 0:005. Computed solution for different time level has been compared with the exact solution and the absolute errors for t ¼ 1; 3; 5; 7; 9 and M ¼ 105 are reported in Table 3. We also computed solution for different M and the absolute errors for time t ¼ 2 with h ¼ 1 and k ¼ 0:01. These results are reported in Table 4. We observe that as we decrease the initial condition (M), the maximum absolute error tends to zero. Example 4. Consider the nonlinear Convection-Diffusion-Reaction equation [43]:
with the following initial condition:
ðx þ 0:5Þ2 /ðxÞ ¼ exp 0:00125
The boundary conditions can be obtained from the exact solution. For this example, we put a ¼ 1 and b ¼ 0:01. We take h ¼ 0:005 and k ¼ 0:001 and the maximum absolute error in the solution for t ¼ 0:6 and t ¼ 0:9 are computed. The results are reported in Table 2 and compared with [2]. In this problem, our numerical results are more accurate in comparison with the method in [2].
!
ut þ aux ¼ buxx þ cu2 ð1 uÞ;
t P 0;
06x61
and the exact solution is given by
The exact solution is
0:025 ðx þ 0:5 tÞ2 uðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 0:00125 þ 0:04t 0:000625 þ 0:02t
uðx; tÞ ¼ ð1 þ expðAðx BtÞ þ CÞÞ1
!
where A ¼
qffiffiffiffi c
, B ¼ 2a þ
4b
pffiffiffiffiffiffi cb and C ¼ AðB 1Þ and the initial con-
dition /ðxÞ ¼ ð1 þ expðAx þ CÞÞ1 .
Table 2 Maximum absolute error for Example 2. x
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
LBT [2]
QBT [2]
Present method
t ¼ 0:6
t ¼ 0:9
t ¼ 0:6
t ¼ 0:9
t ¼ 0:6
t ¼ 0:9
3.1(3) 1.3(3) 8.1(4) 1.7(5) 9.2(5) – – – –
0 1.1(3) 6.8(4) 0.5(6) 6.6(5) – – – –
1.3(4) 6.2(4) 3.8(5) 3.5(4) 4.5(5) – – – –
3.2(4) 1.02(4) 1.2(5) 3.2(4) 2.2(5) – – – –
3.90549(6) 1.76112(5) 5.19360(5) 1.78248(5) 8.85773(5) 2.03339(5) 3.24045(5) 2.06487(5) 4.53367(6)
1.11127(5) 6.67626(5) 1.47407(5) 3.15751(4) 7.0042(6) 5.00361(7) 1.37155(8) 1.54961(10) 7.51235(13)
Table 3 Maximum absolute error for Example 3. x
t¼1
t¼3
t¼5
t¼7
t¼9
1 2 3 4 5 6 7 8 9
2.54383(6) 2.65653(6) 2.75816(6) 2.84930(6) 2.92874(6) 2.99538(6) 3.04828(6) 3.08690(6) 3.11415(6)
1.09811(6) 1.14358(6) 1.19711(6) 1.25080(6) 1.30442(6) 1.35750(6) 1.40931(6) 1.45093(6) 1.49247(6)
6.93872(7) 6.5356(7) 6.81864(7) 7.08063(7) 7.29557(7) 7.41553(7) 7.35247(7) 6.93160(7) 5.98364(7)
4.29387(7) 4.02496(7) 4.12376(7) 4.15886(7) 4.05557(7) 3.68950(7) 2.83034(7) 1.49359(7) 3.35850(7)
4.24303(7) 2.56864(7) 2.53775(7) 2.42229(7) 2.10149(7) 1.39552(7) 2.00475(8) 3.05948(7) 1.52747(6)
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx Table 4 Maximum absolute error for Example 3. x
M ¼ 103
M ¼ 105
M ¼ 107
1 2 3 4 5 6 7 8 9
8.18795(4) 1.98019(3) 3.57075(3) 6.76625(3) 8.77917(4) 3.73394(3) 8.60632(4) 5.83540(4) 3.82148(4)
8.19047(6) 1.98074(5) 3.5717(5) 6.76608(5) 8.77694(6) 3.73414(5) 8.60557(6) 5.83505(6) 3.82128(6)
8.1905(8) 1.98075(7) 3.57118(7) 6.76608(7) 8.77692(8) 3.73414(7) 8.60556(8) 5.83504(8) 3.82128(8)
Table 5 Maximum absolute error for Example 4. x
c¼1 b ¼ 0:001
c ¼ 10 b ¼ 0:01
c ¼ 100
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
2.69954(8) 5.47176(9) 1.12576(9) 2.31615(10) 4.76527(11) 7.14614(12) 1.47025(12) 3.0249(13) 6.22346(14)
1.6958(8) 3.48895(9) 7.17818(10) 1.47684(10) 3.03847(11) 6.25136(12) 1.28617(12) 2.64531(13) 5.5412(14)
1.88386(10) 3.87587(11) 7.92424(12) 1.64063(12) 3.37543(13) 6.94464(14) 1.42879(14) 2.93948(15) 6.00635(16)
c¼1 a ¼ 100
c ¼ 100
c ¼ 1000
a¼5
a¼1
5.41124(8) 3.53667(7) 7.22626(7) 1.20347(6) 1.71858(6) 2.09293(6) 2.14112(6) 1.7592(6) 9.83998(7)
1.38435(3) 2.85513(4) 5.86763(5) 1.02727(5) 2.48386(6) 5.11032(7) 1.0514(7) 2.16311(8) 4.42828(9)
3.05902(7) 2.06115(9) 1.38879(11) 9.35762(14) 6.30512(16) 4.24835(18) 2.86252(20) 7.09547(23) 7.87728(25)
b ¼ 0:1
Table 6 Maximum absolute error for Example 4. x 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
The boundary conditions can be obtained from the exact solution. The nonlinear terms can be approximated non locally in many different ways. In this example, we consider an approximate of the j
j
2
j j j j form ðu2 Þl ¼ 13 ðulþ1 þ ulj þ ul1 Þulj and ðu3 Þl ¼ 12 ðul1 þ ulþ1 Þðulj Þ for nonlinear terms. In our numerical computation, in Table 5, we take a ¼ 10, h ¼ k ¼ 0:01. The results are computed for different c and b in time t ¼ 0:9. In Table 6, we take b ¼ 0:1, h ¼ k ¼ 0:01. The absolute errors for different c and a are reported in Table 6 in time t ¼ 0:9.
5. Conclusion In this study a numerical scheme has been established to approximate solution of Convection-Reaction-Diffusion equation. This scheme is based on non-polynomial tension spline in x direction, finite difference in t direction and nonstandard finite difference in nonlinear terms. Theoretical formulation of the proposed method was explained in detail. The optimal values of spline parameters and truncation error have been determined. The stability and convergence analyses of the present numerical approach by energy method, have been described. The performance of the method for these problems was measured by calculating the abso-
9
lute errors for different time levels which are reported in tabular form. In our proposed method, we developed the unconditionally three level implicit scheme. Our proposed method yields the tridiagonal system, which can be solved by efficient algorithm. The accuracy and applicability of our method have been tested by four examples. The numerical result shows that our scheme performs more accurately than the previous one. References [1] Ismail HNA, Elbarbary EM, Salem GSE. Restrictive Taylor approximation for solving convection-diffusion equation. Appl Math Comput 2004;147:355–63. [2] Dhawan S, Kapoor S, Kumar S. Numerical method for advection diffusion equation using FEM and B-splines. J Comput Sci 2012;3:429–37. [3] Makungu M, Haario H, Mahera WC. A generalised 1-dimensional particle transport method for convection diffusion reaction. Afr Math Un 2012;23:21–39. [4] Morton KW. Numerical solution of convection-diffusion problems. London: Chapman and Hall; 1996. [5] Ben-Yu G, Pascual PJ, Rodriguez MJ, Vazquez L. Numerical solution of the sineGordon equation. Appl Math Comput 1986;18(1):1–14. [6] Richards LA. Capillarity conduction of liquid through porous medium. Physica 1931;1:318. [7] Zeldovich Ya. Theory of flame propagation. Nat Adv Comm Aeronaut Tech Memor 1951;1282:39. [8] Kolmogorov A, Petrovskii I, Piskunoff N. In: Pelce P, editor. Dynamics of Curved Fronts. Boston: Academic Press; 1988. p. 105–30. [9] Cui M. Fourth-order compact scheme for the one-dimensional sine-Gorden equation. Numer Meth Partl Diff Eq 2009;25:685–711. [10] Murray JD. Mathematical biology. Berlin: Springer; 1989. [11] Witelski TP. Segregation and mixing in degenerate diffusion in population dynamics. J Math Biol 1997;35:695–712. [12] Wang YM. A modified accelerated monotone iterative method for finite difference reaction-diffusion-convection equations. J Comput Appl Math 2011;235:3646–60. [13] Aronson DG, Weinberger HF. Nonlinear diffusion in population genetics, combustion, and nerve propagation. In: Partial differential equations and related topics. Lecture notes in mathematics, vol. 446. New York: Springer; 1975. p. 5–49. [14] Aris R. The mathematical theory of diffusion and reaction in permeable catalysts. Clarendon (London): Oxford University Press; 1975. [15] Hetrick DK. Dynamics of nuclear reactors. Chicago: University of Chicago; 1971. [16] Isenberg J, Gutfinger C. Heat transfer to a draining film. Int J Heat Transfer 1972;16:505–12. [17] Wilmott P, Dewynne J, Howison S. Option pricing, mathematical models and computation. Oxford: Oxford Financial Press; 1993. [18] Roache P. Computational fluid dynamics. Albuquerque (NM): Hermosa Press; 1972. [19] Cengel YA, Cimbala JM. Fluid mechanics: fundamental and application. McGraw-Hill Higher Education; 2006. [20] McRea GJ, Goodin WR, Seinfeld JH. Numerical solution of atmospheric diffusion for chemically reacting flows. J Comput Phys 1982;77:1–42. [21] Cunningham AB, Charaklis WG, Abedeen F, Crawford D. influence of biofilm accumulation on porous media hydrodynamics. Environ Sci Tech 1991;25:1305–10. [22] Harten A, Enquist B, Osher S, Chakravarthy S. Uniformly high order essentially non-oscillatory schemes III. J Comput Phys 1987;71:231–303. [23] Liu XD, Osher S, Chan T. Weighted essentially non-oscillatory schemes. J Comput Phys 1994;115:200–12. [24] John V, Novo J. On (essentially) non-oscillatory discretizations of evolutionary convection-diffusion equations. J Comput Phys 2012;231:1570–86. [25] Kailash PC. On the use of nonstandard finite difference methods. J Diff Eq Appl 2005;11:735–58. [26] Kojouharov HV, Chen BM. Nonstandard methods for the convective-dispersive transport equation with nonlinear reactions. John Wiley and Sons; 1999. p. 617–24. [27] Kaya A. Finite difference approximations of multidimensional unsteady convection-diffusion-reaction equations. J Comput Phys 2015;285:331–49. [28] Kaya A. A finite difference scheme for multidimensional convection-diffusionreaction equations. Comput Meth Appl Mech Eng 2014;278:347–60. [29] Brooks AN, Hughes TJR. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [30] Bause M. Stabilized finite element methods with shock-capturing for nonlinear convection-diffusion-reaction models. Int Numer Math Adv Appl 2009:125–34. [31] Bause M, Schwegler K. Analysis of stabilized higher-order finite element approximation of nonstationary and nonlinear convection-diffusion-reaction equations. Comput Meth Appl Mech Eng 2012;209-212:184–96. [32] Bause M, Schwegler K. Higher order finite element approximation of systems of convection diffusion reaction equations with small diffusion. J Comput Appl Math 2013;246:52–64.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
10
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
[33] Yucel H, Stoll M, Benner P. Discontinuous Galerkin finite element methods with shock-capturing for nonlinear convection dominated models. Comp Chem Eng 2013;58:278–87. [34] Uzunca M, Karasozen B, Manguoglu M. Adaptive discontinuous Galerkin methods for non-linear diffusion-convection-reaction equations. Comp Chem Eng 2014;68:24–37. [35] Brezzi F, Hauke G, Marini LD, Sangalli G. Link-cutting bubbles for the stabilization of convection-diffusion-reaction problems. Math Models Meth Appl Sci 2003;13:445–61. [36] Sendur A, Neslitrk AI. Applications of the pseudo residual-free bubbles to the stabilization of convection-diffusion-reaction problems. Calcolo 2011;49:1–19. [37] Quang DA, Ehrhardt M. Adequate numerical solution of air pollution problems by positive difference schemes on unbounded domains. Math Comput Model 2006;44:834–56. [38] Chen-Charpentier BM, Kojouharov HV. An unconditionally positivity preserving scheme for advection-diffusion-reaction equations. Math Comp Model 2013;57:2177–85. [39] Pao CV. Monotone iterative methods for finite difference system of reaction diffusion equations. Numer Math 1985;46:571–86. [40] Rashidinia J, Mohammadi R. Tension spline approach for the numerical solution of nonlinear KleinGordon equation. Comp Phys Commun 2010;181:78–91. [41] Aghamohammadi M, Rashidinia J, Ezzati R. Tension spline method for solution of non-linear Fisher equation. Appl Math Comput 2014;249:399–407. [42] Mickens RE. Nonstandard finite difference model of differential equation. Singapore: World Scientific; 1994. [43] Hundsdorfer W. Numerical solution of advection-diffusion-reaction equations. Lecture notes. Thomas Stieltjes Institute; 2000. [44] Atkinson K, Han W. Theoretical numerical analysis: a functional analysis framework. Springer; 2005. [45] Theeraek P, Phongthanapanich S, Dechaumphai P. Solving convectiondiffusion-reaction equation by adaptive finite volume element method. Math Comp Simul 2011;82:220–33.
[46] Duan H, Hsieh P, Tan R, Yang S. Analysis of a new stabilized finite element method for the reaction-convection-diffusion equations with a large reaction coefficient. Comput Meth Appl Mech Eng 2012;247–248:15–36. [47] Hsieh P, Yang S. A new stabilized linear finite element method for solving reaction-convection-diffusion equations. Comput Meth Appl Mech Eng 2016;307:362–82. [48] Gierea S, Iliescub T, Johna V, Wellsb D. SUPG reduced order models for convection-dominated convection-diffusion-reaction equations. Comput Meth Appl Mech Eng 2015;289:454–74. [49] Cui M. Compact exponential scheme for the time fractional convectiondiffusion-reaction equation with variable coefficients. J Comput Phys 2015;280:143–63.
H.S. Shekarabi is Ph.D. student of College of Basic Science, Karaj Branch, Islamic Azad University.
Dr. J. Rashidinia: For photography and biography see Ain Shams Engineering Journal, Volume 6, Issue 1, March 2015, Pages 381–389.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
Contents lists available at ScienceDirect
Ain Shams Engineering Journal journal homepage: www.sciencedirect.com
Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation H.S. Shekarabi, J. Rashidinia ⇑ Department of Mathematics, Karaj Branch, Islamic Azad University, Karaj, Iran School of Mathematics, Iran University of Science and Technology, Narmak, Tehran, Iran
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 4 June 2016 Revised 22 September 2016 Accepted 30 October 2016 Available online xxxx
In this work, the numerical approximation of Convection-Reaction-Diffusion equation is investigated using the method based on tension spline function and finite difference approximation. For nonlinear term, nonstandard finite difference method by nonlocal approximation is utilized. We describe the mathematical formulation procedure in detail and also analyze the stability and convergence of the method. Numerical results are provided to justify the good performance of the proposed scheme. Ó 2016 Ain Shams University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
MSC: 65M06 12 99 Keywords: Convection-Reaction-Diffusion equation Tension spline Nonstandard finite difference Implicit scheme Convergence analysis
1. Introduction We consider the Convection-Reaction-Diffusion equation:
ut þ aðx; tÞux ¼ bðx; tÞuxx þ gðuÞ;
ðx; tÞ 2 ½d; e ½t 0 ; T;
ð1Þ
with initial conditions
uðx; t 0 Þ ¼ /ðxÞ;
d 6 x 6 e;
ð2Þ
and Dirichlet boundary conditions
uðd; tÞ ¼ p0 ðtÞ;
uðe; tÞ ¼ p1 ðtÞ;
t0 6 t 6 T
ð3Þ
where /ðxÞ, p0 ðtÞ and p1 ðtÞ are known functions while the function uðx; tÞ is unknown. bðx; tÞ and aðx; tÞ are sufficiently smooth and positive functions. The Convection-Reaction-Diffusion equation implies three processes: Peer review under responsibility of Ain Shams University. ⇑ Corresponding author. E-mail addresses: [email protected] (H.S. Shekarabi), [email protected]. ir (J. Rashidinia).
Convection, Diffusion and Reaction. The term Convection is because of movement of molecules from a region (fluid) to another, while the term Diffusion means the spread of particles through random motion from a region with high concentration to a low concentration region. Finally, Reaction is due to decay, adsorption and chemical reactions of substance with another components, or accounting for biodegradation of the contaminants. aðx; tÞ denotes the velocity of the medium in the x direction and bðx; tÞ is the Diffusion coefficient [3]. Convection-Reaction-Diffusion equation (CRD) reveals itself in the mathematical modeling of many important fields of applied science such as physics, mechanics, biology, and hydrology [4]. For example Eq. (1) has been used to describe the following: transport phenomena in porous media [6], theory of combustion and detonation [7], population dynamics and mathematical biology [8,10,11,13], chemical reaction-diffusion [14,15], heat transfer in a draining film [16], finance [17], fluid flow [18,19] as well as transport chemistry in the atmosphere [20]. In hydrology, this type of equation models the transport of adsorbing contaminants and microbe-nutrient systems in groundwater [21]. The solutions of nonlinear Convection-Reaction-Diffusion equation are important in many applications aforesaid. The difficulties
http://dx.doi.org/10.1016/j.asej.2016.10.005 2090-4479/Ó 2016 Ain Shams University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
2
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
in obtaining analytical solutions made many researchers to use numerical methods in order to approximate the solution of these types of equations. There are a lot of studies on the numerical solutions of nonlinear (CRD) equations. Finite difference method, is one of the simplest and the oldest ones to deal with (CRD) equation. In [22,23] essentially non-oscillatory (ENO) and weighted ENO (WENO) schemes have been developed for solutions of partial differential equations. Application and comparison of the above mentioned schemes have been used for solution of ConvectionDiffusion in [24]. Standard finite difference in some problems produces numerical instabilities, like oscillations and numerical dispersion that appear in the solution approximation [3]. Nonstandard finite difference method for overcoming numerical instabilities is used in [25]. Meanwhile, new Eulerian-Lagrangian numerical method that combined the idea of the exact timestepping scheme with the finite difference method has been proposed in [26]. Recently in [27], a new finite difference method based on a special grid has been provided, the approach contains discretization in time, by using Crank-Nicolson and BackwardEuler finite difference schemes, while for the spatial discretization we used the method in [28]. The streamline-upwind PetrovGalerkin (SUPG) method is proposed in [29]. In new research, it is shown that SUPG-SC (streamline-upwind Petrov-Galerkin and shock capturing) technique is capable of reducing the unphysical oscillation in cross-wind direction [30–32,48]. Furthermore, several finite element methods have been developed for solving (CRD) equation [45–47]. In [33] some non-linear Convection dominated problems are solved with Discontinuous Galerkin method with the shock-capturing technique. Adaptive Discontinuous Galerkin method has been extended in [34]. Residual-free bubble method (RFB) which is based on enriching the finite elements space, is described in [35,36]. In [37,38] researchers developed the schemes that preserve some properties of exact solution like positivity. Another approach is the method of upper and a lower solution coupled with its associated monotone iterative techniques [39]. Compact exponential finite difference scheme for solving the time fractional (CRD) equation with variable coefficients is considered in [49]. In [12] a new monotone iterative technique known as modified accelerated monotone iterative method based on improving the rate of convergence, has been suggested. In [40,41] nonpolynomial cubic tension spline scheme has been used for solution of partial differential equation. In this paper we will focus on developing a new method using tension spline function and nonstandard finite difference scheme. We will demonstrate how to construct an implicit scheme for a (CRD) equation. This scheme is unconditionally stable. The arising system is tri-diagonal. Motivation behind the proposed scheme is to present a solution which can be obtained efficiently by the well-known algorithms. The rest of the paper is organized as follows: In Section 2, we describe the formulation of the method and achieve the system of equations. Section 3, is devoted to dealing with stability and convergence of the method. For this purpose, we use energy method. In Section 4, the method has been tested for some problems to demonstrate the effectiveness and accuracy of the proposed scheme. The paper ends with some conclusions in Section 5.
0 such that h ¼ ed , and tj ¼ t 0 þ jk, j ¼ 0; 1; . . . ; m such that j ¼ Tt , in m n which m and n are integers. For each segment ½xl ; xlþ1 , ðl ¼ 0; 1; . . . ; n 1Þ, and point t j , ðj ¼ 0; . . . ; mÞ we define the non-polynomial spline in tension as
j
j
sðx; tj Þ s j ðxÞ ¼ f l þ bl ðx xl Þ þ clj sinh xðx xl Þ þ j
j dl
ð4Þ
cosh xðx xl Þ;
j
j
where the f l , bl , clj , dl are unknown coefficients and x > 0 is a free parameter. The derivatives of non-polynomial tension spline function s j ðxÞ are given by 00
j
j
sj ðxÞ ¼ bl þ xclj cosh xðx xl Þ þ xdl sinh xðx xl Þ; l ¼ 1; 2; . . . ; n ð5Þ 00
j
sj ðxÞ ¼ x2 ðclj sinh xðx xl Þ þ dl cosh xðx xl ÞÞ; l ¼ 1; 2; . . . ; n ð6Þ Adapting the following notation:
s j ðxl Þ ¼ ulj ;
j s j ðxlþ1 Þ ¼ ulþ1 ;
0
0
sj ðxl Þ ¼ mlj ;
j sj ðxlþ1 Þ ¼ mlþ1 ;
00
00
sj ðxl Þ ¼ Mlj ;
j sj ðxlþ1 Þ ¼ Mlþ1 :
ð7Þ
The notation ulj be a grid function which is used in the discrete approximation value of uðlh; jkÞ, l ¼ 0; 1; . . . ; n and j ¼ 0; 1; . . . ; m. Now using (7) and after some algebraic manipulation, we can determine the unknown coefficients in (4) as M
j
j
f l ¼ ulj x2l ; M lj ; 2x2
j
dl ¼
j
bl ¼
u j u j lþ1
h
l
þ
M j M j l
lþ1
xh:
;
clj ¼
2M j ðeh þeh ÞM j lþ1
2x2 ðeh eh Þ
l
;
h ¼ xh:
Noticing the continuity of the first derivative at ðxl ; ul Þ, that is 0
0
jþ j sj ðxj l Þ ¼ s ðxl Þ, we obtain the following equations for l ¼ 1; . . . ; n
j j ulþ1 2ulj þ ul1
h
aþb h
j j ¼ aM lþ1 þ 2bM lj þ aMl1 ;
j j j ðulþ1 ul1 Þ ¼ amlj þ 2bmlj þ aml1 ;
where
a¼
2
1
1
h2
2h ; eh eh
b¼
ð8Þ
ð9Þ
h hðe þ eh Þ 1 : eh eh h2 1
When x ! 0, that h ! 0, then ða; bÞ ¼ 16 ; 13 , and the relations defined by (8) reduced to ordinary cubic spline. In the other words, 0
mlj ¼ sj ðxl Þ ¼ uxjl ¼
j j ulþ1 ul1 ah j j ½Mlþ1 M l1 2 2h
0
j ¼ sj ðxlþ1 Þ ¼ uxjlþ1 ¼ mlþ1
0
j ¼ sj ðxl1 Þ ¼ uxjl1 ¼ ml1
ð10Þ
j ulþ1 ulj j þ aMlj þ h½bM lþ1 h
ð11Þ
j ulj ul1 j þ aMlj h½bM l1 h
ð12Þ
2. Description of the method
Considering finite difference approximation
2.1. Non-polynomial tension spline
xj ¼ u lþ1
In order to develop the numerical method for approximation solution of (1)–(3), the domain of equation,½d; e ½t 0 ; T is divided into an n m mesh. Using the grid points xl ¼ d þ lh, l ¼ 0; 1; . . . ; n
xj ¼ u l
j j 3ulþ1 4ulj þ ul1 2 ¼ uxjlþ1 þ Oðh Þ 2h
j j ulþ1 ul1 2 ¼ uxjl þ Oðh Þ 2h
ð13Þ
ð14Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
3
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx j 3ulþ1 þ 4ulj ulþ1 2 ¼ uxjl1 þ Oðh Þ 2h
xj ¼ u l1
ujþ1 l
tj ¼ u l
2k
uj1 l
ð15Þ
2
¼ utjl þ Oðk Þ
ð16Þ
the spatial derivatives by non-polynomial tension spline would be 2 j xx ¼ s00 ðxl ; t j Þ ¼ M lj þ Oðh Þ u l
ð17Þ
We obtain ðj 1Þth and ðj þ 1Þth time level similar to (8) as follows: j1 uj1 þ uj1 lþ1 2ul l1
h ujþ1 lþ1
2
j1 ¼ aM j1 þ aMj1 lþ1 þ 2bM l l1
ð18Þ
þ
ujþ1 l1
jþ1 ¼ aM jþ1 þ aMjþ1 lþ1 þ 2bM l l1
ð19Þ
By using the relations of (10), (16) and(17), we can obtain a new approximation for the solutions of Eqs. (1)–(3) as follows:
jþ1
ul
j1
ul
2k
l
þ alj mlj g lj
¼
ð20Þ
gMj1 þ ð1 2gÞMlj þ gMjþ1 l l ; j
where 0 6 g 6 1 and bl ¼ bðlh; jkÞ, alj ¼ aðlh; jkÞ and g lj ¼ gðulj Þ. We multiply equations (18) and (19) by g and equation (8) by ð1 2gÞ and add together. We have
1h 2
h
1 h
2
½gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l
! ! jþ1 jþ1 ulþ1 uj1 ul uj1 2b lþ1 l ¼ þ ðA CÞ þ ðB FÞ þ j j 2k 2k bl blþ1 ! ! j jþ1 j j ulþ1 ul1 ul1 uj1 2al b a l1 þ þ þ ðD EÞ þ ðB FÞ j j 2k 2h bl bl1 ! ! j j a j a ulþ1 ulj a j a ulj ul1 þ lþ1 þ l1 j j h h b b
a
lþ1
2ujþ1 l 2 h
1 j b
To simplify Eq. (23), we use (13)–(16); therefore, we achieve
gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l
j þ ðA CÞalþ1
2h j j j 3u a j j l1 þ 4ul ulþ1 þ ðD EÞal1 þ ðA CÞ g lþ1 j 2h blþ1 2b a j j þ ðB FÞ g þ ðD EÞ g l1 ; l j j bl bl1
A¼
D¼
j j þ að1 2gÞM lþ1 þ 2ð1 2gÞbMlj þ að1 2gÞMl1 j1 þ agM j1 þ agMj1 lþ1 þ 2gbM l l1 ;
ð21Þ
We know that
j alþ1 abh j
ðblþ1 Þ
2
alj abh j j bl1 bl
;
;
B¼
j alþ1 a2 h
E¼
j
j
j al1 abh j
ðblþ1 Þ
j j j þ alþ1 mlþ1 g lþ1 jþ1 j1 u u l l þ alj mlj g lj þ 2bj 2k b l jþ1 j1 ul1 ul1 j j j þ a m g þ aj l1 l1 l1 ; 2k
jþ1
j1
ulþ1 ulþ1 2k
X¼h
Z¼h
2
h
Þd2x ulj
½g ¼
a
ð
j
blþ1 þ
þ ð1 2g
ujþ1 lþ1
j alþ1 hb j
2k
uj1 lþ1
þ
j Þ þ alþ1
h
ulj
d1 ¼
j g lþ1
j j ðutjlþ1 þ alþ1 uxjlþ1 g lþ1 Þþ
j alþ1 ha j
ðutjl þ alj uxjl g lj Þ
j
:
Y ¼h
2
ð25Þ
2b j
bl
! þ ðB FÞ ;
! þ ðD EÞ :
a
a
h
j j al1 uxjl1 g l1 Þ
j al1 ha j
bl
d4 ¼
! j h 2al b þ ðB FÞ ; 2 bj l j al1 ha
d2 ¼
j alþ1 ha j
blþ1
;
;
j
bl1
j alþ1 h ðA CÞ; 2
d5 ¼
j al1 h ðD EÞ; 2
and 2
a jh j j ðutjlþ1 þ alþ1 uxjlþ1 g lþ1 Þ þ l j ðutjl1 2bl1 jþ1 j j ul1 uj1 u u j j l1 l1 ð Þ þ al1 l g l1
j 2k bl1 j al1 hb j ðutl1 þ j bl1
j
bl1
d3 ¼
alj h j 2blþ1
a
a
2
þg
j ulþ1
þ ðA CÞ ;
j
blþ1
d2x uj1 l
bl blþ1 jþ1 j j j1 u u 2b ul ul l1 þ j ð Þ þ alj lþ1 g lj 2k 2h bl
j
bl1 bl
If we suppose
where l ¼ 1ð1Þn 1. Now we replace (10)–(12) in (22), and get
d2x ujþ1 l
;
j al1 a2 h
F¼
!
a
ð22Þ
bl1
1
j
j Zg l1 ;
2
j b lþ1
j
bl blþ1
jþ1 j j j j1 k1 ujþ1 þ k3 ujþ1 lþ1 þ k2 ul l1 ¼ k4 ulþ1 þ k5 ul þ k6 ul1 þ k7 ulþ1
where
gd2x ujþ1 þ ð1 2gÞd2x ulj þ gd2x uj1 l l a
;
alj abh
The simplest form of Eq. (24) is
i
¼
2
C¼
j j þ k8 uj1 þ k9 uj1 l l1 Xg lþ1 Yg l
Using Eq. (20), and then eliminating M’s in above equations, lead us to
h
;
bl blþ1
j j d2x ulj ¼ ul1 2ulj þ ulþ1
1 h2
ð24Þ
where
i
jþ1 ¼ agMjþ1 þ agMjþ1 lþ1 þ 2gbM l l1
l1
j j 3ulþ1 4ulj þ ul1
j þ al1 uxjl1
j g l1 Þ
k1 ¼ g h2k ð
a
bj
þ ðA CÞÞ;
lþ1 2
k2 ¼ 2g h2k ð2bj þ ðB FÞÞ; bl
! 2 h a þ ðD EÞ ; k3 ¼ g 2k b j l1
ðutjl þ alj uxjl g lj Þ
then
ð23Þ
k4 ¼ ð1 2gÞ d1 d2 3d4 þ d5 ;
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
4
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
k5 ¼ 2ð1 2gÞ þ d2 d3 þ 4d4 4d5 ;
j j u j ulj 2ulj þ ul1 ujþ1 ulj ju l þ gðulj Þ þ alj lþ1 ¼ bl lþ1 2 k h h
k6 ¼ ð1 2gÞ þ d1 þ d3 d4 3d5 ;
Non-standard finite difference scheme was defined as follows by using two of Mickens rules:
and
k7 ¼ g
2 h a þ ðA CÞ ; 2k blþ1
k8 ¼ 2g
2 h 2b þ ðB FÞ ; 2k bl
k9 ¼ g
Definition 1. A finite difference scheme (27) is called nonstandard if at least one of the following conditions is met: (a) The classical denominator k or h of the discrete derivatives ðut Þlj or ðux Þlj or ðuxx Þlj is replaced by a non-negative function mðkÞ or mðhÞ where
2 h a þ ðD EÞ ; 2k bl1
2
s¼
h ; 2bk
mðzÞ ¼ z þ oðz2 Þ as 0 < z ! 0 mðzÞ is an increasing function of z:
k k¼ : h
mðzÞ may depend on the parameter appearing in the
It is obvious that, if aðx; tÞ and bðx; tÞ are constant functions, we have
A ¼ C;
B ¼ F;
d2 ¼ d3 ;
differential equation:
D ¼ E;
For example the denominator of simple equation ut ¼ ax ah 1
given by Mickens rules is mðkÞ ¼ e
d4 ¼ d5 ¼ 0; k1 ¼ k3 ;
X ¼ Z;
k7 ¼ k9 ;
jþ1 j j j j1 T lj ¼ k1 U jþ1 þ k3 U jþ1 lþ1 þ k2 U l l1 k4 U lþ1 k5 U l k6 U l1 k7 U lþ1
(b) In the expression nonlinear terms that appear in gðuÞ are approximated in a non-local way, by a suitable function of several points of the mesh. The nonlocal terms in the right-hand side of gðuÞ in (1) can be approximated non- locally in many different ways. For example
j j j k8 U j1 k9 U j1 l l1 þ Xg lþ1 þ Yg l þ Zg l1 ;
j
j j ðu2 Þl ’ ul1 ulþ1 ;
By using Taylor series expansion at the grid point ðxl ; tj Þ and
¼
such
aðxlj Þ
j Þ ¼ aðxlj Þ h aðxl1
approximation
þ OðhÞ
aðxlj Þ
j bðxl1 Þ
and
¼
bðxlj Þ
j
@bðx Þ h @xl
j @aðx Þ l
@x
j ulj ulþ1 ;
þ
(
¼
2
h
ð2a þ 2bÞ (
þ
2
h
1 j bl
alj
alj
bl
bl
þ j
j
j
ð2a þ 2bÞ þ
1 j bl
j
ul1 ulþ1 j u l
j
2
;
j j ul1 ulj ulþ1 ;
2
;
ð28Þ
j 1 ðulþ1 2
j þ ul1 Þðulj Þ
ð29Þ
2
This paper focuses on the use of nonlocal approximation strategy (second item in Definition 1) for the construction of nonstandard finite difference. Therefore we choose mðkÞ ¼ k and mðhÞ ¼ h.
@ulj @t
2 j @ ul d2 d3 2 þ h g þ ð1 2gÞ þ þ d4 d5 þ g 1 2 2 @x2 ( !) 3 @ 3 ulj k 2a 2b 2 þ j þ h j 3 @t 3 2bl k 2bl k 4 h 1 1 þ g þ ð1 2gÞ þ d2 d3 þ 2d4 2d5 þ g 12a 2 2 12 @ 4 ulj þ ; @x4
ð26Þ
Consistency condition dictates that 2a þ 2b ¼ 1 (where 2a þ 2b ¼ 1, the error tends to zero as h, k tend to zero simultane2
2
j
ul1 þulþ1
j þ ulj þ ul1 Þulj
2
j
!)
@ulj @x !)
j 1 ðulþ1 3
ulj
j ðu3 Þl ’ ðulj Þ ulþ1 ;
þ ¼ bðxlj Þþ
OðhÞ bðxlj Þ and substituting coefficient defined above, after simplification, we obtain
T lj
.
a
gðulj Þ,
In order to determine the local truncation error in (25), we suppose the solution of the Eqs. (1)–(3) defined by U then, we have
using
ð27Þ
2 2
4
ously). We can obtain the scheme of Oðk þ k h þ h Þ by Simplification of the above equation due to the consistency condition and 1 5 choosing a ¼ 12 and b ¼ 12 correspondingly. 2.2. Nonstandard finite difference Nonstandard finite difference schemes were introduced by Mickens in [42] as powerful numerical methods that preserve significant properties of exact solution of the involved differential equation. It is known that the finite difference scheme for Eq. (1) can be defined as
3. Stability of the method In this section, by following [40], we will investigate the stability and convergence of scheme (25) for numerical solutions of the Convection-Reaction-Diffusion equation using the initial and boundary conditions (2) and (3). For this purpose, we suppose that: first,aðx; tÞ ¼ a and bðx; tÞ ¼ b are constants and second: we introduce the difference operators and some lemmas.
dþt ulj ¼ ujþ1 ulj ; l
d2t ulj ¼ ujþ1 uj1 l l ; dt ulj ¼ ulj uj1 l ;
d2t ulj ¼ ujþ1 2ulj þ uj1 l l ;
d2t ulj ¼ ðdþt þ dt Þulj ; j ulj ; dþx ulj ¼ ulþ1
utjl ¼
1 þ j d u; k t l
d2t ulj ¼ ðdþt dt Þulj ;
j dx ulj ¼ ulj ul1 ;
uxjl ¼
1 þ j d u h x l
If fv i j0 6 i 6 ng and fwi j0 6 i 6 ng are two grid functions on Xh , we can define the inner product as follows: n X hv ; wi ¼ h v i wi ; i¼0
kv k ¼ hv ; v i2 ¼ 1
h
n X
v 2i
!12 :
i¼0
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
5
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
and we have
hdþx
v
; dþx wi
¼h
where n1 X
dþx
v
2
þ i dx wi
k1 ¼
i¼0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kdþx v k ¼ hdþx v ; dþx v i ¼
h
n1 X
ðdþx v l Þ
2
!12
k2 ¼
2
To prove the stability and convergence of the method the following lemmas are required (see Cui[9] and Ben-Yu [5]). Lemma 3.1. For any two mesh function fv i j0 6 i 6 ng and fwi j0 6 i 6 ng, there holds
k01
¼
ah þ 2b þ 10h 20bg þ 2
k02
¼
2
2
8abh þ ð5h þ 8gbÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ah þ 2b þ 10h 20bg
2
2
2
2
2
2
16ab þ ðah þ 2b þ 10h 20bgÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 16ab þ ðah þ 2b þ 10h 20bgÞ : 4ab ð33Þ
W 4;þ1 denotes 4-periodic Sobolev spaces [44].
n1 n1 X X ðd2x v i Þwi ¼ ðdþx v i Þðdþx wi Þ w1 dþx v 0 þ wn dx v n : i¼1
Proof. We rewrite (25) in the following form:
gd2x d2t u j þ ð1 2gÞd2x u j þ sd2t u j 2 þaksd2x u j þ hb ð1 þ ad2x Þg j ¼ 0
It is easy to see that if w0 ¼ wn ¼ 0, then n1 n1 X X ðd2x v i Þwi ¼ ðdþx v i Þðdþx wi Þ: i¼1
5h 8gb þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 8abh þ ð5h þ 8gbÞ ; 4ab qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
i¼0
i¼1
5h 8gb
ð34Þ
Multiplying (34) by ðhd2t ulj Þ and summing up for l from 0 to n yields to
i¼1
this lemma is the discrete Green formula.
hgd2x d2t u j ; d2t u j i þ ð1 2gÞhd2x u j ; d2t u j i þ shd2t u j ; d2t u j i Lemma 3.2. For any mesh function fv i j0 6 i 6 ng and defined on the interval ½L0 ; L1 with v 0 ¼ v n ¼ 0, the inequality holds:
1 þ 2 1 2 kd v k 6 kv k2 6 2 ðb aÞ kdþx v k2 : 4 x 8h
ð35Þ
2
þakshd2x u j ; d2t u j i þ hb hð1 þ ad2x Þg j ; d2t u j i ¼ 0
Using difference operators and Lemma 3.1, the first term of (35) becomes
hgd2x d2t u j ; d2t u j i ¼ ghdþx d2t u j ; dþx d2t u j i 6 gðkdþx d2t u j k2 Þ
ð36Þ
Similarly, for the second term of (35), we obtain Lemma 3.2 is analogous to Poincare inequality.
hd2x u j ; d2t u j i ¼ hdþx u j ; dþx ðujþ1 uj1 Þi
Lemma 3.3. Let xðsÞ and mðsÞ be nonnegative mesh functions. If C P 0 and mðsÞ is nondecreasing and
ð38Þ
Fourth term of (35), by the use of inequalities 2pq 6 p2 þ q2 becomes
i¼0
hd2x u j ; d2t u j i ¼ hdþx u j ; d2t u j i þ hdx u j ; d2t u j i
s
2
this lemma is the discrete Gronwall lemma. We verify the stability analysis of the scheme (25) by using the energy method described in the next theorem. Theorem 3.1. We assume that the exact solution of (1)–(3) satisfy U 2 W 4;þ1 ð0; þ1; L1 ðL0 ; L1 ÞÞ \ L1 ð0; þ1 : W 4;þ1 ðL0 ; L1 Þ. Let u be the solution of the difference scheme (25) and there exist constants ci > 0, i ¼ 1; 2 such that
jgðuÞj 6 c1 ð1 þ jujÞ @gðuÞ @u 6 c2
ð30Þ
Then for sufficiently small k there exits a constant C > 0 (independent of h and k) such that mþ1 kumþ1 k þ kum k6C t k þ kux
ð31Þ
provided that
k1 < k < k2
if g > 2bþah 4b
k01 < k < k02
if g < 2bþah 4b
2
2
2
12 ðkdþx u j k þ kd2t u j k þ kdx u j k þ kd2t u j k Þ
xðsÞ 6 mðsÞeCsk :
(
and for the third term of (35), we have
hd2t u j ; d2t u j i ¼ kd2t u j k2
j1 X xðsÞ 6 mðsÞ þ Ck xðiÞ;
then for all
ð37Þ
¼ ðhdþx u j ; dþx ujþ1 i hdþx u j ; dþx uj1 iÞ
ð32Þ
ð39Þ
For fifth term of (35) the inequalities ðp þ q þ u þ v Þ2 6 4ðp2 þ q2 þ u2 þ v 2 Þ and 2pq 6 p2 þ q2 imply that
1 2 2 hð1 þ ad2x Þg j ; d2t u j i kgðu j Þk þ kd2t u j k2 þ a2 kdþx gðu j Þk 2 1 2 þ kdþx d2t u j k2 kc1 ð1þ j u j jÞk 2 1 þ kd2t u j k2 þ a2 kc2 j dþx u j j k2 2 1 1 þ kdþx d2t u j k2 kd2t u j k2 2 2 1 þ j 2 þ kdx d2t u k þ 4c21 ð1 þ ku j k2 Þ 2 þ 4a2 ðc22 kdþx u j k2
ð40Þ
Substituting (36)–(40) into (35), we get
gkdþx d2t u j k2 þ ð1 2gÞðhdþx u j ; dþx ujþ1 i hdþx u j ; dþx uj1 iÞ þ skd2t u j k2 1 þ aksð ðkdþx u j k2 þ kdx u j k2 Þ þ kd2t u j k2 Þ 2 2 i h h1 1
kd2t u j k2 þ kdþx d2t u j k2 þ 4c21 ð1 þ ku j k2 Þ þ 4a2 ðc22 kdþx u j k2 Þ 2 b 2 ð41Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
6
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
Applying the sum j from 1 to m, we can obtain
g
m X
h ðkdþx d2t u j k2 Þ
m X
2
2b
j¼1
þ ðs þ aks
It can be proved that 2
kdþx d2t u j k2 þ ð1 2
kdþx dþt um k ¼ hdþx ðumþ1 um Þ; dþx ðumþ1 um Þi
Þðhdþx um ; dþx umþ1 iÞ
g
2
2 m m m X h X aks X kd2t u j k2 þ kdþx u j k2 þ kdx u j k2 Þ Þ ð 2b j¼1 2 j¼1 j¼1
ð42Þ Note that
P
kþ1 k i ui1 ui
¼
P
k kþ1 i ui uiþ1 ,
therefore
m h X þ ðs þ aks Þ kd2t u j k2 2b j¼1 2
þ
m m X aks X kdþx ujþ1 k2 kdþx uj1 k2 Þ ð1 2gÞhdþx u1 ; dþx u0 i ð 2 j¼1 j¼1
2 m m i X X h h 2 þ ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð43Þ
by simplifying Eq. (43), then we have m h X ðkdþx d2t u j k2 Þ þ ð1 2gÞðhdþx um ; dþx umþ1 iÞ Þ 2b j¼1 2
ðg
m h X aks kd2t u j k2 þ Þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2b j¼1 2 2
þ ðs þ aks
aks ðkdþx u0 k2 þ kdþx u1 k2 Þ 2 2 m m i X X h h 2 þ ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð1 2gÞhdþx u1 ; dþx u0 i þ
ð44Þ
ð45Þ
and similarly
2hdþx u1 ; dþx u0 i ¼ kdþx u1 k2 þ kdþx u0 k2 kdþx dþt u0 k2
ð46Þ
it follows that m h X ð1 2gÞ aks ðkdþx d2t u j k2 Þ þ ½ þ Þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2b j¼1 2 2
ðg
m ð1 2gÞ þ þ m 2 h X kdx dt u k þ ðs þ aks Þ kd2t u j k2 2 2b j¼1 2
ð1 2gÞ aks ð1 2gÞ þ þ 0 2 þ kdx dt u k ðkdþx u1 k2 þ kdþx u0 k2 Þ 2 2 2 2 m m X X þ 2 h þ ½4c21 ðm þ ku j k2 Þ þ 4a2 c22 kdx u j k ð47Þ b j¼1 j¼1
½
Now, we use Lemma 3.2 m 2h h X ð4g ðkd2t u j k2 Þ þ s þ aks Þ b 2b j¼1 2
2
ð1 2gÞ aks ð1 2gÞ þ þ m 2 þ kdx dt u k ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2 2 ð1 2gÞ aks ð1 2gÞ þ þ 0 2 þ kdx dt u k
½ ðkdþx u1 k2 þ kdþx u0 k2 Þ 2 2 2 2h m m i X X þ 2 h þ ku j k2 Þ þ 4a2 c22 kdx u j k ð48Þ 4c21 ðm þ b j¼1 j¼1 þ½
ð1 2gÞ þ þ m 2 ð1 2gÞ aks kdx dt u k þ ½ þ ðkdþx umþ1 k2 þ kdþx um k2 Þ 2 2 2 ð1 2gÞ aks ½ð1 2gÞ aks þ ð þ Þkdþt um k2 2 2 ð1 2gÞ aks þ þ ð ð50Þ Þðkdþx umþ1 k2 þ kdþx um k2 Þ 4 4
ð1 2gÞ aks þ f1 þ 2g þ aks ð Þgkdþt um k2 2 2 2 m X 5h þ f4g ðkd2t u j k2 þ s þ aksg 2b j¼1 ð1 2gÞ aks þ Þðkdþx umþ1 k2 þ kdþx um k2 Þ 4 4 ð1 2gÞ þ 2ðkdþx u1 k2 þ kdþx u0 k2 Þ
½ 2 aks 2 m m X X h ku j k2 Þ þ 4a2 c22 kdþx u j k2 þ 2ð1 2gÞkdþt u0 k2 þ ½4c21 ðm þ b j¼1 j¼1 þ ð
aks dþx u1 2 dþ u0 dþ u 0 2 1 ð1 2gÞ ¼k f 2ð þ Þðk k þ k x k2 Þ þ 2ð1 2gÞk t k2 2 2 h h k k m m X X 1 2 2 þ þ 2 ½4c21 ðm þ ku j k Þ þ 4a2 c22 kdx u j k g ð51Þ bk j¼1 j¼1 Regarding the arbitrary sufficiently small of , if 1 þ 2g þ aks > 0, i.e. g > 2bþah , it requires 4b 2
4g
With the following identity
2hdþx umþ1 ; dþx um i ¼ kdþx umþ1 k2 þ kdþx um k2 kdþx dþt um k2
Where is a sufficiently small coefficient, the equation may be rewritten as
m h X ðkdþx d2t u j k2 Þ þ ð1 2gÞðhdþx um ; dþx umþ1 iÞ Þ 2b j¼1 2
ðg
ð49Þ
Noticing to (49) and Lemma 3.2, for the second and the third terms on the left-hand side of (48), we have
2 m m i X X h h 2 ku j k2 Þ þ 4a2 c22 kdþx u j k2 4c1 ðm þ b j¼1 j¼1
ð1 2gÞhdþx u1 ; dþx u0 i þ
2
2ðkdþx umþ1 k þ kdþx um k Þ
j¼1
2
5h h þ ak > 0 2b 2bk
ð52Þ
On the other hand, we know that
kd2t u j k2 ¼ kdþt u j þ dt u j k2 6 2ðkdþt u j k2 þ kdt u j k2 Þ
ð53Þ
We need k to satisfy the following condition (for j ¼ m): 2
2
5h h 1 þ 2g þ aks þ 2 4g þ ak 2b 2bk
!
>0
ð54Þ
Then by the use of Eqs. (52) and (54), we proved the stability condition (32). 0 1 þ 1 0 We need to show that 1h dþ and 1k dþ are all x u , h dx u t u bounded. It means that, all terms on the right-hand side of (51) have upper bound. 2 n1 2 n1 X X uðxlþ1 ; t 0 Þ uðxl ; t0 Þ 2 u0 u0l 1 k dþx u0 k ¼ ð lþ1 ð Þ h Þ h h h h l¼0 l¼0 Z b 2 n1 Z xlþ1 1X @/ @/ dxÞ
j dx C ð j ¼ h l¼0 @x xl a @x
ð55Þ
Note that from (25), u1 satisfies 2
h 1 ah þ u þ ðd þ dt Þu0l þ s0l ; 2bk l 2b t s0l ¼ k7 ½/ððl þ 1ÞhÞkðb/xx ððl þ 1ÞhÞ a/x ððl þ 1ÞhÞ þ gð/ðl þ 1ÞhÞÞ þ /ððl 1ÞhÞ kðb/xx ððl 1ÞhÞ a/x ððl 1ÞhÞ þ gð/ðl 1ÞhÞÞ þ k8 ð/ððlÞhÞ kðb/xx ððlÞhÞ a/x ððlÞhÞ þ gð/ðlÞhÞÞÞ
k1 d2x dþt u0l ¼ ð1 þ 2g k1 Þd2x u0l þ
2
h ðaðgðu0lþ1 Þ þ gðu0l1 ÞÞ þ bgðu0l ÞÞ b l ¼ 1ð1Þn 1
ð56Þ
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
7
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx T
Simplifying Eq. (56), we have
ð
2bk1 2 þ 0 2b Þð1 þ 2g k1 Þd2x u0l Þdx dt ul þ dþt u0l ¼ ð ah ah h u1l u0l þ u1 þ s0l l ak
Taking inner product with
0 dþ t u
j where T j ¼ ðT 1j ; T 2j ; . . . ; T n1 Þ is the local truncation error, and we know that
2
ð57Þ
on both sides of (57), we get
2bk1 2 þ 0 þ 0 2b Þð1 þ 2g k1 Þhd2x u0 ; dþt u0 i kdþt u0 k2 ð Þhdx dt u ; dt u i ¼ ð ah ah h h hdþt u0 ; dþt u0 i þ ð1 Þhu0 ; dþt u0 i þ hdþt u0 ; s0 i ð58Þ ak ak By the use of Lemma 3.1, we can show that
2bk1 kdþt u0 k2 þ kdþx dþt u0 k 6 kCkdþt u0 k ah
ð60Þ
Therefore we have
kdþt u0 k þ 4
gd2x d2t E j þ ð1 2gÞd2x E j þ sd2t E j þ 2 2 aksd2x E j þ hb ð1 þ ad2x Þg j ¼ hb T j Multiplying (67) by ðhd2t E j Þ and taking the inner product, we get to
hgd2x d2t E j ; d2t E j i þ ð1 2gÞhd2x E j ; d2t E j i þ shd2t E j ; d2t E j i þ akshd2x E j ; d2t E j i þ
h hð1 þ ad2x Þ½gðU j Þ gðu j Þ; b
2
h j T ; d2t E j i b
ð68Þ
We know the inequalities ðp þ q þ uÞ2 3ðp2 þ q2 þ u2 Þ and pq 6 p2 þ 41 q2 ð8 > 0Þ, by applying these inequalities, imply that
hgd2x d2t E j ; d2t E j i þ ð1 2gÞhd2x E j ; d2t E j i þ shd2t E j ; d2t E j i þ akshd2x E j ; d2t E j i
1 þ 0 1 þ þ 0 d u þ d d u 6C k t k x t Provided that
ð66Þ
Now we subtract (65) from (34), and let E ¼ U u, then we have
d2t E j i ¼ h
g 1 Provided that 2bk > 0, we get k2 < 2b ah ak . Dividing the above inequality þ 0 by kdt u k, we obtain
2bk1 ah
4
2
ð59Þ
kdþt u0 k þ kdþx dþt u0 k 6 kC
2 2
jT j jl1 6 Cðk þ k h þ h Þ
ð61Þ
< 0, by the use of Lemma 3.2, we achieve
bk1 4bk1 kdþt u0 k 6 kC 1 þ kdþt u0 k 6 kC ah ah
ð62Þ
2
1 1 þ 0 1 þ þ 0 know that 1h dþ x u ¼ h dx u þ h dt dx u . It can also be proved that
ð63Þ
Regarding Eq. (63), we proved that all terms on the right-hand side of (51) are bounded. Hence from Lemma 3.2 in both cases we obtain þ mþ1 m mþ1 that kdþ k 6 kC, i.e.,kum k 6 C, under the t k þ kkux t u k þ kkdx u P þ j mþ1 0 ¼u þ m gives condition (32). The relation u j¼1 dt u , kumþ1 k 6 ku0 k þ mkC. So, the proof is completed. h Convergence of the spline difference scheme (25) by using the energy method is the issue of the next theorem.
ðg þ 12Þðkdþx d2t E j k2 Þ þ ak2s ðkdþx E j k2 þ kdx E j k2 Þ
difference scheme (25) converges to U in the discrete L2 -norm and we have
ð72Þ
By the use of Theorem 3.1 and summing up for j from 1 to m and using Lemma 3.3, we obtain (64) and the proof is completed. h 4. Numerical experiment In this section we present the numerical result of the new method on several problems. We can obtain u1 by using Taylor l about u0l and using (1): series expansion of u1 l 2
u1 ¼ u0l kðut Þ0l þ Oðk Þ l 2
¼ u0l kfbl ðuxx Þ0l a0l ðux Þ0l þ gðu0l Þg þ Oðk Þ u1 l By using initial values, we have
n o 0 u1 ¼ /ðlhÞ k bl /xx ðlhÞ a0l /x ðlhÞ þ gð/ðlhÞÞ l Example 1. We consider the following equation: [1]
mþ1 kðU uÞmþ1 k þ kðU uÞm k t k þ kðU uÞx 2 2
ð70Þ
ð71Þ
þð1 2gÞðhdþx E j ; dþx Ejþ1 i hdþx E j ; dþx Ej1 iÞ ! 2 1 h kd2t u j k2 þ s þ aks þ 2 b 2 h b a 1 6 2 3 þ C 22 kE j k2 þ kT j k2 4 b 2 h
0
Theorem 3.2. Suppose that the solution of (1)–(3) is sufficiently regular and defined by U and the assumptions of Theorem 3.1 are also valid. For sufficiently small k, the approximation solution of the spline
2
ð69Þ
Similar to Theorem 3.1, we get the following inequalities
1 Therefore we must have 1 þ 4bk 0, take k, h sufficiently small proah
2 a h ah 0 vided g > 2bk 4b , and then we have 1k dþ 6 C hence t u 1 þ þ 0 d d u 6 4 1 dþ u0 6 C. k x t k t 1 þ þ 0 0 In all cases we obtain 1k dþ þ k dx dt u 6 C. t u 1 is bounded. We For the last step, we must prove that 1h dþ xu
1 þ 1 1 þ 0 d u 6 d u þ k 1 dþ dþ u0 6 C h x h x h t x
2
h h hð1 þ ad2x Þ½gðU j Þ gðu j Þ; d2t E j i ¼ h T j ; d2t E j i b b * + h2 h2 1 kd2t E j k2 þ kT j k2 T j ; d2t E j 6 b 4 b þ
ut þ aux ¼ buxx ;
4
Cðk þ k h þ h Þ
0 6 t 6 T;
06x61
with a ¼ 3:5, b ¼ 0:022 and the following initial condition: Proof. By substituting U by u in (34), we obtain
/ðxÞ ¼ expðcxÞ
gd2x d2t U j þ ð1 2gÞd2x U j
The exact solution of above equation is taken as 2
þ sd2t U j þ aksd2x U j þ
2
h h ð1 þ ad2x Þg j ¼ T j b b
ð65Þ
uðx; tÞ ¼ expðcx þ ntÞ where c ¼ 0:02854797991928 and n ¼ 0:0999.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
8
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
Table 1 Maximum absolute error for Example 1. x
Ismail et al. [1]
0.10 0.50 0.90
t¼1
t¼2
t¼5
2.56(10) 8.93(10) 1.33(9)
2.38(10) 1.38(9) 2.83(9)
5.65(10) 1.91(9) 3.97(9)
Douglas method
0.10 0.50 0.90
t¼1
t¼2
t¼5
2.56(7) 8.37(7) 1.33(6)
2.37(7) 1.38(6) 2.82(6)
5.63(7) 1.90(6) 3.95(6)
Present method
0.10 0.50 0.90
t¼1
t¼2
t¼5
1.32898(11) 6.72048(11) 1.99724(9)
1.20256(11) 6.08158(11) 7.01701(11)
8.91154(12) 4.50668(11) 4.72124(10)
The boundary conditions can be obtained from the exact solution. We solve Example 1 with h ¼ 0:1 and k ¼ 0:001. The result are computed in different time levels. The maximum absolute errors in the solution are reported in Table 1. The computed results by our method are compared with[1]. The table shows our method is more accurate. Example 2. We consider the following equation:
ut þ aux ¼ buxx ;
t P 0;
06x61
Example 3. We consider the following model for the recolonization by oaks of Europe after the last glaciation; therefore, uðx; tÞ is the local density of oaks at time [38]
ut þ aux ¼ buxx þ ru;
t P 0;
0 6 x 6 10
If the population size at time 0 is M and is concentrated at the origin, the exact solution of this equation is
M ðx atÞ2 uðx; tÞ ¼ pffiffiffiffiffiffiffiffi exp rt 4bt 2 pbt
!
The boundary conditions can be obtained from the exact solution. In this example, we take a ¼ b ¼ 1, r ¼ 0:1 with h ¼ 0:1 and k ¼ 0:005. Computed solution for different time level has been compared with the exact solution and the absolute errors for t ¼ 1; 3; 5; 7; 9 and M ¼ 105 are reported in Table 3. We also computed solution for different M and the absolute errors for time t ¼ 2 with h ¼ 1 and k ¼ 0:01. These results are reported in Table 4. We observe that as we decrease the initial condition (M), the maximum absolute error tends to zero. Example 4. Consider the nonlinear Convection-Diffusion-Reaction equation [43]:
with the following initial condition:
ðx þ 0:5Þ2 /ðxÞ ¼ exp 0:00125
The boundary conditions can be obtained from the exact solution. For this example, we put a ¼ 1 and b ¼ 0:01. We take h ¼ 0:005 and k ¼ 0:001 and the maximum absolute error in the solution for t ¼ 0:6 and t ¼ 0:9 are computed. The results are reported in Table 2 and compared with [2]. In this problem, our numerical results are more accurate in comparison with the method in [2].
!
ut þ aux ¼ buxx þ cu2 ð1 uÞ;
t P 0;
06x61
and the exact solution is given by
The exact solution is
0:025 ðx þ 0:5 tÞ2 uðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 0:00125 þ 0:04t 0:000625 þ 0:02t
uðx; tÞ ¼ ð1 þ expðAðx BtÞ þ CÞÞ1
!
where A ¼
qffiffiffiffi c
, B ¼ 2a þ
4b
pffiffiffiffiffiffi cb and C ¼ AðB 1Þ and the initial con-
dition /ðxÞ ¼ ð1 þ expðAx þ CÞÞ1 .
Table 2 Maximum absolute error for Example 2. x
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
LBT [2]
QBT [2]
Present method
t ¼ 0:6
t ¼ 0:9
t ¼ 0:6
t ¼ 0:9
t ¼ 0:6
t ¼ 0:9
3.1(3) 1.3(3) 8.1(4) 1.7(5) 9.2(5) – – – –
0 1.1(3) 6.8(4) 0.5(6) 6.6(5) – – – –
1.3(4) 6.2(4) 3.8(5) 3.5(4) 4.5(5) – – – –
3.2(4) 1.02(4) 1.2(5) 3.2(4) 2.2(5) – – – –
3.90549(6) 1.76112(5) 5.19360(5) 1.78248(5) 8.85773(5) 2.03339(5) 3.24045(5) 2.06487(5) 4.53367(6)
1.11127(5) 6.67626(5) 1.47407(5) 3.15751(4) 7.0042(6) 5.00361(7) 1.37155(8) 1.54961(10) 7.51235(13)
Table 3 Maximum absolute error for Example 3. x
t¼1
t¼3
t¼5
t¼7
t¼9
1 2 3 4 5 6 7 8 9
2.54383(6) 2.65653(6) 2.75816(6) 2.84930(6) 2.92874(6) 2.99538(6) 3.04828(6) 3.08690(6) 3.11415(6)
1.09811(6) 1.14358(6) 1.19711(6) 1.25080(6) 1.30442(6) 1.35750(6) 1.40931(6) 1.45093(6) 1.49247(6)
6.93872(7) 6.5356(7) 6.81864(7) 7.08063(7) 7.29557(7) 7.41553(7) 7.35247(7) 6.93160(7) 5.98364(7)
4.29387(7) 4.02496(7) 4.12376(7) 4.15886(7) 4.05557(7) 3.68950(7) 2.83034(7) 1.49359(7) 3.35850(7)
4.24303(7) 2.56864(7) 2.53775(7) 2.42229(7) 2.10149(7) 1.39552(7) 2.00475(8) 3.05948(7) 1.52747(6)
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx Table 4 Maximum absolute error for Example 3. x
M ¼ 103
M ¼ 105
M ¼ 107
1 2 3 4 5 6 7 8 9
8.18795(4) 1.98019(3) 3.57075(3) 6.76625(3) 8.77917(4) 3.73394(3) 8.60632(4) 5.83540(4) 3.82148(4)
8.19047(6) 1.98074(5) 3.5717(5) 6.76608(5) 8.77694(6) 3.73414(5) 8.60557(6) 5.83505(6) 3.82128(6)
8.1905(8) 1.98075(7) 3.57118(7) 6.76608(7) 8.77692(8) 3.73414(7) 8.60556(8) 5.83504(8) 3.82128(8)
Table 5 Maximum absolute error for Example 4. x
c¼1 b ¼ 0:001
c ¼ 10 b ¼ 0:01
c ¼ 100
0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
2.69954(8) 5.47176(9) 1.12576(9) 2.31615(10) 4.76527(11) 7.14614(12) 1.47025(12) 3.0249(13) 6.22346(14)
1.6958(8) 3.48895(9) 7.17818(10) 1.47684(10) 3.03847(11) 6.25136(12) 1.28617(12) 2.64531(13) 5.5412(14)
1.88386(10) 3.87587(11) 7.92424(12) 1.64063(12) 3.37543(13) 6.94464(14) 1.42879(14) 2.93948(15) 6.00635(16)
c¼1 a ¼ 100
c ¼ 100
c ¼ 1000
a¼5
a¼1
5.41124(8) 3.53667(7) 7.22626(7) 1.20347(6) 1.71858(6) 2.09293(6) 2.14112(6) 1.7592(6) 9.83998(7)
1.38435(3) 2.85513(4) 5.86763(5) 1.02727(5) 2.48386(6) 5.11032(7) 1.0514(7) 2.16311(8) 4.42828(9)
3.05902(7) 2.06115(9) 1.38879(11) 9.35762(14) 6.30512(16) 4.24835(18) 2.86252(20) 7.09547(23) 7.87728(25)
b ¼ 0:1
Table 6 Maximum absolute error for Example 4. x 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90
The boundary conditions can be obtained from the exact solution. The nonlinear terms can be approximated non locally in many different ways. In this example, we consider an approximate of the j
j
2
j j j j form ðu2 Þl ¼ 13 ðulþ1 þ ulj þ ul1 Þulj and ðu3 Þl ¼ 12 ðul1 þ ulþ1 Þðulj Þ for nonlinear terms. In our numerical computation, in Table 5, we take a ¼ 10, h ¼ k ¼ 0:01. The results are computed for different c and b in time t ¼ 0:9. In Table 6, we take b ¼ 0:1, h ¼ k ¼ 0:01. The absolute errors for different c and a are reported in Table 6 in time t ¼ 0:9.
5. Conclusion In this study a numerical scheme has been established to approximate solution of Convection-Reaction-Diffusion equation. This scheme is based on non-polynomial tension spline in x direction, finite difference in t direction and nonstandard finite difference in nonlinear terms. Theoretical formulation of the proposed method was explained in detail. The optimal values of spline parameters and truncation error have been determined. The stability and convergence analyses of the present numerical approach by energy method, have been described. The performance of the method for these problems was measured by calculating the abso-
9
lute errors for different time levels which are reported in tabular form. In our proposed method, we developed the unconditionally three level implicit scheme. Our proposed method yields the tridiagonal system, which can be solved by efficient algorithm. The accuracy and applicability of our method have been tested by four examples. The numerical result shows that our scheme performs more accurately than the previous one. References [1] Ismail HNA, Elbarbary EM, Salem GSE. Restrictive Taylor approximation for solving convection-diffusion equation. Appl Math Comput 2004;147:355–63. [2] Dhawan S, Kapoor S, Kumar S. Numerical method for advection diffusion equation using FEM and B-splines. J Comput Sci 2012;3:429–37. [3] Makungu M, Haario H, Mahera WC. A generalised 1-dimensional particle transport method for convection diffusion reaction. Afr Math Un 2012;23:21–39. [4] Morton KW. Numerical solution of convection-diffusion problems. London: Chapman and Hall; 1996. [5] Ben-Yu G, Pascual PJ, Rodriguez MJ, Vazquez L. Numerical solution of the sineGordon equation. Appl Math Comput 1986;18(1):1–14. [6] Richards LA. Capillarity conduction of liquid through porous medium. Physica 1931;1:318. [7] Zeldovich Ya. Theory of flame propagation. Nat Adv Comm Aeronaut Tech Memor 1951;1282:39. [8] Kolmogorov A, Petrovskii I, Piskunoff N. In: Pelce P, editor. Dynamics of Curved Fronts. Boston: Academic Press; 1988. p. 105–30. [9] Cui M. Fourth-order compact scheme for the one-dimensional sine-Gorden equation. Numer Meth Partl Diff Eq 2009;25:685–711. [10] Murray JD. Mathematical biology. Berlin: Springer; 1989. [11] Witelski TP. Segregation and mixing in degenerate diffusion in population dynamics. J Math Biol 1997;35:695–712. [12] Wang YM. A modified accelerated monotone iterative method for finite difference reaction-diffusion-convection equations. J Comput Appl Math 2011;235:3646–60. [13] Aronson DG, Weinberger HF. Nonlinear diffusion in population genetics, combustion, and nerve propagation. In: Partial differential equations and related topics. Lecture notes in mathematics, vol. 446. New York: Springer; 1975. p. 5–49. [14] Aris R. The mathematical theory of diffusion and reaction in permeable catalysts. Clarendon (London): Oxford University Press; 1975. [15] Hetrick DK. Dynamics of nuclear reactors. Chicago: University of Chicago; 1971. [16] Isenberg J, Gutfinger C. Heat transfer to a draining film. Int J Heat Transfer 1972;16:505–12. [17] Wilmott P, Dewynne J, Howison S. Option pricing, mathematical models and computation. Oxford: Oxford Financial Press; 1993. [18] Roache P. Computational fluid dynamics. Albuquerque (NM): Hermosa Press; 1972. [19] Cengel YA, Cimbala JM. Fluid mechanics: fundamental and application. McGraw-Hill Higher Education; 2006. [20] McRea GJ, Goodin WR, Seinfeld JH. Numerical solution of atmospheric diffusion for chemically reacting flows. J Comput Phys 1982;77:1–42. [21] Cunningham AB, Charaklis WG, Abedeen F, Crawford D. influence of biofilm accumulation on porous media hydrodynamics. Environ Sci Tech 1991;25:1305–10. [22] Harten A, Enquist B, Osher S, Chakravarthy S. Uniformly high order essentially non-oscillatory schemes III. J Comput Phys 1987;71:231–303. [23] Liu XD, Osher S, Chan T. Weighted essentially non-oscillatory schemes. J Comput Phys 1994;115:200–12. [24] John V, Novo J. On (essentially) non-oscillatory discretizations of evolutionary convection-diffusion equations. J Comput Phys 2012;231:1570–86. [25] Kailash PC. On the use of nonstandard finite difference methods. J Diff Eq Appl 2005;11:735–58. [26] Kojouharov HV, Chen BM. Nonstandard methods for the convective-dispersive transport equation with nonlinear reactions. John Wiley and Sons; 1999. p. 617–24. [27] Kaya A. Finite difference approximations of multidimensional unsteady convection-diffusion-reaction equations. J Comput Phys 2015;285:331–49. [28] Kaya A. A finite difference scheme for multidimensional convection-diffusionreaction equations. Comput Meth Appl Mech Eng 2014;278:347–60. [29] Brooks AN, Hughes TJR. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [30] Bause M. Stabilized finite element methods with shock-capturing for nonlinear convection-diffusion-reaction models. Int Numer Math Adv Appl 2009:125–34. [31] Bause M, Schwegler K. Analysis of stabilized higher-order finite element approximation of nonstationary and nonlinear convection-diffusion-reaction equations. Comput Meth Appl Mech Eng 2012;209-212:184–96. [32] Bause M, Schwegler K. Higher order finite element approximation of systems of convection diffusion reaction equations with small diffusion. J Comput Appl Math 2013;246:52–64.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005
10
H.S. Shekarabi, J. Rashidinia / Ain Shams Engineering Journal xxx (2016) xxx–xxx
[33] Yucel H, Stoll M, Benner P. Discontinuous Galerkin finite element methods with shock-capturing for nonlinear convection dominated models. Comp Chem Eng 2013;58:278–87. [34] Uzunca M, Karasozen B, Manguoglu M. Adaptive discontinuous Galerkin methods for non-linear diffusion-convection-reaction equations. Comp Chem Eng 2014;68:24–37. [35] Brezzi F, Hauke G, Marini LD, Sangalli G. Link-cutting bubbles for the stabilization of convection-diffusion-reaction problems. Math Models Meth Appl Sci 2003;13:445–61. [36] Sendur A, Neslitrk AI. Applications of the pseudo residual-free bubbles to the stabilization of convection-diffusion-reaction problems. Calcolo 2011;49:1–19. [37] Quang DA, Ehrhardt M. Adequate numerical solution of air pollution problems by positive difference schemes on unbounded domains. Math Comput Model 2006;44:834–56. [38] Chen-Charpentier BM, Kojouharov HV. An unconditionally positivity preserving scheme for advection-diffusion-reaction equations. Math Comp Model 2013;57:2177–85. [39] Pao CV. Monotone iterative methods for finite difference system of reaction diffusion equations. Numer Math 1985;46:571–86. [40] Rashidinia J, Mohammadi R. Tension spline approach for the numerical solution of nonlinear KleinGordon equation. Comp Phys Commun 2010;181:78–91. [41] Aghamohammadi M, Rashidinia J, Ezzati R. Tension spline method for solution of non-linear Fisher equation. Appl Math Comput 2014;249:399–407. [42] Mickens RE. Nonstandard finite difference model of differential equation. Singapore: World Scientific; 1994. [43] Hundsdorfer W. Numerical solution of advection-diffusion-reaction equations. Lecture notes. Thomas Stieltjes Institute; 2000. [44] Atkinson K, Han W. Theoretical numerical analysis: a functional analysis framework. Springer; 2005. [45] Theeraek P, Phongthanapanich S, Dechaumphai P. Solving convectiondiffusion-reaction equation by adaptive finite volume element method. Math Comp Simul 2011;82:220–33.
[46] Duan H, Hsieh P, Tan R, Yang S. Analysis of a new stabilized finite element method for the reaction-convection-diffusion equations with a large reaction coefficient. Comput Meth Appl Mech Eng 2012;247–248:15–36. [47] Hsieh P, Yang S. A new stabilized linear finite element method for solving reaction-convection-diffusion equations. Comput Meth Appl Mech Eng 2016;307:362–82. [48] Gierea S, Iliescub T, Johna V, Wellsb D. SUPG reduced order models for convection-dominated convection-diffusion-reaction equations. Comput Meth Appl Mech Eng 2015;289:454–74. [49] Cui M. Compact exponential scheme for the time fractional convectiondiffusion-reaction equation with variable coefficients. J Comput Phys 2015;280:143–63.
H.S. Shekarabi is Ph.D. student of College of Basic Science, Karaj Branch, Islamic Azad University.
Dr. J. Rashidinia: For photography and biography see Ain Shams Engineering Journal, Volume 6, Issue 1, March 2015, Pages 381–389.
Please cite this article in press as: Shekarabi HS, Rashidinia J. Three level implicit tension spline scheme for solution of Convection-Reaction-Diffusion equation. Ain Shams Eng J (2016), http://dx.doi.org/10.1016/j.asej.2016.10.005