- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Annual Issues
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Journal of Applied Mathematics
Volume 2013 (2013), Article ID 734374, 14 pages
Numerical Solution of the 1D Advection-Diffusion Equation Using Standard and Nonstandard Finite Difference Schemes
Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria 0002, South Africa
Received 9 October 2012; Accepted 9 January 2013
Academic Editor: Oluwole Daniel Makinde
Copyright © 2013 A. R. Appadu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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.
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 . 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 . 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 . 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 . When the modulus of the amplification factor exceeds one, this indicates an unstable mode .
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 .
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.
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  thus causing phase lead. A ratio less than one implies that the computed waves will move slower than the physical waves, causing phase lag.
In this section, we describe how Takacs  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
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  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 , 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 .
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
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  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 . 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 .
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 .
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).
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).
7. Nonstandard Finite Difference Scheme
In this section, we describe how a nonstandard finite difference scheme (NSFD) is constructed  for the 1D convection-diffusion equation.
The equation has three subequations  which are given by
The square of the modulus of the amplification factor is given by
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:
We next consider the case when . When , and .
Thus, (51) reduces to
We thus require
Since , therefore,
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 .
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).
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 .
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 , Bogey and Bailly  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  to determine the optimal values of for the Lax-Wendroff and NSFD schemes.
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  set as 1.1 and optimise the coefficients in the numerical scheme, such that the integrated error is minimised.
Bogey and Bailly  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 , 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
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.
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 .
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.
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.
|Phase angle in 1D|
|Relative phase error per unit time step|
|Diss. Error:||Dissipation error|
|Disp. Error:||Dispersion error.|
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.
- N. Kumar, “Unsteady flow against dispersion in finite porous media,” Journal of Hydrology, vol. 63, no. 3-4, pp. 345–358, 1983.
- J. Isenberg and C. Gutfinger, “Heat transfer to a draining film,” International Journal of Heat and Mass Transfer, vol. 16, no. 2, pp. 505–512, 1973.
- V. Guvanasen and R. E. Volker, “Numerical solutions for solute transport in unconfined aquifers,” International Journal for Numerical Methods in Fluids, vol. 3, no. 2, pp. 103–123, 1983.
- M. Dehghan, “On the numerical solution of the one-dimensional convection-diffusion equation,” Mathematical Problems in Engineering, vol. 2005, no. 1, pp. 61–74, 2005.
- L. L. Takacs, “A two-step scheme for the advection equation with minimized dissipation and dispersion errors,” Monthly Weather Review, vol. 113, no. 6, pp. 1050–1065, 1985.
- J. E. Fromm, “A method for reducing dispersion in convective difference schemes,” Journal of Computational Physics, vol. 3, no. 2, pp. 176–189, 1968.
- K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, UK, 1994.
- R. J. Babarsky and R. Sharpley, “Expanded stability through higher temporal accuracy for time-centered advection schemes,” Monthly Weather Review, vol. 125, no. 6, pp. 1277–1295, 1997.
- R. Smith and Y. Tang, “Optimal and near-optimal advection-diffusion finite-difference schemes. V. Error propagation,” Proceedings of the The Royal Society of London A, vol. 457, no. 2008, pp. 803–816, 2001.
- C. Hirsch, Numerical Computation of Internal and External Flows, vol. 1, John Wiley & Sons, New York, NY, USA, 1988.
- A. R. Appadu and M. Z. Dauhoo, “The concept of minimized integrated exponential error for low dispersion and low dissipation schemes,” International Journal for Numerical Methods in Fluids, vol. 65, no. 5, pp. 578–601, 2011.
- A. R. Appadu, “Some applications of the concept of minimized integrated exponential error for low dispersion and low dissipation,” International Journal for Numerical Methods in Fluids, vol. 68, no. 2, pp. 244–268, 2012.
- A. Mohammadi, M. Manteghian, and A. Mohammadi, “Numerical solution of the one- dimensional advection-diffusion equation using simultaneously temporal and spatial weighted parameters,” Australian Journal of Basic and Applied Sciences, vol. 5, no. 6, pp. 1536–1543, 2011.
- M. Dehghan, “Weighted finite difference techniques for the one-dimensional advection-diffusion equation,” Applied Mathematics and Computation, vol. 147, no. 2, pp. 307–319, 2004.
- R. E. Mickens, Applications of Nonstandard Finite Difference Schemes, World Scientific Publishing, River Edge, NJ, USA, 2000.
- R. E. Mickens, “Analysis of a new finite-difference scheme for the linear advection-diffusion equation,” Journal of Sound and Vibration, vol. 146, no. 2, pp. 342–344, 1991.
- A. C. Hindmarsh, P. M. Gresho, and D. F. Griffiths, “The stability of explicit euler timeintegration for certain finite difference approximations of the multi-dimensional advection-diffusion equation,” International Journal for Numerical Methods in Fluids, vol. 4, no. 9, pp. 853–897, 1984.
- E. Sousa, “The controversial stability analysis,” Applied Mathematics and Computation, vol. 145, no. 2-3, pp. 777–794, 2003.
- C. K. W. Tam and J. C. Webb, “Dispersion-relation-preserving finite difference schemes for computational acoustics,” Journal of Computational Physics, vol. 107, no. 2, pp. 262–281, 1993.
- C. Bogey and C. Bailly, “A family of low dispersive and low dissipative explicit schemes for flow and noise computations,” Journal of Computational Physics, vol. 194, no. 1, pp. 194–214, 2004.
- A. R. Appadu, “Comparison of some optimisation techniques for numerical schemes discretising equations with advection terms,” International Journal of Innovative Computing and Applications, vol. 4, no. 1, pp. 12–27, 2012.
- C. K. W. Tam and H. Shen, “Direct computation of nonlinear acoustic pulses using high-order finite differences schemes,” AIAA Paper 93-4325, 1993.