Numerical Solution of the 1D Advection-Diffusion Equation Using Standard and Nonstandard Finite Difference Schemes

Three numerical methods have been used to solve the one-dimensional advection-diffusion equation with constant coefficients. This partial differential equation is dissipative but not dispersive. We consider the Lax-Wendroff scheme which is explicit, the Crank-Nicolson scheme which is implicit, and a nonstandard finite difference scheme (Mickens 1991). We solve a 1D numerical experiment with specified initial and boundary conditions, for which the exact solution is known using all these three schemes using some different values for the space and time step sizes denoted by h and k, respectively, for which the Reynolds number is 2 or 4. Some errors are computed, namely, the error rate with respect to the L 1 norm, dispersion, and dissipation errors. We have both dissipative and dispersive errors, and this indicates that the methods generate artificial dispersion, though the partial differential considered is not dispersive. It is seen that the Lax-Wendroff and NSFD are quite good methods to approximate the 1D advectiondiffusion equation at some values of k and h. Two optimisation techniques are then implemented to find the optimal values of k when h = 0.02 for the Lax-Wendroff and NSFD schemes, and this is validated by numerical experiments.


Introduction
The significant applications of advection-diffusion equation lie in fluid dynamics [1], heat transfer [2], and mass transfer [3]. The 3D advection-diffusion equation is given by The coefficient of diffusivity is denoted by and is computed as = / , where , , and denote the pressure, specific heat of the fluid at constant pressure, and thermal conductivity, respectively. Also , , and are the velocity components of the fluid in the directions of , , and , respectively. Equation (1) is also referred to as the convection-diffusion equation. The three terms ( / ), ( / ), and ( / ) are called the advective or convective terms and the terms ( 2 / 2 ), ( 2 / 2 ), and ( 2 / 2 ) are called the diffusive or viscous terms.
In this paper, we consider the one-dimensional convection-diffusion equation given by with = 1, = 0.01, 0 ≤ ≤ 1, and 0 < ≤ . We denote the spatial and temporal step sizes by ℎ and , respectively. The cfl number, , is computed as /ℎ, and the parameter, , is obtained as /ℎ 2 .
The initial condition is ( , 0) = ( ), and boundary conditions are (0, ) = 0 ( ) , 0 < ≤ , where , 0 , and 1 are known functions. There has been little progress in obtaining analytical solution to the 1D advection-diffusion equation when initial and boundary conditions are complicated, even with and being constant [4]. This is the reason why numerical solution of (2) is important.

Journal of Applied Mathematics
The paper is organised as follows. In Section 2, we study the damping and dispersive characteristics of some numerical methods for the 1D advection diffusion equation. In Section 3, we show how to quantify the errors from the numerical results into dissipation and dispersion errors by using a technique devised by Takacs [5]. In Section 4, we describe the numerical experiment that we have considered and show how to choose the parameters and ℎ to run the numerical experiments. Sections 5 and 6 describe some explicit and implicit methods, and we study their dissipative and dispersive properties. Also, we tabulate the errors when the methods are used to solve the numerical experiment described in Section 6. In Section 7, we present a nonstandard finite difference (NSFD) scheme, analyse its spectral properties, and also use it to solve the numerical experiment. In Section 8, we find the optimal value of when ℎ = 0.02 for the Lax-Wendroff and NSFD schemes and validate these using the numerical experiment. Section 9 highlights the salient features of the paper.

