Skip to main content

Theory and Modern Applications

A fourth-order accurate quasi-variable mesh compact finite-difference scheme for two-space dimensional convection-diffusion problems

Abstract

We discuss a new nine-point fourth-order and five-point second-order accurate finite-difference scheme for the numerical solution of two-space dimensional convection-diffusion problems. The compact operators are defined on a quasi-variable mesh network with the same order and accuracy as obtained by the central difference and averaging operators on uniform meshes. Subsequently, a high-order difference scheme is developed to get the numerical accuracy of order four on quasi-variable meshes as well as on uniform meshes. The error analysis of the fourth-order compact scheme is described in detail by means of matrix analysis. Some examples related with convection-diffusion equations are provided to present performance and robustness of the proposed scheme.

1 Introduction

The two-dimensional elliptic equations

$$\begin{aligned} -{\varepsilon} {\nabla}^{{2}} {U}+{a} ( {x},{y} ) {U}^{{x}} +{b} ( {x},{y} ) {U}^{{y}} +{c} ( {x},{y} ) {U}+{d} ( {x},{y} ) ={0},\quad ({x},{y})\in \Omega, \end{aligned}$$
(1.1)

will be considered to develop numerical algorithms for computing the concentration \(U(x,y)\) of mass transfer. Here, ε is viscosity coefficient or diffusion coefficient (\(0<\varepsilon\ll1\)) and \((a(x,y),b(x,y))>( \beta_{1}, \beta_{2} )>(0,0)\) on Ω̅ (closure of \(\Omega = ( 0,1 ) \times ( 0,1 )\)) is velocity vector, \(\beta_{1}, \beta_{2}\) are finite constants. We also assume that \(a ( x,y ),b ( x,y ), c ( x,y ),d(x,y)\) are continuous and \(c ( x,y ) \geq0\) on Ω̅ in order to ensure the existence of a solution. Let the following smooth boundary data be given:

$$\begin{aligned} &U ( 0,y ) = \varphi_{1} ( y ),\qquad U ( 1,y ) = \varphi_{2} ( y ),\quad 0\leq y\leq1, \end{aligned}$$
(1.2)
$$\begin{aligned} & U ( x,0 ) = \psi_{1} ( x ), \qquad U ( x,1 ) = \psi_{2} ( x ),\quad 0 \leq x\leq1. \end{aligned}$$
(1.3)

The second-order partial derivatives in the mathematical model (1.1) describe the diffusion process, while the first-order partial derivatives are associated with the convection phenomenon. When \(\varepsilon\rightarrow0\), convection dominates the diffusion process, and the solution values of (1.1)-(1.3) exhibit boundary layer behavior, that is, solution changes rapidly in a small region, while outside the small region solution behavior is smooth. In two dimensions, the boundary layer may occur at \(x=0\) and \(x=1\), known as normal layer, and/or at \(y=0\) and \(y=1\), known as parallel layer, while the one-dimensional convection-diffusion problems exhibit only normal boundary layer [1]. Such type of a differential equation is said to be singularly perturbed. The solution of singular perturbation problems approaches a discontinuous limit as a small positive quantity (ε), known as perturbation parameter, approaches zero. Thus, the analysis and numerical solution of singular perturbation problems are significant.

The convection-diffusion problems occur in the area of fluid dynamics and several branches of applied mathematics. In the convection-diffusion process, transport phenomena prevail diffusion, whose effects are restricted to a small part of the domain and solution values exhibit multiple characters for small values of diffusion parameter. Thus, a second-order discretization of the Laplace operator may not ensure the consistency and stability of the numerical scheme [2]. The numerical schemes developed by the application of an upwind scheme and a central difference operator result in an unstable solution and are deficient in computational order. Any attempt to obtain solution values up to the desired accuracies may land up in massive computing time despite a well-structured block-tridiagonal matrix of the difference equations received by the application of an upwind scheme and central differences. Thus, acquiring improved finite-difference methods for convection-diffusion problems has a significant impact on numerical approximations of ordinary and partial differential equations [3–6]. The compact scheme amongst the various finite-difference replacements of convection-diffusion problems (1.1)-(1.3) has received more attention due to a minimal width of stencils in the x- and y-coordinate directions and easy computations. In contrast, high-order difference schemes formulated with non-compact stencils yield a higher bandwidth of the iteration matrix, and this involves large arithmetic operations. Numerical simulations with high-order compact difference schemes depict more accurate solution values on quasi-variable meshes as compared to some high-order compact scheme on a uniform mesh network. This happens because a truncation error in a finite-difference approximation depends upon the derivative of the variable as well as mesh spacing. Therefore, to attain uniformly distributed truncation errors, it is essential to employ non-uniform meshes, i.e., finer meshes in the region for largely deviated derivatives and coarse meshes for a smooth function. In this manner, the error disperses almost uniformly over the domain of integration and renders an accurate solution to a greater extent. Thus, high-order finite-difference discretization formulated on a quasi-variable mesh network leads to more precise numerical solutions and brings unconditional numerical stability [7].

Mishra and Sanyasiraju [8] presented an exponential compact scheme of high-order accuracy for convection-dominated equations. Mohanty and Setia [9] described fourth-order discretization of two-space dimensional elliptic equations by means of off-step uniformly spaced meshes. Some fourth-order finite-difference discretizations of two-space dimensional linear and non-linear elliptic problems can be found in [10–13] and references therein. Ghaffar [14] obtained a high accuracy compact scheme using non-uniform meshes for the Helmholtz equation and solved the difference equations with the help of the multigrid method. Jha [15, 16] developed a third-order exponential expanding mesh compact scheme for mildly non-linear elliptic equations. The Galerkin and Petrov-Galerkin finite element method for determining the approximate solution values to the two-dimensional convection-diffusion problems were discussed by Hegarty [17]. A nine-point tailored finite point method for solving convection-diffusion-reaction equation was developed by Shih et al. [18]. In the context of one dimension, the finite-difference replacement of a convection-diffusion equation was extensively discussed in [19, 20], and non-uniform mesh compact finite-difference operators for first- and second-order ordinary derivatives associated with the Numerov fourth-order method were obtained in [21, 22]. Finite-difference methods for convection-diffusion problems showing exponential or parabolic boundary layer behavior were described in detail by Roos et al. [23].

The work presented in this article is organized in the following manner. In Section 2, we describe two-dimensional quasi-variable meshes to deal with parallel and normal layer by means of mesh parameters in the x- and y-directions. Section 3 discusses compact operators and high-order approximations of first- and second-order derivatives on minimum stencils. A new compact scheme of fourth-order accuracy using quasi-variable meshes has been obtained. The suggested scheme is analyzed for the convergence and the bounds of the discretization error are obtained in Section 4. Numerical simulations with some convection-diffusion problems that exhibit normal and/or parallel layers are carried out in Section 5. The paper is concluded at last with remarks and further scope.

2 Quasi-variable meshes

