An optimal time-stepping algorithm for unsteady advection–diffusion problems

Journal of Computational and Applied Mathematics 294 (2016) 57–77

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

An optimal time-stepping algorithm for unsteady advection–diffusion problems✩ Naresh M. Chadha a , Niall Madden b,∗ a

Department of Mathematics, School of Basic Sciences and Research, Sharda University, India

b

School of Mathematics, Statistics and Applied Mathematics, National University of Ireland, Galway, Ireland

article

abstract

info

Article history: Received 1 September 2011 Received in revised form 22 March 2015 Keywords: Time dependent advection–diffusion problem Two-weight finite difference scheme Equivalent differential equation Optimal time-stepping

We consider the numerical solution of a linear time dependent advection–diffusion problem by an implicit two-weight, three-point finite difference scheme. We extend the scheme proposed by Chadha and Madden (2011), to incorporate an optimal time step selection algorithm for the method. The optimal values of the weights involved in the scheme have been obtained by using the notion of an equivalent differential equation (Warming and Hyett, 1974), and, by enforcing certain constraints due to von Neumann stability and a discrete maximum principle, we obtain sharp bounds for optimal time-stepping in terms of mesh width and problem data. Furthermore, we formulate the time step selection as a root-finding problem, and use this to show the optimal time step exists and is unique. The resulting method, based on optimal values of weights and optimal time-stepping, is of fifth-order in space, and third-order in time. The results of extensive numerical experiments are presented for several classic test problems for various values of the advection and diffusion parameters. These show that our proposed method is superior (quantitatively and qualitatively) to other schemes that can be represented within the framework of two-weighted schemes. We also present comparisons with other high-order methods, and an application of the method to a problem with non-smooth data. © 2015 Elsevier B.V. All rights reserved.

1. Introduction We consider the numerical solution of a one-dimensional advection–diffusion problem

∂Φ + LΦ = 0 , ∂t

L := a

∂Φ ∂ 2Φ − ε 2 for (x, t ) ∈ (0, l) × (0, T ], ∂x ∂x

(1.1a)

subject to the boundary and initial conditions

Φ (0, t ) = g0 (t ), Φ (x, 0) = f (x),

Φ (l, t ) = gl (t ),

t ∈ [0, T ],

x ∈ [0, l],

(1.1b) (1.1c)

were f , g0 and gl are known functions and are sufficiently smooth. It is assumed that ε and a, quantifying the diffusion and advection processes respectively, are positive (but arbitrary) constants. It is well known that the advection–diffusion problem (1.1) is quite fundamental to models describing a broad class of processes in physical, chemical, and biological ✩ Based upon works supported by the Science Foundation Ireland under Grant No. 08/RFP/CMS1205.

∗

Corresponding author. E-mail addresses: [email protected] (N.M. Chadha), [email protected] (N. Madden).

http://dx.doi.org/10.1016/j.cam.2015.07.029 0377-0427/© 2015 Elsevier B.V. All rights reserved.

58

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

sciences, see, e.g., [1,2]. A great variety of spatial discretization techniques in combination of various time integration methods have been studied for these type of problems, motivated by the fact that most conventional numerical methods are unsuitable when ε ≪ a, see, e.g., [3–6], and references therein. For the numerical solution of our model problem (1.1), we propose a two-weight, three-point finite difference scheme. The weights involved in the scheme have a well-defined role. One of the weights, θ , governs the weighting between two-levels for the time integration, and the second, φ , controls the weighted average of a two-point upwinding and a three-point central differencing scheme; see [7, Section 3.2] for details. The aim of the present work is to develop an optimal time step selection algorithm, with a view to improving the order of convergence of the scheme, without having to resort to using a wider stencil or nonlinear terms. To do this, we study the interaction between the stability region of a time integration method and the spectrum of the matrix associated with the spatial discretization. We demonstrate that there is a subtle interplay between the weights, θ and φ , and this can be exploited in order to obtain an optimal time-stepping for the scheme. The optimal values of the weights involved in the scheme have been obtained by using the notion of an equivalent differential equation. By further enforcing certain constraints due to von Neumann stability and a discrete maximum principle, we obtain sharp bounds for optimal time-stepping in terms of the spatial resolution, h, and of the advection and diffusion coefficients, a and ε , respectively. Furthermore, we formulate a root finding problem for an optimal time-stepping, which we prove exists and is unique. Combining the optimal values of weights with the optimal time step, yields a linear, compact, three-point, two-level scheme that is shown to be of fifth-order in space and third-order in time, while guaranteeing a stable, monotone solution. The authors do not know of any other linear finite difference scheme that has all these properties. Within the frame work of weighted schemes, several high-order, two-level, implicit finite difference schemes have been studied in recent past e.g., [8,9]. The scheme proposed in [9] involves three weights. However, the advection term is discretized using a two-point upwind operator, which is of first order. The five-point scheme of [8] is most similar to that which we propose; in particular we employ the idea it introduces of using a spatial discretization that is a weighted average of the central and upwinding difference for the advection term. However, that work is mainly concerned with explaining the performance of existing schemes. Also, the new scheme it proposes is explicit in nature and is not analyzed to the level of detail we provide. In [10], Rigal studies a general class of two-level, three-point implicit finite difference schemes for the advectiondominated problem, and shows that certain schemes independently proposed by several other authors, e.g., [11–13], can be analyzed within the framework of the proposed scheme. In fact, the result we derive (for optimal weights and time step) may be regarded as an extension of a result in [13, p188] for the θ -method applied to the diffusion problem, by generalizing it to the advection–diffusion problem while also ensuring stability and monotonicity. Other related work includes that of [14], where the stability of certain two- and three-layer schemes (two and three time-level, in our terminology) is investigated by employing the stability theory for time-dependent operator-difference schemes in Hilbert spaces, as originally suggested by Samarskii [15]. In the same paper [14], the stability boundaries of a two-level scheme applied to a standard heat equation in two dimensions is discussed. However, the scheme involves only one weight, and it is for a self-adjoint problem. In our case, two weights are involved in the scheme, and, due to the presence of the advection term, the operator is not self-adjoint. Thus, applying standard stability analysis developed for the two-layer scheme is not straightforward. As with the approach that we propose, the schemes of [8–10,16] can be considered as weighted finite difference methods. Ours is distinguished by the motivation for the weights, in particular the weighting between central differencing and upwinding for spatial discretization and the weighting between two-levels for the time integration, and how it leads to a way for choosing an optimal time-step for the model problem. We also mention the class of so-called weighted essentially non-oscillatory (WENO) schemes, that are very successful for solving advection dominated partial differential equations, see, e.g., [17,16]. These schemes use a dynamic set of stencils, where a nonlinear convex combination of lower order approximation polynomials adapts either to a high order approximation in smooth regions of the solution, or to a lower-order spatial discretization, in order to achieve high order formal accuracy in smooth regions while maintaining stable, nonoscillatory, computed solution near sharp discontinuities. It is emphasized that the weights involved in WENO are based on certain smoothness indicators. In general, these indicators are scaled sums of the square L2 -norm of all the derivatives of the relevant interpolation polynomial in the corresponding stencil. Our approach bears some resemblance to WENO methods, as it also involves a weighted average of low order difference operators to obtain an operator that gives high-order convergence and non-oscillatory solutions. However, these two approaches are fundamentally different. Firstly, one of weights, namely θ , controls the time-discretization, where as the weights in the WENO method are just for the spacial discretization. Secondly, our weights are chosen based on an analysis involving modified PDEs, whereas as with WENO, the weights are based on the smoothness indicators mentioned above. Notation. Throughout the paper, a, ε denote the advection and diffusion coefficients in (1.1), respectively. We employ uniform N

M

:= {x0 , x1 , . . . , xN } and Ω := {t0 , t1 , . . . , tM } respectively, with h = xj − xj−1 and N N τ = tn − tn−1 . By Ω we mean Ω /{x0 , xN }. The two weights in our finite difference scheme are θ and φ . Their optimal values are denoted by θopt and φopt . To simplify the notation, we further define s := ετ /h2 , c := aτ /h, ψ := 2s + φ c, K := ε/ah and σ := τ /h. For a matrix A = {aij }, by the inequality A ≥ 0, we mean that aij ≥ 0 ∀i, j. meshes in space and time: Ω

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

59

