Time-Splitting Procedures for the Numerical Solution of the 2 D Advection-Diffusion Equation

We perform a spectral analysis of the dispersive and dissipative properties of two time-splitting procedures, namely, locally onedimensional (LOD) Lax-Wendroff and LOD (1, 5) [9] for the numerical solution of the 2D advection-diffusion equation. We solve a 2D numerical experiment described by an advection-diffusion partial differential equation with specified initial and boundary conditions for which the exact solution is known. Some errors are computed, namely, the error rate with respect to the L 1 norm, dispersion and dissipation errors. Lastly, an optimization technique is implemented to find the optimal value of temporal step size that minimizes the dispersion error for both schemes when the spatial step is chosen as 0.025, and this is validated by numerical experiments.


Introduction
The advection-diffusion equation is a parabolic partial differential equation combining the diffusion and advection (convection) equations, which describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: diffusion and advection [1].The numerical solution of advectiondiffusion equation plays an important role in many fields of science and engineering.These include the transport of air and groundwater pollutants, oil reservoir flow [2], heat transfer in draining film, flow through porous media, the dispersion of pollutants in rivers and streams, long range transport of pollutants in the atmosphere, thermal pollution in river systems, and dispersion of dissolved salts in groundwater [3].
The 3D advection-diffusion equation is given by where   ,   , and   are the velocity components of advection in the directions of , , and , respectively, and   ,   , and   are the coefficients of diffusivity in the -, -, and directions, respectively.This study deals with the 2D advection-diffusion equation, where 0 <  ≤ , in the domain 0 ≤  ≤ 1, 0 ≤  ≤ 1, with specified initial and boundary conditions.Dehghan [3] proposed two time-splitting procedures for the solution of the two-dimensional transport equation.The time-splitting procedure used is the locally-one dimensional (LOD) in proceeding from one time step to the next step.LOD replaces the complicated multidimensional partial differential equations by a sequence of solutions of simpler onedimensional partial differential equations.The originality in this work is that we perform a spectral analysis of the dispersion and dissipation properties of the two schemes at some values of the temporal and spatial step sizes.3D plots of the relative phase error per unit time step (RPE) and the modulus of the amplification factor (AFM) versus phase angles in and -directions are obtained.We then use optimization strategies to compute the optimal values of

Numerical Dispersion and Dissipation
Dissipation reduces the amplitude of sinusoids in a Fourier series.This is caused by the presence of derivatives like   and −  in the modified equation [4].On the other hand, the amplitude of sinusoids in a Fourier series is increased by antidissipation.Derivatives like −  and   in the modified equation are generally antidissipative.Dispersion affects the speed of sinusoids in a Fourier series causing phase lag or phase lead and is caused due to the presence of oddorder derivatives in the modified equation [4].
The modulus of the amplification factor (AFM) is a measure of the stability of a scheme and it is also used to measure the dissipative characteristics of the scheme.If the modulus of the amplification factor is equal to one, a disturbance neither grows nor damps [5].If the modulus of the amplification factor is greater than one, then the scheme is unstable [6]; if it is less than one, damping occurs [7].The partial differential equation given by ( 2) is dissipative in nature due to the terms   and   .
The relative phase error (RPE) is a measure of the dispersive characteristics of a scheme.The relative phase error of a scheme approximating the 1D advection-diffusion equation is given by where  is the Courant number,  is phase angle,  is the amplification factor of the numerical scheme approximating the 1D advection-diffusion equation, and R() and I() are the real and imaginary parts of , respectively [8].
We extend the work on the relative phase error in [8] for the case of the 2D advection-diffusion equation.The relative phase error for a numerical scheme approximating (2) is obtained on substituting  by exp(( 1  −  1  −  2 )) [9], where  = √ −1,  is the time,  1 and  2 are wavenumbers in the directions of  and , respectively, and  1 is the dispersion relation.Thus, we get On substituting (4) into (2), we obtain The exact phase velocity is R( 1 )/wavenumber and we obtain The amplification factor can be written as  =  1 +  2 , where  1 and  2 are the real and imaginary parts of , respectively.Also, we can express  as  = exp(−), where  is time step and  is exponential growth rate [9].Thus, we obtain ) .
The relative phase error is the ratio of the numerical phase velocity to the exact phase velocity [10].It is calculated as The phase angles in the directions of  and  are given by   = Δ 1 and   = Δ 2 , respectively.Hence, the relative phase error is where   =   /Δ and   =   /Δ are the Courant numbers in the directions of  and , respectively.