Let L and M be positive integers, and divide the domain \([0,1]\times[0,1]\) into \((L+1)(M+1)\) cells with the coordinates (\(x_{l}, y_{m}\)), where \(x_{0} =0\), \(x_{l} = x_{l-1} + h_{l}\), \(l=1 ( 1 ) L\), \(x_{L+1} =1\) and \(y_{0} =0\), \(y_{m} = y_{m-1} + k_{m}\), \(m=1 ( 1 ) M\), \(y_{M+1} =1\). The mesh step-size is determined by the stretching functions \(h_{l+1} = h_{l} ( 1+\alpha h_{l} )\), \(l=1 ( 1 ) L\) and \(k_{m+1} = k_{m} (1+\beta k_{m} )\), \(m=1 ( 1 ) M\), by suitably chosen normalization \(\tilde{\alpha} =\alpha h_{1}\) and \(\tilde{\beta} =\beta k_{1}\). Since the length of a diffusion space along the x-direction is one, for a given value of α̃, the relation \(\sum_{l=1}^{L+1} h_{l} =1\) easily produces the first mesh step-size \(h_{1}\) in the x-direction (and similarly in the y-direction). As an example, \(h_{1} =1/(1+2 \tilde{\alpha} )\) if \(L=1\). In particular, if \(\alpha=\beta=0\), the meshes are uniformly distributed and \(h=h_{l}, k= k_{m}, \forall l,m\). The mesh step-size is increasing if and only if \(h_{l} < h_{l+1} \forall l=1(1)L \Leftrightarrow h_{l} < h_{l} ( 1+\alpha h_{l} ) \Leftrightarrow \alpha>0\). In a similar manner, we can prove that the mesh step-size is decreasing when \(\alpha<0\). Therefore, the mesh-step sequences \(\{ h_{l} \}_{l=1}^{L+1}\) and, similarly, \(\{ k_{m} \}_{m=1}^{M+1}\) are monotonic.

Let us consider the uniform mesh partition of the domain \(\mathcal{P=} [ 0,1 ] = \{ p_{l} =lh,l=0 ( 1 ) L+1 \}\), \(h= 1/(L+1)\). Since the mesh-step sequence is monotonic, it is possible to define a one-one onto map \(\psi: \mathcal{P\longrightarrow P}\) such that \(\psi ( p_{l} ) = x_{l}, l=0 ( 1 ) L+1\), and the Jacobian \(J ( p ) =d\psi(p)/dp\) is bounded above and below by some positive constants as \(0< n\leq J ( p ) \leq N<\infty, \forall p\in \mathcal{P}\).

Therefore,

$$\begin{aligned} J ( p ) >0 \quad&\Rightarrow\quad \frac{d\psi(p)}{dp} >0 \quad\Rightarrow \quad \frac{\psi ( p_{l+1} ) -\psi ( p_{l} )}{p_{l+1} - p_{l}} >0 \quad\Rightarrow\quad \frac{x_{l+1} - x_{l}}{ ( l+1 ) h-lh} >0\\ &\Rightarrow \quad h_{l+1} >0\quad \forall l. \end{aligned}$$

Also, \(J ( p ) \leq N\Rightarrow \frac{d\psi ( p )}{ dp} \leq N\Rightarrow \frac{\psi ( p_{l+1} ) -\psi ( p_{l} )}{p_{l+1} - p_{l}} \leq N \Rightarrow \frac{x_{l+1} - x_{l}}{ ( l+1 ) h-lh} \leq N\).

That is,

$$ 0< h_{l+1} \leq Nh,\quad \forall l\quad\Rightarrow\quad \max _{l} \vert h_{l+1} \vert \leq Nh \quad\Rightarrow \quad \Vert \boldsymbol {h} \Vert _{\infty} \leq Nh= \frac{N}{L+1} \leq \frac{N}{ L}, $$
(2.1)

where \(\boldsymbol {h} = [ h_{1}, h_{2},\dots, h_{L+1} ]^{t}\).

This implies

$$ \Vert \boldsymbol {h} \Vert _{\infty} =O ( h ) =O \biggl( \frac{1}{L} \biggr)\quad \mbox{and hence}, \quad \Vert \boldsymbol {h} \Vert _{\infty} \longrightarrow0 \quad\mbox{as } L\longrightarrow\infty. $$
(2.2)

Thus, the maximum step-size along the x-direction diminishes with the growth of mesh points. Likewise, the maximum step-size along the y-direction approaches to zero if M is very large. Thus, we observed that the mesh step-size is inversely proportional to the number of mesh points.

Sundqvist and Veronis [24] initially discussed such a mesh network in the context of wind-driven ocean circulation, and later the application to digital electrochemistry was described by Britz [25]. Some compact operators related to a quasi-variable mesh were mentioned in the literature [26, 27].

3 Finite-difference schemes and compact operators

The operators used to obtain first- and second-order partial derivatives with a minimum stencil width are known as compact. The high-order finite-difference replacement of equation (1.1) requires discretization of partial derivatives, and thus we consider the following approximations:

$$\begin{aligned} &\overline{U}_{l,m+\delta}^{xx} =2 h_{l}^{-2} \bigl[ ( 1+\alpha h_{l} )^{-1} \bigl( ( 2+\alpha h_{l} )^{-1} U_{l+1,m+\delta} - U_{l,m+\delta} \bigr) + ( 2+\alpha h_{l} )^{-1} U_{l-1,m+\delta} \bigr], \end{aligned}$$
(3.1)
$$\begin{aligned} & \overline{U}_{l+\delta,m}^{yy} =2 k_{m}^{-2} \bigl[ ( 1+\beta k_{m} )^{-1} \bigl( ( 2+\beta k_{m} )^{-1} U_{l+\delta,m+1} - U_{l+\delta,m} \bigr) + ( 2+\beta k_{m} )^{-1} U_{l+\delta,m-1} \bigr], \end{aligned}$$
(3.2)
$$\begin{aligned} & \left [ \begin{matrix}{} \overline{U}_{l-1,m+\delta}^{x}\\ \overline{U}_{l,m+\delta}^{x}\\ \overline{U}_{l+1,m+\delta}^{x} \end{matrix} \right ] \mathcal{=M(} h_{l},\alpha) \left [ \begin{matrix}{} U_{l-1,m+\delta}\\ U_{l,m+\delta}\\ U_{l+1,m+\delta} \end{matrix} \right ], \qquad \left [ \begin{matrix}{} \overline{U}_{l+\delta,m-1}^{y}\\ \overline{U}_{l+\delta,m}^{y}\\ \overline{U}_{l+\delta,m+1}^{y} \end{matrix} \right ] \mathcal{=M(} k_{m}, \beta) \left [ \begin{matrix}{} U_{l+\delta,m-1}\\ U_{l+\delta,m}\\ U_{l+\delta,m+1} \end{matrix} \right ], \end{aligned}$$
(3.3)

where \(\delta=0,\pm1\), and for \(\tau\in \{ h_{l}, k_{m} \}\), \(\gamma \in \{ \alpha, \beta \}\),

