# Compact schemes for variable coefficient convection-diffusion equations

Anindya Goswami  ,  Kuldip Singh Patel*  and  Pradeep Kumar Sahu

###### Abstract.

Fourth order accurate compact schemes for variable coefficient convection-
diffusion equations are considered. A sufficient condition for stability of
the schemes have been derived using a difference equation based approach. The
constant coefficient problems are considered as a special case, and the
unconditional stability of compact schemes for such case is proved
theoretically. The condition number of the amplification matrix is also
analysed, and an estimate for the same is derived. In order to verify the
derived conditions numerically, MATLAB codes are provided in Appendix of the
manuscript. An example is provided to support the assumption taken to assure
stability.

*Corresponding Author

Anindya Goswami, Department of Mathematics, IISER Pune, India,
<EMAIL_ADDRESS>
Kuldip Singh Patel, Department of Mathematics, IIIT Naya Raipur, India,
<EMAIL_ADDRESS>
Pradeep Kumar Sahu, Department of Mathematics, IIIT Naya Raipur, India,
<EMAIL_ADDRESS>
Corresponding author acknowledge the support under the grant
DST/INTDAAD/P-12/2020, and SERB SIR/2022/000021.

Keywords: Variable coefficient PDEs, Stability, Condition number, Compact
schemes, Convection-diffusion equations.

##  1\. Introduction

The role of convection–diffusion equations in mathematical modelling of
various physical phenomenon is seminal. These physical problems include option
pricing problems in stock market [1], computational fluid dynamics [2], and
various other applications [3, 4]. The development of numerical methods is
inevitable to solve these equations due to the absence of analytical solution
in most of the cases. Various numerical methods are available in literature to
solve these partial differential equations (PDEs) [5, 6, 7].

The second order accurate finite difference method (FDM) is one of them. The
reason for popularity of FDM lies in its effortless implementation and showing
desirable properties (consistency, stability, and convergence etc.)[8, 9, 10].
Indeed, FDM with extended computational stencil may also result with high-
order accuracy, however it complicates the implementation of boundary
conditions. Moreover, proving stability in those cases becomes complicated. To
overcome these complications, high order accurate difference schemes were
developed using compact stencils [11, 12, 13], commonly known as compact
schemes.

Stability of numerical methods is an obvious choice prior to its
implementation for solving the PDEs. The implementation of compact schemes
have already been discussed for various constant-coefficient problems and
their stability have also been proved using von-Neumann method, see [14, 15,
16] and references therein. Specifically, compact schemes have also been
developed for variable coefficient PDEs [17, 18] in literature. However, the
stability analysis in [17] is performed using von-Neumann approach by
considering the coefficients as constants, while [18] have not investigated
the stability of scheme. Note that, the influence of the numerical
implementation of the boundary conditions on the stability of the scheme is
not captured by von-Neumann stability analysis. Moreover, the applicability of
von-Neumann stability analysis is restricted to the linear constant
coefficient finite difference schemes (pp. $166$ [19]). The matrix method for
proving stability of the numerical scheme takes into account the influence of
boundary conditions. Thus, variability in boundary conditions makes it tedious
to derive stability condition using matrix method in most of the cases. The
matrix method stability is discussed in [20] for constant-coefficient
convection-diffusion equations. In [20], the compact schemes are used for
space discretization, however the cubic $C^{1}-$spline collocation method is
applied for temporal discretization.

In this manuscript, we extend the approach discussed in [20] to the variable
coefficient problems. The compact scheme, following [18], is considered to
discretize the spatial dimension of variable coefficient convection diffusion
equation. Further, for temporal discretization, the Backward Euler, and Crank-
Nicolson methods are considered which are different from the discretization
techniques appearing in [20]. In the present paper, the stability of these
fully discrete problems are proved using matrix method.

In addition to the unconditional stability of the scheme, which are shown
theoretically, a highly accurate practical implementation needs good precision
in computation of the amplification matrix too. The related numerical error
can be controlled by controlling the condition number of the matrix which
needs to be inverted for getting the amplification matrix. For this purpose,
an understanding of the upper bound of the condition number becomes important.
In this paper, the relevant condition number is shown to be of order
$\mathcal{O}\left(\frac{\delta v}{\delta z^{2}}\right)$, where $\delta v$, and
$\delta z$ are time and space step sizes respectively. Study of condition
number related to the compact scheme is however less common in the literature.
We refer to [21] where a similar bound has been obtained in a different
context. As per our knowledge, the approach of the derivation, adopted in the
present paper, is a new addition to the literature.

The present manuscript is organized as follows: The compact scheme for
variable coefficient convection-diffusion equation is presented in Sec. 2
followed by its stability analysis. The stability of compact scheme for
constant coefficient problem is proved in Sec. 3. The condition number and
numerical results are given in Sec. 4 and 5 respectively. The conclusion along
with the future scope of proposed work are discussed in Sec.6.

##  2\. Stability analysis for variable coefficient problem

Let $a(z)$ and $b(z)$ be two bounded and twice continuously differentiable
real-valued functions defined on finite open interval $(z_{l},z_{r}),$ where
$b(z)>\epsilon>0$, $\>\forall\>z\in(z_{l},z_{r})$. Then, the variable
coefficient convection-diffusion equation on $(z_{l},z_{r})$ can be written as
follows:

$\frac{\partial u}{\partial v}(v,z)+a(z)\frac{\partial u}{\partial
z}(v,z)-b(z)\frac{\partial^{2}u}{\partial z^{2}}(v,z)=0,$ (1)

where $z\in(z_{l},z_{r})$, and $0\leq v\leq T,$ for some positive constant
$T$. Moreover, we associate the following initial and boundary conditions with
(1)

$\displaystyle u(0,z)$ $\displaystyle=k(z),\quad z\in(z_{l},z_{r}),$ (2)
$\displaystyle u(v,z_{l})$ $\displaystyle=h_{1}(v),\quad
u(v,z_{r})=h_{2}(v),\quad v>0,$ (3)

under the assumptions that $h_{1}$, and $h_{2}$ are smooth functions. Further,
$h_{1}(0)=k(z_{l})$, and $h_{2}(0)=k(z_{r})$. Now, the compact scheme for
variable coefficient convection-diffusion equation (1) is presented following
[22] and considering steady-state form of (1)

$\displaystyle a\frac{\partial u}{\partial z}-b\frac{\partial^{2}u}{\partial
z^{2}}=0.$ (4)

The uniform grids are considered. Indeed we discretize $z$ as
$z_{q}=z_{l}+q\delta z$, where $0\leq q\leq N$ for constant step size $\delta
z$. Standard central difference approximations are considered to discretize
(4) as follows:

$a_{q}\Delta_{z}u_{q}-b_{q}\Delta_{zz}u_{q}-\tau_{q}=0,$ (5)

where $\Delta_{z}u_{q}$, and $\Delta_{zz}u_{q}$ are central difference
operators in space of first and second order respectively. Here
$u_{q}=u(z_{q}),\>a_{q}=a(z_{q})$ and $b_{q}=b(z_{q})$. Moreover, the
truncation error $\tau_{q}$ is given by

$\tau_{q}=\left[\frac{{\delta z}^{2}}{12}\left(2a\frac{\partial^{3}u}{\partial
z^{3}}-b\frac{\partial^{4}u}{\partial
z^{4}}\right)(z_{q})\right]_{q}+\mathcal{O}(\delta z^{4}).$ (6)

In order to accomplish fourth-order accuracy, let us take (4) as an auxiliary
equation and approximate $\frac{\partial^{3}u}{\partial z^{3}}$ and
$\frac{\partial^{4}u}{\partial z^{4}}$ as follows:

