A parameter uniform essentially first order convergent numerical method for a parabolic singularly perturbed differential equation of reaction-diffusion type with initial and Robin boundary conditions

TL;DR

This paper introduces a finite difference method on a Shishkin mesh for reaction-diffusion parabolic equations with boundary layers, achieving uniform first-order convergence in space and time regardless of perturbation parameters.

Contribution

The paper proposes a novel numerical scheme that attains uniform essentially first-order convergence for a class of singularly perturbed parabolic equations with boundary layers.

Findings

Method is uniformly convergent in maximum norm.

Achieves first-order accuracy in time and essentially first-order in space.

Effective for reaction-diffusion equations with boundary layers.

Abstract

In this paper, a class of linear parabolic singularly perturbed second order differential equations of reaction-diffusion type with initial and Robin boundary conditions is considered. The solution u of this equation is smooth, whereas the first derivative in the space variable exhibits parabolic boundary layers. A numerical method composed of a classical finite difference scheme on a piecewise uniform Shishkin mesh is suggested. This method is proved to be first order convergent in time and essentially first order convergent in the space variable in the maximum norm uniformly in the perturbation parameters.

Tables

Table 1. Table 1: Values of    D ε N , D N , p N , p ∗ ​  and  ​ C p ∗ N ​  for  ​ α = 0.9 ​ and ​ N = 64 . superscript subscript 𝐷 𝜀 𝑁 superscript 𝐷 𝑁 superscript 𝑝 𝑁 superscript 𝑝  and  superscript subscript 𝐶 superscript 𝑝 𝑁  for  𝛼 0.9 and 𝑁 64 D_{\varepsilon}^{N},D^{N},p^{N},p^{*}\text{ and }C_{p^{*}}^{N}\text{ for }\;\alpha=0.9\;\text{and}\;N=64.

$\epsilon$	Number of mesh points $N$
	32	64	128	256
$2^{-6}$	0.231E-01	0.122E-01	0.638E-02	0.330E-02
$2^{-8}$	0.253E-01	0.131E-01	0.669E-02	0.338E-02
$2^{-10}$	0.263E-01	0.134E-01	0.675E-02	0.339E-02
$2^{-12}$	0.265E-01	0.134E-01	0.676E-02	0.339E-02
$2^{-14}$	0.266E-01	0.134E-01	0.676E-02	0.339E-02
$D^N$	0.266E-01	0.134E-01	0.676E-02	0.339E-02
$p^N$	0.983E+00	0.991E+00	0.996E+00
$C_{p^{\ast }}^N$	0.162E+01	0.162E+01	0.161E+01	0.160E+01
Computed $t$-order of $\epsilon -$uniform convergence, $p^{\ast }$ = 0.9827155
Computed $\epsilon -$uniform error constant, $C_{p^{\ast }}^{\ast }$ = 1.620163

Table 2. Table 2: Values of    D ε N , D N , p N , p ∗ ​  and  ​ C p ∗ N ​  for  ​ α = 0.9 ​ and ​ M = 256 . superscript subscript 𝐷 𝜀 𝑁 superscript 𝐷 𝑁 superscript 𝑝 𝑁 superscript 𝑝  and  superscript subscript 𝐶 superscript 𝑝 𝑁  for  𝛼 0.9 and 𝑀 256 D_{\varepsilon}^{N},D^{N},p^{N},p^{*}\text{ and }C_{p^{*}}^{N}\text{ for }\;\alpha=0.9\;\text{and}\;M=256.

$\epsilon$	Number of mesh points $N$
	32	64	128	256
$2^{-6}$	0.964E-02	0.399E-02	0.139E-02	0.496E-03
$2^{-8}$	0.119E-01	0.508E-02	0.195E-02	0.560E-03
$2^{-10}$	0.117E-01	0.617E-02	0.258E-02	0.843E-03
$2^{-12}$	0.537E-02	0.298E-02	0.155E-02	0.698E-03
$2^{-14}$	0.272E-02	0.150E-02	0.771E-03	0.290E-03
$D^N$	0.119E-01	0.617E-02	0.258E-02	0.843E-03
$p^N$	0.946E+00	0.126E+01	0.162E+01
$C_{p^{\ast }}^N$	0.655E+00	0.655E+00	0.528E+00	0.332E+00
Computed $x$-order of $\epsilon -$uniform convergence, $p^{\ast }$ = 0.9456793
Computed $\epsilon -$uniform error constant, $C_{p^{\ast }}^{\ast }$ = 0.6552203

Equations

$$\frac{\partial u}{\partial t}(x,t)-\epsilon \frac{{\partial }^{2}u}{\partial x^2}(x,t)+a(x,t)u(x,t)=f(x,t),\text{ on }\Omega ,$$

$$\begin{array}[]{l}u(0,t)-\dfrac{\partial u}{\partial x}(0,t)=\phi_{L}(t),\;\;u(1,t)+\dfrac{\partial u}{\partial x}(1,t)=\phi_{R}(t),\;\;0\leq t\leq T,\\ \hskip 85.35826ptu(x,0)=\phi_{B}(x),\;\;0\leq x\leq 1,\end{array}$$

$$Lu=f\text{ on }\Omega ,$$

$${\beta }_{0}u(0,t)={\varphi }_L(t),{\beta }_{1}u(1,t)={\varphi }_R(t),u(x,0)={\varphi }_B(x),$$

$$L=\frac{\partial }{\partial t}-\epsilon \frac{{\partial }^2}{\partial x^2}+a,{\beta }_0=I-\frac{\partial }{\partial x},{\beta }_1=I+\frac{\partial }{\partial x}$$

$$\frac{\partial u_0}{\partial t}+a{u}_0=f,\text{ on }\Omega ,u_0=u\text{ on }{\Gamma }_B.$$

$$0<\alpha <a(x,t).$$

$${\varphi }_B(0)={\varphi }_L(0)+\frac{d{\varphi }_B}{dx}(0)\text{and}{\varphi }_B(1)={\varphi }_R(0)-\frac{d{\varphi }_B}{dx}(1),$$