$$\mathcal{M} ( \tau,\gamma ) = \frac{1}{\tau(2+3\gamma\tau)} \left [ \begin{matrix}{} -(3+4\gamma\tau) & 4(1+\gamma\tau) & -1\\ -(1+2\gamma\tau) & 2\gamma\tau & 1\\ 1+2\gamma\tau(1+\gamma\tau) & -4-2\gamma\tau(2+\gamma\tau) & 3+2\gamma\tau \end{matrix} \right ]. $$

Now, we define the following operators:

$$\begin{aligned} &\mathcal{P}_{x} U_{l,m} = h_{l} \overline{U}_{l,m}^{x},\qquad \mathcal{P}_{y} U_{l,m} = k_{m} \overline{U}_{l,m}^{y}, \end{aligned}$$
(3.4)
$$\begin{aligned} &\mathcal{Q}_{x} U_{l,m} = h_{l}^{2} \overline{U}_{l,m}^{xx},\qquad \mathcal{Q}_{y} U_{l,m} = k_{m}^{2} \overline{U}_{l,m}^{yy}. \end{aligned}$$
(3.5)

These operators are commutative and derived with a minimum number of stencils essential to discretize the highest-order partial derivatives present in the convection-dominated equation (1.1). In particular, if the meshes are uniformly spaced, that is, \(h= h_{l}, k= k_{m}, \forall l,m\), then \(\mathcal{P}_{x} =2 \mu_{x} \delta_{x}\) and \(\mathcal{Q}_{x} = \delta_{x}^{2}\), where \(\mu_{x}\) and \(\delta_{x}\) are averaging and central differencing operators in the x-direction [28].

By means of operators (3.4)-(3.5) it is easy to approximate partial- and mixed-order derivatives of an analytic function \(G ( x,y )\) at the mesh point (\(x_{l}, y_{m}\)) in the following manner:

$$\begin{aligned} &\overline{G}_{l,m}^{x} = h_{l}^{-1} \mathcal{P}_{x} G_{l,m} +\mathrm{O} \bigl( h_{l}^{2} \bigr),\qquad \overline{G}_{l,m}^{y} = k_{m}^{-1} \mathcal{P}_{y} G_{l,m} + \mathrm{O} \bigl( k_{m}^{2} \bigr), \end{aligned}$$
(3.6)
$$\begin{aligned} &\overline{G}_{l,m}^{xx} = h_{l}^{-2} \mathcal{Q}_{x} G_{l,m} +\mathrm{O} \bigl( h_{l}^{2} \bigr),\qquad \overline{G}_{l,m}^{yy} = k_{m}^{-2} \mathcal{Q}_{y} G_{l,m} + \mathrm{O} \bigl( k_{m}^{2} \bigr), \end{aligned}$$
(3.7)
$$\begin{aligned} &\overline{G}_{l,m}^{xy} = h_{l}^{-1} k_{m}^{-1} \mathcal{P}_{x} \mathcal{P}_{y} G_{l,m} +\mathrm{O} \bigl( h_{l}^{2} + k_{m}^{2} \bigr). \end{aligned}$$
(3.8)

An immediate application of these operators to equation (1.1) yields a five-point second-order accurate discretization scheme. Such kind of a second-order method on a variable mesh network is known as supra-convergent scheme [29].

Now, we describe a new fourth-order scheme for linear Poisson’s equation and then extend it to the elliptic equation (1.1), which involves convection terms \(U^{x} =\partial U/\partial x\) and \(U^{y} =\partial U/\partial y\).

By means of a linear combination of functional values \(G_{\rho, \sigma} =G( x_{\rho}, y_{\sigma} )\), \(( \rho,\sigma ) \in\mathcal{ D}=\{ l-1,l,l+1\}\times\{m-1,m,m+1\}\), the finite-difference replacement for a simple form of two-space dimensional elliptic equations (Poisson’s equation)

$$ -\varepsilon \nabla^{2} U+G ( x,y ) =0 $$
(3.9)

is given by

$$ - \nabla_{h_{l}, k_{m}}^{2} U_{l,m} + h_{l}^{2} k_{m}^{2} \mathcal{L} G_{l,m} = T_{l,m}, \quad l=1 ( 1 ) L, m=1 ( 1 ) M, $$
(3.10)

where

$$\begin{aligned} \nabla_{h_{l}, k_{m}}^{2} \equiv{}& \varepsilon \bigl( k_{m}^{2} \mathcal{Q}_{x} + h_{l}^{2} \mathcal{Q}_{y} \bigr) + \frac{\varepsilon}{3} \bigl[\alpha h_{l}^{3} \mathcal{P}_{x} \mathcal{Q}_{y} +\beta k_{m}^{3} \mathcal{P}_{y} \mathcal{Q}_{x} \bigr] \\ &{} + \frac{\varepsilon}{12} \bigl[ h_{l}^{2} ( 1+\alpha h_{l} ) + k_{m}^{2} ( 1+\beta k_{m} ) \bigr] \mathcal{Q}_{x} \mathcal{Q}_{y} \end{aligned}$$
(3.11)

is a discrete form of the Laplace operator \(\nabla^{2} = \partial_{x}^{2} + \partial_{y}^{2}\),

$$ \mathcal{L} \equiv 1+ \frac{h_{l}}{3} \alpha \mathcal{P}_{x} + \frac{k_{m}}{3} \beta \mathcal{P}_{y} + \frac{h_{l} k_{m}}{9} \alpha \beta \mathcal{P}_{x} \mathcal{P}_{y} + \frac{1}{12} ( 1+\alpha h_{l} ) \mathcal{Q}_{x} + \frac{1}{12} ( 1+ \beta k_{m} ) \mathcal{Q}_{y}, $$
(3.12)

and \(T_{l,m}\) is the local truncation error calculated as

$$ T_{l,m} = h_{l}^{2} k_{m}^{2} O \bigl( h_{l}^{4} + h_{l}^{2} k_{m}^{2} + k_{m}^{4} \bigr). $$
(3.13)

The eighth-order of local truncation error obtained here is irrespective of mesh parameters \(\alpha, \beta\) being chosen zero or non-zero. Since the operator \(\mathcal{L}\) in equation (3.10) is multiplied by \(h_{l}^{2} k_{m}^{2}\), this accomplishes an accuracy of order four employed with quasi-variable meshes (\(\alpha\neq0 \vee \beta\neq0\)) or uniform meshes (\(\alpha=\beta=0\)).

Our main aim is to extend the fourth-order method (3.10) to the convection-dominated equation (1.1) that comprises first-order partial derivatives in the x- and y-directions along with the function \(U(x,y)\).

Let us consider