Time-splitting Procedures and Numerical Experiments
The For simplification we take Δ = Δ = ℎ.Let Δ = ; then the parameters ℎ and  represent the spatial and temporal grid spacing, respectively.We denote the approximated value of  at the grid point (, , ) by   , .

Time-Splitting
Procedures.Since one-dimensional schemes are easier to use than two-dimensional schemes, ( 2) is split into the following two one-dimensional equations: Each of ( 12) and ( 13) can be solved over half of a time step to be used for the complete 2D advection-diffusion equation, using the procedures developed for the 1D advectiondiffusion equation.Some work on time-splitting procedures can be found in [3,11].In this paper, we refer to [3] on how a 2D advectiondiffusion equation is converted into two 1D advectiondiffusion equations using the locally one-dimensional (LOD) time-splitting procedure.Solving ( 12) and ( 13) in each half time step is equivalent to solving the following equations over a full-time step: Then we can use the schemes used for solving the 1D advection-diffusion equation to solve (14) and (15).

Quantification of Errors from Numerical Results
. In this subsection, we describe how errors from numerical results can be quantified into dispersion and dissipation by a technique devised by Takacs [13].
The total mean square error in the 1D case [13] is calculated as where   is the exact solution, V  is the numerical solution at a grid point , and  is the number of grid points.The total mean square error can be expressed as The right hand side of (25) can be rewritten as where  2 () and  2 (V) denote the variance of  and V, respectively, and  and V denote the mean values of  and V, respectively.Then we have Therefore, where Cov(, V) = (1/)(∑  =1   V  −  V).The total mean square error can be expressed as where  = Cov(, V)/(()(V)) is the coefficient of correlation.
We extend the work on quantification of errors in [14,15] for the 2D case.The total mean square error for the 2D case is calculated as where  , and V , are the exact and numerical solutions at a grid point (, ), respectively, and we have Hence, where The error rate with respect to  1 norm for ,  ∈ [0, 1] is calculated as