$\displaystyle\left.\frac{\partial^{3}u}{\partial z^{3}}\right|_{q}$
$\displaystyle=\left[-\frac{1}{b_{q}}\Delta_{z}b_{q}\Delta_{zz}u_{q}+\frac{a_{q}}{b_{q}}\Delta_{zz}u_{q}+\frac{1}{b_{q}}\Delta_{z}a_{q}\Delta_{z}u_{q}\right]+\mathcal{O}({\delta
z}^{2}),$ $\displaystyle\left.-b\frac{\partial^{4}u}{\partial
z^{4}}\right|_{q}$
$\displaystyle=\left[-\frac{2(\Delta_{z}b_{q})^{2}}{b_{q}}+\frac{3a_{q}}{b_{q}}{\Delta_{z}b_{q}}-\frac{{a_{q}}^{2}}{b_{q}}+\Delta_{zz}b_{q}-2\Delta_{z}a_{q}\right]\Delta_{zz}u_{q}$
$\displaystyle+\left[(2\Delta_{z}b_{q}-a_{q})\frac{\Delta_{z}a_{q}}{b_{q}}-\Delta_{zz}a_{q}\right]\Delta_{z}u_{q}+\mathcal{O}({\delta
z}^{2}).$

Using the above expressions and relations (5)-(6), fourth-order accurate
approximation of steady state equation (4) is given as follows:

$\displaystyle\zeta_{q}\Delta_{z}u_{q}-\alpha_{q}\Delta_{zz}u_{q}+\mathcal{O}({\delta
z}^{4})=0,$ (7)

where

$\displaystyle\zeta_{q}$
$\displaystyle=\left[a_{q}-\frac{\delta{z^{2}}}{12}\left(\frac{a_{q}}{b_{q}}\Delta_{z}a_{q}+\frac{2}{b_{q}}\Delta_{z}a_{q}\Delta_{z}b_{q}-\Delta_{zz}a_{q}\right)\right],$
$\displaystyle\alpha_{q}$
$\displaystyle=\left[b_{q}+\frac{\delta{z^{2}}}{12}\left(\frac{a_{q}}{b_{q}}\Delta_{z}b_{q}-2\Delta_{z}a_{q}+\Delta_{zz}b_{q}-\frac{2}{b_{q}}(\Delta_{z}b_{q})^{2}+\frac{a_{q}^{2}}{b_{q}}\right)\right].$

Let $u_{q}^{m}$ denotes the solution of (1) at time grid point $v^{m}=m\delta
v,m=0,1,\dots,M$ and space mesh point $z_{q}$, and
$\Delta_{v}^{+}=\frac{u^{m+1}-u^{m}}{\delta v}$ for constant time step $\delta
v$. Using relation (7), the fully-discrete problem for the unsteady equation
(1) can be written as

$\displaystyle\left[1-\frac{{\delta
z}^{2}}{12}\left(\frac{a_{q}}{b_{q}}+\frac{2\Delta_{z}b_{q}}{b_{q}}\right)\Delta_{z}+\frac{{\delta
z}^{2}}{12}\Delta_{zz}\right]\Delta_{v}^{+}u_{q}^{m}+\zeta_{q}\Delta_{z}u_{q}^{m}-\alpha_{q}\Delta_{zz}u_{q}^{m}+\mathcal{O}({\delta
z}^{4})=0.$ (8)

The expressions for $\Delta_{v}^{+}u_{q}^{m},\>\>\Delta_{z}u_{q}^{m}$ and
$\Delta_{zz}u_{q}^{m}$ are given as

$\Delta_{v}^{+}u_{q}^{m}=\frac{u_{q}^{m+1}-u_{q}^{m}}{\delta
v},\>\>\Delta_{z}u^{m}_{q}=\frac{u^{m}_{q+1}-u^{m}_{q-1}}{2\delta
z},\>\>\Delta_{zz}u^{m}_{q}=\frac{u^{m}_{q+1}-2u^{m}_{q}+u^{m}_{q-1}}{\delta
z^{2}}.$ (9)

Let
$\gamma_{q}=\left(\frac{a_{q}}{b_{q}}+\frac{2}{b_{q}}\Delta_{z}b_{q}\right)$
then fully discrete problem for (1) can be written as follows:

$\displaystyle\Delta_{v}^{+}u_{q}^{m}$ $\displaystyle-\frac{\delta
z^{2}}{12}\left(\gamma_{q}\Delta_{z}\Delta_{v}^{+}u_{q}^{m}-\Delta_{zz}\Delta_{v}^{+}u_{q}^{m}\right)+(1-\theta)\left(\zeta_{q}\Delta_{z}u_{q}^{m}-\alpha_{q}\Delta_{zz}u_{q}^{m}\right)$
$\displaystyle+\theta\left(\zeta_{q}\Delta_{z}u_{q}^{m+1}-\alpha_{q}\Delta_{zz}u_{q}^{m+1}\right)=\mathcal{O}(\\{\delta
v\\}^{p},\\{\delta z\\}^{4}),$ (10)

for $1\leq q\leq N-1,$ and $m=0,1,\dots,M-1$. The order $p$ in the error term
on the right side of (2) depends on $\theta$. The Crank-Nicolson, and backward
Euler time discretization appears for $\theta=\frac{1}{2}$, and $1$
respectively. The order of accuracy in time is $p=2$ for $\theta=1/2$, and
$p=1$ for $\theta=1$. If $U^{m}_{q}$ denotes the approximate value of
$u^{m}_{q}$, then (2) can be rewritten for $1\leq q\leq N-1,$ and $m\geq 0$ as
follows:

$\displaystyle\Delta_{v}^{+}U_{q}^{m}$ $\displaystyle-\frac{\delta
z^{2}}{12}\left(\gamma_{q}\Delta_{z}\Delta_{v}^{+}U_{q}^{m}-\Delta_{zz}\Delta_{v}^{+}U_{q}^{m}\right)+(1-\theta)\left(\zeta_{q}\Delta_{z}U_{q}^{m}-\alpha_{q}\Delta_{zz}U_{q}^{m}\right)$
$\displaystyle+\theta\left(\zeta_{q}\Delta_{z}U_{q}^{m+1}-\alpha_{q}\Delta_{zz}U_{q}^{m+1}\right)=0.$
(11)

with $U^{m}_{0}=h_{1}(v_{m})$, and $U^{m}_{N}=h_{2}(v_{m})$, for all $m\geq
0$. We rewrite (2) by plugging in the expressions of difference operators,
appearing in (9). Let us denote the column vector
$U^{m}=[U^{m}_{1},U^{m}_{2},...,U^{m}_{N-1}],$ $\forall$ $m\geq 0$, and define
the constants depending on $\delta z$ and $\delta v$ as follows:

$\displaystyle p_{i}=\frac{2+\delta z\gamma_{i}}{24\delta
v},\>q_{i}=\frac{5}{6\delta v},\>\text{and}\>r_{i}=\frac{2-\delta
z\gamma_{i}}{24\delta v},\>\forall\>i=1,2,\dots,N-1.\>$ (12) $\displaystyle
l_{i}=-\frac{\zeta_{i}}{2\delta z}-\frac{\alpha_{i}}{\delta
z^{2}},\>m_{i}=\frac{2\alpha_{i}}{\delta
z^{2}},\>\text{and}\>n_{i}=\frac{\zeta_{i}}{2\delta
z}-\frac{\alpha_{i}}{\delta z^{2}},\>\forall\>i=1,2,\dots,N-1.$ (13)

Then Eq. (2) can be written as

$[(1-\theta)X+\theta(X+Y)]U^{m+1}=[\theta X+(1-\theta)(X-Y)]U^{m}+F^{m},$ (14)

where $X$ and $Y$ are the following tridiagonal matrices of order
$(N-1)\times(N-1)$:

$X\vcentcolon=\left[\begin{array}[]{c c c c c c c c
c}q_{1}&r_{1}&&&\cdots&&&&\\\ p_{2}&q_{2}&r_{2}&&\cdots&&&&\\\
&p_{3}&q_{3}&r_{3}&\cdots&&&&\\\ &&&&\vdots&&&&\\\
&&&&\cdots&p_{N-3}&q_{N-3}&r_{N-3}&\\\ &&&&\cdots&&p_{N-2}&q_{N-2}&r_{N-2}\\\
&&&&\cdots&&&p_{N-1}&q_{N-1}\end{array}\right],\quad$ (15)
$Y\vcentcolon=\left[\begin{array}[]{c c c c c c c c
c}m_{1}&n_{1}&&&\cdots&&&&\\\ l_{2}&m_{2}&n_{2}&&\cdots&&&&\\\
&l_{3}&m_{3}&n_{3}&\cdots&&&&\\\ &&&&\vdots&&&&\\\
&&&&\cdots&l_{N-3}&m_{N-3}&n_{N-3}&\\\ &&&&\cdots&&l_{N-2}&m_{N-2}&n_{N-2}\\\
&&&&\cdots&&&l_{N-1}&m_{N-1}\end{array}\right],$ (16)