$$\begin{aligned} &F_{\rho,\sigma} = a_{\rho,\sigma} \overline{U}_{\rho,\sigma}^{x} + b_{\rho,\sigma} \overline{U}_{\rho,\sigma}^{y} + c_{\rho,\sigma} U_{\rho,\sigma} + d_{\rho,\sigma},\quad (\rho,\sigma)\in \hat{\mathcal{D}} \mathcal{=D\sim \bigl\{ (} l,m) \bigr\} , \end{aligned}$$
(3.14)
$$\begin{aligned} &\overline{\overline{U}}_{l,m}^{x} = \overline{U}_{l,m}^{x} + h_{l} \bigl( \alpha_{1} F_{l+1,m} + \alpha_{2} F_{l-1,m} + \alpha_{3} \overline{U}_{l+1,m}^{yy} + \alpha_{4} \overline{U}_{l-1,m}^{yy} \bigr), \end{aligned}$$
(3.15)
$$\begin{aligned} &\overline{\overline{U}}_{l,m}^{y} = \overline{U}_{l,m}^{y} + k_{m} \bigl( \beta_{1} F_{l,m+1} + \beta_{2} F_{l,m-1} + \beta_{3} \overline{U}_{l,m+1}^{xx} + \beta_{4} \overline{U}_{l,m-1}^{xx} \bigr), \end{aligned}$$
(3.16)
$$\begin{aligned} &F_{l,m} = a_{l,m} \overline{\overline{U}}_{l,m}^{x} + b_{l,m} \overline{\overline{U}}_{l,m}^{y} + c_{l,m} U_{l,m} + d_{l,m}, \end{aligned}$$
(3.17)

where \(\alpha_{i}, \beta_{i}, i=1 ( 1 ) 4\) are unknown parameters to be measured in such a way that

$$ \mathcal{L} (F_{l,m} - G_{l,m} )= O \bigl( h_{l}^{4} + h_{l}^{2} k_{m}^{2} + k_{m}^{4} \bigr). $$
(3.18)

By making use of (3.10), (3.14)-(3.17) and (3.18), the algebraic calculations give us

$$\begin{aligned} &\alpha_{1} =-(1+\alpha h_{l} )/ \bigl[8\varepsilon(2+ \alpha h_{l} ) \bigr], \qquad \alpha_{2} =- \alpha_{1} -3 \alpha^{2} h_{l}^{2} / \bigl[4\varepsilon \bigl(2+ \alpha^{2} h_{l}^{2} + \beta^{2} k_{m}^{2} \bigr) \bigr], \\ &\beta_{1} =-(1+\beta k_{m} )/ \bigl[8\varepsilon(2+\beta k_{m} ) \bigr],\qquad \beta_{2} =- \beta_{1} -3 \beta^{2} k_{m}^{2} / \bigl[4\varepsilon \bigl(2+ \alpha^{2} h_{l}^{2} + \beta^{2} k_{m}^{2} \bigr) \bigr], \\ &\alpha_{3} =- \varepsilon\alpha_{1},\qquad \alpha_{4} =- \varepsilon\alpha_{2}, \qquad \beta_{3} =- \varepsilon \beta_{1},\qquad \beta_{4} =- \varepsilon\beta_{2}. \end{aligned}$$

As a result, we obtain a single compact discretization scheme

$$ - \nabla_{h_{l}, k_{m}}^{2} U_{l,m} + h_{l}^{2} k_{m}^{2} \mathcal{L} F_{l,m} = T_{l,m}, \quad l=1 ( 1 ) L, m=1 ( 1 ) M, $$
(3.19)

that numerically approximates equation (1.1) with fourth-order accuracy on a quasi-variable mesh network as well as on a uniform mesh network. The discretization (3.19) yields a non-symmetric matrix after incorporating the boundary values (1.2)-(1.3) as

$$\begin{aligned} &U_{0,m} = \varphi_{1} ( y_{m} ), \qquad U_{L+1,m} = \varphi_{2} ( y_{m} ),\quad m=0 ( 1 ) M+1, \end{aligned}$$
(3.20)
$$\begin{aligned} &U_{l,0} = \psi_{1} ( x_{l} ),\qquad U_{l,M+1} = \psi_{2} ( x_{l} ), \quad l=0 ( 1 ) L+1. \end{aligned}$$
(3.21)

The lexicographical ordering of the unknown values \(U_{l,m}\) in equation (3.19) gives rise to a block tri-diagonal linear system of equations and can be easily computed by means of the Gauss-Seidel iterative formula. For the programming, we must neglect the truncation error \(T_{l,m}\) from equation (3.19) and replace the exact value \(U_{l,m} =U( x_{l}, y_{m} )\) by its approximate value \(u_{l,m}\).

4 Convergence analysis and error bounds

In this section, we discuss the upper bounds of discretization errors and derive the necessary convergence conditions for scheme (3.19). The compact scheme (3.19) determines the exact solution values \(\boldsymbol {U} =[ U_{l,m} ]\), and in terms of the mesh-ratio parameter \(\zeta_{l,m} = k_{m} / h_{l}\), it can be expressed as

$$ - \overline{\nabla}_{h_{l}, k_{m}}^{2} U_{l,m} + h_{l}^{2} \zeta_{l,m}^{2} \mathcal{L} F_{l,m} = \overline{T}_{l,m},\quad l=1 ( 1 ) L, m=1 ( 1 ) M, $$
(4.1)

where \(\overline{T}_{l,m} =O( h_{l}^{6} )\) and

$$\begin{aligned} \overline{\nabla}_{h_{l}, k_{m}}^{2} ={}& h_{l}^{-2} \nabla_{h_{l}, k_{m}}^{2} \equiv \varepsilon \bigl( \zeta_{l,m}^{2} \mathcal{Q}_{x} + \mathcal{Q}_{y} \bigr) + \frac{\varepsilon}{3} h_{l} \bigl[ \alpha \mathcal{P}_{x} \mathcal{Q}_{y} +\beta \zeta_{l,m}^{3} \mathcal{P}_{y} \mathcal{Q}_{x} \bigr]\\ &{} + \frac{\varepsilon}{12} \bigl[1+\alpha h_{l} + \zeta_{l,m}^{2} (1+\beta h_{l} \zeta_{l,m} ) \bigr] \mathcal{Q}_{x} \mathcal{Q}_{y}. \end{aligned}$$

Our aim is to compute approximate solution vector \(\boldsymbol {u} =[ u_{l,m} ]\) that satisfies

$$ - \overline{\nabla}_{h_{l}, k_{m}}^{2} u_{l,m} + h_{l}^{2} \zeta_{l,m}^{2} \mathcal{L} f_{l,m} =0, \quad l=1 ( 1 ) L, m=1 ( 1 ) M, $$
(4.2)

where

$$\begin{aligned} &f_{p,q} = a_{\rho,\sigma} \overline{u}_{\rho,\sigma}^{x} + b_{\rho,\sigma} \overline{u}_{\rho,\sigma}^{y} + c_{\rho,\sigma} u_{\rho,\sigma} + d_{\rho,\sigma} \approx F_{\rho,\sigma},\quad (\rho,\sigma) \in \hat{\mathcal{D}}, \end{aligned}$$
(4.3)
$$\begin{aligned} &f_{l,m} = a_{l,m} \overline{\overline{u}}_{l,m}^{x} + b_{l,m} \overline{\overline{u}}_{l,m}^{y} + c_{l,m} u_{l,m} + d_{l,m} \approx F_{l,m} \end{aligned}$$
(4.4)

