- Email: [email protected]
A second order numerical method for solving advection-diffusion models
Mathematical and Computer Modelling 50 (2009) 806–811
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
A second order numerical method for solving advection-diffusion modelsI R. Company, E. Ponsoda, J.-V. Romero ∗ , M.-D. Roselló Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, Camino de Vera s/n, Edificio 8G, 2◦ , 46022 Valencia, Spain
article
info
Article history: Received 24 September 2008 Accepted 12 December 2008 Keywords: Advection-diffusion equation CE-SE numerical scheme Second order Taylor’s expansion Accuracy
abstract The space–time conservation element and solution element (CE-SE) scheme is a method that improves the well-established methods, like finite differences or finite elements: the integral form of the problem exploits the physical properties of conservation of flow, unlike the differential form. Also, this explicit scheme evaluates the variable and its derivative simultaneously in each knot of the partitioned domain. The CE-SE method can be used for solving the advection-diffusion equation. In this paper, a new numerical method for solving the advection-diffusion equation, based in the CE-SE method is developed. This method increases the spatial precision and it is validated with an analytical solution. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction In [1], the CE-SE numerical method is applied to advection-diffusion problems with constant coefficients
∂ ∂ ∂2 u(x, t ) + a u(x, t ) − µ 2 u(x, t ) = 0, ∂t ∂x ∂x (x, t ) ∈ R × [0, +∞[ ; µ ≥ 0,
(1)
and the resultant a − µ scheme presents many advantages such as behavior and stability if it is compared with other numerical methods, like finite differences or finite elements, for example. Usually, the CE-SE method has been applied to solve problems of conservation laws [2], and also in advection-diffusion problems with constant coefficients [1]. Recently, a modified CE-SE method has been proposed in [3] for advection-diffusion problems with non-autonomous coefficients. In all of these cases, this method presents significants advantages because the integral form of the problem exploits the physical properties of conservation of flow, unlike the differential form. Also, this explicit scheme evaluates the variable and its derivative simultaneously in each knot of the partitioned domain, see [3] for details. Let us consider problem (1) under the initial condition u(x, 0) = f (x).
(2)
Applying the CE-SE numerical scheme, the space-time domain is partitioned in a grid such that the knots (j, n) are obtained for n = 0, 1/2, 1, 3/2, . . . and, for each n, j = n, n ± 1, n ± 2, n ± 3, . . . see [2] for details. Then, we define the solution
I This work has been partially supported by the Generalitat Valenciana grants GV/2007/009 and GVPRE/2008/092 and the Spanish M.C.Y.T. and FEDER grant TRA2007–68006–C02–02. ∗ Corresponding author. E-mail addresses: [email protected] (R. Company), [email protected] (E. Ponsoda), [email protected] (J.-V. Romero), [email protected] (M.-D. Roselló).
0895-7177/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2009.05.009
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
807
Fig. 1. Conservation elements of CE-SE method.
element SE (j, n) as the space-time region enclosed inside the rhombus centered in xj , t n and whose diagonals are ∆t and ∆x, see [3] for details. Using the standard CE-SE method, in each solution element, we define
U (x, t ; j, n) = Ujn + (Ux )nj x − xj + (Ut )nj t − t n , ∀(x, t ) ∈ SE (j, n),
(3)
where Ujn , (Ux )nj , (Ut )nj are constants to be determined into SE (j, n), see [3]. In this way, we define the conservation elements, CE± (j, n), as in [2] to evaluate these constants. In this work we propose a second order Taylor expansion in Eq. (3) in order to obtain a better approximate numerical solution of system (1)–(2) in the knots of the grid. Thus, a new conservation elements must be defined to be able to evaluate the new constants that appear in the larger order development. In this sense, the aim of this paper is to check, testing by an illustrative example, the behavior of this modified CE-SE method. This work is organized as follows. In Section 2, the second order CE-SE method to solve (1)–(2) is presented. We obtain an explicit expression such that the approximate values of the solution u(x, t ) of problem (1)–(2) can be evaluated in a set of knots of a certain partition of space-time domain. Finally, in Section 3, an illustrative example is presented such that the quality on the approximate solution given by the second order CE-SE method is compared with the answer obtained by means of the standard first order CE-SE method. 2. The numerical scheme In a similar way from the a − µ scheme, [1], in each solution element we define U (x, t ; j, n) = Ujn + (Ux )nj x − xj + (Ut )nj t − t n +
∀(x, t ) ∈ SE (j, n),
1 2
(Uxx )nj x − xj
2 ) ,
(4)
where Ujn , (Ux )nj , (Ut )nj and (Uxx )nj are constants to be determined into SE (j, n). To calculate these constants the so called elements of conservation CE (j, n) are in use, which exploit the conservation of the flow in certain space-time regions that we will describe later. It is for it that turns out to be important the integral formulation of problem Eq. (1), whose proof can be found in [3]. In this way, in the solution elements we define G (x, t ; j, n) = (G1 (j, n), G2 (j, n)) = −U (x, t ; j, n), aU (x, t ; j, n) − µ (Ux )nj − µ (Uxx )nj (x − xj ) , ∀(x, t ) ∈ SE (j, n).
(5)
In the usual first order CE-SE method, only two constants are to be determined, then two conservation elements are necessary, see Fig. 1. Now, four constants have to be determined, Ujn , (Ut )nj , (Ux )nj and (Uxx )nj , then we need to introduce two new conservation elements. In this way, we define the conservation elements CEi (j, n), 1 ≤ i ≤ 4. CE1 (j, n) and CE2 (j, n) are the rectangular regions such that the knot (xj , t n ) is located in the right top corner, and whose sides have a length ∆x/2, ∆t /2 for CE1 (j, n) and ∆x/4, ∆t /2 for CE2 (j, n). The elements CE3 (j, n) and CE4 (j, n) are defined in the same way, but the knot (xj , t n ) is located now in the left top corner, see Fig. 2. Note that CE1 (j, n) and CE4 (j, n) are the same elements CE− (j, n) and CE+ (j, n), respectively, from the usual first order CE-SE method.
808
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
Fig. 2. Conservation elements of modified CE-SE method.
Fig. 3. Conservation elements and solution elements of modified CE-SE method.
In order to calculate the constants Ujn , (Ut )nj , (Ux )nj and (Uxx )nj in (4), we use the following approximation of the integral equation Fi (j, n) =
I S (CEi (j,n))
I = S (CEi (j,n))
G (x, t ; j, n) dr G1 (j, n) dx + G2 (j, n) dt = 0,
(6)
where CEi (j, n) is the conservation element and G (x, t ; j, n) is given by (5). Note that for evaluating the integral in the first conservation element, from (xj , t n ) to (xj−1/2 , t n ) the variable t does not change, and the solution element is the SE (j, n), but from (xj−1/2 , t n ) to (xj−1/2 , t n−1/2 ) the solution element is SE (j − 1/2, n − 1/2) and now, the variable x is constant. Thus, the first integral needs to compute four different integrals xj−1/2
Z F1 =
G1 (j, n)dx +
Z
xj
t n−1/2
G2 (j − 1/2, n − 1/2)dt +
tn
Z
xj
G1 (j − 1/2, n − 1/2)dx +
xj−1/2
Z
tn
t n−1/2
G2 (j, n)dt ,
(7)
where Gi (j, n), i = 1, 2, are defined in (5). Similar arguments are necessary to compute the others integrals, see Figs. 2 and 3, xj−1/4
Z F2 =
G1 (j, n)dx +
xj
Z
xj
+ xj−1/4
Z
t n−1/4 tn
G2 (j, n)dt +
G1 (j − 1/2, n − 1/2)dx +
Z
Z
t n−1/2
t n−1/4
G2 (j − 1/2, n − 1/2)dt
tn t n−1/2
G2 (j, n)dt ,
(8)
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811 t n−1/2
Z
G2 (j, n)dt +
F3 = tn
t n−1/2
G2 (j + 1/2, n − 1/2)dt +
t n−1/2
Z
G2 (j, n)dt +
F4 = tn
xj+1/2
Z
tn
Z
t n−1/4
G2 (j, n)dt +
Z
xj
G1 (j, n)dx,
(9)
xj+1/4
G1 (j + 1/2, n − 1/2)dx
xj
tn
Z
G1 (j + 1/2, n − 1/2)dx
xj
t n−1/4
Z +
xj+1/4
Z
809
+ t n−1/2
G2 (j + 1/2, n − 1/2)dt +
Z
xj
G1 (j, n)dx.
(10)
xj+1/2
We define
Ujn
1 ∆x n 2! 2 (Ux )j 2 q(j, n) = 1 ∆x (Uxx )nj 3! 2 3 1 ∆x (Ut )nj 4! 2
,
(11)
and, in this way, the initial condition
u(xj , 0)
∆x ∂ u(x , 0) 4 ∂x j q(j, 0) = ∆x ∂ 2 . u ( x , 0 ) j 24 ∂ x2 ∆x ∂ u(xj , 0) 192 ∂ t From equations Fi (j, n) = 0, 1 ≤ i ≤ 4 , given by Eqs. (7)–(10), and expressing them the matrix form, we have q(j, n + 1) = (Q+ )2 q(j − 1, n) + [Q+ Q− + Q− Q+ ] q(j, n) + (Q− )2 q(j + 1, n),
(12)
(13)
where q(j, n) is given by (11) and Q+ and Q− are the 4 × 4 matrix Q+ =
1 2ν − 3 2 + 5ξ + 2ξ 2
2
q+1
q+2
q+3
q+4
q−1
q−2
q−3
q−4 ,
(14)
and Q− =
1 2ν 2 − 3 2 + 5ξ + 2ξ 2
(15)
where
ν=a
q+1
q+2
∆t , ∆x 1
ξ = 4µ
∆t , (∆x)2
(1 + ν) −4ν + 6ν 2 − 3 2 + 5ξ + 2ξ 2 2 3 − −1 + ν 2 ( 1 + 2 ξ ) , 2 = −2 ν −1 + ν 2 2 2 2 µ −1 + ν −1 + 2ν − 2ξ 2aξ 2 9ξ 2 11ξ −5 − 4ν + ν 2 (7 − 3ξ ) − + + 3ξ 3 2 2 3 15ξ 2 2 ν − ξ ( 1 + 2 ξ ) + ν − 4 − + 3 ξ 2 2 = 2 + 4ν + 2ν 2 (−3 + ξ ) + ξ µ2 −3 + 6ν 2 − 6ξ + ν 3 (−7 + 3ξ ) + ν 1 − 6ξ − 9ξ 2 − 3aξ 2
,
810
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
q+3
q+4
q−1
q−2
q−3
q−4
6ν 2 − 3(2 + ξ ) −
1
ν(4 + 15ξ )
2 3 2 − 1 + ν + 4ξ + 12ξ 2 − 4ν(1 + 3ξ ) 4 , = 15ξ 9ξ 2 2 − 3ξ + ν 2 + 3 − 5ν − 2 2 µ2 5 + 8ν 3 + 7ξ − 12ξ 2 − 12ξ 3 − ν 2 (5 + 11ξ ) + 4ν −1 + 3ξ + 3ξ 2 4aξ 2 2 2 3aξ 2 + 8ν − 6ν + 13ξ + 6ξ 2 − 2µ2 2 3a ξ (− 3 − 6 ξ + ν(2 + 3ξ )) − 2 µ = , 2 2 3aξ 2 − 4ν + 2ν + ξ − 2 µ 2 2 3 aξ 10ν − 6ν + 3ξ (1 + 2ξ ) + ν(−2 + 3ξ ) − 2aξ 2 1 − (−1 + ν) 4ν + 6ν 2 − 3 2 + 5ξ + 2ξ 2 2 3 2 − 1 + ν ( 1 + 2 ξ ) 2 = , 2 2ν −1 + ν µ2 −1 + ν 2 −1 + 2ν 2 − 2ξ − 2 2aξ 1 10 − 8ν + 11ξ − 9ξ 2 − 6ξ 3 + 2ν 2 (−7 + 3ξ ) 2 1 2 2 2ν − 3ξ (1 + 2ξ ) + ν 8 + 15ξ − 6ξ , 2 = 2 −2 + 4ν − 2ν (−3 + ξ ) − ξ µ2 3 − 6ν 2 + 6ξ + ν 3 (−7 + 3ξ ) + ν 1 − 6ξ − 9ξ 2 3aξ 2 15ξ 2 6ν − 3(2 + ξ ) + ν 2 + 2 3 2 2 − −1 + ν + 4ξ + 12ξ + 4ν(1 + 3ξ ) 4 = , 1 2 2 6 − 4ν − 10ν − 9ξ − 15νξ − 6ξ 2 µ2 −5 + 8ν 3 − 7ξ + 12ξ 2 + 12ξ 3 + ν 2 (5 + 11ξ ) + 4ν −1 + 3ξ + 3ξ 2 4aξ 2 2 2 3aξ 2 − 8ν − 6ν + 13ξ + 6ξ 2 2µ2 3aξ 2 (3 + 6ξ + ν(2 + 3ξ )) 2 µ = . 2 2 3aξ 2 + 4ν + 2ν + ξ 2 µ aξ 2 10ν 2 + 6ν 3 + ν(2 − 3ξ ) + 3ξ (1 + 2ξ ) − 2aξ 2
3. Numerical example In order to test the proposed scheme we consider an analytical solution. We choose the one-dimensional advectiondiffusion equation Eq. (1) of a Gaussian pulse of unit height, centered at x0 = 1 in a region bounded by 0 ≤ x ≤ 9 whose analytical solution [4] is 1 (x − x0 − at )2 ua (x, t ) = √ exp − µ(4t + 1) 4t + 1
.
(16)
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
811
Table 1 Errors in t = 5, x = 5.
∆x = 0.04 ∆x = 0.02
CE-SE
Second order CE-SE
0.0027 0.0011
0.0014 0.00033
u 0.25
u 0.25
0.20
0.20
0.15
0.15
0.10
0.10
0.05
0.05
4.0
4.5
5.0
5.5
x 6.0
4.0
4.5
5.0
5.5
x 6.0
Fig. 4. Analytical and numerical solution for ∆x = 0.4, t = 5.
The initial condition is given by
(x − x0 )2 u(x, 0) = exp − µ
and the boundary condition at the two ends at any time t is obtained by substituting x = 0 and x = 9, respectively, in (16). The values of the parameters used are ν = 0.005 m2 /s and a = 0.8 m/s. In this section we compare the numerical method presented in this article with the original CE-SE method. We compare the errors of two methods (Table 1) with ν = 0.8 and ∆x = 0.4. We can observe that the error of the modified CE-SE method improves the one obtained with CE-SE method. In the left side of Fig. 4 is represented the analytical solution and numerical solution for the CE-SE scheme at time t = 5 with ∆x = 0.4. On the right side of Fig. 4, the analytical solution is represented versus the proposed CE-SE scheme at the same time considering ∆x = 0.4. We can see that the approximation solution obtained with the second order CE-SE is better than that obtained with the traditional CE-SE scheme. References [1] X.Y. Wang, C.Y. Chow, S.C. Chang, Application of the space-time conservation element and solution element method to one-dimensional advectiondifussion problems, NASA-TM-1999-209068 (1999). [2] S.C. Chang, The method of space-time conservation element and solution element. A new approach for solving the Navier–Stokes and Euler equations, J. Comput. Phys. 119 (1995) 295–324. [3] E. Ponsoda, E. Defez, M.D. Roselló, J.V. Romero, A stable numerical method for solving for solving variable coefficient advection-diffusion models, Comput. Math. Appl. 56 (2008) 754–768. [4] H. Karahan, A third-order upwind scheme for the advection-diffusion equation using spreadsheets, Adv. Eng. Softw. 38 (2007) 688–697.
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
A second order numerical method for solving advection-diffusion modelsI R. Company, E. Ponsoda, J.-V. Romero ∗ , M.-D. Roselló Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, Camino de Vera s/n, Edificio 8G, 2◦ , 46022 Valencia, Spain
article
info
Article history: Received 24 September 2008 Accepted 12 December 2008 Keywords: Advection-diffusion equation CE-SE numerical scheme Second order Taylor’s expansion Accuracy
abstract The space–time conservation element and solution element (CE-SE) scheme is a method that improves the well-established methods, like finite differences or finite elements: the integral form of the problem exploits the physical properties of conservation of flow, unlike the differential form. Also, this explicit scheme evaluates the variable and its derivative simultaneously in each knot of the partitioned domain. The CE-SE method can be used for solving the advection-diffusion equation. In this paper, a new numerical method for solving the advection-diffusion equation, based in the CE-SE method is developed. This method increases the spatial precision and it is validated with an analytical solution. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction In [1], the CE-SE numerical method is applied to advection-diffusion problems with constant coefficients
∂ ∂ ∂2 u(x, t ) + a u(x, t ) − µ 2 u(x, t ) = 0, ∂t ∂x ∂x (x, t ) ∈ R × [0, +∞[ ; µ ≥ 0,
(1)
and the resultant a − µ scheme presents many advantages such as behavior and stability if it is compared with other numerical methods, like finite differences or finite elements, for example. Usually, the CE-SE method has been applied to solve problems of conservation laws [2], and also in advection-diffusion problems with constant coefficients [1]. Recently, a modified CE-SE method has been proposed in [3] for advection-diffusion problems with non-autonomous coefficients. In all of these cases, this method presents significants advantages because the integral form of the problem exploits the physical properties of conservation of flow, unlike the differential form. Also, this explicit scheme evaluates the variable and its derivative simultaneously in each knot of the partitioned domain, see [3] for details. Let us consider problem (1) under the initial condition u(x, 0) = f (x).
(2)
Applying the CE-SE numerical scheme, the space-time domain is partitioned in a grid such that the knots (j, n) are obtained for n = 0, 1/2, 1, 3/2, . . . and, for each n, j = n, n ± 1, n ± 2, n ± 3, . . . see [2] for details. Then, we define the solution
I This work has been partially supported by the Generalitat Valenciana grants GV/2007/009 and GVPRE/2008/092 and the Spanish M.C.Y.T. and FEDER grant TRA2007–68006–C02–02. ∗ Corresponding author. E-mail addresses: [email protected] (R. Company), [email protected] (E. Ponsoda), [email protected] (J.-V. Romero), [email protected] (M.-D. Roselló).
0895-7177/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2009.05.009
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
807
Fig. 1. Conservation elements of CE-SE method.
element SE (j, n) as the space-time region enclosed inside the rhombus centered in xj , t n and whose diagonals are ∆t and ∆x, see [3] for details. Using the standard CE-SE method, in each solution element, we define
U (x, t ; j, n) = Ujn + (Ux )nj x − xj + (Ut )nj t − t n , ∀(x, t ) ∈ SE (j, n),
(3)
where Ujn , (Ux )nj , (Ut )nj are constants to be determined into SE (j, n), see [3]. In this way, we define the conservation elements, CE± (j, n), as in [2] to evaluate these constants. In this work we propose a second order Taylor expansion in Eq. (3) in order to obtain a better approximate numerical solution of system (1)–(2) in the knots of the grid. Thus, a new conservation elements must be defined to be able to evaluate the new constants that appear in the larger order development. In this sense, the aim of this paper is to check, testing by an illustrative example, the behavior of this modified CE-SE method. This work is organized as follows. In Section 2, the second order CE-SE method to solve (1)–(2) is presented. We obtain an explicit expression such that the approximate values of the solution u(x, t ) of problem (1)–(2) can be evaluated in a set of knots of a certain partition of space-time domain. Finally, in Section 3, an illustrative example is presented such that the quality on the approximate solution given by the second order CE-SE method is compared with the answer obtained by means of the standard first order CE-SE method. 2. The numerical scheme In a similar way from the a − µ scheme, [1], in each solution element we define U (x, t ; j, n) = Ujn + (Ux )nj x − xj + (Ut )nj t − t n +
∀(x, t ) ∈ SE (j, n),
1 2
(Uxx )nj x − xj
2 ) ,
(4)
where Ujn , (Ux )nj , (Ut )nj and (Uxx )nj are constants to be determined into SE (j, n). To calculate these constants the so called elements of conservation CE (j, n) are in use, which exploit the conservation of the flow in certain space-time regions that we will describe later. It is for it that turns out to be important the integral formulation of problem Eq. (1), whose proof can be found in [3]. In this way, in the solution elements we define G (x, t ; j, n) = (G1 (j, n), G2 (j, n)) = −U (x, t ; j, n), aU (x, t ; j, n) − µ (Ux )nj − µ (Uxx )nj (x − xj ) , ∀(x, t ) ∈ SE (j, n).
(5)
In the usual first order CE-SE method, only two constants are to be determined, then two conservation elements are necessary, see Fig. 1. Now, four constants have to be determined, Ujn , (Ut )nj , (Ux )nj and (Uxx )nj , then we need to introduce two new conservation elements. In this way, we define the conservation elements CEi (j, n), 1 ≤ i ≤ 4. CE1 (j, n) and CE2 (j, n) are the rectangular regions such that the knot (xj , t n ) is located in the right top corner, and whose sides have a length ∆x/2, ∆t /2 for CE1 (j, n) and ∆x/4, ∆t /2 for CE2 (j, n). The elements CE3 (j, n) and CE4 (j, n) are defined in the same way, but the knot (xj , t n ) is located now in the left top corner, see Fig. 2. Note that CE1 (j, n) and CE4 (j, n) are the same elements CE− (j, n) and CE+ (j, n), respectively, from the usual first order CE-SE method.
808
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
Fig. 2. Conservation elements of modified CE-SE method.
Fig. 3. Conservation elements and solution elements of modified CE-SE method.
In order to calculate the constants Ujn , (Ut )nj , (Ux )nj and (Uxx )nj in (4), we use the following approximation of the integral equation Fi (j, n) =
I S (CEi (j,n))
I = S (CEi (j,n))
G (x, t ; j, n) dr G1 (j, n) dx + G2 (j, n) dt = 0,
(6)
where CEi (j, n) is the conservation element and G (x, t ; j, n) is given by (5). Note that for evaluating the integral in the first conservation element, from (xj , t n ) to (xj−1/2 , t n ) the variable t does not change, and the solution element is the SE (j, n), but from (xj−1/2 , t n ) to (xj−1/2 , t n−1/2 ) the solution element is SE (j − 1/2, n − 1/2) and now, the variable x is constant. Thus, the first integral needs to compute four different integrals xj−1/2
Z F1 =
G1 (j, n)dx +
Z
xj
t n−1/2
G2 (j − 1/2, n − 1/2)dt +
tn
Z
xj
G1 (j − 1/2, n − 1/2)dx +
xj−1/2
Z
tn
t n−1/2
G2 (j, n)dt ,
(7)
where Gi (j, n), i = 1, 2, are defined in (5). Similar arguments are necessary to compute the others integrals, see Figs. 2 and 3, xj−1/4
Z F2 =
G1 (j, n)dx +
xj
Z
xj
+ xj−1/4
Z
t n−1/4 tn
G2 (j, n)dt +
G1 (j − 1/2, n − 1/2)dx +
Z
Z
t n−1/2
t n−1/4
G2 (j − 1/2, n − 1/2)dt
tn t n−1/2
G2 (j, n)dt ,
(8)
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811 t n−1/2
Z
G2 (j, n)dt +
F3 = tn
t n−1/2
G2 (j + 1/2, n − 1/2)dt +
t n−1/2
Z
G2 (j, n)dt +
F4 = tn
xj+1/2
Z
tn
Z
t n−1/4
G2 (j, n)dt +
Z
xj
G1 (j, n)dx,
(9)
xj+1/4
G1 (j + 1/2, n − 1/2)dx
xj
tn
Z
G1 (j + 1/2, n − 1/2)dx
xj
t n−1/4
Z +
xj+1/4
Z
809
+ t n−1/2
G2 (j + 1/2, n − 1/2)dt +
Z
xj
G1 (j, n)dx.
(10)
xj+1/2
We define
Ujn
1 ∆x n 2! 2 (Ux )j 2 q(j, n) = 1 ∆x (Uxx )nj 3! 2 3 1 ∆x (Ut )nj 4! 2
,
(11)
and, in this way, the initial condition
u(xj , 0)
∆x ∂ u(x , 0) 4 ∂x j q(j, 0) = ∆x ∂ 2 . u ( x , 0 ) j 24 ∂ x2 ∆x ∂ u(xj , 0) 192 ∂ t From equations Fi (j, n) = 0, 1 ≤ i ≤ 4 , given by Eqs. (7)–(10), and expressing them the matrix form, we have q(j, n + 1) = (Q+ )2 q(j − 1, n) + [Q+ Q− + Q− Q+ ] q(j, n) + (Q− )2 q(j + 1, n),
(12)
(13)
where q(j, n) is given by (11) and Q+ and Q− are the 4 × 4 matrix Q+ =
1 2ν − 3 2 + 5ξ + 2ξ 2
2
q+1
q+2
q+3
q+4
q−1
q−2
q−3
q−4 ,
(14)
and Q− =
1 2ν 2 − 3 2 + 5ξ + 2ξ 2
(15)
where
ν=a
q+1
q+2
∆t , ∆x 1
ξ = 4µ
∆t , (∆x)2
(1 + ν) −4ν + 6ν 2 − 3 2 + 5ξ + 2ξ 2 2 3 − −1 + ν 2 ( 1 + 2 ξ ) , 2 = −2 ν −1 + ν 2 2 2 2 µ −1 + ν −1 + 2ν − 2ξ 2aξ 2 9ξ 2 11ξ −5 − 4ν + ν 2 (7 − 3ξ ) − + + 3ξ 3 2 2 3 15ξ 2 2 ν − ξ ( 1 + 2 ξ ) + ν − 4 − + 3 ξ 2 2 = 2 + 4ν + 2ν 2 (−3 + ξ ) + ξ µ2 −3 + 6ν 2 − 6ξ + ν 3 (−7 + 3ξ ) + ν 1 − 6ξ − 9ξ 2 − 3aξ 2
,
810
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
q+3
q+4
q−1
q−2
q−3
q−4
6ν 2 − 3(2 + ξ ) −
1
ν(4 + 15ξ )
2 3 2 − 1 + ν + 4ξ + 12ξ 2 − 4ν(1 + 3ξ ) 4 , = 15ξ 9ξ 2 2 − 3ξ + ν 2 + 3 − 5ν − 2 2 µ2 5 + 8ν 3 + 7ξ − 12ξ 2 − 12ξ 3 − ν 2 (5 + 11ξ ) + 4ν −1 + 3ξ + 3ξ 2 4aξ 2 2 2 3aξ 2 + 8ν − 6ν + 13ξ + 6ξ 2 − 2µ2 2 3a ξ (− 3 − 6 ξ + ν(2 + 3ξ )) − 2 µ = , 2 2 3aξ 2 − 4ν + 2ν + ξ − 2 µ 2 2 3 aξ 10ν − 6ν + 3ξ (1 + 2ξ ) + ν(−2 + 3ξ ) − 2aξ 2 1 − (−1 + ν) 4ν + 6ν 2 − 3 2 + 5ξ + 2ξ 2 2 3 2 − 1 + ν ( 1 + 2 ξ ) 2 = , 2 2ν −1 + ν µ2 −1 + ν 2 −1 + 2ν 2 − 2ξ − 2 2aξ 1 10 − 8ν + 11ξ − 9ξ 2 − 6ξ 3 + 2ν 2 (−7 + 3ξ ) 2 1 2 2 2ν − 3ξ (1 + 2ξ ) + ν 8 + 15ξ − 6ξ , 2 = 2 −2 + 4ν − 2ν (−3 + ξ ) − ξ µ2 3 − 6ν 2 + 6ξ + ν 3 (−7 + 3ξ ) + ν 1 − 6ξ − 9ξ 2 3aξ 2 15ξ 2 6ν − 3(2 + ξ ) + ν 2 + 2 3 2 2 − −1 + ν + 4ξ + 12ξ + 4ν(1 + 3ξ ) 4 = , 1 2 2 6 − 4ν − 10ν − 9ξ − 15νξ − 6ξ 2 µ2 −5 + 8ν 3 − 7ξ + 12ξ 2 + 12ξ 3 + ν 2 (5 + 11ξ ) + 4ν −1 + 3ξ + 3ξ 2 4aξ 2 2 2 3aξ 2 − 8ν − 6ν + 13ξ + 6ξ 2 2µ2 3aξ 2 (3 + 6ξ + ν(2 + 3ξ )) 2 µ = . 2 2 3aξ 2 + 4ν + 2ν + ξ 2 µ aξ 2 10ν 2 + 6ν 3 + ν(2 − 3ξ ) + 3ξ (1 + 2ξ ) − 2aξ 2
3. Numerical example In order to test the proposed scheme we consider an analytical solution. We choose the one-dimensional advectiondiffusion equation Eq. (1) of a Gaussian pulse of unit height, centered at x0 = 1 in a region bounded by 0 ≤ x ≤ 9 whose analytical solution [4] is 1 (x − x0 − at )2 ua (x, t ) = √ exp − µ(4t + 1) 4t + 1
.
(16)
R. Company et al. / Mathematical and Computer Modelling 50 (2009) 806–811
811
Table 1 Errors in t = 5, x = 5.
∆x = 0.04 ∆x = 0.02
CE-SE
Second order CE-SE
0.0027 0.0011
0.0014 0.00033
u 0.25
u 0.25
0.20
0.20
0.15
0.15
0.10
0.10
0.05
0.05
4.0
4.5
5.0
5.5
x 6.0
4.0
4.5
5.0
5.5
x 6.0
Fig. 4. Analytical and numerical solution for ∆x = 0.4, t = 5.
The initial condition is given by
(x − x0 )2 u(x, 0) = exp − µ
and the boundary condition at the two ends at any time t is obtained by substituting x = 0 and x = 9, respectively, in (16). The values of the parameters used are ν = 0.005 m2 /s and a = 0.8 m/s. In this section we compare the numerical method presented in this article with the original CE-SE method. We compare the errors of two methods (Table 1) with ν = 0.8 and ∆x = 0.4. We can observe that the error of the modified CE-SE method improves the one obtained with CE-SE method. In the left side of Fig. 4 is represented the analytical solution and numerical solution for the CE-SE scheme at time t = 5 with ∆x = 0.4. On the right side of Fig. 4, the analytical solution is represented versus the proposed CE-SE scheme at the same time considering ∆x = 0.4. We can see that the approximation solution obtained with the second order CE-SE is better than that obtained with the traditional CE-SE scheme. References [1] X.Y. Wang, C.Y. Chow, S.C. Chang, Application of the space-time conservation element and solution element method to one-dimensional advectiondifussion problems, NASA-TM-1999-209068 (1999). [2] S.C. Chang, The method of space-time conservation element and solution element. A new approach for solving the Navier–Stokes and Euler equations, J. Comput. Phys. 119 (1995) 295–324. [3] E. Ponsoda, E. Defez, M.D. Roselló, J.V. Romero, A stable numerical method for solving for solving variable coefficient advection-diffusion models, Comput. Math. Appl. 56 (2008) 754–768. [4] H. Karahan, A third-order upwind scheme for the advection-diffusion equation using spreadsheets, Adv. Eng. Softw. 38 (2007) 688–697.