Construction of the LOD Lax-Wendroff Procedure
We use the following approximations in the first time step of the LOD procedure [3]: where   is the spatial weighting factor in the direction of .
The following approximations are used for the second step of the LOD procedure: where   is the spatial weighting factor in the direction of .By using the following relationships the finite difference formula for the first half time step of the LOD Lax-Wendroff procedure is given by At the second half time step the finite difference is given by where We obtain a single expression for  +1 , in a complete time step as follows: To find the modified equation of the scheme, we first find Taylor's expansion of each term in (42) about   , .The Taylor series expansion of  +1 , is given by The Taylor series expansions for some grid points about   , are given as follows: Using (2), we convert the temporal derivatives   ,   , and   into spatial derivatives.We thus have Then on substituting the Taylor series expansions of each term of the difference scheme we get the following modified equation for the LOD Lax-Wendroff The scheme is second order accurate in space and the leading error terms are dispersive in nature (presence of odd-order derivatives   and   ).As the time and spatial increments go to zero, the modified equation (46) reduces to its original equation, that is, (2).Hence LOD Lax-Wendroff is consistent.We now study the spectral analysis of the dispersive and dissipative properties of the scheme for the case   =   =  and   =   = .To obtain the amplification factor we use the Von Neumann stability analysis by substituting   , by   exp((  +   )) in (42), where  = √ −1.We thus have The real and imaginary parts of the amplification factor are given by respectively.The modulus of the amplification factor, AFM, is obtained as We find the region of stability using the approach of Hindmarsh et al. [16] and Sousa [17].We consider the case when   =   =  and   → 0 and   → 0. gives From the Von Neumann stability analysis, the scheme is stable if and only if || ≤ 1.Thus, we get On simplification, we get which reduces to When   → 0 and   → 0, we use the following approximations: We consider (49) and use the approximations in (54) to obtain Thus LOD Lax-Wendroff is stable if −2 ≤ 0. Therefore, we have On combining ( 53) and (56), we obtain the region of stability for the LOD Lax-Wendroff procedure as We choose  = 0.01 and  = 0.8 [3].For ℎ = 0.025 from (57), we have On solving for , we get Therefore, for ℎ = 0.025, the stability region for the LOD Lax-Wendroff procedure is 0 ≤  ≤ 0.0193.We choose  such that 0.3/ is an integer as for our numerical experiments,  = 0.3.Then 0.3/ gives the number of time intervals.We choose  = 0.005, 0.01, 3/160, and 1/60.
For ℎ = 0.05, we have  = 16 and  = 4.When   =  and   = , using (51) we get which gives Therefore, for ℎ = 0.05, LOD Lax-Wendroff is stable if 0 ≤  ≤ 0.0487985.We then choose some values of  ∈ (0, 0.048) for the numerical experiments.for our numerical experiments.3D plots of the modulus of the amplification factor versus phase angles in and -directions for some different values of ℎ and  are depicted in Figures 1 and 2. 2D plots of the modulus of the amplification factor versus   , when   = 0, are shown in Figure 3.The scheme is in general less dissipative at ℎ = 0.05 as compared to ℎ = 0.025.Out of the eight combinations of values of ℎ and , the scheme is Mathematical Problems in Engineering least dissipative when ℎ = 0.05 and  = 0.005 and it is most dissipative when ℎ = 0.025 and  = 0.01.Figures 4 and 5 show the 3D plots of the relative phase error versus   versus   for some different values of ℎ and . Figure 6 shows the 2D plots of the relative phase error versus   , for the case   = 0 at ℎ = 0.025 and ℎ = 0.05.For ℎ = 0.025, we observe phase lag behaviour at  = 0.01 and  = 0.005 and phase lead behaviour at  = 1/60 and  = 3/160.For ℎ = 0.05, we have phase lag behaviour at  = 0.005, 0.01, 0.02, and 0.03.The scheme is least dispersive when ℎ = 0.05;  = 0.03 and ℎ = 0.025;  = 0.01.In the following section, we consider the LOD (1, 5) scheme.

LOD (1, 5) Explicit Procedure
In this procedure the following approximations are used in the first half time step [3] ) .
(62)  On substituting (62) into ( 14), we get the following finite difference equation for the first half time step: where The following approximations are used in the second half time step:
We obtain the following modified equation: The scheme is essentially dispersive as the leading error terms are dispersive in nature due to the presence of the oddorder derivative terms   and   .Also, the scheme is consistent and it is fourth order accurate in space.
The amplification factor of the LOD (1, 5) scheme is given by We consider   =   =  and   =   = .We use the Von Neumann stability analysis and the approach of Hindmarsh et al. [16] to obtain the stability region.When   =   = , on simplification of (71), we get For stability, we must have We consider (71) and for   → 0 and   → 0 we use the following approximations: We thus have On neglecting higher order terms, we have Thus, the numerical method is stable if −2 ≤ 0. Therefore, we have When ℎ = 0.025, we have  = 32 and  = 16.Using ( 73) and (77), we get respectively.On combining (78) and (79), we get the stability region when ℎ = 0.025 for the LOD (1, 5) procedure as When ℎ = 0.05, we have  = 16 and  = 4 and the modulus of the amplification factor is given by From ( 77) and (81), we obtain the stability region for the LOD (1, 5) procedure when ℎ = 0.05 as We analyse the spectral analysis; for ℎ = 0.025, we choose  = 0.005, 0.01, 0.02, 0.025 and for ℎ = 0.05, we choose  = 0.01, 0.02, 0.03 and 0.05.3D plots of the modulus of the amplification factor versus phase angle in -direction versus phase angle in -direction at two values of ℎ, namely, 0.025 and 0.05, at some different values of  are shown in Figures 7 and 8. 2D plots of the modulus of the amplification factor versus   when   = 0 are illustrated in Figure 9.The scheme is in general less dissipative at ℎ = 0.05 as compared to ℎ = 0.025.Out of the eight combinations of ℎ and  values, the scheme is least dissipative when ℎ = 0.05 and  = 0.01.
Figures 10 and 11 show the 3D plots of the relative phase error versus   versus   for some different values of ℎ and . Figure 12 shows the 2D plots of the relative phase error versus   , for the case   = 0 at ℎ = 0.025 and ℎ = 0.05, respectively.

