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

Authors: A. R. Appadu
Journal: Journal of Applied Mathematics
 (2013)

---

## 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 byh and k, respectively,
for which the Reynolds number is 2 or 4. Some errors are computed, namely, the error rate with respect to the L1 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 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.


---

## Body

## 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
(1)∂u∂t+A∂u∂x+B∂u∂y+C∂u∂z=α(∂2u∂x2+∂2u∂y2+∂2u∂z2).The coefficient of diffusivity is denoted byα and is computed as α=CT/pDp, where p, Dp, and CT denote the pressure, specific heat of the fluid at constant pressure, and thermal conductivity, respectively. Also A, B, and C are the velocity components of the fluid in the directions of x, y, and z, respectively. Equation (1) is also referred to as the convection-diffusion equation. The three terms A(∂u/∂x), B(∂u/∂y), and C(∂u/∂z) are called the advective or convective terms and the terms α(∂2u/∂x2), α(∂2u/∂y2), and α(∂2u/∂z2) are called the diffusive or viscous terms.In this paper, we consider the one-dimensional convection-diffusion equation given by(2)∂u∂t+a∂u∂x=α∂2u∂x2,
with a=1, α=0.01, 0≤x≤1, and 0<t≤T.We denote the spatial and temporal step sizes byh and k, respectively. The cfl number, c, is computed as ak/h, and the parameter, s, is obtained as αk/h2.The initial condition isu(x,0)=f(x), and boundary conditions are
(3)u(0,t)=g0(t),0<t≤T,u(1,t)=g1(t),0<t≤T,
where f, g0, and g1 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 a being constant [4]. This is the reason why numerical solution of (2) is important.The paper is organised as follows. In Section2, 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 k and h 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 k when h=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.


## 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 byAFM 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 isut+aux=αuxx, we will have dissipation, this is caused because of the term uxx, 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
(4)ξ=ξ1+Iξ2.
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 foru is obtained by substituting u by exp(I(w1t-θx)) where t and x represent time and space, respectively, θ is the wavenumber, and w1 is the dispersion relation [9].We then obtain(5)Iw1-aθI+αθ2=0,
where I=-1, which on simplification gives
(6)w1=aθ+αθ2I.Hence, the dispersion relation is given by(7)w1=aθ+αθ2I.The exact phase velocity is computed asℜ(w1)/θ which simplifies as θ. We next obtain the numerical phase velocity.From (4), we have ξ=ξ1+Iξ2. We can express ξ as exp(-bk) [9] where b is the exponential growth rate.Therefore, we haveexp(-bk)=ξ1+Iξ2 which implies
(8)b=1klog(ξ1-Iξ2ξ12+ξ22).
The numerical phase velocity is computed as ℑ(b)/θ and is equal to
(9)-1kθtan-1(ξ2ξ1).
The phase angle, w, is computed as w=θh where θ is the wavenumber and h 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
(10)RPE=-1kθatan-1(ℑ(ξ)ℜ(ξ)).Sincew=θh and c=ak/h, we can express (10) as
(11)RPE=-1cwtan-1(ℑ(ξ)ℜ(ξ)).If theRPE 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(12)1N∑i=1N‍(ui-vi)2,
where ui represents the analytical solution and vi represents the numerical (discrete) solution at a given grid point, i.The Total Mean Square Error can be expressed as(13)1N∑i=1N‍(ui-vi)2=1N∑i=1N‍(ui)2+1N∑i=1N‍(vi)2-2N∑i=1N‍uivi.
Next, one has
(14)1N∑i=1N‍(ui-ui¯)2=1N∑i=1N‍((ui)2-2uiui¯+(ui¯)2),1N∑i=1N‍(vi-vi¯)2=1N∑i=1N‍((vi)2-2vivi¯+(vi¯)2).
The Total Mean Square Error can be further expressed as
(15)1N∑i=1N‍(ui-ui¯)2+1N∑i=1N‍(vi-vi¯)2+2N∑i=1N‍uiui¯+2N∑i=1N‍vivi¯-1N∑i=1N‍(ui¯)2-1N∑i=1N‍(vi¯)2-2N∑i=1N‍uivi.
The expression in (15) can be rewritten as
(16)σ2(u)+σ2(v)+2(u-)2+2(v-)2-(u-)2-(v-)2-2N∑i=1N‍uivi,
where σ2(u) and σ2(v) denote the variance of u and v, respectively, u- and v- denote the mean values of u and v, respectively.Thus, the Total Mean Square Error is given by(17)σ2(u)+σ2(v)+((u-)2-2u-v-+(v-)2)+(2u-v--2N∑i=1N‍uivi)
which on further simplification yields
(18)σ2(u)+σ2(v)+(u--v-)2-2(1N∑i=1N‍uivi-ui¯vi¯).
Thus, we have
(19)1N∑i=1N‍(ui-vi)2=σ2(u)+σ2(v)+(u--v-)2-2Cov(u,v).
But the correlation coefficient, ρ, is given by Cov(u,v)/σ(u)σ(v). Hence, the Total Mean Square Error can be written as
(20)1N∑i=1N‍(ui-vi)2=σ2(u)+σ2(v)+(u--v-)2-2ρσ(u)σ(v),
which simplifies to
(21)1N∑i=1N‍(ui-vi)2=(σ(u)-σ(v))2+(u--v-)2+2(1-ρ)σ(u)σ(v).
On putting ρ=1, we get 2(1-ρ)σ(u)σ(v)=0. Thus, we define (2(1-ρ)σ(u)σ(v)) as the dispersion error as correlation coefficient in statistics is analogous with phase lag or phase lead in Computational Fluid Dynamics.Consequently,(σ(u)-σ(v))2+(u--v-)2 measures the dissipation error.We also obtain values of the error rate with respect to theL1 norm which is calculated as
(22)Enum=1N∑i=1N‍|ui-vi|,
where ui and vi are the exact and computed values, respectively, and N is the number of spatial grid points.