Outline. The article is organized as follows. In Section 2 we present a general two-weight finite difference scheme. Next, in Section 3, the scheme is analyzed in terms of stability, and certain useful bounds on the weights are obtained. Using the notion of an equivalent differential equation the optimal values of the weights, θopt and φopt , are obtained in Section 4.1. These optimal values of the weights are then combined with the bounds that were obtained in Section 3 and expressed in terms of c and s (defined above): this can be interpreted as a feasible region in the c–s plane. That is, if given values of a, ε , h and τ correspond values of c and s in this feasible region then the resulting scheme is guaranteed to produce a von Neumann stable, monotone numerical solution. We devote Section 5 to a detailed analysis of an approach for determining an optimal value of σ (the ratio of the time to space steps) that yields an optimal time-stepping for our model problem. The approach is based on eliminating higher-order terms in the truncation error. Sharp bounds on σ , in terms of a, ε , and h, are deduced in Section 5.1, and the existence and the uniqueness of an optimal σ are discussed in Section 5.2. In Section 6, we conclude with the results of extensive numerical experiments demonstrating the accuracy and convergence of the two-weight scheme. We show that it is superior to other schemes which can be expressed using the same three-point, two-level stencil. In addition, we provide a comparison with some classic high-order, five-point methods. In all cases it is observed that our proposed method offers significantly better accuracy. Section 6 concludes with a numerical study of how the method performs when applied to problems where the initial data are not smooth. 2. The two-weight finite difference scheme In this section, we present a general two-weight finite difference scheme which has been analyzed for von Neumann stability and monotonicity in [7]. We outline a geometrical interpretation of the role of the weights in the stability of the scheme, where necessary repeating some of the analysis of [7]. Define the standard discrete difference operators: D 0 uj =

δxx uj =

uj + 1 − uj − 1 2h

,

D − uj =

uj+1 − 2uj + uj−1 h2

,

uj − uj−1 h

and δt unj =

, unj +1 − unj

τ

.

We can now define a spatial finite difference operator, weighted with the parameter φ that balances between the standard second-order central difference operator D0 , which may be unstable for small ε , and the two-point upwinding operator D− that is stable, but only first-order accurate: LNφ uj := −εδxx + aδx uj ,





where δx := φ D− + (1 − φ)D0 .

We then introduce the parameter θ that weights the scheme between being implicit and explicit in nature, giving our general method for (1.1) as

  δt Φjn + LNφ θ Φjn+1 + (1 − θ )Φjn = 0,

for j = 1, . . . , N − 1, n = 0, . . . , M − 1.

(2.1)

Many standard schemes are encompassed by the framework of this two-weight scheme. For example, taking θ = 0 and φ = 0 gives the forward (i.e., explicit) Euler method combined with central differencing; θ = 0 and φ = aτ /h give the standard Lax–Wendroff scheme. The implicit Euler method with upwinding is given by taking θ = 1 and φ = 1. Crank–Nicolson type methods correspond to having θ = 1/2. The role of the weights in the stability of the scheme has been explained geometrically in Section 3.2 in [7]. 3. Analysis of the scheme and constraints on the weights We wish to show that our numerical scheme is convergent, and that it yields physically meaningful solutions. A prerequisite for convergence is consistency. This is discussed in Section 3.1. In addition, we require that the scheme is stable. Since this is a time dependent problem, even small numerical errors can be magnified and accumulate to dominate the entire solution. This can be avoided if the scheme is von Neumann stable, as we discuss in Lemma 1. Furthermore, it is well known that classical schemes for advection-dominated problems may also be unstable in the sense that, even for steady problems (and in the absence of any round-off error), computed solutions may exhibit spurious oscillations. This can be avoided if the scheme satisfies a discrete maximum principle. This is discussed in Lemmas 3 and 4. 3.1. Consistency In order to be convergent the numerical scheme must be both consistent and stable. Addressing consistency first, we write the scheme (2.1) in matrix form: T I Φ n +1 = T E Φ n ,

n = 0, . . . , M − 1,

(3.1a)

60

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

where T I (the ‘‘implicit’’ part of the scheme) and T E (the ‘‘explicit’’) are N × N tridiagonal matrices. Written out in detail this is

Φ0n+1 = g0 (tn+1 ), TjI,j−1 Φjn−+11 + TjI,j Φjn+1 + TjI,j+1 Φjn++11 = TjE,j−1 Φjn−1 + TjE,j Φjn + TjE,j+1 Φjn+1 ,

ΦNn+1

j = 2, . . . , N − 1,

(3.1b)

= gl (tn+1 ),

where the only nonzero entries in the first and last rows of T I are T1I ,1 = TN1,N = 1. For all other rows

θ

θ

TjI,j−1 = − (c + ψ), TjI,j = 1 + θ ψ, TjI,j+1 = (c − ψ), 2 2 1 − θ 1−θ TjE,j−1 = (c + ψ), TjE,j = 1 − (1 − θ )ψ, TjE,j+1 = (−c + ψ). 2 2 The scheme is consistent since, for j = 2, . . . , N − 1, we have that

(3.1c) (3.1d)

TjI,j−1 + TjI,j + TjI,j+1 = TjE,j−1 + TjE,j + TjE,j+1 . 3.2. Stability The von Neumann stability of θ -methods is well-studied. In our notation, it may be expressed as follows. Lemma 1 (von Neumann Stability). The two-weight scheme (2.1) is von Neumann stable for any θ ≥ 1/2. For θ < 1/2, it is stable if

(1 − 2θ )c 2 − ψ ≤ 0,

and

Proof. See [7, Section 3.3].

  ψ (1 − 2θ )ψ − 1 ≤ 0.

(3.2)



We next discuss the stability of the scheme in the sense of excluding spurious oscillations. We need two ingredients. The first is that the matrix T I is an M-matrix. The second, which follows directly from this, and that T E ≥ 0, is that the scheme satisfies a discrete maximum principle. We now explain these terms. Definition 1 (M-Matrix). The matrix M is an M-matrix if it is nonsingular, has non-positive off-diagonal entries, and M −1 ≥ 0. There are many characterizations on M-matrices. We use the following, referring to [18, Thm. 7.2.5] for proof. Lemma 2. The nonsingular matrix M is an M-matrix if it is irreducible, if mii > 0,



mij ≥ 0,

mij ≤ 0, (diagonal dominance)

j

and at least one row is strictly diagonally dominant. Lemma 3. The matrix T I as given in (3.1a) is an M-matrix if

θ ≥ 0 and ψ ≥ c ,

where ψ = 2s + φ c .

(3.3a)

Also, if in addition to the conditions in (3.3a), we have that

θ ≤ 1 and ψ ≤

1 1−θ

,

(3.3b)

then T E ≥ 0. Proof. The matrix T I is irreducible since the domain is connected. Also the first and last rows are strictly diagonally dominant. Then for all other rows we need that TjI,j−1 ≤ 0,

TjI,j > 0,

TjI,j+1 ≤ 0,

and TjI,j ≥ −TjI,j−1 − TjI,j+1 .

Note that c = aτ /h > 0, and further that ψ ≥ c > 0. Then if θ ≥ 0, we have TjI,j−1 ≤ 0, TjI,j > 0, and TjI,j+1 ≤ 0. Also,

−TjI,j−1 − TjI,j+1 = θ ψ < 1 + θ ψ = TjI,j , so T I is diagonally dominant. Using ψ ≥ c ≥ 0 and θ ≤ 1 it is easy to verify that TjE,j−1 and TjE,j+1 are non-negative. The final assumption in (3.3b) gives that TjE,j ≥ 0. 

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

61

Note that roughly (i.e., momentarily neglecting boundary conditions) the scheme (2.1) can be written as φ n+1 = (T I )−1 T φ n . Therefore, if the initial and boundary conditions in (1.1) are non-negative, and the weights θ and φ satisfy (3.3), then the solution to (2.1) is always non-negative. This is a clearly as desirable property since the solution to (1.1) represents a physical quantity that cannot be negative. We now make this notion more formal. E

Definition 2 (Discrete Maximum Principle). Let us write the scheme (2.1) as

(Lh,τ un+1 ) := T I unj +1 − T E unj for n = 0, . . . , M − 1, N

M

N

for an arbitrary mesh function u on Ω × Ω . Then if y and z are mesh functions on Ω × Ω (a) |y | ≤ z ; and for n = 0, . . . , M, | | ≤ 0

yn0

0

z0n

ynN

and |

(b) |(Lh,τ yn+1 )| ≤ Lh,τ z n+1 , for n = 0, . . . , M − 1,

|≤

zNn ;

M

such that

and

then |y| ≤ z. Lemma 4. If conditions (3.3) hold, then the scheme (2.1) satisfies the discrete maximum principle. Proof. From Lemma 3 the conditions in (3.3) give that T I is an M-matrix, and that T E ≥ 0. Now apply [6, Lemma 3.12] to get the desired result.  Remark 1. One of the conditions in Lemma 3, namely ψ ≥ c can also be obtained using the notion of monotonicity or using the standard eigenvalue analysis [7]. This gives a necessary condition

φ ≥1−

2s c

,

(3.4)

for the difference scheme (2.1) to produce non-oscillatory computed solution. We finish this section by collecting the conditions in Lemmas 1 and 3 and expressing them as bounds on the weight φ in terms of θ , s and c. Theorem 1. When a ̸= 0, a general two-weight scheme in the form (2.1) is von Neumann stable and satisfies the discrete maximum principle for a given value of θ in the range 0 ≤ θ ≤ 1, if the weight φ ≥ 0 satisfies the following bounds c−