and $F^{m}$ is a $N-1$ dimensional column vector consists of boundary values.
Its expression is given for each $\theta$ in the following subsections.

###  2.1. Backward Euler compact scheme

Taking $\theta=1$ in (14), the fully discrete backward Euler compact scheme
can be written as:

$(X+Y)U^{m+1}=XU^{m}+F^{m},$ (17)

where $X$ and $Y$ are given in (15), and (16) respectively. Moreover

$F^{m}=[p_{1}h_{1}(v^{m})-(p_{1}+l_{1})h_{1}(v^{m+1}),0,\dots,0,r_{N-1}h_{2}(v^{m})-(r_{N-1}+n_{N-1})h_{2}(v^{m+1})].$

$X$ is diagonally dominant for sufficiently small $\delta z$ due to the
boundedness of $a$, $b$ and $b^{-1}$ on the domain. Hence $X$ is invertible.
Using $W=X^{-1}Y$ in (17), we get

$(I+W)U^{m+1}=U^{m}+X^{-1}F^{m}.$ (18)

The stability of (18), is asserted in Theorem 3. To establish that we need to
study the eigenvalues of $W$.

######  Remark 1.

To locate the eigenvalues of matrix $W$, we aim to derive the characteristic
polynomial. Note that the entries of $W$ are not tractable since $W$ involves
$X^{-1}$. Hence the determinant approach for characteristic polynomial is not
applicable. Using the tridiagonal structure of matrices $X$ and $Y$, we
propose a difference equation approach for deriving the characteristic
polynomial of $W$.

Let $\beta=\left[\beta_{1},\beta_{2},\ldots,\beta_{N-1}\right]^{T}$ be an
eigenvector corresponding to $\lambda$. So we have $X^{-1}Y\beta=\lambda\beta$
or $(Y-\lambda X)\beta=0$. This linear system in $\beta$ can be written
explicitly as

$\displaystyle\left.\begin{array}[]{rl}\left(m_{1}-\lambda
q_{1}\right)\beta_{1}+\left(n_{1}-\lambda r_{1}\right)\beta_{2}&=0,\\\
\left(l_{j}-\lambda p_{j}\right)\beta_{j-1}+\left(m_{j}-\lambda
q_{j}\right)\beta_{j}+\left(n_{j}-\lambda r_{j}\right)\beta_{j+1}&=0,\quad
j=2,\ldots,N-2,\\\ \left(l_{N-1}-\lambda
p_{N-1}\right)\beta_{N-2}+\left(m_{N-1}-\lambda
q_{N-1}\right)\beta_{N-1}&=0.\end{array}\right\\}$

Then the above system can be rewritten as

$C_{j}\beta_{j-1}+B_{j}\beta_{j}+A_{j}\beta_{j+1}=0,\>\forall\>j=1,2,\dots,N-1,$
(19)

with constraints $\beta_{0}=\beta_{N}=0$, where

$C_{j}=l_{j}-\lambda p_{j},\>B_{j}=m_{j}-\lambda
q_{j},\>\text{and}\>A_{j}=n_{j}-\lambda r_{j}.$ (20)

The solutions of (19) for the unknown vector
$(\beta_{0},\beta_{1},\dots,\beta_{N})$ form a linear space of dimension $2$
(Theorem 2.21 [23]). If $\beta^{{}^{\prime}}$ and $\beta^{{}^{\prime\prime}}$
denote two independent solutions, without loss of generality, we may take
$\beta_{0}^{{}^{\prime}}=\beta_{0}^{{}^{\prime\prime}}=1$ and
$\beta_{1}^{{}^{\prime}}\neq\beta_{1}^{{}^{\prime\prime}}$. For $2\leq j\leq
N$, (19) gives

$\displaystyle\beta_{j}^{{}^{\prime}}=-\left(\frac{C_{j-1}}{A_{j-1}}\right)\beta_{j-2}^{{}^{\prime}}-\left(\frac{B_{j-1}}{A_{j-1}}\right)\beta_{j-1}^{{}^{\prime}},\quad\beta_{j}^{{}^{\prime\prime}}=-\left(\frac{C_{j-1}}{A_{j-1}}\right)\beta_{j-2}^{{}^{\prime\prime}}-\left(\frac{B_{j-1}}{A_{j-1}}\right)\beta_{j-1}^{{}^{\prime\prime}}.$
(21)

Thus by successive substitution

$\beta_{N}^{{}^{\prime}}=\sum_{s\in\mathscr{A}}s^{\prime}_{0}\prod_{r=1}^{N-1}s_{r},\quad\beta_{N}^{{}^{\prime\prime}}=\sum_{s\in\mathscr{A}}s^{\prime\prime}_{0}\prod_{r=1}^{N-1}s_{r},$
(22)

where $s\in\mathscr{A}$ if $s:=\\{s_{r}\\}_{r=1}^{N-1}$ is such that
$s_{N-1}\in\left\\{-\frac{B_{N-1}}{A_{N-1}},-\frac{C_{N-1}}{A_{N-1}}\right\\}$,
and for all other $r\geq 2$

