Journal of Applied Mathematics

Volume 2013 (2013), Article ID 734374, 14 pages

http://dx.doi.org/10.1155/2013/734374

Research Article

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

A. R. Appadu

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.

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

Figure 1: Plot of and versus phase angle for the Lax-Wendroff scheme.

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 .

Table 1: Errors for Lax-Wendroff scheme.

Error at (0.5, 1.0)
Diss. error
Disp. error

0.005
0.02
0.25

0.01
0.02
0.50

0.0024

0.01
0.04
0.25
0.0021
0.0065
0.0011

0.02
0.04
0.50

(a)

(b)

Figure 2: Comparison of numerical results with exact results using Lax-Wendroff scheme at some values of 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).

(a)  

(b)

Figure 3: Plot of versus phase angle for the Crank-Nicolson scheme.

(a)  

(b)

Figure 4: Plot of versus phase angle for the Crank-Nicolson scheme.

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).

Table 2: Errors for Crank-Nicolson scheme.

Error at (0.5, 1.0)
Diss. error
Disp. error

0.005
0.02
0.25

0.0032

0.01
0.02
0.50
0.0011
0.0035

0.02
0.02
1.0
0.0015
0.0046

0.04
0.02
2.0
0.0029
0.0092
0.0013

0.005
0.04
0.125
0.0037
0.0114
0.0020

0.01
0.04
0.25
0.0038
0.0116
0.0020

0.02
0.04
0.5
0.0042
0.0126
0.0021

0.04
0.04
1.0
0.0055
0.0162
0.0028

(a)

(b)

Figure 5: Comparison of numerical results with exact results using Crank-Nicolson scheme at some values of and .

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

Figure 6: Plot of and versus phase angle for the NSFD scheme.

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)

Figure 7: Comparison of numerical results with exact results using NSFD scheme at some values of and .

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 .

Table 3: Errors for NSFD scheme.

Error at (0.5, 1.0)
Diss. error
Disp. error

0.005
0.02
0.25

0.0026
0.0026

0.01
0.02
0.50
0.0028
0.0085

0.01
0.04
0.25
0.0068
0.0194
0.0192

0.02
0.04
0.50
0.0010
0.0032
0.0032

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].

Figure 8: Plot of RPE versus versus for the Lax-Wendroff scheme at .

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.

Table 4: Errors for the Lax-Wendroff scheme for .

Error at (0.5, 1.0)
Diss. error
Disp. error

0.001

0.0026

0.002

0.0021

1/333

0.0016

0.004

0.0011

0.005

1/164

1/143

1/125

0.0011

1/111

0.0017

0.01
0.0011
0.0035

(a)   IETAM

(b) IEBOGEY

Figure 9: Plots of IETAM versus and EBOGEY versus for the Lax-Wendroff scheme when .

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].

Figure 10: Plot of RPE versus versus for the NSFD scheme at .

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.

Table 5: Errors for NSFD scheme for .

Error at (0.5, 1.0)
Diss. error
Disp. error

0.001
0.0035
0.0100
0.0100

0.002
0.0028
0.0083
0.0083

1/333
0.0022
0.0064
0.0064

0.004
0.0015
0.0045
0.0045

0.005

0.0026
0.0026

1/164

1/143

0.0016

1/125
0.00126
0.0038

1/111
0.0020
0.0061

0.01
0.0028
0.0085

(a)   IETAM

(b) IEBOGEY

Figure 11: Plots of IETAM versus and EBOGEY versus .

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.

References

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.