2s c

≤φ≤

 max 1 −

1− 1−

2s c 2s c

2s c

1 − 2s c

,

for θ = 0;

, (1 − 2θ )c − 

≤φ≤

1

≤ φ,

for θ = 1.

c

1 1−θ

2s



c

(3.5a)

≤φ≤



− 2s ,

for

1 c 1 2



1 1−θ



− 2s ,

for 0 < θ <

1 2

;

(3.5b)

≤ θ < 1;

(3.5c) (3.5d)

Proof. For (3.5a), the bounds in (3.2) give that c − 2s/c ≤ φ ≤ (1 − 2s)/c. Since the scheme is explicit, the only other bound needed is (3.3b), which again gives φ ≤ 1/c − 2s/c.  For (3.5b), we note that the two inequalities in (3.2) can be rearranged to give φ ≥ (1 − 2θ )c − 2s/c, and φ ≤ 1/(1 − 2θ ) − 2s /c, respectively, where in the latter case we have also used that φ ≥ 0, and so ψ ≥ 0. Also, ψ ≥ c in (3.3a) can be rewritten as φ ≥ 1 − 2s/c. Finally, (3.3b) gives that φ ≤ 1/(1 − θ ) − 2s /c, a stricter condition than that coming from (3.2). For the remaining case, the scheme is unconditionally stable, so we need only the conditions (3.3). These are easily rearranged to get (3.5c). For the final case, the scheme is fully implicit, so we need only (3.3a). 





Remark 2. In the case where a = 0, (1.1) reduces to a pure diffusion problem, and φ is not present in the difference scheme. So only conditions on θ are required. In order to satisfy the stability conditions (3.2) and (3.3) in this case, we have 1−

1 2s

≤ θ ≤ 1,

which follows from ψ ≤ 1/(1 − θ ) in (3.3b).

62

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

4. An optimal two-weight scheme The analysis of Section 3 applies to any scheme within the framework of a two-weight scheme (2.1). We now consider particular values of the weights θ and φ which optimize the accuracy of the scheme for given ε , a, h and τ . These values of the weights were derived in [7, Section 4] using the notion of an equivalent differential equation. For the sake of completeness, we present these optimal weights and their derivation in Section 4.1. In Section 4.2 we combine them with the analysis presented in Section 3 to obtain a particular region of interest in the c–s plane. For a given set of values of a, ε and h, any τ that corresponds to a point in this region will produce an oscillation-free, von Neumann stable solution. Thus, we shall refer to this as a feasible region in c–s plane as it satisfies all the constraints due to positivity and the stability of the scheme as given in Theorem 1. 4.1. Optimal values of the parameters In order to calculate optimal values of the weights, we employ the notion of a Modified Equivalent Partial Differential Equation (MEPDE), first introduced by Warming and Hyett [19], and there after used in many articles in the literature, see, e.g., [8,20,9,12,21] and references therein. The idea is to replace the finite differences in (2.1) with Taylor’s series, which yield an infinite-order PDE whose solution agrees exactly with (2.1). Then (1.1) and its derivatives are used to eliminate all but the first time derivative, so that the terms in modified PDE which are different from those in (1.1) involve only spatial derivatives. Applying the technique to the scheme (2.1), yields ∞ ∂Φ ∂Φ ∂ 2Φ  ∂ qΦ +a −ε 2 + Kq (c , s) q = 0, ∂t ∂x ∂x ∂x q =2

(4.1)

where

K2 = K3 = K4 =

ah

κ2 ,

2!

ah2 3! ah3 4!

κ2 = −(φ − c (1 − 2θ ));

κ3 ,

(4.2a)

κ3 = 1 − 6s − c 2 + 6sθ + 3θ ψ + 3c 2 θ;   1 κ4 = (12s2 + 12c 2 s + c 4 ) − 4θ c 2 (1 + 6s + c 2 ) − ψ(1 + 12sθ + 6c 2 θ ) .

κ4 ,

c

(4.2b) (4.2c)

The infinite sum in (4.1) represents the truncation error in the setting of an MEPDE. We require that the function Φ be sufficiently differentiable for the terms K2 , K3 and K4 to exist. (We investigate the ramifications of Φ not being sufficiently smooth in Section 6.4.) The first two terms, (4.2a) and (4.2b), correspond respectively to numerical dissipation and dispersion [3, Section 9.2]. Moreover, the weights involved in the scheme can be chosen such that the leading dominant error terms in (4.1) can be eliminated. This will lead to more accurate methods. We choose θ and φ so that κ2 = κ3 = 0, and we get

φ = c (1 − 2θ ),

θ=

3(c 2 + 2s) ±



3(2c 2 + c 4 + 12s2 ) 6c 2

.

Since the greater of these two values for θ yields a negative value of φ , which is not appropriate in view of the geometrical interpretation given in [7, Section 3.2], the optimal values of the weights are:

φopt = c (1 − 2θ ),

θopt =

3(c 2 + 2s) −



3(2c 2 + c 4 + 12s2 ) 6c 2

.

(4.3)

The expression for θopt can be rewritten as

θopt = Since

1 2

3(2c 2 + c 4 + 12s2 ) − 6s

 −

6c 2

.

3(2c 2 + c 4 + 12s2 ) > 6s, we get θopt < 1/2.



Remark 3. If the Modified Partial Differential Equation [22, Section 7.7] approach is used, the coefficients corresponding to those in (4.2) are given by

κ˜ 2 = −(φ − c (1 − 2θ ));

(4.4a)

κ˜ 3 = 1 − 6s + 2c − 3c φ − 6c θ + 6c θ + 6ψθ ;   1 κ˜ 4 = ψ 36θ c (1 − θ ) + (3ψ(1 − 2θ ) − 1) − 12c − 12θ c 3 (2θ 2 − 3θ + 2) + 2c (2 + 3c 2 − 4θ ). 2

2

c

2 2

(4.4b) (4.4c)

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

63

It is interesting to note that even though the coefficients in (4.2) and (4.4) are different, setting κ˜ 2 = κ˜ 3 = 0 yields the same optimal values for θopt and φopt . Furthermore, for these optimal values, we get that κ˜ 4 = κ4 . This does not hold for higher order terms in the modified and equivalent PDE expansions. Remark 4. For a pure diffusion problem, i.e., a = 0, there is just one weight, θ . In this case, optimal value θ , denoted by θdiff , follows from the fourth order term in equivalent differential equation is and given by

θdiff =

1 2

−

1 12s

.

(4.5)

This is exactly the same weight derived by Richtmyer and Morton in [13, p188] for the heat equation. Remark 5. Careful calculation can show that if, for example, τ = h, then K4 in (4.2c) is O (h2 ). If one takes τ = h2 , then one gets K4 = O (h4 ). Moreover, we devote Section 5 to showing how to optimize the time step so that K4 = 0, thus obtaining a high-order scheme. 4.2. Stability of the scheme with optimal weights Given the optimal values θopt and φopt in (4.3), we now express the bounds due to Lemmas 1 and 3 in terms of c = aτ /h and s = ετ /h2 . This leads to a range of values of s and c for which the scheme yields a von Neumann stable solution that satisfies the discrete maximum principle, which we refer to as the feasible region in the c–s plane. This region is determined by the following inequalities:

√ φ≤1≡ 0≤θ ≡

ψ ≥c≡

−6s +

6c 2 + 3c 4 + 36s2 3c

3c 2 + 6s −

≤ 1,

(4.6)

√

6c 2 + 3c 4 + 36s2 6c 2

≥ 0,

(4.7)

1 2 6c + 3c 4 + 36s2 ≥ c , 3

(4.8)

√

ψ≤

1 1−θ

≡

(c 2 − 2s) 6c 2 + 3c 4 + 36s2 + 12s2 + c 4 − 4c 2 ≤ 0, √ 3c 2 − 6s + 6c 2 + 3c 4 + 36s2

(1 − 2θ )c 2 − ψ ≤ 0 ≡ −2s ≤ 0, √ √   6c 2 + 3c 4 + 36s2 (2s 6c 2 + 3c 4 + 36s2 + c 2 − c 4 − 12s2 ) ψ (1 − 2θ )ψ − 1 ≤ 0 ≡ − ≤ 0. 2 9c

(4.9) (4.10) (4.11)

These inequalities are plotted in Fig. 1. The shaded region is the feasible region in the c–s plane, i.e., any values of a, ε , h and τ that give a point (c , s) in this region will lead to a numerical solution that is von Neumann stable and free from spurious oscillations. 4.3. Further refinement of the feasible region One may use the bounds represented by the shaded region of Fig. 1 to check if a given set of values of ε , a, h and τ yield a scheme guaranteed to be von Neumann stable and to satisfy the discrete maximum principle. However, our goal is to be able to express these bounds in terms of τ for given ε , a, h. In order to achieve this, we further simplify the inequalities (4.6)–(4.11). It can easily be verified that the denominators in each of the inequalities (4.6)–(4.11) are always positive. Thus