$$\begin{array}[]{lcl}\displaystyle\frac{d\phi_{L}}{dt}(0)&=&-\varepsilon\displaystyle\frac{d^{3}\phi_{B}}{dx^{3}}(0)+\varepsilon\displaystyle\frac{d^{2}\phi_{B}}{dx^{2}}(0)+a(0,0)\displaystyle\frac{d\phi_{B}}{dx}(0)-[a(0,0)-\dfrac{\partial a}{\partial x}(0,0)]\phi_{B}(0)\\ &&+f(0,0)-\dfrac{\partial f}{\partial x}(0,0),\\ \displaystyle\frac{d\phi_{R}}{dt}(0)&=&\varepsilon\displaystyle\frac{d^{3}\phi_{B}}{dx^{3}}(1)+\varepsilon\displaystyle\frac{d^{2}\phi_{B}}{dx^{2}}(1)-a(1,0)\displaystyle\frac{d\phi_{B}}{dx}(1)-[a(1,0)+\dfrac{\partial a}{\partial x}(1,0)]\phi_{B}(1)\\ &&+f(1,0)+\dfrac{\partial f}{\partial x}(1,0),\end{array}$$

$$\begin{array}[]{lcl}\dfrac{d^{2}\phi_{L}}{dt^{2}}(0)&=&-\varepsilon^{2}\dfrac{d^{5}\phi_{B}}{dx^{5}}(0)+\varepsilon^{2}\dfrac{d^{4}\phi_{B}}{dx^{4}}(0)+2\varepsilon a(0,0)\dfrac{d^{3}\phi_{B}}{dx^{3}}(0)+[-2\varepsilon a(0,0)\\ &&+4\varepsilon\dfrac{\partial a}{\partial x}(0,0)]\dfrac{d^{2}\phi_{B}}{dx^{2}}(0)+[-2\varepsilon\dfrac{\partial a}{\partial x}(0,0)+3\varepsilon\dfrac{\partial^{2}a}{\partial x^{2}}(0,0)-a^{2}(0,0)\\ &&+\dfrac{\partial a}{\partial t}(0,0)]\dfrac{d\phi_{B}}{dx}(0)+[-\varepsilon\dfrac{\partial^{2}a}{\partial x^{2}}(0,0)+a^{2}(0,0)-\dfrac{\partial a}{\partial t}(0,0)+\varepsilon\dfrac{\partial^{3}a}{\partial x^{3}}(0,0)\\ &&-2a(0,0)\dfrac{\partial a}{\partial x}(0,0)+\dfrac{\partial^{2}a}{\partial x\partial t}(0,0)]\phi_{B}(0)+[-a(0,0)+\dfrac{\partial a}{\partial x}(0,0)]f(0,0)\\ &&+\dfrac{\partial f}{\partial t}(0,0)-\varepsilon\dfrac{\partial^{3}f}{\partial x^{3}}(0,0)+\varepsilon\dfrac{\partial^{2}f}{\partial x^{2}}(0,0)+a(0,0)\dfrac{\partial f}{\partial x}(0,0)-\dfrac{\partial^{2}f}{\partial x\partial t}(0,0),\end{array}$$

$$\begin{array}[]{lcl}\dfrac{d^{2}\phi_{R}}{dt^{2}}(0)&=&\varepsilon^{2}\dfrac{d^{5}\phi_{B}}{dx^{5}}(1)+\varepsilon^{2}\dfrac{d^{4}\phi_{B}}{dx^{4}}(1)-2\varepsilon a(1,0)\dfrac{d^{3}\phi_{B}}{dx^{3}}(1)+[-2\varepsilon a(1,0)\\ &&-4\varepsilon\dfrac{\partial a}{\partial x}(1,0)]\dfrac{d^{2}\phi_{B}}{dx^{2}}(1)+[-2\varepsilon\dfrac{\partial a}{\partial x}(1,0)-3\varepsilon\dfrac{\partial^{2}a}{\partial x^{2}}(1,0)+a^{2}(1,0)\\ &&-\dfrac{\partial a}{\partial t}(1,0)]\dfrac{d\phi_{B}}{dx}(1)+[-\varepsilon\dfrac{\partial^{2}a}{\partial x^{2}}(1,0)+a^{2}(1,0)-\dfrac{\partial a}{\partial t}(1,0)-\varepsilon\dfrac{\partial^{3}a}{\partial x^{3}}(1,0)\\ &&+2a(1,0)\dfrac{\partial a}{\partial x}(1,0)-\dfrac{\partial^{2}a}{\partial x\partial t}(1,0)]\phi_{B}(1)+[-a(1,0)-\dfrac{\partial a}{\partial x}(1,0)]f(1,0)\\ &&+\dfrac{\partial f}{\partial t}(1,0)+\varepsilon\dfrac{\partial^{3}f}{\partial x^{3}}(1,0)+\varepsilon\dfrac{\partial^{2}f}{\partial x^{2}}(1,0)-a(1,0)\dfrac{\partial f}{\partial x}(1,0)+\dfrac{\partial^{2}f}{\partial x\partial t}(1,0).\end{array}$$

$$L\psi (x^{\ast },t^{\ast })=\frac{\partial \psi }{\partial t}(x^{\ast },t^{\ast })-\epsilon \frac{{\partial }^2\psi }{\partial x^2}(x^{\ast },t^{\ast })+a(x^{\ast },t^{\ast })\psi (x^{\ast },t^{\ast })<0,$$

$$|\psi (x,t)|\le \max \left\{\parallel {\beta }_0\psi (0,t)\parallel ,\parallel {\beta }_1\psi (1,t)\parallel ,\parallel \psi (x,0)\parallel ,\frac{1}{\alpha }\parallel L\psi \parallel \right\}.$$

$$\begin{array}[]{rcl}\theta^{\pm}(x,t)=\max\displaystyle\left\{\parallel\beta_{0}\psi(0,t)\parallel,\parallel\beta_{1}\psi(1,t)\parallel,\parallel\psi(x,0)\parallel,\dfrac{1}{\alpha}\parallel L\psi\parallel\right\}\pm\;\;\psi(x,t),\;\;(x,t)\in\bar{\Omega}.\end{array}$$