$s_{r-1}\in\left\\{\begin{array}[]{ll}\\{1\\}&\textrm{ if
}s_{r}=-\frac{C_{r}}{A_{r}},\\\
\left\\{-\frac{B_{r-1}}{A_{r-1}},-\frac{C_{r-1}}{A_{r-1}}\right\\}&\textrm{ if
}s_{r}\in\left\\{1,-\frac{B_{r}}{A_{r}}\right\\}.\\\ \end{array}\right.$ (23)

If $s_{1}\in\left\\{1,-\frac{B_{1}}{A_{1}}\right\\}$, $s^{\prime}_{0}$ and
$s^{\prime\prime}_{0}$ are $\beta_{1}^{\prime}$ and $\beta_{1}^{\prime\prime}$
respectively, otherwise $s^{\prime}_{0}=s^{\prime\prime}_{0}=1$. A
diagrammatic illustration of the above algebraic description of $\mathscr{A}$
is given in Fig 1 for $N=4$ case. For general $N$, we provide with an
Algorithm A.1 in Appendix involving symbolic computation, that computes
$\beta^{{}^{\prime}}_{N}$ by solving the recurrence relation (21). It is
interesting to note that $\\#(\mathscr{A})$ is $F_{N+1}$, where $F_{N}$ is the
$N^{th}$ Fibonacci number. We introduce the following notations

$D_{N}^{1}=\left(\prod_{j=1}^{N-1}A_{j}\right)\sum_{\left\\{s\in\mathscr{A}|s_{1}\neq-\frac{C_{1}}{A_{1}}\right\\}}\prod_{r=1}^{N-1}s_{r},\quad
D_{N}^{2}=\left(\prod_{j=1}^{N-1}A_{j}\right)\sum_{\left\\{s\in\mathscr{A}|s_{1}=-\frac{C_{1}}{A_{1}}\right\\}}\prod_{r=1}^{N-1}s_{r}.$
(24)

Here $D_{N}^{1}$ and $D_{N}^{2}$ are polynomials of degree atmost $N-1$ and
are not dependent on $\beta^{{}^{\prime}}$ or $\beta^{{}^{\prime\prime}}$. The
general solution of (19) is
$\beta_{j}=K_{1}\beta_{j}^{{}^{\prime}}+K_{2}\beta_{j}^{{}^{\prime\prime}}\>\>\forall\>j=1,2,\dots,N-1.$
We are interested in solution which satisfies the boundary conditions
$\beta_{0}=0=\beta_{N}$. The condition $\beta_{0}=0$ gives
$K_{1}=-K_{2},\>\>K_{1}\neq 0,K_{2}\neq 0$. Additionally using $\beta_{N}=0$,
we get $\beta_{N}^{{}^{\prime}}=\beta_{N}^{{}^{\prime\prime}}$. Multiplying
both sides of this equality by $\prod_{j=1}^{N-1}A_{j}$ and using expressions
in (22) and (24), we get
$D_{N}^{1}\beta_{1}^{{}^{\prime}}+D_{N}^{2}=D_{N}^{1}\beta_{1}^{{}^{\prime\prime}}+D_{N}^{2}.$
Hence $D_{N}^{1}(\beta_{1}^{{}^{\prime}}-\beta_{1}^{{}^{\prime\prime}})=0.$
Since $\beta_{1}^{{}^{\prime}}\neq\beta_{1}^{{}^{\prime\prime}}$, hence
$D_{N}^{1}=0$. Thus, we have proved the following result.

Figure 1. Diagram for (23) with $N=4$. Each path gives a finite sequence
$s\in\mathscr{A}$.

######  Theorem 2.

Let $\lambda$ be an eigenvalue of matrix $W$, then $\lambda$ solves
$D_{N}^{1}=0$, where $D_{N}^{1}$ is given in (24).

######  Theorem 3.

If the roots of polynomial $D^{1}_{N}$ has positive real part, then backward
Euler compact scheme (18) is stable.

###### Proof.

From Theorem 2, the roots of $D_{N}^{1}$ are the eigenvalues of $W$. Rewriting
the equation (18) as

$U^{m+1}=HU^{m}+(X+Y)^{-1}F^{m},$

where $H=(I+W)^{-1}$. If $\lambda$ and $\beta$ is a pair of eigenvalue and
eigenvector of $W$, we have $H\beta=\frac{1}{1+\lambda}\beta$. Consequently,
$\frac{1}{1+\lambda}$ is an eigenvalue of $H$. Hence, backward Euler compact
scheme (36) is stable if we have $\frac{1}{|1+\lambda|}<1$ for each $\lambda$.
The above inequality follows as real part of each $\lambda$ is assumed to be
positive. ∎

######  Remark 4.

Note that, the assertion in Theorem 3 holds under certain assumption. An
Algorithm A.1 is provided in Appendix to compute $D^{1}_{N}$ symbolically. In
Example 1, we showed with specific choices of parameters that all the roots of
$D^{1}_{N}$ have positive real part, hence it is not an unrealistic
assumption.

######  Remark 5.

Using a rough upper bound on the growth of Fibonacci series, the growth of
complexity of computing $D^{1}_{N}$ can be estimated. That estimated
complexity turns out to be dominated by $\mathcal{O}(N2^{N})$. On the other
hand the other approach of deriving characteristic polynomial involves
computation of determinant having a complexity of order factorial of $N$,
according to the Laplace expansion. As $N!$ grows faster than $N2^{N}$,
computation of $D^{1}_{N}$ appears advantageous. However, there are smarter
algorithms for computing determinant with cubic complexity, whereas it is not
yet clear how to further reduce the complexity for computing $D^{1}_{N}$.

###  2.2. Crank-Nicolson compact scheme

Let us define $\tilde{Y}=\frac{1}{2}Y$, and thus
$\tilde{W}:=X^{-1}\tilde{Y}=\frac{1}{2}W.$ Taking $\theta=1/2,$ in (14) the
fully discrete Crank-Nicolson compact scheme can be written as

$(X+\tilde{Y})U^{m+1}=(X-\tilde{Y})U^{m}+\tilde{F}^{m},$ (25)

where

$\displaystyle\tilde{F}^{m}=\left[\left(p_{1}-\frac{l_{1}}{2}\right)h_{1}(v_{m})-\left(p_{1}+\frac{l_{1}}{2}\right)h_{1}(v_{m+1}),0,...,0,\left(r_{N-1}-\frac{n_{N-1}}{2}\right)h_{2}(v_{m})\right.$
$\displaystyle\left.-\left(r_{N-1}+\frac{n_{N-1}}{2}\right)h_{2}(v_{m+1})\right].$

Furthermore, (25) can be written as:

$(I+\tilde{W})U^{m+1}=(I-\tilde{W})U^{m}+X^{-1}\tilde{F}^{m}.$ (26)

######  Theorem 6.

If the roots of polynomial $D^{1}_{N}$ has positive real part, then Crank-
Nicolson compact scheme (26) is stable.

###### Proof.

Let $\tilde{\lambda}$ be an eigenvalue of matrix $\tilde{W}$, then using
$2\tilde{W}=W$, and Theorem 2, $2\tilde{\lambda}$ is a root of the polynomial
$D^{1}_{N}$. Hence, if every root of $D^{1}_{N}$ has positive real part, the
real part of $\tilde{\lambda}$ is also positive. Therefore $\tilde{W}$ is
invertible, and we can rewrite the Eq. (26) as

$U^{m+1}=\tilde{H}U^{m}+(X+\tilde{Y})^{-1}\tilde{F}^{m},$

where $\tilde{H}=(I+\tilde{W})^{-1}(I-\tilde{W})$. Moreover,
$\frac{1-\tilde{\lambda}}{1+\tilde{\lambda}}$ is an eigenvalue of $\tilde{H}$
iff $\tilde{\lambda}$ is an eigenvalue of $\tilde{W}$. Thus, Crank-Nicolson
compact scheme (26) is stable if we have
$\left\lvert\frac{1-\tilde{\lambda}}{1+\tilde{\lambda}}\right\rvert<1$ for
each $\tilde{\lambda}$. Since real part of $\tilde{\lambda}$ is positive,
therefore $\left\lvert 1-\tilde{\lambda}\right\rvert<\left\lvert
1+\tilde{\lambda}\right\rvert$. Thus
$\left\lvert\frac{1-\tilde{\lambda}}{1+\tilde{\lambda}}\right\rvert<1,$ as
desired. ∎

##  3\. Special Case: Constant coefficient problem

In this section we consider constant coefficient version of convection
diffusion equation (1). We prove below that the assumption on the
characteristic polynomial $D^{1}_{N}$ in Theorems 3 and 4 are true for this
special case. Without loss of generality, we take $a(z)=c=b(z),$ then we have

$\frac{\partial u}{\partial v}(v,z)+c\frac{\partial u}{\partial
z}(v,z)-c\frac{\partial^{2}u}{\partial z^{2}}(v,z)=0.$ (27)

Using (2), the compact scheme for equation (27) is as follows:

$\displaystyle\Delta_{v}^{+}U_{q}^{m}$ $\displaystyle-\frac{\delta
z^{2}}{12}\left(\Delta_{z}\Delta_{v}^{+}U_{q}^{m}-\Delta_{zz}\Delta_{v}^{+}U_{q}^{m}\right)+(1-\theta)c\left(\Delta_{z}U_{q}^{m}-\left(1+\frac{\delta
z^{2}}{12}\right)\Delta_{zz}U_{q}^{m}\right)$ $\displaystyle+\theta
c\left(\Delta_{z}U_{q}^{m+1}-\left(1+\frac{\delta
z^{2}}{12}\right)\Delta_{zz}U_{q}^{m+1}\right)=0.$ (28)

We define the constants depending on $\delta z$ and $\delta v$ as

$\displaystyle c_{1}=\frac{2+\delta z}{24\delta v},\>c_{2}$
$\displaystyle=\frac{5}{6\delta v},\>\text{and}\>c_{3}=\frac{2-\delta
z}{24\delta v},$ (29) $\displaystyle y_{1}=-\frac{c}{2\delta
z}-\frac{\left(c+\frac{c\delta z^{2}}{12}\right)}{\delta z^{2}},\>y_{2}$
$\displaystyle=2\frac{\left(c+\frac{c\delta z^{2}}{12}\right)}{\delta
z^{2}},\>\text{and}\>y_{3}=\frac{c}{2\delta z}-\frac{\left(c+\frac{c\delta
z^{2}}{12}\right)}{\delta z^{2}}.$ (30)

Further, denote $\frac{\delta v}{\delta z^{2}}$ as $d$ and we write

$y_{1}=\frac{-c}{2\delta v}\left[\left(2+\delta z\right)d+\frac{\delta
v}{6}\right],\quad y_{2}=\frac{c}{\delta v}\left(2d+\frac{\delta
v}{6}\right),\quad\mbox{and}\quad y_{3}=\frac{-c}{2\delta
v}\left[\left(2-\delta z\right)d+\frac{\delta v}{6}\right].$

Then Eq. (3) can be written as

$[(1-\theta)X_{1}+\theta(X_{1}+Y_{1})]U^{m+1}=[\theta
X_{1}+(1-\theta)(X_{1}-Y_{1})]U^{m}+F^{m},$ (31)

where $X_{1}$ and $Y_{1}$ are the following tridiagonal matrices of order
$(N-1)\times(N-1)$:

$X_{1}\vcentcolon=\left[\begin{array}[]{c c c c c c c c
c}c_{2}&c_{3}&&&\cdots&&&&\\\ c_{1}&c_{2}&c_{3}&&\cdots&&&&\\\
&c_{1}&c_{2}&c_{3}&\cdots&&&&\\\ &&&&\vdots&&&&\\\
&&&&\cdots&c_{1}&c_{2}&c_{3}&\\\ &&&&\cdots&&c_{1}&c_{2}&c_{3}\\\
&&&&\cdots&&&c_{1}&c_{2}\end{array}\right],\quad$ (32)
$Y_{1}\vcentcolon=\left[\begin{array}[]{c c c c c c c c
c}y_{2}&y_{3}&&&\cdots&&&&\\\ y_{1}&y_{2}&y_{3}&&\cdots&&&&\\\
&y_{1}&y_{2}&y_{3}&\cdots&&&&\\\ &&&&\vdots&&&&\\\
&&&&\cdots&y_{1}&y_{2}&y_{3}&\\\ &&&&\cdots&&y_{1}&y_{2}&y_{3}\\\
&&&&\cdots&&&y_{1}&y_{2}\end{array}\right],$ (33)

and $F_{1}^{m}$ is a $N-1$ dimensional coloum vector consists of boundary
values. Its expression is given for each $\theta$ in the following section. In
the next section, the stability of the above fully discrete problem (31) for
Backward Euler ($\theta=1$) and Crank-Nicolson ($\theta=\frac{1}{2}$) schemes
are proved with no additional assumption on the coefficient.

If the real constants $A,B$ and $C$ are such that $\frac{B^{2}}{4AC}\in[0,1]$,
then clearly there is a unique $\phi\in[0,\pi]$ so that
$\cos^{2}{\frac{\phi}{2}}=\frac{B^{2}}{4AC}$. A direct calculation gives that
this identity is equivalent to the following relation

$\frac{(2B^{2}-4AC+2B\sqrt{B^{2}-4AC})}{4AC}=2\cos^{2}{\frac{\phi}{2}}-1+2\cos{\frac{\phi}{2}}\sqrt{-\sin^{2}{\frac{\phi}{2}}}=\cos{\phi}+i\sin{\phi},$
(34)

where $i=\sqrt{-1}$. The above result has also appeared in [20]. We will use
this relation in the following sections.

###  3.1. Backward Euler compact scheme

Taking $\theta=1$ in (31), the fully discrete backward Euler compact scheme
can be written as:

$(X_{1}+Y_{1})U^{m+1}=X_{1}U^{m}+F_{1}^{m},$ (35)

where $X_{1}$ and $Y_{1}$ are given in (32), and (33) respectively. Moreover

$F_{1}^{m}=[c_{1}h_{1}(v^{m})-(c_{1}+y_{1})h_{1}(v^{m+1}),0,\dots,0,c_{3}h_{2}(v^{m})-(c_{3}+y_{3})h_{2}(v^{m+1})].$

Evidently, $X_{1}$ is diagonally dominant, and hence invertible. Using
$W_{1}=X_{1}^{-1}Y_{1}$ in (35), we get

$(I+W_{1})U^{m+1}=U^{m}+X_{1}^{-1}F_{1}^{m}.$ (36)

The stability of (36), is asserted in Lemma 8. To establish that we prove the
following results.

######  Theorem 7.

Let $\lambda$ be a eigenvalue of matrix $W_{1}$, then the real part of
$\lambda$ is positive.

###### Proof.

Let $\beta=\left[\beta_{1},\beta_{2},\ldots,\beta_{N-1}\right]^{T}$ be an
eigenvector corresponding to $\lambda$. So we have
$X_{1}^{-1}Y_{1}\beta=\lambda\beta$ or $(Y_{1}-\lambda X_{1})\beta=0$. This
linear system in $\beta$ can be written explicitly as

$\displaystyle\left.\begin{array}[]{rl}\left(y_{2}-\lambda
c_{2}\right)\beta_{1}+\left(y_{3}-\lambda c_{3}\right)\beta_{2}&=0,\\\
\left(y_{1}-\lambda c_{1}\right)\beta_{j-1}+\left(y_{2}-\lambda
c_{2}\right)\beta_{j}+\left(y_{3}-\lambda c_{3}\right)\beta_{j+1}&=0,\quad
j=2,\ldots,N-2,\\\ \left(y_{1}-\lambda
c_{1}\right)\beta_{N-2}+\left(y_{2}-\lambda
c_{2}\right)\beta_{N-1}&=0.\end{array}\right\\}$