Numerical Results
6.1.LOD Lax-Wendroff Procedure.The results of the numerical experiment at some values of  using the LOD Lax-Wendroff procedure at  = 0.3 at ℎ = 0.025 and ℎ = 0.05 are shown in Tables 1 and 2, respectively.We observe that the dispersion error is significantly greater than the dissipation error.Out of the five combinations of values of  and ℎ, we observe that the dispersion error and total mean square error are both least when ℎ = 0.05 and  = 0.03 and are both greatest when ℎ = 0.05 and  = 0.0025.(1,5) Procedure.The results of the numerical experiment using the LOD (1, 5) procedure at  = 0.3 for ℎ = 0.025 and ℎ = 0.05 are shown in Tables 3 and 4, respectively.Out of the 12 combinations of values of  and ℎ, we observe that the dissipation, dispersion error, and error rate are all least when ℎ = 0.05 and  = 0.05.Also, the dissipation and dispersion errors and error rate are greatest when ℎ = 0.05 and  = 0.0625.We observe that at ℎ = 0.025 the total mean square error, dispersion error, and dissipation error are not much affected by the values of .Also, the dispersion error is greater than the at all values of ℎ and  considered.

Optimization
In this section, we obtain the optimal value of  at ℎ = 0.025 that minimizes the dispersion error for the two time-splitting procedures.
Since the partial differential equation we consider is slightly dissipative and also we observe from the numerical experiments carried out that the dissipative errors are much less than the dispersive errors, we choose to minimize the square of the dispersion error of the two splitting schemes.
7.1.Proposed Techniques of Optimization.Tam and Webb [18], Bogey and Bailly [19], and Hixon [20] among others have implemented techniques which enable coefficients to be determined in numerical schemes specifically designed for computational aeroacoustics.We now describe briefly how Tam and Webb [18] and Bogey and Bailly [19] define their Mathematical Problems in Engineering      measures and consequently their technique of optimization in computational aeroacoustics.The dispersion-relation-preserving (DRP) scheme was designed so that the dispersion relation of the finite difference scheme is formally the same as that of the original partial differential equations.The integrated error is defined as where the quantities  * ℎ and ℎ represent the numerical and exact wavenumbers, respectively.The dispersion error and dissipation error are calculated as |R( * ℎ) − ℎ| and |( * ℎ)|, respectively.Tam and Shen [21] set  as 1.1 and optimize the coefficients in the numerical scheme such that the integrated error is minimized.
Bogey and Bailly minimize the relative difference between the exact wavenumber ℎ and the effective/numerical wavenumber  * ℎ and define their integrated errors as  = ∫ In computational fluid dynamics for a particular method under consideration, the dispersion error is calculated as We have modified the measures used by Tam and Webb and Bogey and Bailly in a computational aeroacoustics framework to suit them in a computational fluid dynamics framework [22] such that the optimal parameter can be obtained.We have defined the following integrated errors integrated error from Tam and Webb (IETAM) integrated error from Bogey and Bailly (IEBOGEY) [22] as follows: In [8], the integrated error for a scheme discretising the 1D advection-diffusion equation is obtained as The value of ℎ was fixed as 0.02 and the range of values of , was determined.Then the integrated error, which is a function of , was minimized and the optimal value of  is determined using NLPSolve function in maple.
We extend the work on optimization of parameters in [8] which is for the 1D advection-diffusion equation for the case of the 2D advection-diffusion equation.We define the integrated error as We first obtain an expression for the RPE of the LOD Lax-Wendroff when discretizing the equation Since  = 0.8/ℎ and  = 0.8/ℎ 2 and we choose ℎ = 0.025, we have  = 32 and  = 16.Hence, the RPE is a function of ,   , and   .Since we can have phase wrapping, we make use of Taylor's expansion to obtain an approximation for the RPE up to the terms (  ) 3 (  ) 2 , (  ) 5 , (  ) 2 (  ) 3 , (  ) 5 ,   (  ) 4 , (  ) 4   . (92) The integrated error ∫ 1.1 0 ∫ 1.1 0 (RPE − 1) 2     is obtained by using Simpson's method and it is a function of  only.A plot of the integrated error versus  is shown in Figure 13(a) for  ∈ [0, 0.0193].Using the NLPSolve function in maple, this optimal value of  is found to be 0.009593 and the minimum value of the integrated error is 1.883960 × 10 −7 .
To validate our results, we perform the same numerical experiment described in Section 3.2 at ℎ = 0.025 and use the optimal value of  or a value of  close to this optimal value, in that case,  = 3/310 ≈ 0.0096, and compute the errors.The error rate, total mean square error, and dispersion error are 3.9967 × 10 −4 , 4.1603 × 10 −6 , 3.9878 × 10 −6 .These three errors are all least as compared to when other values of  are used as shown in Table 1.
We adopt the same procedure to compute the optimal value of  for the LOD (1, 5) scheme when ℎ = 0.025.We obtain an approximate expression for the RPE of the LOD (1, 5) when  = 32 and  = 16.We use Simpson's rule to approximate the integral given by (90) which is a function of .A plot of the integrated error versus  for  ∈ [0, 0.026288] is shown in Figure 13(b) and using NLPSolve function in maple the optimal value of  is 0.013782 and also the minimum value of the integral is 1.139313 × 10 −6 .