and the formula for \(\overline{u}_{\rho,\sigma}^{x}\), \(\overline{u}_{\rho,\sigma}^{y}, (\rho,\sigma)\in \hat{\mathcal{D}}\) is obtained from (3.3), and \(\overline{\overline{u}}_{l,m}^{x}\), \(\overline{\overline{u}}_{l,m}^{y}\) are obtained from (3.15)-(3.16) upon replacing the symbol U by u.

Let \(\boldsymbol {\epsilon} = \boldsymbol {U} - \boldsymbol {u}\) be the discretization error vector and \(\epsilon_{l,m} = U_{l,m} - u_{l,m}\), \(l=1 ( 1 ) L, m=1 ( 1 ) M\) be the point-wise error.

By subtracting (4.2) from (4.1), the discretization error satisfies the relation

$$ - \overline{\nabla}_{h_{l}, k_{m}}^{2} \epsilon_{l,m} + h_{l}^{2} \zeta_{l,m}^{2} \mathcal{L} E_{l,m} = \overline{T}_{l,m},\quad l=1 ( 1 ) L, m=1 ( 1 ) M, $$
(4.5)

where \(E_{\rho,\sigma} =F_{\rho,\sigma} - f_{\rho,\sigma}\), \((\rho,\sigma)\in \mathcal{D}\) and to be explicit

$$\begin{aligned} &E_{\rho,\sigma} = a_{\rho,\sigma} \overline{\epsilon}_{\rho,\sigma}^{x} + b_{\rho,\sigma} \overline{\epsilon}_{\rho,\sigma}^{y} + c_{\rho,\sigma} \epsilon_{\rho,\sigma},\quad (\rho,\sigma)\in \hat{ \mathcal{D}}, \end{aligned}$$
(4.6)
$$\begin{aligned} &E_{l,m} = a_{l,m} \overline{\overline{\epsilon}}_{l,m}^{x} + b_{l,m} \overline{\overline{\epsilon}}_{l,m}^{y} + c_{l,m} \epsilon_{l,m}. \end{aligned}$$
(4.7)

Here, the expressions for \(\overline{\epsilon}_{\rho,\sigma}^{x}\), \(\overline{\epsilon}_{\rho,\sigma}^{y}, (\rho,\sigma)\in \hat{\mathcal{D}}\) and \(\overline{\overline{\epsilon}}_{l,m}^{x}\), \(\overline{\overline{\epsilon}}_{l,m}^{y}\) are obtained from equations (3.3) and (3.15)-(3.16), respectively, upon interchange of U by ϵ.

Representing the system of linear equations (4.5) in a matrix-vector form, one obtains

$$ \boldsymbol {\mathcal{M} \epsilon} + \overline{\boldsymbol {T}} = \boldsymbol {0}, $$
(4.8)

where \(\boldsymbol {T} ( h_{l} ) = [ \overline{T}_{l,m} ]_{l=1 ( 1 ) L,m=1 ( 1 ) M}^{\mathrm{t}}\) is the sixth-order error vector and \(\boldsymbol {\mathcal{M}} = [ \mathcal{M}_{i,j} ]_{i,j=1 ( 1 ) LM} = [ \boldsymbol {\mathcal{R}} \ \boldsymbol {\mathcal{S}} \ \boldsymbol {\mathcal{R}} ]\) is the block tri-diagonal coefficient matrix and, in particular, when \(\max_{l} h_{l} \rightarrow0\), that is, for a sufficiently small value of mesh spacing, we find

$$\boldsymbol {\mathcal{R}} = \frac{\varepsilon}{12} \left [ -( \begin{matrix}{} 1+ \zeta_{i,j}^{2} ) & -2(5- \zeta_{i,j}^{2} ) & -(1+ \zeta_{i,j}^{2} ) \end{matrix} \right ] $$

and

$$\boldsymbol {\mathcal{S}} = \frac{\varepsilon}{6} \left [ \begin{matrix}{} -(5 \zeta_{i,j}^{2} -1) & 10(1+ \zeta_{i,j}^{2} ) & -(5 \zeta_{i,j}^{2} \end{matrix} -1) \right ]\quad \mbox{as tri-diagonal matrices}. $$

Theorem 4.1

The block tri-diagonal matrix \(\boldsymbol {\mathcal{M}}\) is irreducible provided the mesh-ratio parameter \(\zeta_{i,j} \in ( 1/ \sqrt{5}, \sqrt{5} )\) for all \(i,j=1 ( 1 ) LM\).

Proof

It is evident that the matrix \(\boldsymbol {\mathcal{M}}\) has positive diagonal values since \(\zeta_{i,j}\) is positive for all i and j. Further, \(\boldsymbol {\mathcal{M}}\) has non-positive off-diagonal values provided \(1/ \sqrt{5} <\zeta_{i,j} < \sqrt{5}\). Now, if we label LM distinct points in the xy-plane as \(1,2,\ldots,LM\) and draw an arrow from i to j such that \(\mathcal{M}_{i,j} \neq0\), then, for any two distinct points i and j, there exists a directed path that joins the ordered pair of nodes i and j. Therefore, the graph \(\boldsymbol {\mathcal{G}} ( \boldsymbol {\mathcal{M}} )\) of the matrix \(\boldsymbol {\mathcal{M}}\) is connected, and hence \(\boldsymbol {\mathcal{M}}\) is an irreducible matrix [30, 31]. This completes the proof. □

Theorem 4.2

For a sufficiently small value of mesh step-sizes and \(c ( x,y ) \geq0\), the block tri-diagonal matrix \(\boldsymbol {\mathcal{M}}\) is monotone.

Proof

Let \(\zeta = \min_{l,m} \zeta_{l,m}\), \(c = \min_{l,m} c_{l,m}\), \(l=1 ( 1 ) L\), \(m=1 ( 1 ) M\), and \(\vartheta_{j}\), \((j=1,\ldots,LM)\) denote the \(jth\) row element sum of the matrix \(\boldsymbol {\mathcal{M}}\). Since \(c ( x,y ) \geq0\), therefore \(c = \min_{l,m} c ( x_{l}, y_{m} ) \geq0\). Thus, for a sufficiently small mesh step-size (i.e., \(\max_{l} h_{l} \rightarrow0\)), the following asymptotic bounds on the weak row sums can be computed:

$$\begin{aligned} &\vartheta_{j} \geq \frac{11}{12} \varepsilon \bigl( 1+ \zeta^{2} \bigr) >0,\quad j=1,L, \\ &\vartheta_{j} \geq \varepsilon>0,\quad j=2 ( 1 ) L-1, \\ &\vartheta_{ ( r-1 ) L+j} \geq \varepsilon \zeta^{2} >0,\quad r=2 ( 1 ) M-1, j=1,L, \\ &\vartheta_{ ( r-1 ) L+j} \geq \textstyle\begin{cases} c \zeta^{2} h_{1}^{2}, &\alpha>0,\\ c \zeta^{2} h_{L+1}^{2}, &\alpha< 0, \end{cases}\displaystyle \quad r=2 ( 1 ) M-1,j=2 ( 1 ) L-1, c\geq0 \\ &\vartheta_{ ( M-1 ) L+j} \geq \frac{11}{12} \varepsilon \bigl( 1+ \zeta^{2} \bigr) >0,\quad j=1,L, \\ &\vartheta_{ ( M-1 ) L+j} \geq \varepsilon>0,\quad j=2 ( 1 ) L-1. \end{aligned}$$