Let us introduce $\beta_{0}=\beta_{N}=0$. Then the above system can be
rewritten as

$\left(y_{1}-\lambda c_{1}\right)\beta_{j-1}+\left(y_{2}-\lambda
c_{2}\right)\beta_{j}+\left(y_{3}-\lambda c_{3}\right)\beta_{j+1}=0,$ (37)

for all $j=1,2,\dots,N-1$. Let $A=\left(y_{3}-\lambda
c_{3}\right),B=\left(y_{2}-\lambda c_{2}\right)$, and $C=\left(y_{1}-\lambda
c_{1}\right)$. We first notice that if either of $A$ or $C$ is zero, then by
using (37) iteratively, one obtains $\beta_{j}=0\>\forall\>j=1,2,...,N-1$.
Thus solutions of $A=0$ or $C=0$ for $\lambda$ do not give eigenvalues. Hence,
we consider the cases where $AC\neq 0$. We complete the argument by
considering two complementary cases.

Case I: Assume that $B^{2}-4AC\neq 0$. Since $\lambda$ is such that $C\neq 0$,
the solution-space of (37) is known to have dimension $2$ (Theorem $2.21$,
[23]). Since $C\neq 0$ and $B^{2}-4AC\neq 0$, the pair of roots $\nu_{1}$ and
$\nu_{2}$ of $Aq^{2}+Bq+C$ are distinct and nonzero. Using these roots one can
directly verify that the general solution of difference scheme
$(\ref{linear_system_FB})$ is

$\beta_{j}=K_{1}{\nu_{1}}^{j}+K_{2}{\nu_{2}}^{j},\quad j=0,1,\ldots,N.$ (38)

In order to satisfy $\beta_{0}=0$, we must have $K_{1}=-K_{2}\neq 0$. Note
that only the nonzero coefficients are relevant as $\beta$ is an eigenvector.
Again, to satisfy $\beta_{N}=0$, i.e.
$K_{1}\left(\nu_{1}^{N}-\nu_{2}^{N}\right)=0$ we must have
$\left(\frac{\nu_{1}}{\nu_{2}}\right)^{N}=1$. Therefor, every eigenvalue
$\lambda$ of $W$ for which $B^{2}-4AC\neq 0$, is such that the roots $\nu_{1}$
and $\nu_{2}$ of $Aq^{2}+Bq+C$ satisfy

$\frac{\nu_{1}}{\nu_{2}}=e^{i\frac{2k\pi}{N}}=\cos\left(\frac{2k\pi}{N}\right)+i\sin\left(\frac{2k\pi}{N}\right),\quad\textrm{
for some }k=1,2,\ldots,N-1.$

Again since $\nu_{1}$ and $\nu_{2}$ are the solutions of quadratic equation,
we also get

$\frac{\nu_{1}}{\nu_{2}}=\frac{-B-\sqrt{B^{2}-4AC}}{-B+\sqrt{B^{2}-4AC}}=\frac{2B^{2}-4AC+2B\sqrt{B^{2}-4AC}}{4AC}.$

From above two expressions of $\frac{\nu_{1}}{\nu_{2}}$, and (34), we have a
simpler relation

$\frac{B^{2}}{4AC}=\cos^{2}\left(\frac{k\pi}{N}\right),\quad\textrm{ for some
}k=1,2,\ldots,N-1.$

Thus the eigenvalues can be obtained by solving the following family of
equations

