Abstract
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 and , respectively, for which the Reynolds number is 2 or 4. Some errors are computed, namely, the error rate with respect to the 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 advection-diffusion equation at some values of and . Two optimisation techniques are then implemented to find the optimal values of when for the Lax-Wendroff and NSFD schemes, and this is validated by numerical experiments.
1. 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 , , and are called the diffusive or viscous terms.
In this paper, we consider the one-dimensional convection-diffusion equation given by with , , , and .
We denote the spatial and temporal step sizes by and , respectively. The cfl number, , is computed as , and the parameter, , is obtained as .
The initial condition is , and boundary conditions are where , , and 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.
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 for the Lax-Wendroff and NSFD schemes and validate these using the numerical experiment. Section 9 highlights the salient features of the paper.
2. 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 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 , is computed as . We now show how the relative phase error of a given numerical scheme approximating (2) is obtained.
A perturbation for is obtained by substituting by where and represent time and space, respectively, is the wavenumber, and is the dispersion relation [9].
We then obtain where , which on simplification gives
Hence, the dispersion relation is given by
The exact phase velocity is computed as which simplifies as . We next obtain the numerical phase velocity.
From (4), we have . We can express as [9] where is the exponential growth rate.
Therefore, we have which implies The numerical phase velocity is computed as and is equal to The phase angle, , is computed as where is the wavenumber and is the spatial step. The relative phase error 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 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.
3. 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 represents the numerical (discrete) solution at a given grid point, .
The Total Mean Square Error can be expressed as Next, one has The Total Mean Square Error can be further expressed as The expression in (15) can be rewritten as where and denote the variance of and , respectively, and denote the mean values of and , 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 . Hence, the Total Mean Square Error can be written as which simplifies to On putting , we get . Thus, we define as the dispersion error as correlation coefficient in statistics is analogous with phase lag or phase lead in Computational Fluid Dynamics.
Consequently, measures the dissipation error.
We also obtain values of the error rate with respect to the norm which is calculated as where and are the exact and computed values, respectively, and is the number of spatial grid points.
4. 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, [4].
Since and , we can express in terms of , in that case we have . Thus, for , the corresponding values of are 0.02, 0.04, and 0.08, respectively.
Since and , we have the following relationships between and , namely, , , and .
Then three values of were chosen as , , and , 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 , 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 and .
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 and [3], with , and with the following initial and boundary conditions:
The exact solution is given by
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 subject to boundary conditions given by (25).
We consider two values for , say 2 and 4. Thus, we have as and . For and , we have and , respectively. Hence, and , and therefore we have and .
We consider the case when . If we choose , 0.50, and 1.0, then the values taken by are 0.01, 0.02, and 0.04, respectively.
Next we consider . If we choose , 0.5, and 1.0, we have , 0.02, and 0.04, respectively.
Hence, for , the values taken by are 0.005, 0.01, and 0.02. For , 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 instead of only the three cases, namely, , 0.5, and 1.0.
5. 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 is or
Hence, a discretization for is
On plugging approximations for and as given by (30) and (33) into (2), we obtain a family of explicit and implicit numerical schemes given by where where and .
6. Standard Schemes
6.1. Lax-Wendroff Scheme
The Lax-Wendroff scheme is given by and is obtained on replacing by zero and by , 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 region of stability is [14]. Plots of the and , both versus the phase angle, , for four combinations of values of and are shown in Figures 1(a) and 1(b). The combination , is the least dissipative one. The scheme is not dispersive when , . Phase lag behaviour is observed when , and , . Phase lead phenomenon occurs when and .
(a) AFM
(b) RPE
The results of our numerical experiment for the four combination of values of and are shown in Figures 2(a) and 2(b). We tabulate the errors in Table 1. The errors are the least when and and greatest when and .
(a)
(b)
6.2. Crank-Nicolson Scheme
The Crank-Nicolson method is obtained if we plug and into (34). A single expression for the scheme is
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 , , and .
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).
| (a) |
| (b) |
| (a) |
| (b) |
In the case of Crank-Nicolson, the scheme is less dissipative at as compared to for all the four values of , namely, 0.005, 0.01, 0.02, and 0.04. The combination , 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 . However, if we choose , the dispersion character is much affected by the value of . In general for , the case is in general the least dispersive one.
We tabulate the errors for the eight combinations of and in Table 2, and we observe that the errors are the least when and . The results of the numerical experiment are shown in Figures 5(a) and 5(b).
(a)
(b)
7. 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
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
On rearranging the terms in (48), we get the NSFD method which is [15, 16] where
The square of the modulus of the amplification factor is given by
For stability, and this implies that . 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 and are positive, is the trivial inequality. Hence, we consider the inequality
Since, and , 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 , and therefore (59) reduces to
We next consider the case when . When , and .
Thus, (51) reduces to
We thus require
Using (50), (62) becomes From (63), we deduce , which on expansion and simplification gives .
Since , therefore,
Combining (60) and (64), we obtain (60), and therefore the region of stability is described by
Case 1. For and , using (60), we have . Hence, for and , we have . Also, for and , we have . However, if and , we have , but this combination will give rise to an unstable method.
Case 2. For and , using (60), the region of stability is given by . Therefore, for , we consider and .
Plots of the and versus phase angle are shown in Figures 6(a) and 6(b), respectively.
(a) AFM
(b) RPE
The NSFD scheme considered is an explicit one, and we have four combinations of and , namely,(i), ;(ii), ;(iii), ;(iv), .
The scheme is least dissipative when , and , . The scheme is least dispersive when , . 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).
(a)
(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 and and the greatest when and .
Based on Tables 1, 2, and 3, we can see that the Lax-Wendroff and the NSFD schemes are most effective when and . The errors are smaller for the Lax-Wendroff as compared to NSFD scheme when and .
8. 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 . 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.
8.1. Proposed Techniques of Optimisation
Tam and Webb [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.
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 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 or
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:
8.2. Optimisation Procedure
Lax-Wendroff. We consider the Lax-Wendroff scheme given by (36), with . The amplification factor of the resulting method is and therefore the is computed as
A plot of the exact versus 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 which minimise each of these two integrals. In the case of IETAM, we obtain while in the case of IEBOGEY, we are out with We next validate whether this value of computed does indeed minimise the errors by performing the numerical experiment using Lax-Wendroff with at some different values of and then compare the errors. The errors are tabulated in Table 4, and we can see that indeed for , all the five types of errors are least.
(a) IETAM
(b) IEBOGEY
NSFD. We consider the NSFD scheme given by (49), with . The amplification factor of the resulting method is and therefore the is computed as where .
A plot of the exact versus 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 while in the case of IEBOGEY, we are out with We next validate whether this value of computed does indeed minimise the errors by performing the numerical experiment using NSFD with at some different values of and then compare the errors. The errors are tabulated in Table 5, and we can see that indeed for , all the five types of errors are least.
(a) IETAM
(b) IEBOGEY
9. 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 advection-diffusion 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 for the Lax-Wendroff and NSFD schemes. We observe that when , 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 convection-diffusion problems, and we can use appropriate optimisation techniques to choose parameters and for minimal numerical dispersion and numerical dissipation.
Nomenclature
| Time step | |
| Spatial step | |
| Time level | |
| Advection velocity | |
| cfl/Courant number | |
| Phase angle in 1D | |
| Relative phase error per unit time step | |
| Amplification factor | |
| Exact results | |
| Computed results | |
| Diss. Error: | Dissipation error |
| Disp. Error: | Dispersion error. |
Acknowledgments
This work was funded through the Research Development Programme of the University of Pretoria, and the period of funding is from January 2012 to January 2013. The author would like to thank Professor Jean Lubuma for some discussion on nonstandard finite difference schemes. Lastly, the author is grateful to the two anonymous reviewers for their comments which were useful in clarifying and focusing the presentation.