Note that, except corresponding to the main diagonal, all of the weak row element sums are positive and off-diagonal values of the matrix \(\boldsymbol {\mathcal{M}}\) are non-positive for sufficiently small values of mesh step-size. Since \(\boldsymbol {\mathcal{M}}\) is irreducible (Theorem 4.1), it follows that the matrix \(\boldsymbol {\mathcal{M}}\) is monotone, and this completes the proof. □

Theorem 4.3

The matrix \(\boldsymbol {\mathcal{M}}\) is monotone if and only if the elements of the inverse matrix \(\boldsymbol {\mathcal{M}}^{-1}\) are non-negative.

Proof

See Henrici [32]. □

Theorem 4.4

If \(h= \max_{l} h_{l}\), then \(\Vert \boldsymbol {\epsilon} \Vert _{\infty} \leq O ( h^{4} )\).

Proof

By means of Theorems 4.2 and 4.3, one obtains that the matrix \(\boldsymbol {\mathcal{M}}\) is invertible and \(\boldsymbol {\mathcal{M}}^{-1} \geq0\). For notational convenience, we denote \(\boldsymbol {\mathcal{M}}^{-1} = [ \mathcal{M}_{i,j}^{-1} ]_{i,j=1 ( 1 ) LM}\) , and I is a matrix of order \(LM\times1\) having all entries as one. Then the matrix identity \(\boldsymbol {\mathcal{M}}^{-1} ( \boldsymbol {\mathcal{M} I} ) =\boldsymbol {I}\) gives rise to

$$ \sum_{j=1}^{LM} \mathcal{M}_{i,j}^{-1} \vartheta_{j} =1,\quad i=1 ( 1 ) LM. $$
(4.9)

If \(h= \max_{l} h_{l}\), then \(h= h_{1}\) for \(\alpha<0\) and \(h=h_{L+1}\) for \(\alpha>0\). Thus, in view of each \(\mathcal{M}_{i,j}^{-1} \geq0\) and \(\vartheta_{j} >0\), relation (4.9) gives rise to the following inequalities for \(i=1 ( 1 ) LM\):

$$\begin{aligned} &\mathcal{M}_{i,j}^{-1} \leq \frac{1}{\vartheta_{j}} \leq \frac{12}{ 11\varepsilon ( 1+ \zeta^{2} )} +O ( h ), \quad j=1,L, \\ &\sum_{j=2}^{L-1} \mathcal{M}_{i,j}^{-1} \leq \frac{1}{\min_{j=2 ( 1 ) L-1} \vartheta_{j}} \leq \frac{1}{\varepsilon} +O(h), \\ &\sum_{r=2}^{M-1} \mathcal{M}_{i, ( r-1 ) L+j}^{-1} \leq \frac{1}{ \min_{r=2 ( 1 ) M-1} \vartheta_{ ( r-1 ) L+j}} \leq \frac{1}{\varepsilon \zeta^{2}} +O ( h ),\quad j=1,L, \\ &\sum_{r=2}^{M-1} \sum _{j=2}^{L-1} \mathcal{M}_{i, ( r-1 ) L+j}^{-1} \leq \textstyle\begin{cases} \sum_{j=1}^{LM} \mathcal{M}_{i,j}^{-1} \vartheta_{j} =1,& c=0,\\ \frac{1}{\mathop{\min\limits_{j=2 ( 1 ) L-1}}\limits_{r=2 ( 1 ) M-1} \vartheta_{ ( r-1 ) L+j}} \leq \frac{1}{c\zeta^{2} h^{2}} + O ( h ),& c>0, \end{cases}\displaystyle \\ &\mathcal{M}_{i, ( M-1 ) L+j}^{-1} \leq \frac{1}{ \vartheta_{ ( M-1 ) L+j}} \leq \frac{12}{11\varepsilon ( 1+ \zeta^{2} )} +O ( h ),\quad j=1,L, \\ &\sum_{j=2}^{L-1} \mathcal{M}_{i, ( M-1 ) L+j}^{-1} \leq \frac{1}{ \min_{j=2 ( 1 ) L-1} \vartheta_{ ( M-1 ) L+j}} \leq \frac{1}{\varepsilon} +O(h). \end{aligned}$$

The matrix-vector norm in the following analysis is taken to be

$$\begin{aligned} &\bigl\Vert \boldsymbol {\mathcal{M}}^{-1} \bigr\Vert _{\boldsymbol {\infty}} \\ &\quad= \max_{i=1 ( 1 ) LM} \left [ \begin{matrix}{} \vert \mathcal{M}_{i,1}^{-1} \vert + \sum_{j=2}^{L-1} \vert \mathcal{M}_{i,j}^{-1} \vert + \vert \mathcal{M}_{i,L}^{-1} \vert + \vert \mathcal{M}_{i, ( M-1 ) L+1}^{-1} \vert + \sum_{j=2}^{L-1} \vert \mathcal{M}_{i, ( M-1 ) L+j}^{-1} \vert \\ {}+ \vert \mathcal{M}_{i,LM}^{-1} \vert + \sum_{r=2}^{M-1} ( \vert \mathcal{M}_{i, ( r-1 ) L+1}^{-1} \vert + \sum_{j=2}^{L-1} \vert \mathcal{M}_{i, ( r-1 ) L+j}^{-1} \vert + \vert \mathcal{M}_{i,rL}^{-1} \vert ) \end{matrix} \right ] \end{aligned}$$

and \(\Vert \overline{\boldsymbol {T}} \Vert _{\boldsymbol {\infty}} = \max_{l=1 ( 1 ) L} \vert O ( h_{l}^{6} ) \vert \approx O ( h^{4} )\).

As a result of combining the above inequalities and equation (4.8), we obtain

$$ \Vert \boldsymbol {\epsilon} \Vert _{\infty} \leq \bigl\Vert \boldsymbol { \mathcal{M}}^{-1} \bigr\Vert _{\infty} \boldsymbol {\cdot} \Vert \overline{\boldsymbol {T}} \Vert _{\infty} \leq \frac{1}{ c\zeta^{2}} h^{4} + \frac{2}{11} \frac{11+46\zeta^{2} +11 \zeta^{4}}{ \varepsilon(1+\zeta^{2} ) \zeta^{2}} h^{6} +O \bigl( h^{7} \bigr). $$
(4.10)

That is, \(\Vert \boldsymbol {\epsilon} \Vert _{\infty} \leq O ( h^{4} )\) and \(\Vert \boldsymbol {\epsilon} \Vert _{\infty} \rightarrow0\) as \(h= \max h_{l} \rightarrow0\). This completes the proof. □