$$\begin{array}[]{lcl}|u(x,t)|&\leq&C(\parallel\phi_{L}(t)\parallel+\parallel\phi_{R}(t)\parallel+\parallel\phi_{B}(x)\parallel+\parallel f\parallel),\\ |\dfrac{\partial^{l}u}{\partial t^{l}}(x,t)|&\leq&C(\parallel u\parallel+\displaystyle\sum_{q=0}^{l}\parallel\dfrac{\partial^{q}f}{\partial t^{q}}\parallel),\;\;l=1,2,\\ |\dfrac{\partial u}{\partial x}(x,t)|&\leq&C(\parallel u\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel),\\ |\dfrac{\partial^{l}u}{\partial x^{l}}(x,t)|&\leq&C\varepsilon^{\frac{-(l-1)}{2}}(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel),\;\;l=2,3,\\ |\dfrac{\partial^{4}u}{\partial x^{4}}(x,t)|&\leq&C\varepsilon^{\frac{-3}{2}}(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial^{2}f}{\partial t^{2}}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel+\parallel\dfrac{\partial^{3}f}{\partial x\partial t^{2}}\parallel\\ &&+\varepsilon^{\frac{1}{2}}\parallel\dfrac{\partial^{2}f}{\partial x^{2}}\parallel),\\ |\dfrac{\partial^{2}u}{\partial x\partial t}(x,t)|&\leq&C(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel),\\ |\dfrac{\partial^{3}u}{\partial x^{2}\partial t}(x,t)|&\leq&C\varepsilon^{\frac{-1}{2}}(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial^{2}f}{\partial t^{2}}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel+\parallel\dfrac{\partial^{3}f}{\partial x\partial t^{2}}\parallel).\end{array}$$

$$\frac{{\partial }^{2}u}{\partial x\partial t}-\epsilon \frac{{\partial }^{3}u}{\partial x^3}+a\frac{\partial u}{\partial x}=\frac{\partial f}{\partial x}-\frac{\partial a}{\partial x}u,$$

$$\begin{array}[]{lcl}\dfrac{\partial u}{\partial x}(0,t)=u(0,t)-\phi_{L}(t),\;\;\dfrac{\partial u}{\partial x}(1,t)=\phi_{R}(t)-u(1,t),\;\;\dfrac{\partial u}{\partial x}(x,0)=\dfrac{\partial\phi_{B}(x)}{\partial x}=\eta(x).\end{array}$$

$$\frac{\partial z}{\partial t}-\epsilon \frac{{\partial }^{2}z}{\partial x^2}+az=h,$$

$$z(0,t)=u(0,t)-{\varphi }_L(t),z(1,t)={\varphi }_R(t)-u(1,t),z(x,0)=\eta (x),$$

$$\begin{array}[]{lcl}|\dfrac{\partial z}{\partial t}(x,t)|\leq C(\parallel z\parallel+\parallel h\parallel+\parallel\dfrac{\partial h}{\partial t}\parallel),\;\;|\dfrac{\partial^{2}z}{\partial t^{2}}(x,t)|&\leq C(\parallel z\parallel+\parallel h\parallel+\parallel\dfrac{\partial h}{\partial t}\parallel+\parallel\dfrac{\partial^{2}h}{\partial t^{2}}\parallel).\end{array}$$

$$|\frac{\partial z}{\partial x}(y,t)|\le C{\epsilon }^{\frac{-1}{2}}\parallel z\parallel .$$

$$u=v+w,$$

$$Lv=f,\text{ on }\Omega$$

$${\beta }_{0}v(0,t)={\beta }_0{u}_0(0,t),{\beta }_{1}v(1,t)={\beta }_1{u}_0(1,t),v(x,0)=u_0(x,0)$$

$$Lw=0,\text{ on }\Omega$$

$$\begin{array}[]{rcl}\beta_{0}{w}(0,t)=\beta_{0}({u}-{v})(0,t),\;\;\beta_{1}{w}(1,t)=\beta_{1}({u}-{v})(1,t),\;\;w(x,0)=0.\end{array}$$

$$\begin{array}[]{rclc}|\dfrac{\partial^{l}v}{\partial t^{l}}(x,t)|&\leq&C,\;\;l=0,1,2,&|\dfrac{\partial^{l}v}{\partial x^{l}}(x,t)|\leq C,\;\;l=1,2,3\\ |\dfrac{\partial^{4}v}{\partial x^{4}}(x,t)|&\leq&C\varepsilon^{-1/2},&|\dfrac{\partial^{l+1}v}{\partial x^{l}\partial t}(x,t)|\leq C,\;\;l=1,2.\end{array}$$

$$v=u_0+\epsilon v_1.$$

$$\frac{\partial v_1}{\partial t}-\epsilon \frac{{\partial }^2{v}_1}{\partial x^2}+a{v}_1=\frac{{\partial }^2{u}_0}{\partial x^2}$$

$$(v_1-\frac{\partial v_1}{\partial x})(0,t)=0,(v_1+\frac{\partial v_1}{\partial x})(1,t)=0,v_1(x,0)=0,$$

Full text

1 Introduction

A differential equation in which small parameters multiply the highest order derivative and some or none of the lower order derivatives is known as a singularly perturbed differential equation. In this paper, a class of linear parabolic singularly perturbed second order differential equation of reaction-diffusion type with initial and Robin boundary conditions is considered.

For a general introduction to parameter-uniform numerical methods for singular perturbation problems, see [1], [2], [8] and [9]. In [3], a Dirichlet boundary value problem for a linear parabolic singularly perturbed differential equation is studied and a numerical method comprising of a standard finite difference operator on a fitted piecewise uniform mesh is considered and it is proved to be uniform with respect to the small parameter in the maximum norm. In [4], a boundary-value problem for a singularly perturbed parabolic PDE with convection is considered on an interval in the case of the singularly perturbed Robin boundary condition is considered and using a defect correction technique, an ε-uniformly convergent schemes of high-order time-accuracy is constructed. The efficiency of the new defect-correction schemes is confirmed by numerical experiments. In [5], a one-dimensional steady-state convection dominated convection-diffusion problem with Robin boundary conditions is considered and the numerical solutions obtained using an upwind finite difference scheme on Shishkin meshes are uniformly convergent with respect to the diffusion cofficient.