## 4. Choice of the Parametersh and k

We refer to [4] where three explicit methods are used to solve the partial differential equation
(23)ut+0.8ux=0.008uxx,
where
(24)u(x,t=0)=exp(-(x-2)28),g0(t)=2020+texp[-(5+4t)210(t+20)],g1(t)=2020+texp[-2(5+2t)25(t+20)].
Tests were carried out for three values of the cell Reynolds number, RΔ=c/s, namely, RΔ=2,4,8 [4].Sincec=0.8k/h and s=0.008k/h2, we can express RΔ in terms of h, in that case we have RΔ=100h. Thus, for RΔ=2,4,8, the corresponding values of h are 0.02, 0.04, and 0.08, respectively.Sincec=0.8k/h and h=0.02,0.04,0.08, we have the following relationships between c and k, namely, c=40k, c=20k, and c=10k.Then three values ofc were chosen as 0.16, 0.32, and 0.64, and then the corresponding values of k determined as 0.004, 0.008, 0.016, 0.032, 0.064. For these values of k, the number of time steps, M, are calculated as M=1/k, and hence M take the following values, namely, 250, 125, 62.5, 31.25, and 16.625, respectively. However, we note that M can only be an integer. Hence, an improvement can be made when choosing c and k while keeping RΔ=2,4,8 and h=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≤x≤1 and 0≤t≤1 [3], with a=1, α=0.01 and with the following initial and boundary conditions:
(25)u(x,0)=exp(-(x+0.5)20.00125),u(0,t)=0.0250.000625+0.02texp(-(0.5-t)2(0.00125+0.04t)),u(1,t)=0.0250.000625+0.02texp(-(1.5-t)2(0.00125+0.04t)).The exact solution is given by(26)u(x,t)=0.0250.000625+0.02texp(-(x+0.5-t)2(0.00125+0.04t)).The values ofh and k 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(27)ut+1.0ux=0.01uxx,
subject to boundary conditions given by (25).We consider two values forRΔ, say 2 and 4. Thus, we have RΔ=100h as a=1 and α=0.01. For RΔ=2 and 4, we have h=0.02 and 0.04, respectively. Hence, c=k/0.02 and c=k/0.04, and therefore we have c=50k and c=25k.We consider the case whenc=50k. If we choose c=0.25, 0.50, and 1.0, then the values taken by k are 0.01, 0.02, and 0.04, respectively.Next we considerc=25k. If we choose c=0.25, 0.5, and 1.0, we have k=0.01, 0.02, and 0.04, respectively.Hence, forh=0.02, the values taken by k are 0.005, 0.01, and 0.02. For h=0.04, k 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 ofk and h are possible, and we can also consider the case when c=2.0 instead of only the three cases, namely, c=0.25, 0.5, and 1.0.


## 5. Construction of Explicit and Implicit Finite Difference Methods

We can approximate∂u/∂x as
(28)(1-γ)(uin-ui-1n)+γ(ui+1n-uin)h,
or
(29)(1-γ)(uin+1-ui-1n+1)+γ(ui+1n+1-uin+1)h.Hence, an approximation for∂u/∂x is
(30)(1-ϕ)[(1-γ)(uin-ui-1n)+γ(ui+1n-uin)h]+ϕ[(1-γ)(uin+1-ui-1n+1)+γ(ui+1n+1-uin+1)h],
where h represents the spatial step size, ϕ and γ are the temporal and spatial weighting factors, respectively.An approximation for∂2u/∂x2 is
(31)ui+1n-2uin+ui-1nh2,
or
(32)ui+1n+1-2uin+1+ui-1n+1h2.Hence, a discretization for∂2u/∂x2 is
(33)(1-ϕ)[ui+1n-2uin+ui-1nh2]+ϕ[ui+1n+1-2uin+1+ui-1n+1h2].On plugging approximations for∂u/∂x and ∂2u/∂x2 as given by (30) and (33) into (2), we obtain a family of explicit and implicit numerical schemes given by
(34)uin+1=1A0×(A1ui-1n+A2uin+A3ui+1n+A4ui-1n+1+A5ui+1n+1),
where
(35)A0=1-ϕ[c(2γ-1)-2s],A1=(ϕ-1)[c(γ-1)-s],A2=1+(ϕ-1)[c(1-2γ)+2s],A3=(1-ϕ)[s-cγ],A4=ϕ[s+c(1-γ)],A5=ϕ[s-γc],
where c=ak/h and s=αk/h2.


## 6. Standard Schemes

### 6.1. Lax-Wendroff Scheme