φ ≤ 1 yields c 3 ≤ c + 12s, 0≤θ

(4.12)

yields c 2 + 6s ≥ 1,

(4.13)

ψ ≥ c yields c + 12s ≥ c . 4

2

2

(4.14)

The inequality due to stability criteria (4.11) can be rearranged as



2s 6c 2 + 3c 4 + 36s2 ≥ c 4 + 12s2 − c 2 .

(4.15)

The right hand side in (4.15) is positive, provided the condition (4.14) is satisfied. Assuming (4.14) is satisfied, then one can square both the sides in (4.15) and simplify to get 48s2 + 2c 4 ≥ c 2 (1 + c 4 + 12s2 ).

(4.16)

64

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

Fig. 1. A feasible region in the c–s plane, satisfying all the constraints (4.6)–(4.11). Any combination of (h, τ , ε, a) that falls into the shaded region produces a non-oscillatory, von Neumann stable solution.

Finally, consider the inequality (4.9). The denominator is always positive, thus the inequality reduces to

 (c 2 − 2s) 6c 2 + 3c 4 + 36s2 + 12s2 + c 4 − 4c 2 ≤ 0. This can be rearranged as





c 2 ( 6c 2 + 3c 4 + 36s2 − 3) − (2s 6c 2 + 3c 4 + 36s2 − c 4 − 12s2 + c 2 ) ≤ 0.

(4.17)

Note that, in view of stability condition given by (4.15), the second part of the above inequality is positive. If √ 6c 2 + 3c 4 + 36s2 − 3 ≤ 0, then inequality (4.17) is readily satisfied. This implies that the scheme satisfies the discrete maximum principle if the stability constraint (4.15) is satisfied along with 6c 2 + 3c 4 + 36s2 − 9 ≤ 0.

(4.18)

√

A feasible region which satisfies all the desired bounds along with 6c 2 + 3c 4 + 36s2 − 3 ≤ 0 is shown in Fig. 2 as a shaded region. Note that (4.16) follows from (4.15) and it may be interesting to compare the corresponding graphs. The former gives another subregion inside a region due to ψ ≥ c (≡ c 4 + 12s2 ≥ c 2 ). However, if c 4 + 12s2 ≥ c 2 is satisfied, then this apparent subregion is not going to affect the boundaries of the feasible region shown in Fig. 2. Note that the factors involved in (4.13), (4.14), and (4.18), namely c 2 + 6s − 1, c 4 + 12s2 − c 2 , and 6c 2 + 3c 4 + 36s2 − 9 are all polynomial functions in c√and s. Thus, they may be rearranged to get an expression for τ in terms of a, ε , and h. However, this is not the case when 6c 2 + 3c 4 + 36s2 − 3 ≥ 0. For instance one may rearrange the inequality (4.17) as follows:

(c + 2c − 12s )(c + 2c + 12s ) + 4

2

2

4

2

2



6c 2

 +

3c 4

+

36s2

 −3c + 2s(c + 12s − c ) ≤ 0, 4

4

2

2

by squaring both the factors in (4.17), which is not of much practical use. Thus, for our further analysis we shall only consider the inequalities which construct the feasible region shown in Fig. 2. 5. Towards optimal time-stepping We now present an approach which, for given values of a, ε , and h, gives an optimal value of τ that one may use in the difference scheme (2.1) taking the optimal values of θ and φ from (4.3). In Section 5.1 we express the inequalities for constructing the ‘‘feasible’’ region shown in Fig. 2 in terms of the ratio of the time step to the space mesh width: σ := τ /h. This gives us a range of σ with sharp bounds. Any value of σ in this range is guaranteed to give a scheme that is von Neumann stable solution and satisfies the discrete maximum principle. Next, we express κ4 from (4.2c) in terms of σ . We aim to find a value of σ that solves κ4 = 0 within the range of σ obtained in Section 5.1. We shall refer to this problem as our root-finding problem. The existence and the uniqueness of the solution of this problem is discussed in Section 5.2.

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

65

Fig. 2. A subregion bounded by (4.12)–(4.16), and 6c 2 + 3c 4 + 36s2 ≤ 9.

5.1. Sharp bounds on σ for optimal time-stepping We now express the inequalities of Section 4.3 in terms of σ = τ /h, in which we might seek an optimal τ for given value of a, ε , and h. To simplify the notation we define K := ε/(ah). Lemma 5. The optimal weights given by (4.3) satisfy the conditions θ ≥ 0, φ ≤ 1, and ψ ≥ c if

σ ∈ σ ∈



1 a

(−3K +



  1

1 − 12K 2 ,

a

1√ 1 + 12K , a



1 + 9K 2 ),

1√ 1 + 12K , a



for K ≥ 1/4,

(5.1a)

for 0 ≤ K < 1/4.

(5.1b)

Proof. The upper bound follows from (4.12) and corresponds to φ ≤ 1. The inequality (4.13), which corresponds to θopt ≥ 0 yields

σ ≥

1 a

−3K +





1 + 9K 2 .

√

The inequality (4.14) comes from ψ ≥ c, and can be expressed as a2 σ 2 + 12K 2 − 1 ≥ 0. It holds for any K ≥ 1/ 12, and otherwise can be expressed as

σ ≥

1 1 − 12K 2 . a

√

2 It is easily checked that, √ √ if K ≥ 1/16, then −3K + 2 1 − 12K > −3K + 1 + 9K 2 . 

1 + 9K 2 ≥

√

1 − 12K 2 . Otherwise, if 0 < K 2 < 1/16, then

Theorem 2. The scheme (2.1) with θ = θopt , and φ = φopt is von Neumann stable and has θ ≥ 0 if

σ ∈ σ ∈



1 a

(−3K +

  1 a

1−



1 + 9K 2 ),

12K 2

,

1

1



a

1 − 6K 2 + 6K







 1−

a

6K 2

+ 6K

K2



K2 + 1 ,

+1 ,

for K ≥ 1/4,

for 0 ≤ K < 1/4.

(5.2a)

(5.2b)

Proof. We need to show that if σ lies in the region specified by (5.2), then the conditions (4.7), (4.10), and (4.11) all hold. The second of these, (4.10), is trivial. As we have shown in Section 4.3, the condition (4.11) can be simplified as (4.16), providing that (4.14) also holds. Considering first (4.16), that expression can be rearranged as



σ − 2

1 − 6K 2 + 6K a2

√

K2 + 1



σ − 2

1 − 6K 2 − 6K a2

√

K2 + 1



≤ 0.

66

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

Requiring the first factor to be negative gives the upper bound:

σ ≤

1

 1 − 6K 2 + 6K

a



K 2 + 1.

(5.3)

For the second to be positive, we need that 1

σ2 ≥ If K 

2

a2

1 − 6K 2 − 6K



K2 + 1 .



(5.4)

≥ 1/48, this just corresponds to the trivial condition σ ≥ 0. Otherwise we must explicitly require that σ ≥ √ 

1 − 6K 2 − 6K K 2 + 1 /a. It remains to show that (4.14) and (4.7) both hold, though the  first of these also requires (4.12). These three conditions

will hold if σ satisfies the bounds in (5.1). It is easy to see that follows that

 1 − 6K 2 + 6K





1 + 9K 2 ,



K 2 + 1 ≤ max −3K +



1 − 6K 2 + 6K

√

K2 + 1 ≤

√

1 − 12K 2 for any K , and so it

1 − 12K 2 .



So we have now established the lower bounds in (5.2). Furthermore, for any K ,

 1 − 6K 2 + 6K



1 + K2 <

√

giving the upper bounds in (5.2).

1 + 12K ,



We now collect the above results. Theorem 3. The two-weight scheme (2.1) with the optimal weights (4.3) satisfies all the desired conditions for von Neumann stability and the discrete maximum principle, given by (3.2), and (3.3), respectively, if

σ ∈ σ ∈



1 a

(−3K +

  1 a

1−



1 + 9K 2 ),

12K 2

,

1

1 a



  −1 − 6K 2 + 2 1 + 9K 4 + 3K 2 ,



a

−1 −

6K 2





+2 1+

9K 4

+

3K 2

,

for K ≥ 1/4,

for K < 1/4.

(5.5a)

(5.5b)

Proof. We need conditions (5.1) and (5.2) to hold, in addition to having ψ ≤ 1/(1 − θ ). Proving that (5.1) hold, this last condition is expressed in (4.18). Rewriting this inequality in terms of σ , we get

√

 σ + 2

1 + 6K 2 + 2 1 + 9K 4 + 3K 2

√

 σ + 2

a2

1 + 6K 2 − 2 1 + 9K 4 + 3K 2 a2

 ≤ 0.