$B^{2}-4ACt=0,\quad
t=\cos^{2}\left(\frac{k\pi}{N}\right)\in[0,1),\quad\textrm{ for some
}k=1,2,\ldots,N-1.$ (39)

If $t$ happens to be zero, to satisfy (39), $B$ must be zero. That corresponds
to an eigenvalue $\lambda=\frac{6c}{5}\left(2d+\frac{\delta v}{6}\right)$,
which is real and positive. For $t\in(0,1)$, substituting the expressions of
$A$, $B$, and $C$ in $(\ref{eq:FB_cos_eq_B})$, we obtain after rearrangement

$(c_{2}^{2}-4c_{1}c_{3}t)\lambda^{2}-(2c_{2}y_{2}-4c_{1}y_{3}t-4c_{3}y_{1}t)\lambda+(y_{2}^{2}-4y_{1}y_{3}t)=0.$
(40)

Using (29)-(30), the coefficient of $\lambda^{2}$ is
$(c_{2}^{2}-4c_{1}c_{3}t)=(25-t+\delta z^{2}t/4)/(36\delta v^{2})$ which is
positive as $t\in(0,1)$. Similarly the constant term
$(y_{2}^{2}-4y_{1}y_{3}t)=c^{2}[4(1-t)(\frac{1}{\delta
z^{2}}+\frac{1}{12})^{2}+\frac{t}{\delta z^{2}}]$ is also positive. On the
other hand the coefficient of $\lambda$, that is
$-(2c_{2}y_{2}-4c_{1}y_{3}t-4c_{3}y_{1}t)=-\frac{2c}{3\delta
v}\left(\frac{1}{\delta z^{2}}(5+t)+\frac{5}{12}-\frac{t}{6}\right)$, is a
negative quantity. Therefore, both the sum and product of the roots of (40)
are positive. Hence, the real parts of the roots are positive.

Case II: Assume $B^{2}-4AC=0$, which gives

$(c_{2}^{2}-4c_{1}c_{3})\lambda^{2}-(2c_{2}y_{2}-4c_{1}y_{3}-4c_{3}y_{1})\lambda+(y_{2}^{2}-4y_{1}y_{3})=0.$
(41)

Again by simplifying the coefficients using (29)-(30), as in Case I, we
observe that both the sum and product of the roots of (41) are positive.
Hence, the real parts of the roots are positive. ∎

######  Theorem 8.

The backward Euler compact scheme (36) is stable.

###### Proof.

The proof follows using the assertion of Theorem 7. The argument is analogous
to that in the proof of Theorem 3. We omit the details. ∎

###  3.2. Crank Nicolson compact scheme

Let us define $\tilde{Y}_{1}=\frac{1}{2}Y_{1}$. Hence its sub-diagonal,
diagonal, and super-diagonal entries are $\tilde{y}_{1}=\frac{1}{2}y_{1}$,
$\tilde{y}_{2}=\frac{1}{2}y_{2}$, and $\tilde{y}_{3}=\frac{1}{2}y_{3}$
respectively. Taking $\theta=1/2$ in (31), we get the following fully discrete
Crank-Nicolson compact scheme:

$(X_{1}+\tilde{Y}_{1})U^{m+1}=(X_{1}-\tilde{Y}_{1})U^{m}+\tilde{F}_{1}^{m},$
(42)

where
$\tilde{F}_{1}^{m}=[(c_{1}-\tilde{y}_{1})h_{1}(v_{m})-(c_{1}+\tilde{y}_{1})h_{1}(v_{m+1}),0,...,0,(c_{3}-\tilde{y}_{3})h_{2}(v_{m})-(c_{3}+\tilde{y}_{3})h_{2}(v_{m+1})]$.
Take $\tilde{W}_{1}=X_{1}^{-1}\tilde{Y}_{1}$, then (42) can be written as:

$(I+\tilde{W}_{1})U^{m+1}=(I-\tilde{W}_{1})U^{m}+X_{1}^{-1}\tilde{F}_{1}^{m}.$
(43)

######  Theorem 9.

The Crank-Nicolson compact scheme (43) for convection diffusion equation (27)
is stable.

###### Proof.

Since $2\tilde{W}_{1}=W_{1}$, if $\tilde{\lambda}$ is an eigenvalue of
$\tilde{W}_{1}$, then $2\tilde{\lambda}$ is an eigenvalue of $W_{1}$. Now
using Theorem 7, we get that real part of $2\tilde{\lambda}$ is positive.
Since, $\tilde{\lambda}$ has positive real part, by following the line of
arguments in the proof of Theorem 6, we can show that the spectral radius of
the expansion matrix for (43) is less than 1. Thus the Crank-Nicolson compact
scheme (43) is stable. ∎

##  4\. Condition Number

In this section we obtain upper bound for the condition number of the matrix
$(I+W).$ First we bound $\left\lVert W\right\rVert_{2},$ where
$\left\lVert.\right\rVert_{2}$ is spectral norm of a matrix.

######  Theorem 10.

If $\delta z<2,$ then we have

$\left\lVert W\right\rVert_{2}=\mathcal{O}\left(\frac{\delta v}{{\delta
z}^{2}}\right).$

###### Proof.

It is enough to show that $\left\lVert X^{-1}\right\rVert_{2}\left\lVert
Y\right\rVert_{2}=\mathcal{O}\left(\frac{\delta v}{{\delta z}^{2}}\right)$.
The upper bound for $\left\lVert X^{-1}\right\rVert_{2}$ will be obtained
using the notable Gerschgorin’s circle theorem (GCT) [see for example, pp. 61,
[24]]. Let $X^{*}$ be the transpose of $X$ which
implies$(X^{-1})^{*}X^{-1}=(XX^{*})^{-1}$. Hence, following the definition of
singular value, $h^{2}$ is an eigenvalue of $(XX^{*})^{-1}$ if and only if $h$
is a singular value of $X^{-1}$. From (15), let us write the matrix $P=XX^{*}$
as follows:

$P_{l,l^{\prime}}=\left\\{\begin{array}[]{ll}p_{l}^{2}+q_{l}^{2}+r_{l}^{2}&\textrm{
if }1<l^{\prime}=l<N-1\\\ q_{1}^{2}+r_{1}^{2}&\textrm{ if }l^{\prime}=l=1\\\
q_{N-1}^{2}+r_{N-1}^{2}&\textrm{ if }l^{\prime}=l=N-1\\\ q_{l\wedge
l^{\prime}}p_{(l\wedge l^{\prime})+1}+r_{l\wedge l^{\prime}}q_{(l\wedge
l^{\prime})+1}&\textrm{ if }|l^{\prime}-l|=1\\\ r_{l\wedge
l^{\prime}}p_{(l\wedge l^{\prime})+2}&\textrm{ if }|l^{\prime}-l|=2\\\
0&\textrm{ else},\end{array}\right.$ (44)

for all $1\leq l,l^{\prime}\leq N-1$. Then $h^{-2}$ is an eigenvalue of the
positive definite matrix $P$. Thus

$\left\lVert X^{-1}\right\rVert_{2}=\max\\{h\mid h\textrm{ is a singular value
of }X^{-1}\\}=\frac{1}{\sqrt{\sigma_{min}(P)}},$

where $\sigma_{min}(P)$ is the smallest eigenvalue of $P$. Note that each
eigenvalue of $P$ is positive. We obtain the lower bound of eigenvalues of $P$
using GCT. To this end, we consider the following $N-1$ different
Gerschgorin’s discs to estimate the eigenvalues of matrix $P$:

$R_{1}\left(q_{1}^{2}+r_{1}^{2},\quad
q_{1}p_{2}+r_{1}(p_{3}+q_{2})\right),\quad
R_{2}\left(p_{2}^{2}+q_{2}^{2}+r_{2}^{2},\quad
p_{2}q_{1}+q_{2}(p_{3}+r_{1})+r_{2}(p_{4}+q_{3})\right),$
$R_{N-2}\left(p_{N-2}^{2}+q_{N-2}^{2}+r_{N-2}^{2},\quad
p_{N-2}(q_{N-3}+r_{N-4})+q_{N-2}(p_{N-1}+r_{N-3})+r_{N-2}q_{N-1}\right),$
$R_{N-1}\left(p_{N-1}^{2}+q_{N-1}^{2},\quad
p_{N-1}(q_{N-2}+r_{N-3})+q_{N-1}r_{N-2})\right),$
$R_{l}\left(p_{l}^{2}+q_{l}^{2}+r_{l}^{2},\quad
p_{l}(q_{l-1}+r_{l-2})+q_{l}(p_{l+1}+r_{l-1})+r_{l}(p_{l+2}+q_{l+1})\right),\>\>for\>\>3\leq
l\leq N-3.$