The Lax-Wendroff scheme is given by(36)uin+1=12(2s+c+c2)ui-1n+(1-2s-c2)uin+12(2s-c+c2)ui+1n,
and is obtained on replacing ϕ by zero and γ by (1-c)/2, in (34).The modified equation is given by [14]
(37)ut+aux-αuxx+16ah2(1-c2-6s)uxxx+⋯=0,
and this indicates that the leading error terms are dispersive in nature.The amplification factor and the relative phase error are obtained as(38)ξ=1+(2s+c2)(cos(w)-1)-Icsin(w),RPE=1cwtan-1(csin(w)1+(2s+c2)(cos(w)-1)).The region of stability is0<s≤(1-c2)/2 [14]. Plots of the AFM and RPE, both versus the phase angle, w, for four combinations of values of k and h are shown in Figures 1(a) and 1(b). The combination k=0.01, h=0.04 is the least dissipative one. The scheme is not dispersive when k=0.01, h=0.04. Phase lag behaviour is observed when k=0.005, h=0.02 and k=0.01, h=0.02. Phase lead phenomenon occurs when k=0.02 and h=0.04.Plot ofAFM and RPE versus phase angle for the Lax-Wendroff scheme.
(a)
AFM
(b)
RPEThe results of our numerical experiment for the four combination of values ofh and k  are shown in Figures 2(a) and 2(b). We tabulate the errors in Table 1. The errors are the least when k=0.005 and h=0.02 and greatest when k=0.01 and h=0.04.Table 1
Errors for Lax-Wendroff scheme.
k
h
cfl
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.005
0.02
0.25
1.8166
×
1
0
-
4
5.8157
×
1
0
-
4
1.6348
×
1
0
-
4
6.3582
×
1
0
-
9
5.4502
×
1
0
-
8
0.01
0.02
0.50
7.3296
×
1
0
-
4
0.0024
1.6348
×
1
0
-
4
9.1960
×
1
0
-
8
8.9741
×
1
0
-
7
0.01
0.04
0.25
0.0021
0.0065
0.0011
9.1500
×
1
0
-
7
7.1622
×
1
0
-
6
0.02
0.04
0.50
1.2252
×
1
0
-
4
3.7946
×
1
0
-
4
3.7946
×
1
0
-
4
4.0896
×
1
0
-
9
2.4477
×
1
0
-
8Comparison of numerical results with exact results using Lax-Wendroff scheme at some values ofk and h.
(a)
(b)


### 6.2. Crank-Nicolson Scheme

The Crank-Nicolson method is obtained if we plugγ=1/2 and ϕ=1/2 into (34). A single expression for the scheme is
(39)uin+1=14(1+s)×(∫(c+2s)ui-1n+1-(c-2s)ui+1n+1+(c+2s)ui-1n-(c-2s)ui+1n+(4-4s)uin∫).The modified equation is given by(40)ut+aux-αuxx+112ah2(2+c2)uxxx+⋯=0,
and this indicates that the leading error terms are dispersive in nature.The amplification factor is given by(41)ξ=(B2D2-C22B22+C22)-I(B2C2+C2D2B22+C22),
and the RPE is computed as
(42)1cwtan-1(B2C2+C2D2B2D2-C22),
where B2=1+s-scos(w), C2=2csin(w), and D2=4+4scos(w)-4s.The scheme is unconditionally stable. We next plot the variation of the AFM versus phase angle for some values ofk and h in Figures 3(a) and 3(b). Plots of the RPE versus phase angle are depicted in Figures 4(a) and 4(b).Plot ofAFM versus phase angle for the Crank-Nicolson scheme.
(a)
h=0.02
(b)
h
=
0.04Plot ofRPE versus phase angle for the Crank-Nicolson scheme.
(a)
h=0.02
(b)
h
=
0.04In the case of Crank-Nicolson, the scheme is less dissipative ath=0.04 as compared to h=0.02 for all the four values of k, namely, 0.005, 0.01, 0.02, and 0.04. The combination h=0.04, k=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 k used when h=0.04. However, if we choose h=0.02, the dispersion character is much affected by the value of k. In general for h=0.02, the case k=0.02 is in general the least dispersive one.We tabulate the errors for the eight combinations ofh and k in Table 2, and we observe that the errors are the least when k=0.005 and h=0.02. The results of the numerical experiment are shown in Figures 5(a) and 5(b).Table 2
Errors for Crank-Nicolson scheme.
k
h
cfl
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.005
0.02
0.25
9.9859
×
1
0
-
4
0.0032
7.3954
×
1
0
-
4
1.4704
×
1
0
-
7
1.6929
×
1
0
-
6
0.01
0.02
0.50
0.0011
0.0035
7.3475
×
1
0
-
4
1.6161
×
1
0
-
7
2.0307
×
1
0
-
6
0.02
0.02
1.0
0.0015
0.0046
7.4486
×
1
0
-
4
2.2664
×
1
0
-
7
3.7049
×
1
0
-
6
0.04
0.02
2.0
0.0029
0.0092
0.0013
5.9447
×
1
0
-
7
1.5389
×
1
0
-
5
0.005
0.04
0.125
0.0037
0.0114
0.0020
2.2440
×
1
0
-
6
2.3758
×
1
0
-
5
0.01
0.04
0.25
0.0038
0.0116
0.0020
2.2981
×
1
0
-
6
2.4852
×
1
0
-
5
0.02
0.04
0.5
0.0042
0.0126
0.0021
2.5212
×
1
0
-
6
2.9471
×
1
0
-
5
0.04
0.04
1.0
0.0055
0.0162
0.0028
3.5178
×
1
0
-
6
5.1581
×
1
0
-
5Comparison of numerical results with exact results using Crank-Nicolson scheme at some values ofk and h.
(a)
(b)