Theorem 4.4 corroborates the fourth-order convergence of the new compact scheme (3.19) for obtaining the numerical solution values of convection-dominated problems (1.1)-(1.3) in the quasi-variable mesh network. Note that the restriction \(c= \min_{l,m} c_{l,m} \geq0\) in (4.10) pertains the assumption \(c ( x,y ) \geq0\) on Ω̅ in Section 1.

5 Comparison of numerical and exact solution values

In order to examine the performance of the quasi-variable mesh compact finite difference method discussed in Section 4, we consider some examples having parallel and/or normal boundary layers and investigate the nature of solution values. The maximum absolute errors (\(\mathcal{E}_{\infty}\)) and root mean squared errors (\(\mathcal{E}_{2}\)) of analytical solution values \(U_{l,m}\) and approximate solution values \(u_{l,m}\) at the mesh point (\(x_{l}, y_{m}\)) and corresponding numerical convergence order are computed by means of the formulas

$$\begin{aligned} &\mathcal{E}_{2} = \Biggl[ \frac{1}{LM} \sum _{l=1}^{L} \sum_{m=1}^{M} \vert U_{l,m} - u_{l,m} \vert ^{2} \Biggr]^{1/2}, \\ &\mathcal{E}_{2}^{*} = \Biggl[ \frac{1}{(2L+1)(2M+1)} \sum _{l=1}^{2L+1} \sum _{m=1}^{2M+1} \vert U_{l,m} - u_{l,m} \vert ^{2} \Biggr]^{1/2},\qquad \Theta_{2} = \log_{2} \biggl[ \frac{\mathcal{E}_{2}}{ \mathcal{E}_{2}^{*}} \biggr], \\ &\mathcal{E}_{\infty} =\mathop{ \max_{ l=1 ( 1 ) L}}_{ m=1 ( 1 ) M } |U_{l,m} - u_{l,m} |, \\ &\mathcal{E}_{\infty}^{*} = \mathop{\max_{l=1 ( 1 ) 2L+1}}_{ m=1 ( 1 ) 2M+1 } |U_{l,m} - u_{l,m} |,\qquad \Theta_{\infty} = \log_{2} \biggl[ \frac{\mathcal{E}_{\infty}}{\mathcal{E}_{\infty}^{*}} \biggr]. \end{aligned}$$

The above estimates are obtained using quasi-variable meshes (\(\alpha\neq0\) or \(\beta\neq0\)) as well as for constant mesh step-size (\(\alpha=0\) and \(\beta=0\)). To ease the computational work, we have taken \(L=M\), and Dirichlet’s boundary values are received from the known analytic solutions. The Gauss-Seidel iterative method for the solution of linear difference equations uses the tolerance of error as 10−10 [33]. Simulations with second-order accurate supra-convergent scheme show slow converging results, and thus they are not mentioned in the tabulated results. The Maple’s CodeGeneration tool symbolically obtained the compact schemes, and C programming on the Macintosh operating system performed numerical computing.

Example 5.1

[3, 16]

Consider the singularly perturbed problems

$$\varepsilon \nabla^{2} U ( x,y ) = \frac{1}{1+y} U^{x} ( x,y ) - \frac{1-2\varepsilon(1+y)}{1+y} e^{y-x},\quad 0< x,y< 1 $$

that possess an analytic solution as \(U ( x,y ) = e^{y-x} + (1+y)^{1+1/\varepsilon} / 2^{1/\varepsilon}\). Given \(\varepsilon= 10^{-2}\), the behavior of the solution changes sharply and can be easily captured by means of the quasi-variable mesh compact scheme as presented in Table 1.

Table 1 Accuracies in solution values at \(\pmb{\varepsilon= 10^{-2}}\) in Example 5.1

Example 5.2

[34]

Consider the constant coefficients convection-diffusion problem

$$\nabla^{2} U ( x,y ) =\eta U^{x} (x,y)+ \eta U^{y} (x,y)-\delta U(x,y)+g(x,y),\quad 0< x,y< 1. $$

The analytic solution is taken as \(U ( x,y ) = xye^{xy} \sin (\pi x) \sin (\pi y)\), and thus \(g(x,y)\) can be calculated accordingly. The convergence order with \(\eta=100\pi\), \(\delta=30 \pi^{2}\) is reasonable with coarse quasi-variable meshes, whereas the uniform mesh solution values deteriorate in both order and accuracies as shown in Table 2.

Table 2 Accuracies in solution values in Example 5.2 with \(\pmb{\eta=100\pi,\delta=30 \pi^{2}}\)

Example 5.3

[13]

Consider the variable coefficient two-space dimensional convection-diffusion equation on a rectangular domain \(0< x\), \(y<1\):

$$\varepsilon \nabla^{2} U ( x,y ) = \bigl( x^{2} -1 \bigr) ( 2y-1 ) U^{x} ( x,y ) +2xy ( 1-y ) U^{y} ( x,y ) +g ( x,y ). $$

The analytical solution is \(U ( x,y ) = \sin ( \pi x ) + \sin ( 13\pi x ) + \cos ( \pi x ) + \cos (13\pi x)\). Here, \(\varepsilon=1/ R_{e}\) and \(R_{e}\) is the Reynolds number. Table 3 presents the error in solution values from lower to higher Reynolds number (\(1\leq R_{e} \leq 10^{4}\)) by using the uniform mesh compact scheme, and the expected computational order of accuracy is not achieved. At the same time, the fourth-order compact method on quasi-variable meshes yields a satisfactory result as shown in Table 4.

Table 3 Accuracies in solution values with uniform meshes in Example 5.3
Table 4 Accuracies in solution values with quasi-variable meshes in Example 5.3

6 Concluding remarks

We have proposed a compact finite-difference algorithm for finding fourth-order accurate numerical solution values to the two-dimensional convection-diffusion problems on a rectangular network. Based on the idea of supra-convergence that provides second-order convergence on variable meshes, a new fourth-order convergent method is proposed on quasi-variable meshes, and the application of the new scheme is illustrated with several convection-diffusion problems exhibiting parallel and/or normal boundary layers. The coarse quasi-variable meshes make the algorithm more efficient in comparison with uniformly spaced fine meshes despite the fact that both of them are fourth-order accurate methods. The new method improves maximum absolute and root mean squared errors of solution values as well as their computational order of convergence. It is possible to extend such a technique to three-dimensional convection-diffusion problems.