Here $R_{l}(g_{l},s_{l})$ $1\leq l\leq N-1$ denotes the $l^{th}$ disc with
center $g_{l}$ and radius $s_{l}$. If we set $r_{0}=p_{1}=r_{N-1}=p_{N}=0$
then for each $1\leq l\leq N-1$,

$g_{l}=p_{l}^{2}+q_{l}^{2}+r_{l}^{2},\quad\textrm{and }\quad
s_{l}=p_{l}(q_{l-1}+r_{l-2})+q_{l}(p_{l+1}+r_{l-1})+r_{l}(p_{l+2}+q_{l+1}).$

Clearly from (12), for sufficiently small $\delta z$, $(g_{l}-s_{l})$ is
positive and of order $\mathcal{O}({\delta v}^{-2})$. Then an application of
GCT gives

$\displaystyle\sigma_{min}(P)\geq\min_{1\leq l\leq N-1}(g_{l}-s_{l}).$

Therefore, $\left\lVert X^{-1}\right\rVert_{2}<\displaystyle\left(\min_{1\leq
l\leq N-1}(g_{l}-s_{l})\right)^{-1/2}=\mathcal{O}({\delta v}).$ Let
$n_{0}=l_{1}=n_{N-1}=l_{N}=0$ then we have

$\displaystyle\left\lVert Y\right\rVert_{\infty}$ $\displaystyle=\max_{1\leq
i\leq N-1}\sum_{j=1}^{N-1}\left\lvert Y_{i,j}\right\rvert=\max_{1\leq i\leq
N-1}\left(\left\lvert l_{i}\right\rvert+\left\lvert
m_{i}\right\rvert+\left\lvert n_{i}\right\rvert\right)$
$\displaystyle=\mathcal{O}\left(\frac{1}{\delta z^{2}}\right),$

and similarly,

$\displaystyle\left\lVert Y\right\rVert_{1}$ $\displaystyle=\max_{1\leq j\leq
N-1}\sum_{i=1}^{N-1}\left\lvert Y_{i,j}\right\rvert=\max_{1\leq j\leq
N-1}\left(\left\lvert n_{j-1}\right\rvert+\left\lvert
m_{j}\right\rvert+\left\lvert l_{j+1}\right\rvert\right)$
$\displaystyle=\mathcal{O}\left(\frac{1}{\delta z^{2}}\right).$

From the properties of subordinate matrix norm, $\left\lVert
Y\right\rVert_{2}\leq\sqrt{\left\lVert Y\right\rVert_{1}\left\lVert
Y\right\rVert_{\infty}}$. Then from above we obtain

$\left\lVert Y\right\rVert_{2}=\mathcal{O}\left(\frac{1}{\delta
z^{2}}\right).$

Therefore, we have

$\displaystyle\left\lVert X^{-1}\right\rVert_{2}\left\lVert Y\right\rVert_{2}$
$\displaystyle=\mathcal{O}\left(\frac{\delta v}{\delta z^{2}}\right).$

Hence the proof is complete. ∎

######  Theorem 11.

The condition number of the matrix $(I+W)$ is $\mathcal{O}\left(\frac{\delta
v}{{\delta z}^{2}}\right)$.

###### Proof.

Note that if $\lambda$ is an eigenvalue of $W$ then $\frac{1}{1+\lambda}$ is
an eigenvalue of $(I+W)^{-1}.$ Given that real part of $\lambda$ is positive,
we have $\left|\frac{1}{1+\lambda}\right|<1.$ It gives
$\left\lVert(I+W)^{-1}\right\rVert_{2}<1.$ Thus condition number of matrix
$(I+W)=\left\lVert(I+W)\right\rVert_{2}\left\lVert(I+W)^{-1}\right\rVert_{2}$
is upper bounded by $1+\left\lVert
W\right\rVert_{2}=\mathcal{O}\left(\frac{\delta v}{{\delta z}^{2}}\right)$,
using Lemma 10. Thus the condition number of $(I+W)$ is
$\mathcal{O}\left(\frac{\delta v}{{\delta z}^{2}}\right).$ ∎

######  Remark 12.

Similarly, the bound on the condition number of matrix $(I+\tilde{W})$ can be
derived. Moreover, for the constant coefficient case, we have already derived
the bound on condition number for the matrix $(I+\tilde{W}_{1})$ in [21].

##  5\. Numerical Illustrations

In this section, the stability of Backward Euler compact scheme (18) discussed
in Sec. 2.1 is numerically verified for special choices of $a(z)$, $b(z)$, and
$N$. Moreover, the general algorithm to compute $D^{1}_{N}$ is also presented.

######  Example 1.

(Stability of Backward Euler compact scheme (18) discussed in Sec. 2.1 for
$N=4$ case)

From Fig. 1, and Eq. (24), the expression of $D^{1}_{N}$ for $N=4$ is written
as follows:

$D_{4}^{1}=A_{2}B_{1}C_{3}+A_{1}B_{3}C_{2}-B_{1}B_{2}B_{3}.$ (45)

Particularly, we consider $a(z)=z+1$, $b(z)=(z+1)^{2}$, $z_{l}=0$, $z_{r}=1$,
and $\delta v=0.1$. Then using (12), (13), and (20), the zeros of Eq. (45) are
$28.3037$, $8.4060$, and $1.6272$. Therefore from Theorem 3, the backward
Euler compact scheme (18) for variable coefficient equation (1) is stable
under these specific choices.

######  Example 2.

(General algorithm to check the stability of backward Euler scheme (18))

Three MATLAB algorithms are provided in Appendix A. Algorithm A.1 computes the
symbolic expressions for $\beta^{{}^{\prime}}_{N}$ and $D^{1}_{N}$ given in
(22), and (24) respectively. Further, the function A.2 gives the
symbolic/numerical values for the roots of $D^{1}_{N}$. Finally, both of these
functions are called in another algorithm A.3 to conclude the stability of
Backward Euler compact scheme (18) under the specific choices of parameters
entered by user.

######  Remark 13.

A similar computation can be performed for the Crank-Nicolson compact scheme
(26) discussed in Sec. 2.2 in an analogous way, we have omitted the details
for the same. Moreover, for constant coefficient case discussed in Sec. 3, the
unconditional stability of Backward Euler and Crank-Nicolson compact schemes
are theoretically established in Theorems 8, and 9 respectively. Therefore,
the numerical verification for the constant coefficient case is also omitted.

##  6\. Conclusions and future directions

A novel difference equation based approach has been developed to investigate
the stability of variable coefficient convection diffusion equations using
compact schemes. It has been observed that the compact schemes for variable
coefficient problems are stable under certain restrictions. Further, the
unconditional stability of compact schemes for the constant coefficient case
has been proved theoretically. An estimate on the condition number of
amplification matrix has also been derived. The MATLAB algorithms have been
developed, and a supporting numerical example has been provided. As a future
scope, the proposed approach can be extended to the multi-dimensional
problems, system of PDEs, etc. This work finds an application to investigate
the stability of the variable coefficient PDEs which can not be transformed to
constant coefficient problems.

##  Appendix A MATLAB codes

### A.1.

It gives symbolic expressions of $\beta^{\prime}_{N}$ and $D^{1}_{N}$.

⬇

1function [BetaN,D1N]=betaandD(N)

2syms(sym(’A’,[1 N-1]));A=(sym(’A’,[1 N-1]));

3syms(sym(’B’,[1 N-1]));B=(sym(’B’,[1 N-1]));

4syms(sym(’C’,[1 N-1]));C=(sym(’C’,[1 N-1]));

5syms(sym(’P’,[1 N+1]));P=(sym(’P’,[1 N+1]));P(1)=1;

6 j=3;

7 for i=1:N-1

8 P(j)=(-C(j-2)/A(j-2))*P(j-2)+(-B(j-2)/A(j-2))*P(j-1);