Since the first factor is positive, this reduces to requiring that

σ ≤

1 a



 −1 − 6K 2 + 2 1 + 9K 4 + 3K 2 . 

√

It can be easily verified that −1 − 6K 2 + 2 1 + 9K 4 + 3K 2 < combining (5.1), (5.2), and (5.6). 

(5.6)



1 − 6K 2 + 6K

√

K 2 + 1. The desired result follows from

It is instructive to visualize the bounds in (5.1), (5.2), and (5.5) as presented in Fig. 1 and Fig. 2. For example, there are two possible lower bounds in each of these that correspond to the conditions θ ≥ 0 and ψ ≥ c; as we can see from Fig. 2, the lower bound in the feasible region is defined in a piecewise manner. On the other hand, the two possible upper bounds expressed in (5.1) and (5.2) are subsumed by the stricter bound coming from (5.6). This is because, as demonstrated in Fig. 2, the condition (4.18) is stronger than the other conditions which produce potential upper bounds for σ for any value of c and s. Remark 6. For a pure diffusion problem, one may combine the bounds in Remark 2 with the optimal value of θ in Remark 4 to obtain the bound σ ≤ 5h/(6ε) to satisfy the desired properties, given by (3.2), and (3.3).

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

67

5.2. The existence and uniqueness of an optimal time step As mentioned above, we seek a time step τ that eliminates the fourth order term in (4.2c). That is, we wish to find σ that satisfies the bounds (5.5) and solves κ4 = 0, which expressed in terms of σ is



F (σ ) := 6K (a2 σ 2 − 12K 2 ) + (−1 + a2 σ 2 + 12K 2 ) 6 + 3a2 σ 2 + 36K 2 = 0.

(5.7)

This is a relatively simple nonlinear problem that can be easily solved using a Newton iteration method. However, it is instructive to analyze the existence and the uniqueness of the solution (5.7) in the desired range(s) of σ , see, (5.5). To simplify the notation, we define

  1 −3K + 1 + 9K 2 , σl2 := 1 − 12K 2 ; a a   1 σr := −1 − 6K 2 + 2 1 + 9K 4 + 3K 2 .

σl1 :=

1

(5.8)

a

Theorem 4. There exists a unique solution to our root finding problem, i.e., the equation F (σ ) = 0, given in (5.7), in the following intervals:

[σl1 , σr ] if K 2 ≥ 1/16, [σl2 , σr ] if 1/24 ≤ K 2 ≤ 1/16. Proof. The existence of the solution. Note that the expressions F (σl1 ), F (σl2 ) and F (σr ) depend only on K , and not on a. Rearranging these expressions, it can be easily seen that for any K > 0, F (σl1 ) < 0 < F (σr ), and so (5.7) has a solution for K 2 ≥ 16. For the second interval (σl2 , σr ), we evaluate F (σl2 ) := 6K (1 − 24K 2 ); and F (σr ) > 0. This implies that we have a solution of F (σ ) = 0 in this interval for K 2 ≥ 1/24. The uniqueness of the solution.

 3a2 σ (−1 + a2 σ 2 + 12K 2 ) = 12Ka2 σ + 2a2 σ 6 + 3a2 σ 2 + 36K 2 + √ . dσ 6 + 3a2 σ 2 + 36K 2 dF

The factor (−1 + a2 σ 2 + 12K 2 ) ≥ 0 which follows from (4.14) with the φopt . This implies that F (σ ) is an increasing function in the specified intervals.  Remark 7. Since any value of σ outside the interval [σl2 , σr ] specified in Theorem 4, in particular if the lower bound is violated, may produce oscillations in the computed solution. In this case, we recommend that if K 2 < 1/24 then, rather than using the σ that is optimal – in the sense of solving (5.7) – one should instead take some value in the range [σl2 , σr ] so as to satisfy all the desired properties for the scheme. For example, one may take σ = σl2 , which should minimize the term   F (σ ), or take σ = σl2 + σr /2, which will ensure that the scheme is von Neumann stable and produce computed solution which are free from spurious oscillations. In practise, we recommend taking the value σ = σl2 + σr /2 when K 2 < 1/24. This is because the bound associated with σl2 is quite sharp; violating it even by a small amount can lead to oscillations in the computed solution. We also point out that the discussion of Section 4.3 shows the upper bound associated with σr is not sharp. However, where the computed σopt is outside the region [σl2 , σr ], it is because σopt < σl2 , and never that σopt > σr . Therefore, taking the simpler bound (4.18) instead of (4.9) is completely reasonable.





Remark 8. When K (= ε/(ah)) ≫ 0, for example as in the diffusion-dominated case, some care may be required in the evaluation of expressions for σl1 and σr due to sensitivity of the expression to the numerical error associated with machine precision. In that case, they should be approximated using, for example, a standard Maclaurin series, giving expressions

σl1 ≈

3h



2ε

1 9

−

1

a2 h 2

22 92

ε

2

+

1

a4 h4

23 93

ε

4



,

σr ≈

1h 2ε

 1−

1 a2 h2 6 ε

2

−

31 a4 h4 54 ε 4

.

Remark 9. In the case of a pure diffusion problem, i.e., when a = 0 in (1.1), the optimal time-stepping problem is comparatively simple. As described in Remark 4, choosing θopt eliminates the fourth-order term in (4.1). Also, the fifth-order term in (4.1) vanishes when a = 0. Therefore, we can choose σ to eliminate the sixth-order term, giving

−

1 180

+

s2 3

− θs



1 6

 + s = 0.

68

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

Recall that the optimal value of √ θ denoted by θdiff is given by θdiff = 1/2 − 1/(12s), see (4.5). Solving the above equation with this value of θ gives s = 1/ 20. This yields the optimal value of σ , denoted by σdiff as 1

σdiff = √

h

2 5ε

.

(5.9)

Again we note that this corresponds with the value given in [13, p188]. 6. Numerical experiments In this section we examine the implementation of the proposed scheme, and its accuracy for a selection of test problems. In Section 6.1 we gather the theoretical results of the above sections to summarize the proposed scheme. Our main test problem is introduced in Section 6.2, and we report on extensive investigations of the accuracy and convergence of the scheme for various values of the discretization parameters, and the diffusion parameter, ε . This is followed in Section 6.2.2 with a study of the dependence of the optimal weights on the adjective term, a. We conclude with a comparison of various schemes that fall into the framework of the general two-weight scheme (2.1). In Section 6.3, we compare the proposed scheme with several high-order methods from the literature. These use a fivepoint stencil, and so are somewhat more expensive to implement than the proposed one. However, they share certain commonalities with our scheme: they are all linear, and they may be analyzed in the setting of a modified equivalent PDE. In this section we also include a comparison of the efficiency of the schemes for advection- and diffusion-dominated problems. It was noted in Section 4.1 that the MEPDE approach is valid only if the problem data are sufficiently differentiable. In Section 6.4 we study problems where the initial data lack smoothness, including the transport of a rectangular pulse. We observe that the accuracy of the method is compromised if the initial condition, Φ (x, 0), is not in at least C 2 (0, l). However, stability of the solutions is always maintained. In particular, no spurious oscillations are observed for any of the test cases. 6.1. The optimal time-stepping algorithm Given problem data for the model problem (1.1), i.e., values for ε and a, and for the boundary and initial data, select the desired mesh width h. Step 1: Calculate K := ε/ah, and

 1 (−3K + 1 + 9K 2 ), σl2 := 1 − 12K 2 , a a   1 σr := −1 − 6K 2 + 2 1 + 9K 4 + 3K 2 . 1

σl1 :=

a

Step 2: Calculate τopt

• If K 2 ≥ 1/16, solve     F (σ ) ≡ 6K a2 σ 2 − 12K 2 + −1 + a2 σ 2 + 12K 2 6 + 3a2 σ 2 + 36K 2 = 0, for σopt in [σl1 , σr ]; • else if 1/24 ≤ K 2 ≤ 1/16, solve F (σ ) = 0 for σopt , in [σl2 , σr ]; • else (i.e., K 2 < 1/24), take σopt = (1/2)(σl2 + σr ). Calculate τopt = σopt h. Step 3: Calculate θopt and φopt Calculate the optimal values of the weights θ and φ denoted by θopt and φopt , respectively, and given as follows:

θopt =

3a2 σopt h + 6ε −



2 6a2 h2 + 3a4 σopt h2 + 36ε 2

3ha

,

φopt = aσopt (1 − 2θopt ).

Step 4: Compute solution using the difference scheme Use τopt , θopt , φopt in the difference scheme given as follows:

  δt Φj + LNφ θopt Φjn+1 + (1 − θopt )Φj = 0,

where LNφ := −εδxx + aδx , δx = φopt D− + (1 − φopt )D0 ,

and the operators D0 , D− , δxx and δt are as defined in Section 2.

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

69