## 6.1. Lax-Wendroff Scheme

The Lax-Wendroff scheme is given by(36)uin+1=12(2s+c+c2)ui-1n+(1-2s-c2)uin+12(2s-c+c2)ui+1n,
and is obtained on replacing ϕ by zero and γ by (1-c)/2, in (34).The modified equation is given by [14]
(37)ut+aux-αuxx+16ah2(1-c2-6s)uxxx+⋯=0,
and this indicates that the leading error terms are dispersive in nature.The amplification factor and the relative phase error are obtained as(38)ξ=1+(2s+c2)(cos(w)-1)-Icsin(w),RPE=1cwtan-1(csin(w)1+(2s+c2)(cos(w)-1)).The region of stability is0<s≤(1-c2)/2 [14]. Plots of the AFM and RPE, both versus the phase angle, w, for four combinations of values of k and h are shown in Figures 1(a) and 1(b). The combination k=0.01, h=0.04 is the least dissipative one. The scheme is not dispersive when k=0.01, h=0.04. Phase lag behaviour is observed when k=0.005, h=0.02 and k=0.01, h=0.02. Phase lead phenomenon occurs when k=0.02 and h=0.04.Plot ofAFM and RPE versus phase angle for the Lax-Wendroff scheme.
(a)
AFM
(b)
RPEThe results of our numerical experiment for the four combination of values ofh and k  are shown in Figures 2(a) and 2(b). We tabulate the errors in Table 1. The errors are the least when k=0.005 and h=0.02 and greatest when k=0.01 and h=0.04.Table 1
Errors for Lax-Wendroff scheme.
k
h
cfl
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.005
0.02
0.25
1.8166
×
1
0
-
4
5.8157
×
1
0
-
4
1.6348
×
1
0
-
4
6.3582
×
1
0
-
9
5.4502
×
1
0
-
8
0.01
0.02
0.50
7.3296
×
1
0
-
4
0.0024
1.6348
×
1
0
-
4
9.1960
×
1
0
-
8
8.9741
×
1
0
-
7
0.01
0.04
0.25
0.0021
0.0065
0.0011
9.1500
×
1
0
-
7
7.1622
×
1
0
-
6
0.02
0.04
0.50
1.2252
×
1
0
-
4
3.7946
×
1
0
-
4
3.7946
×
1
0
-
4
4.0896
×
1
0
-
9
2.4477
×
1
0
-
8Comparison of numerical results with exact results using Lax-Wendroff scheme at some values ofk and h.
(a)
(b)


## 6.2. Crank-Nicolson Scheme