Consider the following parabolic initial-boundary value problem for a singularly perturbed linear second order differential equation

[TABLE]

with

[TABLE]

where $\Omega=\{(x,t):0<x<1,\;0<t\leq T\},\;\;\bar{\Omega}=\Omega\cup\Gamma,\;\;\Gamma=\Gamma_{L}\cup\Gamma_{B}\cup\Gamma_{R}$ with $\Gamma_{L}=\{(0,t):0\leq t\leq T\},\;\Gamma_{R}=\{(1,t):0\leq t\leq T\}\text{ and }\Gamma_{B}=\{(x,0):0<x<1\}.$

The problem (1), (2) can also be written in the operator form

[TABLE]

[TABLE]

where the operators $L,\beta_{0},\beta_{1}$ are defined by

[TABLE]

where $I$ is the identity operator. The reduced problem corresponding to (1), (2) is defined by

[TABLE]

The problem (1), (2) is said to be singularly perturbed in the following sense.

The solution $u$ of (1), (2) is expected to exhibit weak twin layers of width $O(\sqrt{\varepsilon})$ at $x=0$ and $x=1.$

2 Solution of the continuous problem

Standard theoretical results on the existence of the solution of (1), (2) are stated, without proof, in this section. See [6] and [7] for more details. For all $(x,t)\in\bar{\Omega},$ it is assumed that $a(x,t)$ satisfies the condition

[TABLE]

Sufficient conditions for the existence, uniqueness and regularity of solution of (1), (2) are given in the following theorem.

Theorem 2.1

Assume that $a\text{ and }f$ are sufficiently smooth. Also assume that $\phi_{L}\in C^{2}(\Gamma_{L}),\;\phi_{B}\in C^{5}(\Gamma_{B})$,  $\phi_{R}\in C^{2}(\Gamma_{R})$ and that the following compatibility conditions are fulfilled at the corners $(0,0)$ and $(1,0)$ of  $\Gamma.$

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

Then there exists a unique solution $u$ of (1), (2) satisfying $u\in C_{\lambda}^{(4)}(\bar{\Omega})$.

3 Analytical results

The operator $L$ satisfies the following maximum principle:

Lemma 3.1

Let the assumptions (4) - (8) hold. Let $\psi$ be any function in the domain of $L$ such that $\beta_{0}\psi(0,t)\geq 0,\;\beta_{1}\psi(1,t)\geq 0,\;\psi(x,0)\geq 0.$ Then  $L\psi(x,t)\geq 0$  on  $\Omega$  implies that $\psi(x,t)\geq 0$  on  $\bar{\Omega}.$

Proof. Let $x^{*},t^{*}$ be such that $\psi(x^{*},t^{*})=\displaystyle\min_{\bar{\Omega}}\psi(x,t)$ and assume that the lemma is false. Then $\psi(x^{*},t^{*})<0.$ For $x^{*}=0,$ $\beta_{0}\psi(0,t^{*})=\psi(0,t^{*})-\dfrac{\partial\psi}{\partial x}(0,t^{*})<0,\;\;$for $x^{*}=1,\beta_{1}\psi(1,t^{*})=\psi(1,t^{*})+\dfrac{\partial\psi}{\partial x}(1,t^{*})<0\;\;$and for $t^{*}=0,\;\;\psi(x^{*},0)<0,$ contradicting the hypotheses. Therefore, $(x^{*},t^{*})\notin\Gamma.$ Let $(x^{*},t^{*})\in\Omega.$ Then $\dfrac{\partial\psi}{\partial t}(x^{*},t^{*})=0\text{ and }\dfrac{\partial^{2}\psi}{\partial x^{2}}(x^{*},t^{*})\geq 0$ lead to

[TABLE]

which contradicts the assumption and proves the result for $L.$

Lemma 3.2

Let the assumptions (4) - (8) hold. If $\psi$ is any function in the domain of $L,$ then,  for each  $(x,t)\in\bar{\Omega},$

[TABLE]

Proof. Define the two functions

[TABLE]

It is not hard to verify that $\beta_{0}{\theta}^{\pm}(0,t)\geq 0,\;\beta_{1}{\theta}^{\pm}(1,t)\geq 0,\;{\theta}^{\pm}(x,0)\geq 0$ and $L{\theta}^{\pm}\geq 0$ on $\Omega.$ It follows from Lemma 3.1 that ${\theta}^{\pm}\geq 0$ on $\bar{\Omega}$ as required.

A standard estimate of the solution $u$ of the problem (1), (2) and its derivatives is contained in the following lemma.

Lemma 3.3

Let the assumptions (4) - (8) hold and let $u$ be the solution of the problem (1), (2). Then,  for all $(x,t)\in\bar{\Omega}$

[TABLE]

Proof. The bound on $u$ is an immediate consequence of Lemma 3.2.

Differentiating (1) partially with respect to $`t`$ once and twice respectively, and applying Lemma 3.1, the bounds on $\dfrac{\partial u}{\partial t}$ and $\dfrac{\partial^{2}{u}}{\partial t^{2}}$ respectively are derived. Now, differentiating (1) partially with respect to $`x`$ once, gives

[TABLE]

and from the initial and boundary conditions, we derive

[TABLE]

Denoting $\dfrac{\partial u}{\partial x}$ by $z$ in (9) and (10), we get,

[TABLE]

[TABLE]

where $h=\dfrac{\partial f}{\partial x}-\dfrac{\partial a}{\partial x}u.$

This problem (11), (12) is similar to the problem in [3]. Now, using the stability result in [3], the bound on $z$ or $\dfrac{\partial u}{\partial x}$ is determined. Thus,

$|\dfrac{\partial u}{\partial x}(x,t)|\leq C(\parallel u\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel).$