References

  1. Segal, A: Aspects of numerical methods for elliptic singular perturbation problems. SIAM J. Sci. Comput. 3, 327-349 (1982)

    Article  MathSciNet  MATH  Google Scholar 

  2. Stynes, M: Steady-state convection-diffusion problems. Acta Numer. 14, 445-508 (2005)

    Article  MathSciNet  MATH  Google Scholar 

  3. Tian, ZF, Dai, SQ: High-order compact exponential finite difference methods for convection-diffusion type problems. J. Comput. Phys. 220, 952-974 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  4. Thomas, JW: Numerical Partial Differential Equations: Conservation Laws and Elliptic Equations. Springer, New York (1999)

    Book  MATH  Google Scholar 

  5. Boisvert, RF: Families of high order accurate discretization of some elliptic problems. SIAM J. Sci. Stat. Comput. 2, 268-285 (1981)

    Article  MathSciNet  MATH  Google Scholar 

  6. Strikwerda, JC: Finite Difference Schemes and Partial Differential Equations. SIAM, Philadelphia (2004)

    Book  MATH  Google Scholar 

  7. Ferziger, JH, Peric, M: Computational Methods for Fluid Dynamics. Springer, Berlin (2002)

    Book  MATH  Google Scholar 

  8. Mishra, N, Sanyasiraju, VSSY: Efficient exponential compact higher order difference scheme for convection dominated problems. Math. Comput. Simul. 82, 617-628 (2011)

    Article  MathSciNet  MATH  Google Scholar 

  9. Mohanty, RK, Setia, N: A new compact high order off-step discretization for the system of 2D quasi-linear elliptic partial differential equations. Adv. Differ. Equ. 2013, 223 (2013)

    Article  MathSciNet  Google Scholar 

  10. Saldanha, G, Ananthakrishnaiah, U: A fourth-order finite difference scheme for two-dimensional nonlinear elliptic partial differential equations. Numer. Methods Partial Differ. Equ. 11, 33-40 (1995)

    Article  MathSciNet  MATH  Google Scholar 

  11. Zhai, S, Feng, X, He, Y: A new method to deduce high-order compact difference schemes for two-dimensional Poisson equation. Appl. Math. Comput. 230, 9-26 (2014)

    MathSciNet  Google Scholar 

  12. Mohanty, RK, Singh, S: A new fourth order discretization for singularly perturbed two dimensional non-linear elliptic boundary value problems. Appl. Comput. Math. 175, 1400-1414 (2006)

    MathSciNet  MATH  Google Scholar 

  13. Zhang, J, Kouatchou, J, Ge, L: A family of fourth-order difference schemes on rotated grid for two-dimensional convection-diffusion equation. Math. Comput. Simul. 59, 413-429 (2002)

    Article  MathSciNet  MATH  Google Scholar 

  14. Ghaffar, F, Badshah, N, Islam, S, Khan, MA: Multigrid method based on transformation-free high-order scheme for solving 2D Helmholtz equation on nonuniform grids. Adv. Differ. Equ. 2016 19 (2016)

    Article  MathSciNet  Google Scholar 

  15. Jha, N, Kumar, N: An exponential expanding meshes sequence and finite difference method adopted for two-dimensional elliptic equations. Int. J. Model. Simul. Sci. Comput. 7, 1-17 (2016)

    Article  Google Scholar 

  16. Jha, N, Kumar, N, Sharma, KK: A third (four) order accurate, nine-point compact scheme for mildly-nonlinear elliptic equations in two space variables. Differ. Equ. Dyn. Syst. (2015). doi:10.1007/s12591-015-0263-9

    Google Scholar 

  17. Hegarty, AF, O’Riordan, E, Stynes, M: A comparison of uniformly convergent difference schemes for two-dimensional convection-diffusion problems. J. Comput. Phys. 105, 24-32 (1993)

    Article  MathSciNet  MATH  Google Scholar 

  18. Shih, Y, Kellogg, RB, Tsai, P: A tailored finite point method for convection-diffusion- reaction problems. J. Sci. Comput. 43, 239-260 (2010)

    Article  MathSciNet  MATH  Google Scholar 

  19. Kadalbajoo, MK, Sharma, KK: A numerical method based on finite difference for boundary value problems for singularly perturbed delay differential equations. Appl. Math. Comput. 197, 692-707 (2008)

    MathSciNet  MATH  Google Scholar 

  20. Natesan, S, Ramanujam, N: Booster method for singularly perturbed one-dimensional convection-diffusion Neumann problems. J. Optim. Theory Appl. 99, 53-72 (1998)

    Article  MathSciNet  MATH  Google Scholar 

  21. Bieniasz, LK: A set of compact finite-difference approximations to first and second derivatives, related to the extended Numerov method of Chawla on non-uniform grids. Computing. 81, 77-89 (2007)

    Article  MathSciNet  MATH  Google Scholar 

  22. Manteuffel, TA, White, AB: The numerical solution of second-order boundary value problems on nonuniform meshes. Math. Comput. 47, 511-535 (1986)

    Article  MathSciNet  MATH  Google Scholar 

  23. Roos, HG, Stynes, M, Tobiska, L: Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems, vol. 24. Springer, Berlin (2008)

    MATH  Google Scholar 

  24. Sundqvist, H, Veronis, G: A simple finite-difference grid with non-constant intervals. Tellus A 22, 26-31 (1970)

    MathSciNet  Google Scholar 

  25. Britz, D: Digital Simulation in Electrochemistry. Springer, Berlin (2005)

    Book  MATH  Google Scholar 

  26. Samarskii, AA, Matus, PP, Vabishchevich, PN: Difference Schemes with Operator Factors. Springer, Berlin (2002)

    Book  MATH  Google Scholar 

  27. Saul’yev, VK: Integration of Equations of Parabolic Type by the Method of Nets, vol. 54. Elsevier, Amsterdam (2014)

    MATH  Google Scholar 

  28. Jain, MK, Iyengar, SRK, Jain, RK: Computational Methods for Partial Differential Equations. Wiley, New York (1994)

    Google Scholar 

  29. Kreiss, HO, Manteuffel, TA: Supra-convergent schemes on irregular grids. Math. Comput. 47(176), 537-554 (1986)

    Article  MathSciNet  MATH  Google Scholar 

  30. Young, DM: Iterative Solution of Large Linear Systems. Elsevier, Amsterdam (2014)

    Google Scholar 

  31. Varga, RS: Matrix Iterative Analysis. Springer, Berlin (2000)

    Book  MATH  Google Scholar 

  32. Henrici, P: Discrete Variable Methods in Ordinary Differential Equations. Wiley, New York (1962)

    MATH  Google Scholar 

  33. Saad, Y: Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia (2003)

    Book  MATH  Google Scholar 

  34. Gopaul, A, Sunhaloo, MS, Boojhawon, R, Bhuruth, M: Analysis of incomplete factorizations for a nine-point approximation to a convection-diffusion model problem. J. Comput. Appl. Math. 224, 719-733 (2009)

    Article  MathSciNet  MATH  Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Navnit Jha.

Additional information

Competing interests

The authors declare that they have no competing interest.

Authors’ contributions

NJ discussed the quasi-variable meshes and supra-convergent scheme. NK described the high-order compact scheme and matrix analysis of the difference equations. NJ and NK jointly obtained numerical results. All authors read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jha, N., Kumar, N. A fourth-order accurate quasi-variable mesh compact finite-difference scheme for two-space dimensional convection-diffusion problems. Adv Differ Equ 2017, 64 (2017). https://doi.org/10.1186/s13662-017-1115-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13662-017-1115-4

MSC

Keywords