Fig. 3. Plots of the solution (6.1) to the test problem with ε = 10−3 , ε = 10−4 , and ε = 10−5 . Table 1 Table of errors for various N, demonstrating that for σ = σopt , the scheme is highly accurate. N

M

σ

θopt

φopt

EN

125 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000

32 64 130 206 301 486 646 832 1044 1283 1550 1844 2166

3.991 3.962 3.847 3.644 3.329 2.573⋆ 2.324⋆ 2.106⋆ 1.917⋆ 1.754⋆ 1.614⋆ 1.492⋆ 1.385⋆

0.024 0.045 0.080 0.107 0.123 0.090 0.098 0.103 0.108 0.111 0.114 0.116 0.117

0.950 0.902 0.807 0.716 0.628 0.527 0.467 0.418 0.376 0.341 0.312 0.287 0.265

1.98e−02 5.48e−03 1.32e−03 5.11e−04 2.11e−04 3.22e−06 1.39e−06 6.57e−07 3.35e−07 1.81e−07 1.03e−07 6.14e−08 3.80e−08

6.2. Test problem 1 We consider an example of a model predicting the transport of a Gaussian pulse with unit amplitude centered at x = a, given by

−(x − a)2 Φ (x, 0) = exp 4ε 



,

0 ≤ x ≤ 2.

This problem is widely used for comparison of different numerical schemes for the advection–diffusion problem; see, e.g., [9,23]. The exact solution to this problem is

Φ (x, t ) = √

1

1+t

 exp

−(x − (1 + t )a)2 4ε(1 + t )



,

0 ≤ x ≤ 2, t ≥ 0.

(6.1)

The boundary conditions are taken from the exact solution. In Fig. 3 we plot Φ (x, 2) for ε = 10−3 , ε = 10−4 , and ε = 10−5 . 6.2.1. Accuracy and convergence of the optimal scheme We wish to show that the proposed scheme is convergent, and to estimate the observed rate of convergence. To do that we have fixed a = 1/4 and ε = 10−4 and solved (1.1) for various values of N. In Table 1 we give these results. Here M is the number of time steps computed √ and σ is the value returned by the algorithm described above. If this corresponds to the optimal value of σ (i.e., K ≥ 1/ 24) we append the ⋆-symbol to its value. We also give the values of θopt and φopt . The last column contains error, reported here is the maximum error at any mesh point at the final time step:

EN = max |Φ (xj , tM ) − ΦjM |. j=0,...,N

Recall that when φ = 1 the scheme corresponds to the stable, but low-order, upwinding scheme, whereas φ = 0 gives the central differencing scheme that is suitable for a diffusion-dominated problem. We see in Table 1 that indeed when N is small, we get a value of φopt ≈ 1. By contrast, for larger N, φopt tends to zero.

70

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

(a) Convergence with respect to N.

(b) Convergence with respect to M.

Fig. 4. A log–log plot for the convergence of scheme with optimal weights and optimal time-stepping for ε = 10−4 . For smaller N and M, when σ ̸= σopt , the error appears to be O (N −2 + M −2 ). When σ = σopt , it appears to be O (N −5 + M −3 ). Table 2 Table of errors for various N with ε = 10−3 (left) and ε = 10−5 (right). N 125 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000

M 50 156 576 1 275 2 253 3 511 5 048 6 865 8 961 11 337 13 992 16 927 20 142

σ

2.573⋆ 1.614⋆ 0.870⋆ 0.589⋆ 0.444⋆ 0.356⋆ 0.297⋆ 0.255⋆ 0.223⋆ 0.198⋆ 0.179⋆ 0.162⋆ 0.149⋆

θopt

φopt

EN

0.090 0.114 0.123 0.126 0.126 0.127 0.127 0.127 0.127 0.127 0.127 0.127 0.127

0.527 0.312 0.164 0.110 0.083 0.066 0.055 0.048 0.042 0.037 0.033 0.030 0.028

9.95e−05 3.29e−06 6.55e−08 6.04e−09 1.09e−09 2.89e−10 9.76e−11 3.88e−11 1.83e−11 9.14e−12 5.33e−12 2.83e−12 1.78e−12

N 125 250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000

M 33 64 127 189 252 315 378 441 505 568 632 697 762

σ

θopt

φopt

EN

4.000 4.000 3.998 3.997 3.994 3.991 3.986 3.982 3.976 3.969 3.962 3.954 3.946

0.002 0.005 0.010 0.015 0.019 0.024 0.028 0.032 0.037 0.041 0.045 0.049 0.053

0.995 0.990 0.980 0.970 0.960 0.950 0.941 0.931 0.921 0.912 0.902 0.892 0.883

1.54e−01 6.04e−02 1.21e−02 5.82e−03 3.32e−03 2.12e−03 1.50e−03 1.11e−03 8.52e−04 6.75e−04 5.46e−04 4.52e−04 3.80e−04

Most importantly, though, we see that in all cases the error decreases as N and M increase. For smaller N and M, when

σ ̸= σopt , the scheme appears to be second-order accurate in both space and time. However, when σ = σopt , the scheme

appears to be at least fifth-order accurate in space and third-order accurate in time; see Fig. 4. Next, in Table 2 we give results corresponding to those in Table 1 but with ε = 10−3 (left) and ε = 10−5 (right). We notice that in the former case, all N lead to an optimal value of σ . Furthermore, the errors are extremely small, though with very large values of M. By contrast, in the right of Table 2 when ε = 10−5 none of the values of σ are ‘‘optimal’’. However, the scheme is very efficient, requiring relatively few time steps in order to ensure a stable computed solution. 6.2.2. Dependence of the scheme on the advection coefficient As mentioned earlier one of our main goals has been to design a scheme that adjusts automatically as the problem changes from being advection-dominated problems to diffusion-dominated. To verify that the proposed scheme achieves this we have fixed N = 1000 and ε = 10−4 , and solved (1.1) for a set of values of a in the interval [0, 1/2]. Consider first Fig. 5 where we show the values of θopt and φopt as a function of a. As one would expect, the value of φopt varies from φopt = 0 to φopt = 1 as a varies from a = 0 (diffusion-dominated) to a = 1/2 (advection-dominated). The value of θopt also changes, though not as significantly. Notice in Fig. 5 that the value of θopt , in particular, changes suddenly near a = 1/4. This is because, for this test problem,

√

when a ≤ 1/4, the optimal θ is chosen to be the solution to (5.7), whereas for a > 1/4 we have that K < 1/ 24 and it is taken as described in Step 3 of the algorithm above. This is demonstrated further in Table 3 where we present as summary of our results, where we see that smaller values of a require smaller values of M in order to compute an accurate solution. 6.2.3. Comparison of the optimal scheme with other methods As noted in Section 2, for particular values of θ and φ in the two-weight scheme (2.1) one obtains a variety of standard numerical schemes. We now consider a representative selection of these: Scheme 1: (θ , φ) = (0, 0), which is the forward Euler method, with central differencing in space; Scheme 2: (θ , φ) = (0, aτ /h), which is an explicit Lax–Wendroff type scheme [8];

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

71

Fig. 5. Values of θopt and φopt as a function of a. Table 3 Table of algorithm parameters and computed errors for N = 1000, ε = 10−4 and for various a. a

M

σ

θopt

φopt

EN

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50

224 230 248 276 311 301 338 381 426 473 520

4.472⋆ 4.352⋆ 4.035⋆ 3.625⋆ 3.216⋆ 3.329 2.961 2.628 2.348 2.116 1.923

0.127 0.123 0.114 0.102 0.090 0.123 0.113 0.103 0.095 0.087 0.080

0.000 0.164 0.312 0.433 0.527 0.628 0.687 0.729 0.762 0.787 0.807

1.69e−09 2.62e−07 1.62e−06 3.88e−06 6.30e−06 2.11e−04 2.66e−04 2.96e−04 3.13e−04 3.24e−04 3.31e−04

Table 4 Comparison of various schemes for test problem with ε = 10−4 and with σ = τ /h = N /M = 1. Errors are underlined if the corresponding solution exhibits spurious oscillations. N

500 1000 1500 2000 2500 3000

M

500 1000 1500 2000 2500 3000

Scheme 1

2

3

4

5

4.93e−01 1.78e−01 1.02e−01 7.16e−02 5.51e−02 4.48e−02

5.72e−02 1.17e−02 3.94e−03 1.53e−03 5.43e−04 7.55e−05

3.22e−01 2.49e−01 2.04e−01 1.73e−01 1.51e−01 1.34e−01

7.63e−02 1.92e−02 8.39e−03 4.70e−03 3.00e−03 2.08e−03

4.32e−03 4.70e−04 1.10e−04 3.55e−05 1.33e−05 5.21e−06