Dissipative and Dispersive Characteristics of Numerical Methods
Dissipation is defined as the constant decrease with time of the amplitude of plane waves, as they propagate in time. If the modulus of the amplification factor, denoted by AFM is equal to one, a disturbance neither grows nor damps [6]. The modulus of the amplification factor is also a measure of the stability of a scheme. If this value is greater than one, this creates instability, while damping is present whenever the value is less than one [7]. When the modulus of the amplification factor exceeds one, this indicates an unstable mode [8]. Since our partial differential equation is + = , we will have dissipation, this is caused because of the term , and such dissipation is called implicit dissipation. We can also have artificial dissipation which is caused due to the numerical method.
We let the amplification factor of the scheme approximating (2) be Then the modulus of the amplification factor, denoted by AFM, is computed as | |. We now show how the relative phase error (RPE) of a given numerical scheme approximating (2) is obtained. A perturbation for is obtained by substituting by exp( ( 1 − )) where and represent time and space, respectively, is the wavenumber, and 1 is the dispersion relation [9].
We then obtain where = √ −1, which on simplification gives Hence, the dispersion relation is given by The exact phase velocity is computed as R( 1 )/ which simplifies as . We next obtain the numerical phase velocity.
From (4), we have = 1 + 2 . We can express as exp(− ) [9] where is the exponential growth rate. Therefore, we have exp(− ) = 1 + 2 which implies The numerical phase velocity is computed as I( )/ and is equal to The phase angle, , is computed as = ℎ where is the wavenumber and ℎ is the spatial step. The relative phase error (RPE) is a measure of the dispersive character of a scheme. This quantity is a ratio and measures the velocity of the computed waves to that of the physical waves. Hence, we have Since = ℎ and = /ℎ, we can express (10) as If the RPE is greater than one, the computed waves appear to move faster than the physical waves [10] thus causing phase lead. A ratio less than one implies that the computed waves will move slower than the physical waves, causing phase lag.

Quantification of Errors from Numerical
Results [5,11,12] In this section, we describe how Takacs [5] quantifies errors from numerical results into dispersion and dissipation errors. The Total Mean Square Error is calculated as where represents the analytical solution and V represents the numerical (discrete) solution at a given grid point, .
The Total Mean Square Error can be expressed as Next, one has Journal of Applied Mathematics 3 The Total Mean Square Error can be further expressed as The expression in (15) can be rewritten as where 2 ( ) and 2 (V) denote the variance of and V, respectively, and V denote the mean values of and V, respectively.
Thus, the Total Mean Square Error is given by which on further simplification yields Thus, we have But the correlation coefficient, , is given by Cov( , V)/ ( ) (V). Hence, the Total Mean Square Error can be written as which simplifies to On putting = 1, we get 2(1 − ) ( ) (V) = 0. Thus, we define (2(1− ) ( ) (V)) as the dispersion error as correlation coefficient in statistics is analogous with phase lag or phase lead in Computational Fluid Dynamics. Consequently, ( ( ) − (V)) 2 + ( − V) 2 measures the dissipation error.
We also obtain values of the error rate with respect to the 1 norm which is calculated as where and V are the exact and computed values, respectively, and is the number of spatial grid points.

Choice of the Parameters ℎ and
We refer to [4] where three explicit methods are used to solve the partial differential equation where Tests were carried out for three values of the cell Reynolds number, Δ = / , namely, Δ = 2, 4, 8 [4]. Since = 0.8 /ℎ and = 0.008 /ℎ 2 , we can express Δ in terms of ℎ, in that case we have Δ = 100ℎ. Thus, for Δ = 2, 4, 8, the corresponding values of ℎ are 0.02, 0.04, and 0.08, respectively.
Then three values of were chosen as 0.16, 0.32, and 0.64, and then the corresponding values of determined as 0.004, 0.008, 0.016, 0.032, 0.064. For these values of , the number of time steps, , are calculated as = 1/ , and hence take the following values, namely, 250, 125, 62.5, 31.25, and 16.625, respectively. However, we note that can only be an integer. Hence, an improvement can be made when choosing and while keeping Δ = 2, 4, 8 and ℎ = 0.02, 0.04, 0.08.
We next refer to [13] where both explicit and implicit methods have been used for numerical solution of the one-dimensional advection-diffusion equation in a region bounded by 0 ≤ ≤ 1 and 0 ≤ ≤ 1 [3], with = 1, = 0.01 and with the following initial and boundary conditions: ) .

Journal of Applied Mathematics
The exact solution is given by ) .

(26)
The values of ℎ and used were 0.02 and 0.004, respectively, for all the numerical methods considered in [13].
In our work, we consider both implicit and explicit schemes to solve Hence, for ℎ = 0.02, the values taken by are 0.005, 0.01, and 0.02. For ℎ = 0.04, can take the values 0.01, 0.02, and 0.04. Some of these possibilities might give rise to an unstable scheme and must be ignored.
However, for implicit methods, all the 6 combinations of and ℎ are possible, and we can also consider the case when = 2.0 instead of only the three cases, namely, = 0.25, 0.5, and 1.0.