The Crank-Nicolson method is obtained if we plugγ=1/2 and ϕ=1/2 into (34). A single expression for the scheme is
(39)uin+1=14(1+s)×(∫(c+2s)ui-1n+1-(c-2s)ui+1n+1+(c+2s)ui-1n-(c-2s)ui+1n+(4-4s)uin∫).The modified equation is given by(40)ut+aux-αuxx+112ah2(2+c2)uxxx+⋯=0,
and this indicates that the leading error terms are dispersive in nature.The amplification factor is given by(41)ξ=(B2D2-C22B22+C22)-I(B2C2+C2D2B22+C22),
and the RPE is computed as
(42)1cwtan-1(B2C2+C2D2B2D2-C22),
where B2=1+s-scos(w), C2=2csin(w), and D2=4+4scos(w)-4s.The scheme is unconditionally stable. We next plot the variation of the AFM versus phase angle for some values ofk and h in Figures 3(a) and 3(b). Plots of the RPE versus phase angle are depicted in Figures 4(a) and 4(b).Plot ofAFM versus phase angle for the Crank-Nicolson scheme.
(a)
h=0.02
(b)
h
=
0.04Plot ofRPE versus phase angle for the Crank-Nicolson scheme.
(a)
h=0.02
(b)
h
=
0.04In the case of Crank-Nicolson, the scheme is less dissipative ath=0.04 as compared to h=0.02 for all the four values of k, namely, 0.005, 0.01, 0.02, and 0.04. The combination h=0.04, k=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 k used when h=0.04. However, if we choose h=0.02, the dispersion character is much affected by the value of k. In general for h=0.02, the case k=0.02 is in general the least dispersive one.We tabulate the errors for the eight combinations ofh and k in Table 2, and we observe that the errors are the least when k=0.005 and h=0.02. The results of the numerical experiment are shown in Figures 5(a) and 5(b).Table 2
Errors for Crank-Nicolson scheme.
k
h
cfl
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.005
0.02
0.25
9.9859
×
1
0
-
4
0.0032
7.3954
×
1
0
-
4
1.4704
×
1
0
-
7
1.6929
×
1
0
-
6
0.01
0.02
0.50
0.0011
0.0035
7.3475
×
1
0
-
4
1.6161
×
1
0
-
7
2.0307
×
1
0
-
6
0.02
0.02
1.0
0.0015
0.0046
7.4486
×
1
0
-
4
2.2664
×
1
0
-
7
3.7049
×
1
0
-
6
0.04
0.02
2.0
0.0029
0.0092
0.0013
5.9447
×
1
0
-
7
1.5389
×
1
0
-
5
0.005
0.04
0.125
0.0037
0.0114
0.0020
2.2440
×
1
0
-
6
2.3758
×
1
0
-
5
0.01
0.04
0.25
0.0038
0.0116
0.0020
2.2981
×
1
0
-
6
2.4852
×
1
0
-
5
0.02
0.04
0.5
0.0042
0.0126
0.0021
2.5212
×
1
0
-
6
2.9471
×
1
0
-
5
0.04
0.04
1.0
0.0055
0.0162
0.0028
3.5178
×
1
0
-
6
5.1581
×
1
0
-
5Comparison of numerical results with exact results using Crank-Nicolson scheme at some values ofk and h.
(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 equationut+ux=αuxx has three subequations [16] which are given by
(43)ut+ux=0,(44)ux=αuxx,(45)ut=αuxx.Equations (43) and (44) have known exact finite difference scheme which are
(46)uin+1-uink+uin-ui-1nh=0,
with k=h and
(47)ui-ui-1h=α(ui+1-2ui+ui-1αh(exp(h/α)-1)),
respectively.A finite difference scheme that englobes the features of the two equations, namely, (43) and (44) is
(48)uin+1-uink+uin-ui-1nh=α(ui+1n-2uin+ui-1nαh(exp(h/α)-1)).On rearranging the terms in (48), we get the NSFD method which is [15, 16]
(49)uin+1=βui+1n+(1-α1-2β1)uin+(α1+β1)ui-1n,
where
(50)α1=kh,β1=α1exp(h/α)-1.The square of the modulus of the amplification factor is given by(51)|ξ|2=((1-α1-2β1)+(α1+2β1)cos(w))2+(α1sin(w))2.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 whenw=π. The square of the modulus of the amplification factor is given by
(52)|ξ|2=(1-2α1-4β1)2.We thus need,(53)(1-2α1-4β1)2≤1,
which implies that
(54)|1-2α1-4β1|≤1.Thus, for stability, we have the following inequality:(55)-1≤1-2α1-4β1≤1,
Which was simplified to
(56)0≤α1+2β1≤1.Sinceα1 and β1 are positive, α1+2β1≥0 is the trivial inequality. Hence, we consider the inequality
(57)α1+2β1≤1.Since,α1=k/h and β1=α1/(exp(h/α)-1), we have
(58)kh+2kh(exp(h/α)-1)≤1.For stability, we need the following condition:(59)k≤(exp(h/α)-1exp(h/α)+1)h.However, the stability condition in (59) is difficult to achieve. We use a Maclaurin’s series for exp(h/α), and therefore (59) reduces to
(60)kh+2αkh2≤1.We next consider the case whenw→0. When w→0, cos(w)≈1-(1/2)w2 and sin(w)≈w.Thus, (51) reduces to
(61)|ξ|2≈1+(-2β1+α12-α1)w2.We thus require(62)-2β1+α12-α1≤0.Using (50), (62) becomes
(63)-2kh+k(k-h)(exp(h/α)-1)h2(exp(h/α)-1)≤0.
From (63), we deduce -2kh+(k2-kh)(exp(h/α)-1)≤0, which on expansion and simplification gives (exp(h/α)-1)(k-h)≤0.Sinceexp(h/α)-1≥0, therefore,
(64)k≤h.Combining (60) and (64), we obtain (60), and therefore the region of stability is described by
(65)kh+2αkh2≤1.Case 1.
Forh=0.02 and α=0.01, using (60), we have k≤0.01. Hence, for h=0.02 and c=0.25, we have k=0.005. Also, for h=0.02 and c=0.5, we have k=0.010. However, if h=0.02 and c=1.0, we have k=0.020, but this combination will give rise to an unstable method.Case 2.
Forh=0.04 and α=0.01, using (60), the region of stability is given by k≤0.027. Therefore, for h=0.04, we consider k=0.01 and 0.02.Plots of theAFM and RPE versus phase angle are shown in Figures 6(a) and 6(b), respectively.Plot ofAFM and RPE versus phase angle for the NSFD scheme.
(a)
AFM
(b)
RPEThe NSFD scheme considered is an explicit one, and we have four combinations ofk and h, namely,(i)
k=0.005, h=0.02;(ii)
k=0.01, h=0.02;(iii)
k=0.01, h=0.04;(iv)
k=0.02, h=0.04.The scheme is least dissipative whenk=0.01, h=0.04 and k=0.005, h=0.02. The scheme is least dispersive when k=0.02, h=0.04. The scheme experiences both phase lead and phase lag behaviour, depending on the values of k and h. The results of our numerical experiment are shown in Figures 7(a) and 7(b).Comparison of numerical results with exact results using NSFD scheme at some values ofk and h.
(a)
(b)The modified equation is given by(66)ut+ux+12h(ac-1-1exp(h/α)-1)uxx+16h2(1-c2-6s)uxxx+⋯=0,
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 k=0.005 and h=0.02 and the greatest when k=0.02 and h=0.04.Table 3
Errors for NSFD scheme.
k
h
cfl
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.005
0.02
0.25
8.7288
×
1
0
-
4
0.0026
0.0026
8.0435
×
1
0
-
7
5.5063
×
1
0
-
7
0.01
0.02
0.50
0.0028
0.0085
-
0.0084
8.3500
×
1
0
-
6
5.8963
×
1
0
-
6
0.01
0.04
0.25
0.0068
0.0194
0.0192
4.8616
×
1
0
-
5
3.2639
×
1
0
-
5
0.02
0.04
0.50
0.0010
0.0032
0.0032
8.7192
×
1
0
-
7
1.0934
×
1
0
-
6Based on Tables1, 2, and 3, we can see that the Lax-Wendroff and the NSFD schemes are most effective when k=0.005 and h=0.02. The errors are smaller for the Lax-Wendroff as compared to NSFD scheme when k=0.005 and h=0.02.


## 8. Optimising Parameters in the Lax-Wendroff and NSFD Schemes

Our aim in this section is to compute an optimal value ofk for a given value of h, say h=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.


### 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 k 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(67)E=∫-ηη‍|θ*h-θh|2d(θh),
where the quantities θ*h and θh represent the numerical and exact wavenumbers, respectively. The dispersion error and dissipation error are calculated as |ℜ(θ*h)-θh| and |Im(θ*h)|, 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, θh, the effective/numerical wavenumber, θ*h, and define their integrated errors as
(68)E=∫(θh)l(θh)h‍|θ*h-θh|θhd(θh),
or
(69)E=∫ln(θh)lln(θh)h‍|θ*h-θh|d(ln(θh)).In Computational Fluid Dynamics for a particular method under consideration, the dispersion error is calculated as(70)|1-RPE|.
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:
(71)IETAM=∫0w1‍|1-RPE|2dw,IEBOGEY=∫0w1‍|1-RPE|dw.


### 8.2. Optimisation Procedure

Lax-Wendroff.   We consider the Lax-Wendroff scheme given by (36), with h=0.02. The amplification factor of the resulting method is
(72)ξLW=1-50k-2500k2+(50k+2500k2)cos(w)-50ksin(w)I,
and therefore the RPE is computed as
(73)RPELW=0.02kwtan-1(ℑ(ξLW)ℜ(ξLW)).A plot of the exactRPE versus w∈[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].Figure 8
Plot of RPE versusk versus w for the Lax-Wendroff scheme at h=0.02.We compute the following:(74)IETAM=∫01.1‍(RPELW-1)2dw,IEBOGEY=∫01.1‍|RPELW-1|dw.We plot the integrated errors versusk in Figures 9(a) and 9(b) and obtain the optimal value of k. We can also use the function NLPSolve from Maple to determine the value of k which minimise each of these two integrals. In the case of IETAM, we obtain
(75)k=0.00615029705055891978
while in the case of IEBOGEY, we are out with
(76)k=0.006112886302132816582.
We next validate whether this value of k computed does indeed minimise the errors by performing the numerical experiment using Lax-Wendroff with h=0.02 at some different values of k∈(0,0.01236) and then compare the errors. The errors are tabulated in Table 4, and we can see that indeed for k=1/164≈0.00601, all the five types of errors are least.Table 4
Errors for the Lax-Wendroff scheme forh=0.02.
k
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.001
8.2049
×
1
0
-
4
0.0026
6.3090
×
1
0
-
4
1.0506
×
1
0
-
7
1.1375
×
1
0
-
6
0.002
6.6803
×
1
0
-
4
0.0021
5.1728
×
1
0
-
4
7.1763
×
1
0
-
8
7.5203
×
1
0
-
7
1/333
5.1018
×
1
0
-
4
0.0016
4.0089
×
1
0
-
4
4.3376
×
1
0
-
8
4.3728
×
1
0
-
7
0.004
3.4840
×
1
0
-
4
0.0011
2.8318
×
1
0
-
4
2.1247
×
1
0
-
8
2.0307
×
1
0
-
7
0.005
1.8166
×
1
0
-
4
5.8157
×
1
0
-
4
1.6348
×
1
0
-
4
6.3582
×
1
0
-
9
5.4502
×
1
0
-
8
1/164
1.3952
×
1
0
-
5
4.3926
×
1
0
-
5
3.0697
×
1
0
-
5
1.1388
×
1
0
-
11
3.6197
×
1
0
-
10
1/143
1.6764
×
1
0
-
4
5.3721
×
1
0
-
4
-
7.8352
×
1
0
-
5
3.8582
×
1
0
-
9
4.7968
×
1
0
-
8
1/125
3.5162
×
1
0
-
4
0.0011
-
2.0134
×
1
0
-
4
1.9362
×
1
0
-
8
2.0827
×
1
0
-
7
1/111
5.4133
×
1
0
-
4
0.0017
-
3.2451
×
1
0
-
4
4.8426
×
1
0
-
8
4.9113
×
1
0
-
7
0.01
0.0011
0.0035
-
6.5674
×
1
0
-
4
2.1234
×
1
0
-
7
1.9665
×
1
0
-
6Plots of IETAM versusk and EBOGEY versus k for the Lax-Wendroff scheme when h=0.02.
(a)
IETAM
(b)
IEBOGEYNSFD.  We consider the NSFD scheme given by (49), with h=0.02. The amplification factor of the resulting method is
(77)ξNSFD=1+65.65176427k(cos(w)-1)-50(ksin(w))I,
and therefore the RPE is computed as
(78)RPENSFD=1α1wtan-1(ℑ(ξNSFD)ℜ(ξNSFD)),
where α1=k/h.A plot of the exactRPE versus e∈[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].Figure 10
Plot of RPE versusk versus w for the NSFD scheme at h=0.02.We compute the following:(79)IETAM=∫01.1‍(RPENSFD-1)2dw,IEBOGEY=∫01.1‍|RPENSFD-1|dw.We plot the integrated errors versusk in Figures 11(a) and 11(b) and obtain the optimal value of k. We can also use the function NLPSolve from Maple to determine the value of k which minimise each of these two integrals. In the case of IETAM, we obtain
(80)k=0.00611388415557632438
while in the case of IEBOGEY, we are out with
(81)k=0.00611348537281972832.
We next validate whether this value of k computed does indeed minimise the errors by performing the numerical experiment using NSFD with h=0.02 at some different values of k∈(0,0.01] and then compare the errors. The errors are tabulated in Table 5, and we can see that indeed for k=1/164≈0.00601, all the five types of errors are least.Table 5
Errors for NSFD scheme forh=0.02.
k
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.001
0.0035
0.0100
0.0100
1.2802
×
1
0
-
5
8.0987
×
1
0
-
6
0.002
0.0028
0.0083
0.0083
8.6039
×
1
0
-
6
5.5211
×
1
0
-
6
1/333
0.0022
0.0064
0.0064
5.1527
×
1
0
-
6
3.3593
×
1
0
-
6
0.004
0.0015
0.0045
0.0045
2.5390
×
1
0
-
6
1.6870
×
1
0
-
6
0.005
8.7288
×
1
0
-
4
0.0026
0.0026
8.0435
×
1
0
-
7
5.5063
×
1
0
-
7
1/164
1.1717
×
1
0
-
4
3.5591
×
1
0
-
4
3.4999
×
1
0
-
4
1.2220
×
1
0
-
8
1.2981
×
1
0
-
8
1/143
5.2775
×
1
0
-
4
0.0016
-
0.0015
3.0107
×
1
0
-
7
1.9762
×
1
0
-
7
1/125
0.00126
0.0038
-
0.0037
1.7163
×
1
0
-
6
1.1639
×
1
0
-
6
1/111
0.0020
0.0061
-
0.0060
4.3904
×
1
0
-
6
3.0429
×
1
0
-
6
0.01
0.0028
0.0085
-
0.0084
8.3500
×
1
0
-
6
5.8963
×
1
0
-
6Plots of IETAM versusk and EBOGEY versus k.
(a)
IETAM
(b)
IEBOGEY


## 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 k 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(67)E=∫-ηη‍|θ*h-θh|2d(θh),
where the quantities θ*h and θh represent the numerical and exact wavenumbers, respectively. The dispersion error and dissipation error are calculated as |ℜ(θ*h)-θh| and |Im(θ*h)|, 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, θh, the effective/numerical wavenumber, θ*h, and define their integrated errors as
(68)E=∫(θh)l(θh)h‍|θ*h-θh|θhd(θh),
or
(69)E=∫ln(θh)lln(θh)h‍|θ*h-θh|d(ln(θh)).In Computational Fluid Dynamics for a particular method under consideration, the dispersion error is calculated as(70)|1-RPE|.
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:
(71)IETAM=∫0w1‍|1-RPE|2dw,IEBOGEY=∫0w1‍|1-RPE|dw.


## 8.2. Optimisation Procedure

Lax-Wendroff.   We consider the Lax-Wendroff scheme given by (36), with h=0.02. The amplification factor of the resulting method is
(72)ξLW=1-50k-2500k2+(50k+2500k2)cos(w)-50ksin(w)I,
and therefore the RPE is computed as
(73)RPELW=0.02kwtan-1(ℑ(ξLW)ℜ(ξLW)).A plot of the exactRPE versus w∈[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].Figure 8
Plot of RPE versusk versus w for the Lax-Wendroff scheme at h=0.02.We compute the following:(74)IETAM=∫01.1‍(RPELW-1)2dw,IEBOGEY=∫01.1‍|RPELW-1|dw.We plot the integrated errors versusk in Figures 9(a) and 9(b) and obtain the optimal value of k. We can also use the function NLPSolve from Maple to determine the value of k which minimise each of these two integrals. In the case of IETAM, we obtain
(75)k=0.00615029705055891978
while in the case of IEBOGEY, we are out with
(76)k=0.006112886302132816582.
We next validate whether this value of k computed does indeed minimise the errors by performing the numerical experiment using Lax-Wendroff with h=0.02 at some different values of k∈(0,0.01236) and then compare the errors. The errors are tabulated in Table 4, and we can see that indeed for k=1/164≈0.00601, all the five types of errors are least.Table 4
Errors for the Lax-Wendroff scheme forh=0.02.
k
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.001
8.2049
×
1
0
-
4
0.0026
6.3090
×
1
0
-
4
1.0506
×
1
0
-
7
1.1375
×
1
0
-
6
0.002
6.6803
×
1
0
-
4
0.0021
5.1728
×
1
0
-
4
7.1763
×
1
0
-
8
7.5203
×
1
0
-
7
1/333
5.1018
×
1
0
-
4
0.0016
4.0089
×
1
0
-
4
4.3376
×
1
0
-
8
4.3728
×
1
0
-
7
0.004
3.4840
×
1
0
-
4
0.0011
2.8318
×
1
0
-
4
2.1247
×
1
0
-
8
2.0307
×
1
0
-
7
0.005
1.8166
×
1
0
-
4
5.8157
×
1
0
-
4
1.6348
×
1
0
-
4
6.3582
×
1
0
-
9
5.4502
×
1
0
-
8
1/164
1.3952
×
1
0
-
5
4.3926
×
1
0
-
5
3.0697
×
1
0
-
5
1.1388
×
1
0
-
11
3.6197
×
1
0
-
10
1/143
1.6764
×
1
0
-
4
5.3721
×
1
0
-
4
-
7.8352
×
1
0
-
5
3.8582
×
1
0
-
9
4.7968
×
1
0
-
8
1/125
3.5162
×
1
0
-
4
0.0011
-
2.0134
×
1
0
-
4
1.9362
×
1
0
-
8
2.0827
×
1
0
-
7
1/111
5.4133
×
1
0
-
4
0.0017
-
3.2451
×
1
0
-
4
4.8426
×
1
0
-
8
4.9113
×
1
0
-
7
0.01
0.0011
0.0035
-
6.5674
×
1
0
-
4
2.1234
×
1
0
-
7
1.9665
×
1
0
-
6Plots of IETAM versusk and EBOGEY versus k for the Lax-Wendroff scheme when h=0.02.
(a)
IETAM
(b)
IEBOGEYNSFD.  We consider the NSFD scheme given by (49), with h=0.02. The amplification factor of the resulting method is
(77)ξNSFD=1+65.65176427k(cos(w)-1)-50(ksin(w))I,
and therefore the RPE is computed as
(78)RPENSFD=1α1wtan-1(ℑ(ξNSFD)ℜ(ξNSFD)),
where α1=k/h.A plot of the exactRPE versus e∈[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].Figure 10
Plot of RPE versusk versus w for the NSFD scheme at h=0.02.We compute the following:(79)IETAM=∫01.1‍(RPENSFD-1)2dw,IEBOGEY=∫01.1‍|RPENSFD-1|dw.We plot the integrated errors versusk in Figures 11(a) and 11(b) and obtain the optimal value of k. We can also use the function NLPSolve from Maple to determine the value of k which minimise each of these two integrals. In the case of IETAM, we obtain
(80)k=0.00611388415557632438
while in the case of IEBOGEY, we are out with
(81)k=0.00611348537281972832.
We next validate whether this value of k computed does indeed minimise the errors by performing the numerical experiment using NSFD with h=0.02 at some different values of k∈(0,0.01] and then compare the errors. The errors are tabulated in Table 5, and we can see that indeed for k=1/164≈0.00601, all the five types of errors are least.Table 5
Errors for NSFD scheme forh=0.02.
k
E
num
max
|
u
e
-
u
c
|
Error at (0.5, 1.0)
Diss. error
Disp. error
0.001
0.0035
0.0100
0.0100
1.2802
×
1
0
-
5
8.0987
×
1
0
-
6
0.002
0.0028
0.0083
0.0083
8.6039
×
1
0
-
6
5.5211
×
1
0
-
6
1/333
0.0022
0.0064
0.0064
5.1527
×
1
0
-
6
3.3593
×
1
0
-
6
0.004
0.0015
0.0045
0.0045
2.5390
×
1
0
-
6
1.6870
×
1
0
-
6
0.005
8.7288
×
1
0
-
4
0.0026
0.0026
8.0435
×
1
0
-
7
5.5063
×
1
0
-
7
1/164
1.1717
×
1
0
-
4
3.5591
×
1
0
-
4
3.4999
×
1
0
-
4
1.2220
×
1
0
-
8
1.2981
×
1
0
-
8
1/143
5.2775
×
1
0
-
4
0.0016
-
0.0015
3.0107
×
1
0
-
7
1.9762
×
1
0
-
7
1/125
0.00126
0.0038
-
0.0037
1.7163
×
1
0
-
6
1.1639
×
1
0
-
6
1/111
0.0020
0.0061
-
0.0060
4.3904
×
1
0
-
6
3.0429
×
1
0
-
6
0.01
0.0028
0.0085
-
0.0084
8.3500
×
1
0
-
6
5.8963
×
1
0
-
6Plots of IETAM versusk and EBOGEY versus k.
(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 ofk and h. 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 k when h=0.02 for the Lax-Wendroff and NSFD schemes. We observe that when k≈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 convection-diffusion problems, and we can use appropriate optimisation techniques to choose parameters h and k for minimal numerical dispersion and numerical dissipation.



---

*Source: 734374-2013-03-07.xml*