Scheme 3: (θ , φ) = (1, 1), which the backward Euler method with upwinding; Scheme 4: (θ , φ) = (1/2, 0), gives the Crank–Nicolson scheme with central differencing; Scheme 5: The two-weight scheme with the optimal weights (4.3). In the first instance, we consider these schemes where the value of τ is not chosen in the optimal fashion as described in Section 5. We take ε = 10−4 , a = 1/4, a range of values of N, and M = N (i.e., σ = 1). The results are shown in Table 4. In all cases, Scheme 5 yields the smallest error. In Fig. 6 we give a log–log plot for the errors, but for a larger range of N. This emphasizes that the optimal two-weight scheme (shown using triangles), is superior to the others, in terms of accuracy. Next best is Scheme 3, but that is not stable for larger N. Our main goal in this paper has not been just to design a highly accurate scheme, but also one that can generate solutions which are stable and free from spurious oscillations. Note that (6.1), the solution to the test problem, should be non-negative for all t (it represents the concentration of a solute in a flow), so we can detect spurious oscillations by checking for negative values in the computed solution. Where these are present, we have underlined the corresponding error value in Table 5. We see that for N = 500, only the low-order Backward Euler/Upwinding scheme is free from oscillations. In this case, the value of σ = τ /h = N /M is not within the bounds given in (5.8). For further comparisons between these schemes for different ε and N we refer the reader to [7, Section 5]. Next, we compare the same schemes but when σ is chosen in the optimal fashion. Again we have taken ε = 10−4 , a = 1/4 and a range of values of N. The results are shown in Table 5. In all cases, Scheme 5, i.e., the optimal one, yields

72

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

Fig. 6. A log–log plot of the errors for each of the five schemes taking σ = 1.

Fig. 7. A log–log of the errors for each of the five schemes taking σ = σopt .

Table 5 Comparison of various schemes for test problem with ε = 10−4 and with σ = σopt . Errors are underlined if the corresponding solution exhibits spurious oscillations. N

M

Scheme 1

2

3

4

5

500 1000 1500 2000 2500 3000

130 301 646 1044 1550 2166

5.55e+08 1.37e+03 3.87e−01 1.67e−01 9.79e−02 6.51e−02

3.93e−02 1.28e−02 3.09e−03 1.72e−03 1.08e−03 7.38e−04

3.66e−01 2.90e−01 2.30e−01 1.91e−01 1.62e−01 1.41e−01

1.05e−01 2.51e−02 9.54e−03 5.08e−03 3.14e−03 2.14e−03

1.32e−03 2.11e−04 1.39e−06 3.35e−07 1.03e−07 3.80e−08

a smaller error than the other methods, particularly for larger values of N. Again we have underlined the errors for which the corresponding solution exhibited detectable oscillations. Now we see that Scheme 5 is stable for all N. The only other scheme that has that property is the low-order Backward Euler/Upwinding scheme. A log–log plot of the errors is given in Fig. 7. Again we see the superiority of the optimal scheme. Finally, in Fig. 8 we plot solutions obtained using Schemes 3, 4 and 5 with N = 200, along with the true solution. Even for such a small value of N it is difficult to discern the difference between the true solution and that obtained using the optimal scheme. By contrast, significant damping of the solution is obvious with Scheme 4 (backward Euler, with upwinding). Scheme 4, which employs central differencing in space, exhibits large oscillations.

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

73

Fig. 8. Solutions computed by Schemes 3, 4, and 5 with N = 200 for 0.55 ≤ t ≤ 0.95.

6.3. Comparison with other linear high-order schemes In this section, we compare our scheme against four higher-order schemes available in the literature, considering both diffusion-dominated and advection dominated cases. The schemes are related to the proposed scheme in the sense that they may be analyzed in the framework of a modified equivalent PDE. The first method is explicit, and we have implemented the others using the Crank–Nicolson method for time integration. They all use a five-point stencil for the spatial discretization. Thus, these schemes are similar to the framework of the proposed weighted scheme, but with a wider stencil. The methods are given as follows: Scheme 6 (D4): The fourth-order explicit scheme given by Dehghan in [8, (41)–(43)] which uses the spatial discretizations

δx uj =

1



12

(12s + 2c 2 − 3c − 2)

− 4(c 2 + 6s − 4)

uj + 2 − uj 2h

uj + 1 − uj − 1 2h



+ (12s + 2c 2 + 3c − 2)

uj − uj − 2 2h

,

and

δxx uj =

1



4s

(−c 4 + 4c 2 − 12s2 − 12sc 2 + 8s)

+ (c − 4c + 12s + 12sc − 2s) 4

2

2

2

uj−1 − 2uj + uj+1 h2

uj−2 − 2uj + uj+2 4h2

 ,

where, as usual, s = ετ /h2 , and c = aτ /h. Scheme 7 (HV2): The second-order upwind scheme given by Hundsdorfer and Verwer in [24, (3.43)]

δx uj =

 1  uj−2 − 4uj−1 + 3uj .

2h

Scheme 8 (HV3): The third-order upwind scheme [24, (3.27)]

δx uj =

 1  uj−2 − 6uj−1 + 3uj + 2uj+1 .

6h

Scheme 9 (HV4): The fourth-order central scheme [24, (3.29)]

δx uj =

 1  uj−2 − 8uj−1 + 8uj+1 + uj+2 .

12h

With each of the advection discretizations in Schemes 7, 8 and 9, we use a fourth-order central diffusion scheme [24, (3.39)]

δxx uj =

1  12h

uj−2 − 16uj−1 + 30uj − 16uj+1 + uj+2 .



74

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

Fig. 9. Comparison of several high-order schemes applied to test problem (6.2) with ε = 1/10. Table 6 Comparison of various high-order schemes for test problem (6.2) with ε = 10−4 . N

M

Scheme 5

6 (D4)

7 (HV2)

8 (HV3)

9 (HV4)

600 1200 1800 2400 3000 3600 4200

318 914 1743 2879 4331 6104 8198

7.38e−04 6.06e−07 8.97e−08 2.02e−08 6.00e−09 2.10e−09 8.20e−10

7.75e−06 2.74e−07 8.24e−08 2.99e−08 1.29e−08 6.43e−09 3.49e−09

3.73e−01 1.08e−01 5.03e−02 2.91e−02 1.89e−02 1.32e−02 9.78e−03

9.98e−02 1.23e−02 3.41e−03 1.24e−03 5.43e−04 2.70e−04 1.48e−04

1.06e−01 1.31e−02 3.63e−03 1.34e−03 5.91e−04 2.98e−04 1.66e−04

We apply these methods to a variant on [23, Problem 3]. The exact solution is

Φ (x, t ) = A exp(Bt − λx),

where λ = (−a +



a2 − 4ε B)/(2ε),

(6.2)

from which one can obtain the initial and boundary conditions. The domain is [0, 1] × [0, 1]. First, to compare the methods for a diffusion-dominated problem, we take ε = 1/10, a = 1/2, A = 1/2000 and B = 15. The results are summarized in Fig. 9 below. The largest value of N taken is 200, for which the error for our proposed scheme is already O (10−9 ). Taking larger N results in errors so small that the effects of finite precision begin to dominate. The timestep is selected according to the algorithm of Section 6.1. These experiments show that proposed algorithm yields more accurate solutions than any of others. Each of Schemes 6–9 achieve their theoretical order of convergence in N, with the fourth-order scheme of [8] being slightly superior to the fourth-order scheme from [24]. We mentioned that our proposed scheme appears, in fact, to be converging with order O (N −6 ), although we do not have a theoretical explanation for this phenomenon at this time. It should be noted from Fig. 9 that, for this diffusion-dominated case, the algorithm of Section 6.1 yields M ≈ N 2 /2. Therefore, although our method is successful for these problems, it is expensive. A comparison of Schemes 5–9 when ε = 10−4 is given in Fig. 10, and also in Table 6. First note that Schemes 7 and 8 give essentially to same order of accuracy. This is likely due to the fact that Scheme 8 is not upwind-biased, and so more suited to diffusion dominated problems. For small N, when the optimal time-step is not available, the scheme of [8] is clearly best. However, for larger N our proposed scheme yields the most accurate solution and has accuracy that is O (N −5 + M −3 ). Furthermore, M ≈ ε N 2 , making the time step selected according to our algorithm quite efficient. 6.4. Problems with non-smooth initial conditions It is noted in Section 4.1 that the MEPDE approach is valid only if the problem data are sufficiently differentiable. We now investigate how the smoothness of the initial data influence the accuracy and stability of the computed solution. We are primarily interested in the case where the initial condition is discontinuous. However, to allow for a more comprehensive study, we consider (1.1) with the initial condition

Φ (x, 0) =



(1 − (5x − 3)2 )q 0

2/5 ≤ x ≤ 4/5 otherwise,

for q = 0, 1, 2, 3,

(6.3)

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

75

Fig. 10. Comparison of several high-order schemes applied to test problem (6.2) with ε = 10−4 .

Fig. 11. The initial condition (6.3) for q = 0 (left) and q = 2 (right), with corresponding solutions at t = 2.