Construction of Explicit and Implicit Finite Difference Methods
We can approximate / as or Hence, an approximation for / is where ℎ represents the spatial step size, and are the temporal and spatial weighting factors, respectively. An approximation for 2 / 2 is Hence, a discretization for 2 / 2 is On plugging approximations for / and 2 / 2 as given by (30) and (33) into (2), we obtain a family of explicit and implicit numerical schemes given by where where = /ℎ and = /ℎ 2 .

Standard Schemes
The Lax-Wendroff scheme is given by and is obtained on replacing by zero and by (1 − )/2, in (34). The modified equation is given by [14] and this indicates that the leading error terms are dispersive in nature. The amplification factor and the relative phase error are obtained as ) .
The modified equation is given by  and this indicates that the leading error terms are dispersive in nature.
The amplification factor is given by and the RPE is computed as where 2 = 1 + − cos( ), 2 = 2 sin( ), and 2 = 4 + 4 cos( ) − 4 . The scheme is unconditionally stable. We next plot the variation of the AFM versus phase angle for some values of and ℎ in Figures 3(a) and 3(b). Plots of the RPE versus phase angle are depicted in Figures 4(a) and 4(b).
In the case of Crank-Nicolson, the scheme is less dissipative at ℎ = 0.04 as compared to ℎ = 0.02 for all the four values of , namely, 0.005, 0.01, 0.02, and 0.04. The combination ℎ = 0.04, = 0.005 is the least dissipative one. Based on Figure 4(b), we can observe that dispersion character is slightly affected by the value of used when ℎ = 0.04. However, if we choose ℎ = 0.02, the dispersion character is much affected by the value of . In general for ℎ = 0.02, the case = 0.02 is in general the least dispersive one.
We tabulate the errors for the eight combinations of ℎ and in Table 2

Nonstandard Finite Difference Scheme
In this section, we describe how a nonstandard finite difference scheme (NSFD) is constructed [15] for the 1D convection-diffusion equation.
The equation + = has three subequations [16] which are given by = . (45) Equations (43) and (44) have known exact finite difference scheme which are with = ℎ and respectively. A finite difference scheme that englobes the features of the two equations, namely, (43) and (44) is

Journal of Applied Mathematics
On rearranging the terms in (48), we get the NSFD method which is [15,16] The square of the modulus of the amplification factor is given by For stability, 0 < | | ≤ 1 and this implies that 0 < | | 2 ≤ 1. We now obtain the region of stability by using the approach used by Hindmarsh et al. [17] and Sousa [18].
We consider the case when = . The square of the modulus of the amplification factor is given by We thus need, which implies that Thus, for stability, we have the following inequality: Which was simplified to Since 1 and 1 are positive, 1 + 2 1 ≥ 0 is the trivial inequality. Hence, we consider the inequality Since, 1 = /ℎ and 1 = 1 /(exp(ℎ/ ) − 1), we have For stability, we need the following condition: However, the stability condition in (59) is difficult to achieve. We use a Maclaurin's series for exp(ℎ/ ), and therefore (59) reduces to ℎ + 2 ℎ 2 ≤ 1.
The scheme is least dissipative when = 0.01, ℎ = 0.04 and = 0.005, ℎ = 0.02. The scheme is least dispersive when = 0.02, ℎ = 0.04. The scheme experiences both phase lead and phase lag behaviour, depending on the values of and ℎ. The results of our numerical experiment are shown in Figures  7(a) and 7(b).
The modified equation is given by and this indicates that the leading error terms are dissipative. We tabulate the errors in Table 3, and we observe that the errors are the least when = 0.005 and ℎ = 0.02 and the greatest when = 0.02 and ℎ = 0.04. Based on Tables 1, 2, and 3, we can see that the Lax-Wendroff and the NSFD schemes are most effective when = 0.005 and ℎ = 0.02. The errors are smaller for the Lax-Wendroff as compared to NSFD scheme when = 0.005 and ℎ = 0.02.