On differentiating (11) partially with respect to $`t`$ once and twice respectively, and applying the stability in [3], the following bounds on $\dfrac{\partial z}{\partial t}$ or $\dfrac{\partial^{2}{u}}{\partial x\partial t}$ and $\dfrac{\partial^{2}{z}}{\partial t^{2}}$ respectively are derived

[TABLE]

To bound $\dfrac{\partial z}{\partial x},$ for each $(x,t)$, consider an interval $I=[a,a+\sqrt{\varepsilon}],\;a\geq 0$ such that $x\in I.$ Then for some $y$ such that $a<y<a+\sqrt{\varepsilon}$  and  $t\in(0,T],$

$\dfrac{\partial z}{\partial x}(y,t)=\dfrac{z(a+\sqrt{\varepsilon},t)-z(a,t)}{\sqrt{\varepsilon}}.$

Therefore,

[TABLE]

Then, for any $x\in I,$

$\dfrac{\partial z}{\partial x}(x,t)=\dfrac{\partial z}{\partial x}(y,t)+\displaystyle\int^{x}_{y}\dfrac{\partial^{2}z(s,t)}{\partial x^{2}}ds$

$\dfrac{\partial z}{\partial x}(x,t)=\dfrac{\partial z}{\partial x}(y,t)+\varepsilon^{-1}\displaystyle\int^{x}_{y}\left(\dfrac{\partial z(s,t)}{\partial t}-h(s,t)+a(s,t)z(s,t)\right)ds.$

Therefore,

$|\dfrac{\partial z}{\partial x}(x,t)|\leq|\dfrac{\partial z}{\partial x}(y,t)|+C\varepsilon^{-1}\displaystyle\int^{x}_{y}(\parallel z\parallel+\parallel h\parallel+\parallel\dfrac{\partial h}{\partial t}\parallel)ds.$

Using (13) in the above equation, yields

$|\dfrac{\partial z}{\partial x}(x,t)|\leq C\varepsilon^{\frac{-1}{2}}(\parallel z\parallel+\parallel h\parallel+\parallel\dfrac{\partial h}{\partial t}\parallel).$

i.e. $|\dfrac{\partial^{2}u}{\partial x^{2}}(x,t)|\leq C\varepsilon^{\frac{-1}{2}}(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel).$

Rearranging the terms in (11), we get

$|\dfrac{\partial^{2}z}{\partial x^{2}}(x,t)|\leq C\varepsilon^{-1}(\parallel z\parallel+\parallel h\parallel+\parallel\dfrac{\partial h}{\partial t}\parallel).$

i.e. $|\dfrac{\partial^{3}u}{\partial x^{3}}(x,t)|\leq C\varepsilon^{-1}(\parallel u\parallel+\parallel f\parallel+\parallel\dfrac{\partial f}{\partial t}\parallel+\parallel\dfrac{\partial f}{\partial x}\parallel+\parallel\dfrac{\partial^{2}f}{\partial x\partial t}\parallel).$

Following the steps similar to those used to bound $\dfrac{\partial z}{\partial x},$ the bound of the mixed derivative $\dfrac{\partial^{2}z}{\partial x\partial t}$ or $\dfrac{\partial^{3}u}{\partial x^{2}\partial t}$ is also derived.

Differentiating (11) once partially with respect to $`x`$ and rearranging the equation, the bound on $\dfrac{\partial^{3}z}{\partial x^{3}}$ or $\dfrac{\partial^{4}u}{\partial x^{4}}$ follows.

The Shishkin decomposition of the solution $u$ of the problem (1), (2) is

[TABLE]

where the smooth component $v$ of the solution $u$ satisfies

[TABLE]

with

[TABLE]

and the singular component $w$ of the solution $u$ satisfies

[TABLE]

with

[TABLE]

Bounds on the smooth component $v$ of $u$ and its derivatives are contained in

Lemma 3.4

Let the assumptions (4) - (8) hold. Then there exists a constant $C,$ such that, for each $(x,t)\in\bar{\Omega},$

[TABLE]

Proof. The smooth component $v$ is subjected to further decomposition

[TABLE]

The component $v_{1}$ satisfies the following equation:

[TABLE]

with

[TABLE]

where $u_{0}$ is the solution of the reduced problem (3).

From the expressions (19), (20) and using Lemma (3.3), it is found that for $k=1,2,3,4,\;\;l=0,1,2,\;\;m=1,2$

[TABLE]

From (18) and (21), the following bounds hold:

[TABLE]

The layer functions $B^{L},\;B^{R},\;B,$ associated with the solution $u,$ are defined on $\bar{\Omega}$ by $B^{L}(x)=e^{-x\sqrt{\alpha/\varepsilon}},\;B^{R}(x)=B^{L}(1-x),\;B(x)=B^{L}(x)+B^{R}(x).$

The following elementary properties of these layer functions, for all $0\leq x\leq y\leq 1,$ should be noted:

$B(x)=B(1-x),\;B^{L}(x)>B^{L}(y),\;0\;<\;B^{L}(x)\leq\;1,\;B^{R}(x)<B^{R}(y),\;0\;<\;B^{R}(x)\leq\;1,\;B(x)\;\text{\;is monotonically decreasing for increasing}\;\;x\in[0,\frac{1}{2}],\\ B(x)\;\text{\;is monotonically increasing for increasing}\;\;x\in[\frac{1}{2},1],\\ B(x)\leq 2B^{L}(x)\;\text{for}\;x\in[0,\frac{1}{2}],\;\;B(x)\leq 2B^{R}(x)\;\text{for}\;x\in[\frac{1}{2},1],\;B^{L}(2\frac{\sqrt{\varepsilon}}{\sqrt{\alpha}}\ln N)=N^{-2}.$

Bounds on the singular component $w$ of $u$ and its derivatives are contained in

Lemma 3.5

Let the assumptions (4) - (8) hold. Then there exists a constant $C,$ such that, for each $(x,t)\in\bar{\Omega},$

[TABLE]

Proof. To derive the bound of $w$, define two functions

${\psi}^{\pm}(x,t)=Ce^{\alpha t}B(x)\;\pm\;w(x,t)$.