Conclusion
In this paper, two time-splitting procedures are used to solve a 2D advection-diffusion equation with constant coefficients when the advection velocity in both and -directions is 0.8 and also when the coefficient of diffusivity in both and -directions is 0.01.We perform a stability analysis and spectral analysis of the dispersion and dissipation properties of the two schemes at some values of ℎ and .Numerical experiments are carried out and various errors are computed.These errors are dependent on the values of ℎ and .It is observed that in general the dispersion error is more affected by the values of  and ℎ for the LOD Lax-Wendroff scheme as compared to that of the LOD (1, 5) scheme at a given value of ℎ.We then use an optimization technique based on minimisation of the square of the dispersion error to find the optimal value of  when ℎ is chosen as 0.025 and this is validated by numerical experiments.
Future extension of this work to consider other types of advection-diffusion equations when dissipation dominates and to find out which optimization techniques are suitable in these cases.Also, the work can be extended to 2D nonlinear convection-diffusion problems.

Figure 1 :
Figure 1: Plots of modulus of amplification factor versus phase angle in -direction,   , versus phase angle in -direction,   , at ℎ = 0.025 with some different values of  for the LOD Lax-Wendroff scheme.

Figure 2 :
Figure 2: Plots of modulus of amplification factor versus phase angle in -direction,   , versus phase angle in -direction,   , at ℎ = 0.05 with some different values of  for the LOD Lax-Wendroff scheme.

Figure 3 :Figure 4 :
Figure 3: Plots of modulus of amplification factor versus   when   = 0 at ℎ = 0.025 and ℎ = 0.05 with some different values of  for the LOD Lax-Wendroff scheme.

)Figure 5 :Figure 6 :
Figure 5: Plots of relative phase error versus   versus   , at ℎ = 0.05 with some different values of  for the LOD Lax-Wendroff scheme.

Figure 10 :
Figure 10: 3D plots of relative phase error versus   versus   at ℎ = 0.025 at some different values of  for the LOD (1, 5) procedure.

Figure 13 :
Figure 13: Plots of integrated errors versus  at ℎ = 0.025 for the two time-splitting schemes.