and homogeneous boundary conditions, on the domain [0, l] × [0, T ] = [0, 2]2 . If q = 0, then Φ (x, 0) is discontinuous: the problem corresponds to the advection of a rectangular pulse. Otherwise, Φ (x, 0) ∈ C q−1 (0, 2). In Fig. 11 we show the initial data, and solutions at time t = 2 (generated using Chebfun, which should be exact to machine precision [25]) for the cases q = 0 (left) and q = 2 (right), with ε = 10−4 and a = 1/4. If Φ is not continuous, then the analysis of Section 4.1, which underpins the two-weight method (2.1), does not apply. As a result, high-order convergence cannot be obtained. Moreover, it is well known that a linear scheme, such as the one we propose, cannot be monotone across discontinuities and yet have better than first order convergence. This is discussed at length in, e.g., [26, Section 9.2] and [24, Section 1.7] for the advection problem. To circumvent this, one needs nonlinear schemes, see, e.g., [27]. However, studying this challenging problem allows us to verify that solution produced is indeed free from oscillations and other undesirable artifacts. To that end, in Fig. 12(a) we show the computed solutions to our test problem, with q ∈ {0, 1, 2, 3}. This demonstrates that, even if the problem data are not smooth, solutions produced by (2.1) do not exhibit oscillations. As expected, however, the rate of convergence of the method is adversely affected if Φ (x, 0) is not smooth. As shown in (12(b)), the method is only first-order when Φ (x, 0) is discontinuous. When q = 1, so Φ (x, 0) is in C 1 (0, l), but not in C 2 (0, l), the scheme is second order for all N. However, when Φ (x, 0) ∈ C 2 (0, l), we obtain fifth-order convergence, for sufficiently large N. See also Table 7. 7. Conclusions In the present study we have proposed a linear two-weight, three-point finite difference scheme for a model time dependent one-dimensional advection–diffusion problem. Our scheme involves two weights namely φ and θ ; φ is involved in a weighted average between central differencing and upwinding in the spatial discretization, and θ is involved in the time integration. We demonstrate that there is a subtle interplay between these weights, and this can be exploited in to give an optimal time step selection strategy that increases the order of convergence of the scheme.

76

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

(a) Computed solutions with N = 125.

(b) Convergence of the scheme.

Fig. 12. Φ (x, 0) as given in (6.3), with q = 0, 1, 2, 3. Table 7 Convergence of the two-weight scheme with optimal time-stepping when Φ (x, 0) is as given in (6.3), with q = 0, 1, 2, 3. N

Smoothness q=0

q=1

q=2

q=3

250 500 750 1000 1250⋆ 1500⋆ 1750⋆ 2000⋆ 2250⋆ 2500⋆

8.07e−02 4.02e−02 2.66e−02 1.99e−02 1.60e−02 1.33e−02 1.14e−02 9.97e−03 8.86e−03 7.98e−03

1.63e−03 4.06e−04 1.72e−04 8.83e−05 4.26e−05 2.95e−05 2.17e−05 1.66e−05 1.31e−05 1.06e−05

1.64e−04 4.06e−05 1.57e−05 6.46e−06 8.48e−08 3.67e−08 1.74e−08 8.89e−09 4.83e−09 2.76e−09

9.21e−05 2.22e−05 8.61e−06 3.55e−06 3.97e−08 1.72e−08 8.15e−09 4.17e−09 2.26e−09 1.29e−09

The scheme is analyzed for stability and other desirable properties in order to obtain certain bounds on the weights involved. The optimal values of the weights are obtained using the notion of a modified equivalent partial differential equation. These optimal values of the weights are used, in conjunction with the bounds obtained on them, to derive sharp bounds for the time step which guarantees that the finite difference scheme is stable and the resulting computed solution is free from spurious oscillations. Furthermore, a root finding problem in one variable is formulated in order to obtain an optimal time-stepping. The constituents of this root finding problem are the coefficient of the fourth order term in the equivalent differential equation, and the bounds obtained due to stability and other desirable properties. All the theoretical results are gathered in the form of an algorithm, which has been applied to various test problems. For comparative study we have considered various high-order schemes available in the literature. The numerical results validate the theory developed here. References [1] L. Debnath, Nonlinear Partial Differential Equations for Scientists and Engineers, Birkhäuser, Boston, MA, 1997. [2] K.W. Morton, Numerical Solution of Convection–Diffusion Problems, in: Applied Mathematics and Mathematical Computation, vol. 12, Chapman and Hall, London, 1996. [3] Clive A.J. Fletcher, Computational Techniques for Fluid Dynamics 1. Fundamental and General Techniques, second ed., in: Springer Series in Computational Physics, Springer-Verlag, Berlin, 1991. [4] C. Hirsch, Numerical Computation of Internal and External Flows. Volume 2: Computational Methods for Inviscid and Viscous Flows, in: Wiley Series in Numerical Methods in Engineering, Wiley, 1990. [5] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer, Berlin, 1994. [6] Hans-Gorg Roos, Martin Stynes, Lutz Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008. [7] N.M. Chadha, N. Madden, A two-weight scheme for a time-dependent advection–diffusion problem, in: C. Clavero, J.L. Gracia, F.J. Lisbona (Eds.), BAIL 2010—Boundary and Interior Layers, Computational and Asymptotic Methods, in: Lect. Notes Comput. Sci. Eng., vol. 81, Springer, 2011, pp. 99–108. [8] M. Dehghan, Weighted finite difference techniques for the one-dimensional advection–diffusion equation, Appl. Math. Comput. 147 (2) (2004) 307–319. [9] M. Morandi Cecchi, M.A. Pirozzi, High order finite difference numerical methods for time-dependent convection-dominated problems, Appl. Numer. Math. 55 (3) (2005) 334–356. [10] A. Rigal, High order difference schemes for unsteady one-dimensional diffusion–convection problems, J. Comput. Phys. 114 (1) (1994) 59–76.

N.M. Chadha, N. Madden / Journal of Computational and Applied Mathematics 294 (2016) 57–77

77

[11] R. Manohar, S.R.K. Iyengar, U.A. Krishnaiah, High order difference methods for linear variable coefficient parabolic equation, J. Comput. Phys. 77 (2) (1988) 513–523. [12] B.J. Noye, A new third-order finite-difference method for transient one-dimensional advection–diffusion, Commun. Appl. Numer. Methods 6 (4) (1990) 279–288. [13] R.D. Richtmyer, K.W. Morton, Difference Methods for Initial-Value Problems, second ed., Wiley, New York, 1967. [14] A.V. Goolin, The stability boundary of certain two-layer and three-layer difference schemes, in: L.G. Vulkov (Ed.), Numerical analysis and its applications, in: Lect. Notes Comput. Sci., vol. 1998, Springer, 2001, pp. 341–349. [15] A.A. Samarskii, Regularisation of difference schemes, Comput. Math. Math. Phys. 7 (1) (1967) 79–120. [16] Chi-Wang Shu, High order weighted essentially nonoscillatory schemes for convection dominated problems, SIAM Rev. 51 (1) (2009) 82–126. [17] X.Y. Wu, Q. Wang, N.A. Adams, An adaptive central-upwind weighted essentially non-oscillatory scheme, J. Comput. Phys. 229 (2010) 8952–8965. [18] P. Wesseling, Principles of Computational Fluid Dynamics, in: Springer Series in Computational Mathematics, vol. 29, Springer, Berlin, 2000. [19] R.F. Warming, B.J. Hyett, The modified equation approach to the stability and accuracy analysis of finite-difference methods, J. Comput. Phys. 14 (1974) 159–179. [20] M. Dehghan, Quasi-implicit and two-level explicit finite-difference procedures for solving the one-dimensional advection equation, Appl. Math. Comput. 167 (2005) 46–67. [21] John Noye, Some explicit three-level finite-difference simulations of advection, Math. Comput. Simul. 32 (1990) 359–372. [22] J.W. Thomas, Numerical Partial Differential Equations: Finite Difference Methods, in: Texts in Applied Mathematics, vol. 22, Springer-Verlag, New York, 1995. [23] A. Mohebbi, M. Dehghan, High-order compact solution of the one-dimensional heat and advection–diffusion equations, Appl. Math. Model. 34 (10) (2010) 3071–3084. [24] W. Hundsdorfer, J. Verwer, Numerical Solution of Time-Dependent Advection–Diffusion–Reaction Equations, in: Springer Series in Computational Mathematics, vol. 33, Springer-Verlag, Berlin, 2003. [25] T.A. Driscoll, N. Hale, L.N. Trefethen (Eds.), Chebfun Guide, Pafnuty Publications, Oxford, 2014. [26] P. Wesseling, von Neumann stability conditions for the convection–diffusion equation, IMA J. Numer. Anal. 16 (4) (1996) 583–598. [27] Bram van Leer, Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme, J. Comput. Phys. 14 (4) (1974) 361–370.