For a proper choice of $C,$  $\beta_{0}\psi^{\pm}(0,t)\geq 0,\;\beta_{1}\psi^{\pm}(1,t)\geq 0$ and ${\psi^{\pm}}(x,0)\geq 0$ and for $(x,t)\in\Omega$,

[TABLE]

By Lemma 3.1, $\psi^{\pm}\geq 0$ on $\bar{\Omega}$ and it follows that

[TABLE]

Differentiating (16) partially with respect to $`t`$ once and twice, and using Lemma 3.1, it is not hard to see that

[TABLE]

Differentiating (16) with respect to $`x`$ once,

[TABLE]

And from the initial and boundary conditions,

[TABLE]

Denoting $\dfrac{\partial w}{\partial x}$ by $\tilde{z}$ in (22) and (23), yields

[TABLE]

[TABLE]

This problem (24), (25) is similar to the problem in [3]. Now, using the stability result in [3], the bound on $\tilde{z}$ or $\dfrac{\partial w}{\partial x}$ is determined. Thus,

[TABLE]

Differentiating (24) partially with respect to $`t`$ once and twice respectively, and applying the stability in [3], the following bounds on $\dfrac{\partial\tilde{z}}{\partial t}$ or $\dfrac{\partial^{2}w}{\partial x\partial t}$ and $\dfrac{\partial^{2}\tilde{z}}{\partial t^{2}}$ respectively are derived.

[TABLE]

To bound $\dfrac{\partial^{2}w}{\partial x^{2}}$ or $\dfrac{\partial\tilde{z}}{\partial x},$ consider an interval $I=[c,c+\sqrt{\varepsilon}],\;c\geq 0$ such that $x\in I.$ Then for some $y$ such that $c<y<c+\sqrt{\varepsilon}$  and  $t\in(0,T],$

[TABLE]

Therefore,

[TABLE]

Hence,

[TABLE]

Then, for any $x\in I,$ such that $y<\eta<x$

[TABLE]

By using the bounds for $\tilde{z},\;\dfrac{\partial\tilde{z}}{\partial t},\;w$ and (26) in the above equation yields

[TABLE]

Therefore,

[TABLE]

Rearranging the equation (24) yields

[TABLE]

Using the bounds of $w,$ $\tilde{z}$ and $\dfrac{\partial\tilde{z}}{\partial t}$ in the above equation, the following bound holds.

[TABLE]

Differentiating (24) with respect to $`t`$ once and following a similar procedure to bound $\dfrac{\partial\tilde{z}}{\partial x},$ the bound of the mixed derivative $\dfrac{\partial^{2}\tilde{z}}{\partial x\partial t}$ or $\dfrac{\partial^{3}w}{\partial x^{2}\partial t}$ is derived.

[TABLE]

On differentiating (24) with respect to $x$ and rearranging yields

[TABLE]

4 The Shishkin mesh

A piecewise uniform Shishkin mesh is now constructed. Let $\Omega^{M}_{t}=\{t_{k}\}_{k=1}^{M},\;\Omega^{N}_{x}=\{x_{j}\}_{j=1}^{N-1},\;\bar{\Omega}^{M}_{t}=\{t_{k}\}_{k=0}^{M},\;\bar{\Omega}^{N}_{x}=\{x_{j}\}_{j=0}^{N},\;\Omega^{N,M}=\Omega^{N}_{x}\times\Omega^{M}_{t},\;\bar{\Omega}^{N,M}=\bar{\Omega}^{N}_{x}\times\bar{\Omega}^{M}_{t}\;\text{ and }\;\Gamma^{N,M}=\Gamma\cap\bar{\Omega}^{N,M}.$

The mesh $\bar{\Omega}^{M}_{t}$ is chosen to be a uniform mesh with $M$ sub-intervals on $[0,T]$. The mesh $\bar{\Omega}^{N}_{x}$ is a piecewise-uniform Shishkin mesh with $N$ mesh intervals. The interval $[0,1]$ is subdivided into $3$ sub-intervals

[TABLE]

where

[TABLE]

Thus, on the sub-interval $(\sigma,1-\sigma],$ a uniform mesh with $\;\frac{N}{2}$ mesh-points is placed and on the subintervals $(0,\sigma]$ and $(1-\sigma,1],\;$a uniform mesh with $\;\frac{N}{4}$ mesh-points is placed.

In particular, when the parameter $\sigma$ takes on its left-hand value, the Shishkin mesh $\bar{\Omega}_{x}^{N}$ becomes a classical uniform mesh. In practice, it is convenient to take $N=4k,\;k\geq 2.$

From the above construction of $\bar{\Omega}_{x}^{N}$ it is clear that the transition points $\{\sigma,1-\sigma\}$ are the only points at which the mesh size can change and that it does not necessarily change at each of these points. The following notations are introduced: $h_{j}=x_{j}-x_{j-1};\;h_{j+1}=x_{j+1}-x_{j},\;J=\{x_{j}=\sigma,1-\sigma:h_{j+1}\neq h_{j}\}.$ For each point $x_{j}$ in the sub-intervals $(0,\sigma]$ and $(1-\sigma,1],\;x_{j}-x_{j-1}=4N^{-1}\sigma$ and for $x_{j}$ in the sub-interval $(\sigma,1-\sigma],$ $x_{j}-x_{j-1}=2N^{-1}(1-2\sigma).$

The construction of $\bar{\Omega}_{x}^{N}$ as a piecewise uniform Shishkin mesh on $[0,1]$ leads to a piecewise uniform Shishkin mesh $\bar{\Omega}^{N,M}$ on $[0,1]\times[0,1]$ by considering the cartesian product of the discrete set $\bar{\Omega}_{x}^{N}=\{x_{j}\}_{j=0}^{N}$ with $\bar{\Omega}_{t}^{M}=\{t_{j}\}_{j=0}^{M}.$

Thus $\bar{\Omega}^{N,M}$ is a piecewise uniform Shishkin grid with $NM$ mesh elements.

5 The discrete problem

In this section, a classical finite difference operator with an appropriate Shishkin mesh is used to construct a numerical method for the problem (1), (2), which is shown later to be first order parameter-uniform convergent in time and essentially first order parameter-uniform convergent in the space variable.