Optimising Parameters in the Lax-Wendroff and NSFD Schemes
Our aim in this section is to compute an optimal value of for a given value of ℎ, say ℎ = 0.02. By optimal, we mean a value which reduces the errors. Since the partial differential equation considered is slightly dissipative and not dispersive, we aim to minimize the dispersion error of the scheme. [19], Bogey and Bailly [20] among others have implemented techniques which enable coefficients to be determined in numerical schemes specifically designed for Computational Aeroacoustics. We develop these techniques into respective equivalent forms [21] to determine the optimal values of for the Lax-Wendroff and NSFD schemes. We now describe briefly how Tam and Webb [19], Bogey and Bailly [20] define their measures and consequently their technique of optimisation in Computational Aeroacoustics.

Proposed Techniques of Optimisation. Tam and Webb
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 [22] set as 1.1 and optimise the coefficients in the numerical scheme, such that the integrated error is minimised.  Bogey and Bailly [20] minimise the relative difference between the exact wavenumber, ℎ, 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, Bogey and Bailly in a Computational Aeroacoustics framework to suit them in a Computational Fluid Dynamics framework [21], such that the optimal parameter can be obtained. Thus, we define the following integrals: integrated Error from Tam and Webb, (IETAM), integrated error from Bogey and Bailly (IEBOGEY) as follows:

Optimisation Procedure
Lax-Wendroff. We consider the Lax-Wendroff scheme given by (36), with ℎ = 0.02. The amplification factor of the resulting method is LW = 1 − 50 − 2500 2 + (50 + 2500 2 ) cos ( ) and therefore the RPE is computed as A plot of the exact RPE versus ∈ [0, 1.1] is shown in Figure 8, and we do not have phase wrapping phenomenon. We propose two measures, one adapted from Tam and Webb [19] and the other from Bogey and Bailly [20].
We compute the following: We plot the integrated errors versus in Figures 9(a) and 9(b) and obtain the optimal value of . We can also use the function NLPSolve from Maple to determine the value of We next validate whether this value of computed does indeed minimise the errors by performing the numerical experiment using Lax-Wendroff with ℎ = 0.02 at some different values of ∈ (0, 0.01236) and then compare the errors. The errors are tabulated in Table 4, and we can see that and therefore the RPE is computed as where 1 = /ℎ. A plot of the exact RPE versus ∈ [0, 1.1] is shown in Figure 10, and we do not have phase wrapping phenomenon. We propose two measures, one adapted from Tam and Webb [19] and the other from Bogey and Bailly [20].
We compute the following: We plot the integrated errors versus in Figures 11(a) and 11(b) and obtain the optimal value of . We can also use the function NLPSolve from Maple to determine the value of which minimise each of these two integrals. In the case of IETAM, we obtain We next validate whether this value of computed does indeed minimise the errors by performing the numerical experiment using NSFD with ℎ = 0.02 at some different values of ∈ (0, 0.01] and then compare the errors. The errors are tabulated in Table 5, and we can see that indeed for = 1/164 ≈ 0.00601, all the five types of errors are least.

Conclusion
In this paper, three numerical methods have been used to solve a 1D advection-diffusion equation with specified initial and boundary conditions. Both explicit and implicit finite difference methods as well as a nonstandard finite difference scheme have been used. When the 1D linear advection equation is approximated by a numerical method, the amplification factor and relative phase error depend on only the cfl number. However, in the case of the 1D advectiondiffusion equation the modulus of the amplification factor and relative phase error depends on the spatial and temporal step sizes. The results of our numerical experiment are much affected by the choice of and ℎ. In general, we observe that the Lax-Wendroff scheme is the most efficient method followed by the nonstandard finite difference scheme. We perform two optimisation procedures by computing the optimal values of when ℎ = 0.02 for the Lax-Wendroff and NSFD schemes. We observe that when ≈ 0.006, the errors are reduced further for both methods.
This work can be extended to the case when is large. Also, we can consider numerical solution of 1D nonlinear as well as 2D linear and 2D nonlinear convectiondiffusion problems, and we can use appropriate optimisation techniques to choose parameters ℎ and for minimal numerical dispersion and numerical dissipation.