9 j=j+1;

10 end

11 BetaN=P(N+1); % the expression of polynomial BetaN’

12 K=coeffs(P(N+1),P2);Pterm=1;

13 for s=1:N-1

14 Pterm=Pterm*A(s); % defining product term

15 end

16 D1N=simplifyFraction(Pterm*K(2)); % expression for D1N

17end

### A.2.

It gives the symbolic computation of the roots of $D^{1}_{N}$.

⬇

1

2function[u]=roots_value(N,p,q,r,l,m,n)

3syms lambda

4syms(sym(’A’,[1 N-1]));A=(sym(’A’,[1 N-1]));

5syms(sym(’B’,[1 N-1]));B=(sym(’B’,[1 N-1]));

6syms(sym(’C’,[1 N-1]));C=(sym(’C’,[1 N-1]));

7 for h=1:N-1

8 A(h)=n(h)-lambda*r(h);B(h)=m(h)-lambda*q(h);C(h)=l(h)-lambda*p(h);

9 end

10 syms(sym(’P’,[1 N+1]));P=(sym(’P’,[1 N+1]));P(1)=1;j=3;

11 for i=1:N-1

12 P(j)=(-C(j-2)/A(j-2))*P(j-2)+(-B(j-2)/A(j-2))*P(j-1);j=j+1;

13 end

14 K=coeffs(P(N+1),P2);Pterm=1; % product term to be multiplied

15 for s=1:N-1

16 Pterm=Pterm*A(s); % defining product term

17 end

18 D1N=simplifyFraction(Pterm*K(2)); % expression for D1N

19 M=coeffs(D1N,lambda);

20 V=zeros(N,1);

21 for b=1:N

22 V(b)=M(N-b+1); % Reverse order to apply MATLAB function “roots”

23 end

24 u=roots(V); % find roots of D1N

25end

### A.3.

It tests that roots of $D^{1}_{N}$ are positive or not?

⬇

1clc;clear;N=input(’Enter value of N=’);

2z_l=input(’Enter value of left space boundary point z_l=’);

3z_r=input(’Enter value of right space boundary point z_r=’);

4dv=input(’Enter value of time step dv=’);

5syms a(z) b(z) w1(z) y1(z) w2(z) y2(z)

6a(z)=input(’Enter the convection coefficient a(z)=’);

7b(z)=input(’Enter the diffusion coefficient b(z)=’);

8x=linspace(z_l,z_r,N+1);dz=x(2)-x(1); % discretizing space domain

9w=a(x);y=b(x);w1(z)=diff(a(z));y1(z)=diff(b(z));

10w2(z)=diff(w1(z));y2(z)=diff(y1(z)); % defining derivatives of a(z) and b(z)

11wd=w1(x);yd=y1(x);wdd=w2(x);ydd=y2(x);

12 gamma=zeros(N-1,1);xeta=zeros(N-1,1);alpha=zeros(N-1,1);

13 p=zeros(N-1,1);q=zeros(N-1,1);r=zeros(N-1,1); % defining p, q, r, l, m, and
n as per manuscript.

14 l=zeros(N-1,1);m=zeros(N-1,1);n=zeros(N-1,1);

15 for i=1:N-1

16 gamma(i)=(w(i+1)/y(i+1))+(2/y(i+1))*yd(i+1);

17
xeta(i)=w(i+1)-((dz^2)/12)*((w(i+1)/y(i+1))*wd(i+1)+(2/y(i+1))*wd(i+1)*yd(i+1)-wdd(i+1));

18
alpha(i)=y(i+1)+((dz^2)/12)*(((w(i+1)/y(i+1))*yd(i+1))-2*wd(i+1)+ydd(i+1)-((2/y(i+1))*(yd(i+1)^2))+((w(i+1)^2)/y(i+1)));

19 p(i)=((2+dz)*gamma(i))/(24*dv);

20 q(i)=5/(6*dv);

21 r(i)=((2-dz)*gamma(i))/(24*dv);

22 l(i)=(-xeta(i)/(2*dz))-(alpha(i)/(dz^2));

23 m(i)=(2*alpha(i))/(dz^2);

24 n(i)=(xeta(i)/(2*dz))-(alpha(i)/(dz^2));

25 end

26 [BetaN,D1N]=betaandD(N); % calling function to get expressions for
\beta^{’}_N and D^1_N.

27 fprintf(’Take \beta^{’}_1=P2, then expression for \beta^{’}_N is= %s. \n
’,BetaN);

28 fprintf(’The expression for D^1_N is= %s .\n ’,D1N);

29 u=roots_value(N,p,q,r,l,m,n); % to check the roots

30if min(u) > 0

31disp(”Result: Compact Scheme is stable for above parameter”)

32else

33 disp(”Result: The stability condition is not met for above parameter”)

34end

## References

  * [1] F. Black, M. Scholes, Pricing of options and corporate liabilities, The Journal of Political Economy 81 (1973) 637–654.
  * [2] P. J. Roache, Computational fluid dynamics, Hermosa, Albuquerque, NM (1976).
  * [3] J. Isenberg, C. Gutfinger, Heat transfer to a draining film, International Journal of Heat and Mass Transfer 16 (2) (1973) 505–512.
  * [4] J. Parlange, Water transport in soils, Annual review of fluid mechanics 12 (1) (1980) 77–102.
  * [5] Y. Achdou, O. Pironneau, Computational Methods for Options Pricing, SIAM, 2005.
  * [6] C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, 1998.
  * [7] M. Mehra, K. S. Patel, A. Shukla, Wavelet-optimized compact finite difference method for convection–diffusion equations, Int. J. Nonlinear Sci. Numer. 22 (2021) 353–372.
  * [8] J. W. Thomas, Numerical Partical Differential Equation: Finite Difference Methods, Springer–Verlog, 1995.
  * [9] D. Tavella, C. Randall, Pricing Financial Instruments: The Finite Difference Method, John Wiley & Sons, 2000.
  * [10] J. C. Strikewerda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2004.
  * [11] A. Rigal, High-order difference scheme for unsteady one dimensional diffusion convection problems, J. Comput. Phys. 114 (1994) 59–76.
  * [12] W. F. Spotz, G. F. Carey, High-order compact scheme for the steady stream-function vorticity equation, International journal for numerical methods in Engineering 38 (1995) 3497–3512.
  * [13] V. Kumar, High-order compact finite-difference scheme for singularly-perturbed reaction-diffusion problems on a new mesh of Shishkin type, J. Optim. Theory. Appl. 143 (2009) 123–147.
  * [14] K. S. Patel, M. Mehra, A numerical study of Asian option with high-order compact finite difference scheme, J. Appl. Math. Comput. 57 (2018) 467–491.
  * [15] K. S. Patel, M. Mehra, High-order compact finite difference scheme for pricing Asian option with moving boundary condition, Differ. Equ. Dyn. Syst. 27 (2019) 39–56.
  * [16] B. During, C. Heuer, High-order compact schemes for parabolic problems with mixed derivatives in multiple space dimensions., SIAM J. Numer. Anal. 53 (2015) 2113–2124.
  * [17] J. C. Kalita, D. Dalal, A. K. Dass, A class of higher order compact schemes for the unsteady two-dimensional convection–diffusion equation with variable convection coefficients, International Journal for Numerical Methods in Fluids 38 (12) (2002) 1111–1131.
  * [18] W. F. Spotz, G. F. Carey, Extension of high‐-order compact schemes to time‐-dependent problems, Numer. Methods Partial Differential Equations 17 (2001) 657–672.
  * [19] L. N. Trefthen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, Cornell University, 1996.
  * [20] A. Mohebbi, M. Dehghan, High-order compact solution of the one-dimensional heat and advection–diffusion equations, Applied mathematical modelling 34 (10) (2010) 3071–3084.
  * [21] A. Goswami, K. S. Patel, Stability analysis and condition number of compact scheme for convection-diffusion equations, arXiv:2201.05854.
  * [22] W. F. Spotz, High-order compact finite difference schemes for computational mechanics, The University of Texas at Austin, 1995.
  * [23] S. N. Elaydi, An introduction to difference equations, Springer, 2005.
  * [24] G. D. Smith, Numerical Solution of Partial Differential Equations, Oxford University Press,, 1978.