The discrete initial-boundary value problem is now defined by the finite difference scheme on the Shishkin mesh $\bar{\Omega}^{N,M},$ defined in the previous section.

[TABLE]

with

[TABLE]

The problem (27), (28) can also be written in the operator form

[TABLE]

[TABLE]

where $L^{N,M}=D^{-}_{t}-\varepsilon\delta^{2}_{x}+a,\;\;\beta_{0}^{N,M}=I-D_{x}^{+},\;\;\beta_{1}^{N,M}=I+D_{x}^{-}$.

The following discrete results are analogous to those for the continuous case.

Lemma 5.1

Let the assumptions (4) - (8) hold. Then, for any mesh function $\Psi$, the inequalities $\beta_{0}^{N,M}\Psi(0,t_{k})\geq 0,\;\beta_{1}^{N,M}\Psi(1,t_{k})\geq 0,\;\Psi(x_{j},0)\geq 0$ and $L^{N,M}\Psi\;\geq\;0$ on  $\Omega^{N,M}$ imply that ${\Psi}\geq 0$ on $\bar{\Omega}^{N,M}.$

Proof. Let $j^{*},k^{*}$ be such that $\Psi(x_{j^{*}},t_{k^{*}})=\displaystyle\min_{j,k}\Psi(x_{j},t_{k})$ and assume that the lemma is false. Then $\Psi(x_{j^{*}},t_{k^{*}})<0$ . From the hypotheses,  $\Psi(x_{j^{*}},t_{k^{*}})-\Psi(x_{j^{*}},t_{k^{*}-1})<0$, $\Psi(x_{j^{*}},t_{k^{*}})-\Psi(x_{j^{*}-1},t_{k^{*}})<0,\;\Psi(x_{j^{*}+1},t_{k^{*}})-\Psi(x_{j^{*}},t_{k^{*}})>0.$ Hence, $\;D_{t}^{-}\Psi(x_{j^{*}},t_{k^{*}})<0$ and $\delta^{2}_{x}\Psi(x_{j^{*}},t_{k^{*}})\;\geq\;0.$ It follows that, for $(x_{j^{*}},t_{k^{*}})\in\Omega^{N,M},$

[TABLE]

which is a contradiction.

If $x_{j^{*}}=0,$ then $\beta_{0}^{N,M}\Psi(0,t_{k^{*}})=\Psi(0,t_{k^{*}})\;-\;D^{+}_{x}\Psi(0,t_{k^{*}})<0,$ a contradiction. Therefore, $x_{j^{*}}\neq 0,$ and for a similar reason $x_{j^{*}}\neq 1.$ For $t_{k^{*}}=0,\;\Psi(x_{j^{*}},0)<0,$ which is a contradiction. Therefore, $t_{k^{*}}\neq 0.$ Hence the result.

An immediate consequence of this is the following discrete stability result.

Lemma 5.2

Let the assumptions (4) - (8) hold. Then, for any mesh function $\Psi$ defined on $\bar{\Omega}^{N,M},$

[TABLE]

Proof. Define the two mesh functions

[TABLE]

It is not hard to verify that $\beta_{0}^{N,M}{\Theta}^{\pm}(0,t_{k})\geq 0,\;\beta_{1}^{N,M}{\Theta}^{\pm}(1,t_{k})\geq 0,\;{\Theta}^{\pm}(x_{j},0)\geq 0$ and $L^{N,M}{\Theta}^{\pm}\geq 0$ on $\Omega^{N,M}$. It follows from Lemma 5.1 that ${\Theta}^{\pm}\geq 0$ on $\bar{\Omega}^{N,M}$.

The following comparison principle will be used in the proof of the error estimate.

Lemma 5.3

Assume that the mesh functions ${\Phi}$ and ${Z}$ satisfy   $|\beta_{0}^{N,M}Z(0,t_{k})|\leq\\ \beta_{0}^{N,M}\Phi(0,t_{k}),\;|\beta_{1}^{N,M}Z(1,t_{k})|\leq\beta_{1}^{N,M}\Phi(1,t_{k}),\;|Z(x_{j},0)|\leq\Phi(x_{j},0)\text{ and }|L^{N,M}{Z}|\leq L^{N,M}{\Phi}\\ \text{ on }\Omega^{N,M}.$ Then $|Z|\leq\Phi\text{ on }\bar{\Omega}^{N,M}.$

Proof. Define the two mesh functions $\Psi^{\pm}$ by   $\Psi^{\pm}={\Phi}\pm{Z}.$

Then, $\Psi^{\pm}$ satisfy $\beta_{0}^{N,M}\Psi^{\pm}(0,t_{k})\geq 0,\;\beta_{1}^{N,M}\Psi^{\pm}(1,t_{k})\geq 0,\;\Psi^{\pm}(x_{j},0)\geq 0,\text{ and }L^{N,M}\Psi^{\pm}\geq 0\text{ on }\Omega^{N,M}.$ The required result follows from the Lemma 5.1.

6 The local truncation error

From Lemma 5.2, it is seen that in order to bound the error $U-u$, it suffices to bound $\beta_{0}^{N,M}(U-u)(0,t_{k}),\;\beta_{1}^{N,M}(U-u)(1,t_{k}),\;(U-u)(x_{j},0),\text{ and }L^{N,M}(U-u)$. Note that, for $(x_{j},t_{k})\in\Omega^{N,M},$ $L^{N,M}(U-u)=L^{N,M}U-L^{N,M}u=f-L^{N,M}u=Lu-L^{N,M}u=(L-L^{N,M})u.$

It follows that $L^{N,M}(U-u)=(\dfrac{\partial}{\partial t}-D^{-}_{t})u-\varepsilon(\dfrac{\partial^{2}}{\partial x^{2}}-\delta^{2}_{x})u.$

Let $V,W$ be the discrete analogues of $v,w$ respectively, given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $v$ and $w$ are the solutions of (14), (15) and (16), (17) respectively. Further,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The local truncation error of the smooth and singular components can be treated separately. Note that, for any smooth function $\psi$ and for each $(x_{j},t_{k})\in\Omega^{N,M}$, the following distinct estimates of the local truncation error hold:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here $I_{j}=[x_{j-1},x_{j+1}]$.

7 Error estimate

The proof of the theorem on the error estimate is broken into two parts. First, a theorem concerning the error in the smooth component is established. Then the error in the singular component is estimated.

Define the barrier function through

[TABLE]

where $C$ is sufficiently large and $\theta$ is a piecewise linear polynomial on $\bar{\Omega}$ for each $x_{j}=\sigma\in J$ defined by

[TABLE]

Also note that,

[TABLE]

Then, on $\Omega^{N,M}$, ${\Phi}$ satisfies  $0\leq\Phi(x_{j},t_{k})\leq C(M^{-1}+N^{-1}\ln N).$

Also,

[TABLE]

For $x_{j}\notin J$, using (34), it is not hard to see that,

[TABLE]

and, for $x_{j}\in J$, using (34), it is not hard to see that,

[TABLE]

The following theorem gives the estimate of the error in the smooth component $V$.

Theorem 7.1

Let the assumptions (4) - (8) hold. Let $v$ denote the smooth component of the solution of the problem (1), (2) and $V$ denote the smooth component of the solution of the problem (27), (28). Then

[TABLE]

Proof. From the expression (31),

[TABLE]

From the expression (30),

[TABLE]

Also note that, the expressions (7), (35) and (36) yield

[TABLE]

For each mesh point $x_{j},$ there are two possibilities: either $x_{j}\notin J$ or $x_{j}\in J$.

For $x_{j}\notin J$, from the expressions (29), (33) and Lemma 3.4, it follows that

[TABLE]

For $x_{j}\in J$, then $x_{j}=\sigma\text{ or }1-\sigma.$ Here the argument for $x_{j}=\sigma$ is given and for $x_{j}=1-\sigma$ it is analogous.

For $x_{j}=\sigma,$ from the expressions (29), (33) and Lemma 3.4, it follows that

[TABLE]

From the expressions (37), (38), (39) and the comparison principle, the required result follows.

The following theorem gives the estimate of the error in the singular component $\vec{W}$.

Theorem 7.2

Let the assumptions (4) - (8) hold. Let $w$ be the singular component of the solution of the problem (1), (2) and $W$ be the singular component of the solution of the problem (27), (28). Then

[TABLE]

Proof. From the expression (31),

[TABLE]

From the expression (30),

[TABLE]

Also note that,

[TABLE]

For $x_{j}\notin J,$ from the expressions (29), (33) and Lemma 3.5, it follows that

[TABLE]

For $x_{j}\in J,$ from the expressions (29), (33) and Lemma 3.5, it follows that

[TABLE]

From the expressions (42), (43), (44) and the comparison principle, the required result follows.

The following theorem gives a parameter uniform bound which is first order in time and essentially first order in space for the convergence of the discrete solution.

Theorem 7.3

Let the assumptions (4) - (8) hold. Let $u$ denote the solution of the problem (1), (2) and $U$ denote the solution of the problem (27), (28). Then

[TABLE]

Proof. An application of the triangular inequality and the results of Theorem 7.1 and Theorem 7.2 lead to the required result.

8 Numerical Illustration

The numerical method proposed above is illustrated through the example presented in this section. The method proposed above is applied to solve the problem and the parameter-uniform order of convergence and the parameter-uniform error constants are computed. To get the order of convergence in the variable $t$ seperately, a Shishkin mesh is considered for $x$ and the resulting problem is solved for various uniform meshes with respect to $t$. In order to get the order of convergence in the variable $x$ seperately, a uniform mesh is considered for $t$ and the resulting problem is solved for various piecewise uniform Shishkin meshes with respect to $x$. The same two-mesh algorithm found in [2] is applied to get parameter-uniform order of convergence and the error constants. The numerical results are presented in Table 1 and Table 2.

**Example ** Consider the problem

[TABLE]

[TABLE]

For various values of $\varepsilon,$ the maximum errors, the $\varepsilon$- uniform order of convergence and the $\varepsilon$-uniform error constant are computed. Fixing a Shishkin mesh on $[0,1]$ with $64$ points horizontally, the problem is solved by the method suggested above. The order of convergence and the error constant for $u$ are calculated for $t$ using two-mesh algorithm and the results are presented in Table 1. A uniform mesh on $[0,1]$ with $256$ points vertically is considered and the order of convergence and the error constant for $u$ in the variable $x$ using two-mesh algorithm are calculated and the results are presented in Table 2.

The notations $D^{N},\;p^{N}$ and $C_{p^{*}}^{N}$ denote the $\varepsilon$-uniform maximum pointwise two-mesh differences, the $\varepsilon$-uniform order of convergence and the $\varepsilon$-uniform error constant respectively and are given by $D^{N}=\displaystyle\max_{\varepsilon}\;D^{N}_{\varepsilon}$ where $D^{N}_{\varepsilon}=\;\parallel U^{N}_{\varepsilon}-U^{2N}_{\varepsilon}\parallel_{\Omega^{N}}$, $p^{N}=log_{2}\displaystyle\frac{D^{N}}{D^{2N}}$ and $C_{p^{*}}^{N}=\displaystyle\frac{D^{N}N^{p^{*}}}{1-2^{-p^{*}}}.$ Then the parameter-uniform order of convergence and the error constant are given by $p^{*}=\displaystyle\min_{N}p^{N}$ and $C^{*}_{p^{*}}=\displaystyle\max_{N}C^{N}_{p^{*}}$ respectively.

It is evident from the Figures 4-4 that the solution $u$ exhibits no layer whereas the derivative $\dfrac{\partial u}{\partial x}$ exhibits parabolic twin boundary layers at $(0,t)$ and $(1,t),\;0\leq t\leq 1.$ Further, the $t$- order of convergence and the $x$- order of convergence of the numerical method presented in Table 1 and Table 2 agree with the theoretical result.

Acknowledgment

The first author sincerely thanks the University Grants Commission, New Delhi, India, for the financial support through the Rajiv Gandhi National Fellowship to carry out this work.