Robust numerical schemes for singularly perturbed delay parabolic convection diffusion problems with degenerate coefficient

TL;DR

This paper develops and analyzes a robust numerical scheme combining backward Euler and specialized spatial discretization on a modified Shishkin mesh to effectively solve singularly perturbed delay convection diffusion problems with degenerate coefficients.

Contribution

It introduces a new numerical method with stability and convergence analysis for challenging delay convection diffusion equations with degenerate coefficients.

Findings

The proposed scheme is stable and convergent as predicted theoretically.

Numerical results demonstrate the scheme's effectiveness over traditional methods.

Comparison shows improved accuracy on Shishkin mesh.

Abstract

This article studies a dirichlet boundary value problem for singularly perturbed time delay convection diffusion equation with degenerate coefficient. A priori explicit bounds are established on the solution and its derivatives. For asymptotic analysis of the spatial derivatives the solution is decomposed into regular and singular parts. To approximate the solution a numerical method is considered which consists of backward Euler scheme for time discretization on uniform mesh and a combination of midpoint upwind and central difference scheme for the spatial discretization on modified Shishkin mesh. Stability analysis is carried out, numerical results are presented and comparison is done with upwind scheme on uniform mesh as well as upwind scheme on Shishkin mesh to demonstrate the effectiveness of the proposed method. The convergence obtained in practical satisfies the theoretical predictions.

Tables

Table 1. Table 1 : The maximum pointwise errors E ε N , M subscript superscript 𝐸 𝑁 𝑀 𝜀 E^{N,M}_{\varepsilon} and the corresponding order of convergence q ε N , M subscript superscript 𝑞 𝑁 𝑀 𝜀 q^{N,M}_{\varepsilon} for the Problem 1 using different schemes with p = 1 𝑝 1 p=1 and M = N 𝑀 𝑁 M=N .

$\epsilon \downarrow$	Scheme	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$2^{-8}$	Upwind scheme on	7.061e-02	3.726e-02	1.964e-02	1.009e-02	5.122e-03	2.581e-03	1.296e-03
	uniform mesh	9.224e-01	9.236e-01	9.617e-01	9.774e-01	9.887e-01	9.942e-01
	Upwind scheme on	6.350e-02	3.726e-02	1.964e-02	1.009e-02	5.122e-03	2.581e-03	1.296e-03
	Shishkin mesh	7.693e-01	9.236e-01	9.617e-01	9.774e-01	9.887e-01	9.942e-01
	Hybrid scheme on	4.308e-03	2.056e-03	1.031e-03	5.177e-04	2.647e-04	1.352e-04	6.834e-05
	Shishkin mesh	1.068e+00	9.957e-01	9.936e-01	9.676e-01	9.694e-01	9.843e-01
$2^{-12}$	Upwind scheme	1.481e-01	1.390e-01	6.916e-02	3.640e-02	1.927e-02	9.902e-03	5.028e-03
	uniform mesh	9.131e-02	1.007e+00	9.259e-01	9.179e-01	9.604e-01	9.777e-01
	Upwind scheme on	6.438e-02	3.929e-02	2.365e-02	1.385e-02	7.897e-03	4.424e-03	2.443e-03
	Shishkin mesh	7.125e-01	7.320e-01	7.720e-01	8.107e-01	8.360e-01	8.565e-01
	Hybrid scheme on	4.124e-03	2.097e-03	1.054e-03	5.295e-04	2.728e-04	1.395e-04	7.084e-05
	Shishkin mesh	9.760e-01	9.927e-01	9.925e-01	9.571e-01	9.670e-01	9.780e-01
$2^{-16}$	Upwind scheme	2.176e-02	6.828e-02	1.466e-01	1.379e-01	6.858e-02	3.610e-02	1.912e-02
	uniform mesh	-1.649e+00	-1.102e+00	8.809e-02	1.008e+00	9.258e-01	9.169e-01
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.386e-02	7.905e-03	4.428e-03	2.445e-03
	Shishkin mesh	7.128e-01	7.318e-01	7.720e-01	8.107e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.303e-04	2.739e-04	1.402e-04	7.119e-05
	Shishkin mesh	9.752e-01	9.923e-01	9.924e-01	9.530e-01	9.664e-01	9.777e-01
$2^{-20}$	Upwind scheme	3.961e-03	6.240e-03	2.076e-02	6.750e-02	1.461e-01	1.376e-01	6.842e-02
	uniform mesh	-6.558e-01	-1.734e+00	-1.701e+00	-1.114e+00	8.623e-02	1.008e+00
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.905e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.124e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.523e-01	9.662e-01	9.776e-01
$2^{-24}$	Upwind scheme	3.961e-03	2.062e-03	1.843e-03	5.574e-03	2.037e-02	6.727e-02	1.460e-01
	uniform mesh	9.417e-01	1.618e-01	-1.596e+00	-1.870e+00	-1.724e+00	-1.118e+00
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.905e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.125e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01
$2^{-28}$	Upwind scheme	3.961e-03	2.062e-03	1.057e-03	5.943e-04	1.476e-03	5.380e-03	2.027e-02
	uniform mesh	9.417e-01	9.635e-01	8.312e-01	-1.312e+00	-1.866e+00	-1.913e+00
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.906e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.126e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01
$2^{-32}$	Upwind scheme	3.961e-03	2.062e-03	1.057e-03	5.367e-04	2.702e-04	4.030e-04	1.378e-03
	uniform mesh	9.417e-01	9.635e-01	9.782e-01	9.900e-01	-5.767e-01	-1.774e+00
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.906e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.126e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01
$2^{-36}$	Upwind scheme	3.961e-03	2.062e-03	1.057e-03	5.367e-04	2.702e-04	1.356e-04	1.173e-04
	uniform mesh	9.417e-01	9.635e-01	9.782e-01	9.900e-01	9.953e-01	2.082e-01
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.906e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.126e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01
$2^{-40}$	Upwind scheme	3.961e-03	2.062e-03	1.057e-03	5.367e-04	2.702e-04	1.356e-04	6.789e-05
	uniform mesh	9.417e-01	9.635e-01	9.782e-01	9.900e-01	9.953e-01	9.976e-01
	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.906e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
	Hybrid scheme on	4.126e-03	2.099e-03	1.055e-03	5.304e-04	2.741e-04	1.403e-04	7.126e-05
	Shishkin mesh	9.750e-01	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01
${𝐄}^{𝐍,𝐌}$	Upwind scheme on	1.481e-01	1.390e-01	1.466e-01	1.379e-01	1.461e-01	1.376e-01	1.460e-01
	uniform mesh	9.131e-02	-7.659e-02	8.809e-02	-8.307e-02	8.623e-02	-8.492e-02
${𝐄}^{𝐍,𝐌}$	Upwind scheme on	6.444e-02	3.932e-02	2.368e-02	1.387e-02	7.906e-03	4.429e-03	2.446e-03
	Shishkin mesh	7.128e-01	7.317e-01	7.719e-01	8.106e-01	8.360e-01	8.566e-01
${𝐄}^{𝐍,𝐌}$	Hybrid scheme	4.308e-03	2.099e-03	1.055e-03	5.304e-03	2.741e-04	1.403e-04	7.126e-05
	on Shishkin mesh (3.1)	1.0373	9.922e-01	9.923e-01	9.521e-01	9.662e-01	9.776e-01

Table 2. Table 2 : The maximum pointwise errors E ε N , M subscript superscript 𝐸 𝑁 𝑀 𝜀 E^{N,M}_{\varepsilon} and the corresponding order of convergence q ε N , M subscript superscript 𝑞 𝑁 𝑀 𝜀 q^{N,M}_{\varepsilon} for the Problem 1 using the Richardson extrapolation on the hybrid scheme ( 3.1 ) with p = 1 𝑝 1 p=1 and M = N 𝑀 𝑁 M=N .

$\epsilon \downarrow$	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$2^{-8}$	6.397e-03	2.258e-03	7.740e-04	2.414e-04	6.028e-05	1.507e-05	4.443e-06
	1.502e+00	1.545e+00	1.681e+00	2.001e+00	2.000e+00	1.762e+00
$2^{-12}$	6.352e-03	2.246e-03	7.705e-04	2.561e-04	8.279e-05	2.612e-05	8.072e-06
	1.500e+00	1.543e+00	1.589e+00	1.629e+00	1.664e+00	1.694e+00
$2^{-16}$	6.387e-03	2.248e-03	7.695e-04	2.558e-04	8.274e-05	2.611e-05	8.069e-06
	1.506e+00	1.547e+00	1.589e+00	1.629e+00	1.664e+00	1.694e+00
$2^{-20}$	6.407e-03	2.258e-03	7.726e-04	2.563e-04	8.275e-05	2.610e-05	8.068e-06
	1.505e+00	1.547e+00	1.592e+00	1.631e+00	1.665e+00	1.694e+00
$2^{-24}$	6.409e-03	2.261e-03	7.741e-04	2.568e-04	8.292e-05	2.613e-05	8.069e-06
	1.503e+00	1.546e+00	1.592e+00	1.631e+00	1.666e+00	1.695e+00
$2^{-28}$	6.410e-03	2.261e-03	7.745e-04	2.570e-04	8.299e-05	2.616e-05	8.077e-06
	1.503e+00	1.546e+00	1.591e+00	1.631e+00	1.666e+00	1.695e+00
$2^{-32}$	6.410e-03	2.262e-03	7.745e-04	2.570e-04	8.301e-05	2.616e-05	8.080e-06
	1.503e+00	1.546e+00	1.591e+00	1.631e+00	1.666e+00	1.695e+00
$2^{-36}$	6.410e-03	2.262e-03	7.746e-04	2.571e-04	8.301e-05	2.617e-05	8.081e-06
	1.503e+00	1.546e+00	1.591e+00	1.631e+00	1.666e+00	1.695e+00
$2^{-40}$	6.410e-03	2.262e-03	7.746e-04	2.571e-04	8.301e-05	2.617e-05	8.081e-06
	1.503e+00	1.546e+00	1.591e+00	1.631e+00	1.666e+00	1.695e+00
${𝐄}^{𝐍,𝐌}$ (hybrid scheme (3.1)
with extrapolation)	6.410e-03	2.262e-03	7.746e-04	2.571e-04	8.301e-05	2.617e-05	8.081e-06
${𝐪}^{𝐍,𝐌}$ (hybrid scheme (3.1)
with extrapolation)	1.503e+00	1.546e+00	1.591e+00	1.631e+00	1.666e+00	1.695e+00

Table 3. Table 3 : The maximum pointwise errors E ε N , M subscript superscript 𝐸 𝑁 𝑀 𝜀 E^{N,M}_{\varepsilon} and the corresponding order of convergence q ε N , M subscript superscript 𝑞 𝑁 𝑀 𝜀 q^{N,M}_{\varepsilon} for the Problem 2 using different schemes with p = 1 𝑝 1 p=1 and M = N 𝑀 𝑁 M=N .

$\epsilon \downarrow$	Scheme	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$2^{-8}$	Upwind scheme on	6.278e-02	3.256e-02	1.716e-02	8.789e-03	4.454e-03	2.243e-03	1.126e-03
	uniform mesh	9.474e-01	9.243e-01	9.649e-01	9.805e-01	9.896e-01	9.947e-01
	Upwind scheme on	6.278e-02	3.256e-02	1.716e-02	8.789e-03	4.454e-03	2.243e-03	1.126e-03
	Shishkin mesh	9.474e-01	9.243e-01	9.649e-01	9.805e-01	9.896e-01	9.947e-01
	Hybrid scheme on	1.468e-02	3.811e-03	1.817e-03	9.210e-04	4.637e-04	2.326e-04	1.165e-04
	Shishkin mesh	1.946e+00	1.068e+00	9.806e-01	9.901e-01	9.950e-01	9.975e-01
$2^{-12}$	Upwind scheme on	1.203e-01	1.131e-01	5.492e-02	2.879e-02	1.516e-02	7.794e-03	3.959e-03
	Uniform mesh	8.849e-02	1.043e+00	9.315e-01	9.251e-01	9.601e-01	9.773e-01
	Upwind scheme on	1.039e-01	5.739e-02	3.595e-02	2.148e-02	1.253e-02	7.108e-03	3.959e-03
	Shishkin mesh	8.570e-01	6.746e-01	7.429e-01	7.779e-01	8.177e-01	8.443e-01
	Hybrid scheme on	2.615e-02	9.482e-03	3.466e-03	1.329e-03	5.531e-04	2.725e-04	1.365e-04
	Shishkin mesh	1.464e+00	1.452e+00	1.382e+00	1.265e+00	1.022e+00	9.976e-01
$2^{-16}$	Upwind scheme on	1.663e-02	5.109e-02	1.165e-01	1.119e-01	5.427e-02	2.849e-02	1.502e-02
	Uniform mesh	-1.619e+00	-1.189e+00	5.864e-02	1.044e+00	9.297e-01	9.232e-01
	Upwind scheme on	1.023e-01	5.616e-02	3.534e-02	2.115e-02	1.236e-02	7.013e-03	3.910e-03
	Shishkin mesh	8.648e-01	6.681e-01	7.405e-01	7.747e-01	8.182e-01	8.428e-01
	Hybrid scheme on	2.646e-02	9.559e-03	3.666e-03	1.354e-03	5.732e-04	2.800e-04	1.401e-04
	Shishkin mesh	1.469e+00	1.383e+00	1.437e+00	1.240e+00	1.034e+00	9.986e-01
$2^{-20}$	Upwind scheme on	1.121e-02	6.017e-03	1.416e-02	5.130e-02	1.169e-01	1.119e-01	5.430e-02
	Uniform mesh	8.972e-01	-1.235e+00	-1.857e+00	-1.188e+00	6.255e-02	1.044e+00
	Upwind scheme on	1.022e-01	5.612e-02	3.533e-02	2.115e-02	1.236e-02	7.010e-03	3.909e-03
	Shishkin mesh	8.645e-01	6.675e-01	7.405e-01	7.747e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.647e-02	9.530e-03	3.476e-03	1.590e-03	7.209e-04	3.157e-04	1.420e-04
	Shishkin mesh	1.474e+00	1.455e+00	1.129e+00	1.141e+00	1.191e+00	1.153e+00
$2^{-24}$	Upwind scheme on	1.121e-02	6.017e-03	3.162e-03	3.208e-03	1.496e-02	5.180e-02	1.171e-01
	Uniform mesh	8.972e-01	9.285e-01	-2.084e-02	-2.221e+00	-1.792e+00	-1.177e+00
	Upwind scheme on	1.022e-01	5.613e-02	3.534e-02	2.115e-02	1.236e-02	7.012e-03	3.909e-03
	Shishkin mesh	8.643e-01	6.675e-01	7.405e-01	7.747e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.646e-02	9.519e-03	3.480e-03	1.335e-03	6.125e-04	3.029e-04	1.450e-04
	Shishkin mesh	1.475e+00	1.452e+00	1.382e+00	1.124e+00	1.016e+00	1.063e+00
$2^{-28}$	Upwind scheme on	1.121e-02	6.017e-03	3.162e-03	1.637e-03	8.373e-04	3.785e-03	1.528e-02
	Uniform mesh	8.972e-01	9.285e-01	9.499e-01	9.670e-01	-2.177e+00	-2.013e+00
	Upwind scheme on	1.022e-01	5.614e-02	3.534e-02	2.115e-02	1.236e-02	7.012e-03	3.910e-03
	Shishkin mesh	8.642e-01	6.675e-01	7.405e-01	7.748e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.646e-02	9.515e-03	3.654e-03	1.528e-03	5.937e-04	2.816e-04	1.411e-04
	Shishkin mesh	1.475e+00	1.381e+00	1.258e+00	1.364e+00	1.076e+00	9.968e-01
$2^{-32}$	Upwind scheme on	1.121e-02	6.017e-03	3.162e-03	1.637e-03	8.373e-04	4.247e-04	8.944e-04
	Uniform mesh	8.972e-01	9.285e-01	9.499e-01	9.670e-01	9.791e-01	-1.074e+00
	Upwind scheme on	1.022e-01	5.614e-02	3.535e-02	2.115e-02	1.236e-02	7.012e-03	3.910e-03
	Shishkin mesh	8.642e-01	6.675e-01	7.405e-01	7.748e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.645e-02	9.514e-03	3.671e-03	1.472e-03	6.550e-04	2.910e-04	1.406e-04
	Shishkin mesh	1.475e+00	1.374e+00	1.319e+00	1.168e+00	1.171e+00	1.049e+00
$2^{-36}$	Upwind scheme on	1.121e-02	6.017e-03	3.162e-03	1.637e-03	8.373e-04	4.247e-04	2.143e-04
	Uniform mesh	8.972e-01	9.285e-01	9.499e-01	9.670e-01	9.791e-01	9.872e-01
	Upwind scheme on	1.022e-01	5.614e-02	3.535e-02	2.115e-02	1.236e-02	7.012e-03	3.910e-03
	Shishkin mesh	8.642e-01	6.675e-01	7.405e-01	7.748e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.645e-02	9.514e-03	3.675e-03	1.473e-03	6.270e-04	2.924e-04	1.416e-04
	Shishkin mesh	1.475e+00	1.372e+00	1.319e+00	1.232e+00	1.101e+00	1.046e+00
$2^{-40}$	Upwind scheme on	1.121e-02	6.017e-03	3.162e-03	1.637e-03	8.373e-04	4.247e-04	2.143e-04
	Uniform mesh	8.972e-01	9.285e-01	9.499e-01	9.670e-01	9.791e-01	9.872e-01
	Upwind scheme on	1.022e-01	5.614e-02	3.535e-02	2.115e-02	1.236e-02	7.012e-03	3.910e-03
	Shishkin mesh	8.642e-01	6.675e-01	7.405e-01	7.748e-01	8.183e-01	8.428e-01
	Hybrid scheme on	2.645e-02	9.514e-03	3.677e-03	1.473e-03	6.270e-04	2.865e-04	1.417e-04
	Shishkin mesh	1.475e+00	1.372e+00	1.320e+00	1.232e+00	1.130e+00	1.016e+00
${𝐄}^{𝐍,𝐌}$	Upwind scheme on	1.203e-01	1.131e-01	1.165e-01	1.119e-01	1.169e-01	1.119e-01	1.171e-01
	Uniform mesh	8.849e-02	-4.264e-02	5.864e-02	-6.339e-02	6.255e-02	-6.519e-02
${𝐄}^{𝐍,𝐌}$	Upwind scheme	1.039e-01	5.739e-02	3.595e-02	2.148e-02	1.253e-02	7.108e-03	3.959e-03
	Shishkin mesh	8.570e-01	6.746e-01	7.429e-01	7.779e-01	8.177e-01	8.443e-01
${𝐄}^{𝐍,𝐌}$	hybrid scheme on	2.647e-02	9.559e-03	3.677e-03	1.590e-03	7.209e-04	3.157e-04	1.450e-04
	Shishkin mesh (3.1)	1.470e+00	1.378e+00	1.209e+00	1.141e+00	1.191e+00	1.122e+00

Table 4. Table 4 : The maximum pointwise errors E ε N , M subscript superscript 𝐸 𝑁 𝑀 𝜀 E^{N,M}_{\varepsilon} and the corresponding order of convergence q ε N , M subscript superscript 𝑞 𝑁 𝑀 𝜀 q^{N,M}_{\varepsilon} for the Problem 2 using the Richardson extrapolation on the hybrid scheme ( 3.1 ) with p = 1 𝑝 1 p=1 and M = N 𝑀 𝑁 M=N .

$\epsilon \downarrow$	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$2^{-8}$	1.347e-02	3.252e-03	8.099e-04	2.023e-04	5.054e-05	1.999e-05	9.089e-06
	2.051e+00	2.005e+00	2.001e+00	2.001e+00	1.338e+00	1.137e+00
$2^{-12}$	2.301e-02	8.141e-03	2.716e-03	9.060e-04	2.917e-04	9.194e-05	2.838e-05
	1.499e+00	1.584e+00	1.584e+00	1.635e+00	1.666e+00	1.696e+00
$2^{-16}$	2.298e-02	8.123e-03	2.711e-03	9.035e-04	2.910e-04	9.170e-05	2.831e-05
	1.500e+00	1.583e+00	1.585e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-20}$	2.299e-02	8.126e-03	2.712e-03	9.031e-04	2.909e-04	9.166e-05	2.830e-05
	1.500e+00	1.583e+00	1.586e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-24}$	2.299e-02	8.127e-03	2.713e-03	9.031e-04	2.909e-04	9.165e-05	2.829e-05
	1.500e+00	1.583e+00	1.587e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-28}$	2.299e-02	8.127e-03	2.713e-03	9.031e-04	2.909e-04	9.165e-05	2.829e-05
	1.500e+00	1.583e+00	1.587e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-32}$	2.299e-02	8.127e-03	2.713e-03	9.031e-04	2.909e-04	9.165e-05	2.829e-05
	1.500e+00	1.583e+00	1.587e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-36}$	2.299e-02	8.127e-03	2.713e-03	9.031e-04	2.909e-04	9.165e-05	2.829e-05
	1.500e+00	1.583e+00	1.587e+00	1.634e+00	1.666e+00	1.696e+00
$2^{-40}$	2.299e-02	8.127e-03	2.713e-03	9.031e-04	2.909e-04	9.165e-05	2.829e-05
	1.500e+00	1.583e+00	1.587e+00	1.634e+00	1.666e+00	1.696e+00
${𝐄}^{𝐍,𝐌}$ (hybrid scheme (3.1)
with extrapolation)	2.301e-02	8.141e-03	2.716e-03	9.060e-04	2.917e-04	9.194e-05	2.838e-05
${𝐪}^{𝐍,𝐌}$ (hybrid scheme (3.1)
with extrapolation)	1.499e+00	1.584e+00	1.584e+00	1.635e+00	1.666e+00	1.696e+00

Table 5. Table 5 : The ε 𝜀 \varepsilon -uniform maximum pointwise errors E N , M superscript 𝐸 𝑁 𝑀 E^{N,M} and the corresponding order of convergence q N , M superscript 𝑞 𝑁 𝑀 q^{N,M} for the Problem 1 using the Richardson extrapolation on the hybrid scheme ( 3.1 ) for different values of p 𝑝 p and M = N 𝑀 𝑁 M=N .

$p\downarrow$	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$3$	1.014e-02	3.724e-03	1.293e-03	4.318e-04	1.396e-04	4.402e-05	1.360e-05
	1.445e+00	1.526e+00	1.582e+00	1.629e+00	1.665e+00	1.695e+00
$6$	1.183e-02	4.284e-03	1.515e-03	6.023e-04	2.183e-04	7.572e-05	2.527e-05
	1.465e+00	1.500e+00	1.331e+00	1.464e+00	1.528e+00	1.583e+00

Table 6. Table 6 : The ε 𝜀 \varepsilon -uniform maximum pointwise errors E N , M superscript 𝐸 𝑁 𝑀 E^{N,M} and the corresponding order of convergence q N , M superscript 𝑞 𝑁 𝑀 q^{N,M} for the Problem 2 using the Richardson extrapolation on the hybrid scheme ( 3.1 ) for different values of p 𝑝 p and M = N 𝑀 𝑁 M=N .

$p\downarrow$	$N=32$	$N=64$	$N=128$	$N=256$	$N=512$	$N=1024$	$N=2048$
$2$	1.347e-02	4.914e-03	1.714e-03	5.728e-04	1.852e-04	5.843e-05	1.805e-05
	1.454e+00	1.519e+00	1.582e+00	1.629e+00	1.664e+00	1.695e+00
$5$	2.752e-02	1.070e-02	3.766e-03	1.259e-03	4.080e-04	1.287e-04	3.977e-05
	1.363e+00	1.506e+00	1.581e+00	1.626e+00	1.664e+00	1.695e+00

Equations

$$({\partial }_{t}u-\epsilon {\partial }_x^{2}u)(x,t)=v(g(u(x,t-\tau ))){\partial }_{x}u(x,t)+c[f(u(x,t-\tau ))-u(x,t)],(x,t)\in D.$$

$$L_{\epsilon }u(x,t)=\left(\epsilon \frac{{\partial }^{2}u}{\partial x^2}+a\frac{\partial u}{\partial x}-b\frac{\partial u}{\partial t}-cu\right)(x,t)=e(x,t)u(x,t-\tau )+f(x,t),$$

$$u(x,t)$$

$$u(0,t)$$

$$u(1,t)$$

$$a(x,t)$$

$$a_0(x,t)$$

$$b(x,t)$$

$$c(x,t)$$

$$e(x,t)$$

$$\Vert f{\Vert }_{\overline{Q}}=\underset{x\in \overline{Q}}{\sup }|f(x)|.$$

$$\Vert u{\Vert }_{\overline{Q}}\le \Vert u{\Vert }_{\Gamma }+\frac{T}{\beta }\Vert f{\Vert }_{\overline{Q}}.$$

$${\psi }^{\pm }(x,t)=\Vert u{\Vert }_{\Gamma }+\frac{t}{\beta }\Vert f{\Vert }_{\overline{Q}}\pm u(x,t),$$

$$a,b,c,f\in C_{\lambda }^{l,l/2}(\overline{Q}),s(x,t)\in C_{\lambda }^{l+2,l/2+1}({\Gamma }_b),$$

$$q_0(t)\in C_{\lambda }^{l/2+1}({\overline{\Gamma }}_0),q_1(t)\in C_{\lambda }^{l/2+1}({\overline{\Gamma }}_1),l\ge 0,\lambda \in (0,1).$$

$$\frac{{\partial }^{k+k_0}s(x,t)}{\partial x^k\partial t^{k_0}}=0,\frac{{\partial }^{k_0}{q}_0(t)}{\partial t^{k_0}}=0,\frac{{\partial }^{k_0}{q}_1(t)}{\partial t^{k_0}}=0,0\le k+2k_0\le l+2,$$

$$\frac{{\partial }^{k+k_0}f(x,t)}{\partial x^k\partial t_0^k}=0,0\le k+2k_0\le l,(x,t)\in S^c,$$

$${\Vert \frac{{\partial }^{i+j}u}{\partial x^i\partial t^j}\Vert }_{\overline{Q}}\le C{\epsilon }^{-i/2},\forall 0\le i+2j\le 6,$$

$$L_{\epsilon }u(x,t)$$

$$u(x,0)$$

$$u(0,t)$$

$$u(1,t)$$

$$\left({\tilde{u}}_{\tilde{x}\tilde{x}}+\tilde{a}{\epsilon }^{\frac{p-1}{2}}{\tilde{u}}_{\tilde{x}}-\tilde{b}{\tilde{u}}_t-\tilde{c}\tilde{u}\right)(\tilde{x},t)$$

$$\tilde{u}(\tilde{x},t)$$

$$\tilde{u}(0,t)$$

$$\tilde{u}(1/\sqrt{\epsilon },t)$$

$${\Vert \frac{{\partial }^{i+j}\tilde{u}}{\partial x^i\partial t^j}\Vert }_{{\tilde{N}}_{\delta }}\le C(1+\Vert \tilde{u}{\Vert }_{{\tilde{N}}_{2\delta }}).$$

$${\Vert \frac{{\partial }^{i+j}u}{\partial x^i\partial t^j}\Vert }_{{\overline{Q}}_1}\le C{\epsilon }^{-i/2}(1+\Vert u{\Vert }_{{\overline{Q}}_1}).$$

$$L_{\epsilon }u(x,t)$$

$$u(x,\tau )$$

$$u(0,t)$$

$$u(1,t)$$

$$u(x,t)=y(x,t)+z(x,t).$$

$$y=(y_0+\sqrt{\epsilon }{y}_1+\epsilon y_2+{\epsilon }^{3/2}{y}_3)(x,t),$$

$$a(x,t)\frac{\partial y_0}{\partial x}-b(x,t)\frac{\partial y_0}{\partial t}-c(x,t)y_0$$

$$y_0(x,t)$$

$$a(x,t)\frac{\partial y_1}{\partial x}-b(x,t)\frac{\partial y_1}{\partial t}-c(x,t)y_1$$

$$y_1(x,t)$$

$$a(x,t)\frac{\partial y_2}{\partial x}-b(x,t)\frac{\partial y_2}{\partial t}-c(x,t)y_2$$

$$y_2(x,t)$$

Full text

Keywords : Singular perturbation, parabolic delay differential equations, degenerate coefficient, hybrid scheme, Shishkin mesh, extrapolation.

MSC classification 2010: 65N12, 65N30, 65N06, 65N15.

1 Introduction

Singularly perturbed parabolic delay differential equations (SPPDDEs) plays a crucial role in mathematical modeling of various real life phenomena which takes into consideration the past history of the system along with its present state. The delay or lag represent incubation period, gestation time, transport delays etc. The solution and dynamics of singularly perturbed delay partial differential equations are completely different from those of the partial differential equations without time delay.

A characteristic example of SPPDDEs is the following equation arising in numerical control modeling a furnace used to process metal sheets [1]

[TABLE]

Here, $u$ is the temperature distribution of the metal sheet, which is moving with velocity $v$ and heated by a source given by function $f$; both $v$ and $f$ are adapted dynamically by a controlling device monitoring the current temperature distribution. Since speed of the controller is finite it induces a fixed delay of length $\tau$.

Last one decade has witnessed a growing interest in the numerical study of singularly perturbed parabolic delay differential equations (SPPDDEs). However, uniformly convergent numerical methods are not much developed for SPPDDEs. Numerical study of SPPDDEs for the class of reaction diffusion equations was initiated by Ansari et. al. [1]. The authors used classical finite difference scheme on piecewise uniform Shishkin mesh. Gowrisankar and Natesan  [5] used layer adapted meshes obtained via equidistributing a monitor function for the numerical solution of singularly perturbed parabolic delay differential reaction diffusion problems. Sunil and Mukesh [8] constructed a hybrid scheme consistng of HODIE type on generalized Shishkin mesh in spatial direction and implicit Euler scheme on uniform mesh in time direction for the numerical approximation of singularly perturbed parabolic delay differential reaction diffusion problems. Joginder et al. [12] designed and analyzed a domain decomposition method for the numerical solution of SPPDDEs.

For work on convection diffusion problem for SPPDDEs one can refer to [2, 3, 4, 6, 7, 11]. Aditya and Manju [6] analyzed the weighted difference approximations on piecewise uniform mesh for singularly perturbed delay differential convection diffusion problems and established that the proposed scheme is $L_{2}^{h}$ stable. Gowrisankar and Natesan [4] used layer adapted meshes based on equi-distribution of a monitor function for the numerical solution of SPPDDEs of convection diffusion type. Abhishek and Natesan [2] proposed a hybrid scheme on Shishkin mesh for the numerical solution of convection diffusion problem for SPPDDEs which is almost second order accurate in space and first order in time direction. In [3] the authors applied Richardson extrapolation on simple upwind scheme to obtain almost second order of convergence in space direction and second order of convergence in time direction for SPPDDEs of convection diffusion type. The authors [11] derived a higher order uniformly convergent method which is second order accurate in time and fourth order accurate in space for the numerical solution of singularly perturbed parabolic delay convection diffusion problems.

To the best of our knowledge, all the literature on the numerical solution of SPPDDEs of convection type is restricted to the case when the convection coefficient has same sign throughout the domain. Hence, a very first attempt has been made here to construct a parameter uniform numerical scheme for such a class of problem with degenerating convection coefficient.

We consider the following problem on a rectangular domain:

[TABLE]

where $0<\varepsilon\leq 1$, $\tau>0$, $(x,t)\in Q=\Omega\times(0,T]=(0,1)\times(0,T]$, $\overline{Q}=[0,1]\times[0,T]$, $T$ is some finite time such that $T=k\tau$ for some integer $k>1$, $\Gamma=\Gamma_{b}\cup\Gamma_{l}\cup\Gamma_{r}$, with the interval and boundary conditions given by

[TABLE]

The coefficients $a(x,t)$, $b(x,t)$, $d(x,t)$ and $f(x,t)$ are sufficiently smooth functions such that

[TABLE]

Problem (1.2)-(1) covers the multiple turning points for $p>1$. The solution of the problem $(\ref{1.1})$-$(\ref{1.3})$ exhibit a parabolic boundary layer of width $O(\sqrt{\varepsilon})$ in the neighbourhood of the left boundary $\Gamma_{l}$ as all the characteristic curves of the reduced problem are parallel to the boundary $\Gamma_{r}$.

The existence of the unique solution of the Problem (1.2)-(1) is guaranteed under the assumption that the problem data is H$\ddot{o}$lder continuous, sufficiently smooth and satisfy appropriate compatibility conditions at the corner points (0,0), (1,0), $(0,-\tau)$ and $(1,-\tau)$.

In this article our focus is to develop a higher order robust numerical scheme for the solution of SPPDDEs with multiple degeneracy. The article is designed as follows. In Section 2, analytical aspects of the continuous problem are discussed and a priori estimates are established on the exact solution and its derivatives. In an attempt to design a higher order scheme, in Section 3, the considered problem is discretized by the hybrid scheme on a piecewise-uniform modified Shishkin mesh in the space direction and the implicit Euler method on uniform mesh in the time direction. Stability and error analysis have been carried out for the proposed scheme to establish $\varepsilon$-uniform convergence of $O(N^{-2}L^{2}+\Delta t)$. In Section 4, we combine the hybrid scheme with the Richardson extrapolation to increase the order of convergence from $O(N^{-2}L^{2}+\Delta t)$ to $O(N^{-2}L^{2}+(\Delta t)^{2})$. Numerical experiments are conducted in Section 5 to verify the theoretical results and illustrate the efficiency of the proposed schemes as compared to upwind scheme on uniform mesh as well as upwind scheme on Shishkin mesh.

Notations: Throughout this article, we use C as a generic positive constant independent of $\varepsilon$ and the mesh parameters. All the functions defined on a domain Q are measured in supremum norm, denoted by

[TABLE]

2 A Priori Bounds

In this section, a priori bounds for the solution $u(x,t)$ of the problem (1.2)-(1) and its derivatives are estimated on the domain $\overline{Q}$. We derive some a priori bounds using the method of steps and the minimum principle for the opertaor $L_{\varepsilon}$. The delay term $u(x,t-\tau)$ is a known function $s(x,t-\tau)$ for $(x,t-\tau)\in[0,1]\times[0,\tau]$ and hence the RHS of (1.2) becomes $e(x,t)s(x,t-\tau)+f(x,t)$. This gives us the solution $u(x,t)$ for $(x,t)\in[0,1]\times[0,\tau]$. Using this we can compute the solution $u(x,t)$ for $(x,t)\in[0,1]\times[\tau,2\tau]$ and so on. Hence, using the method of steps the existence and uniqueness results can be established for all $(x,t)\in\overline{Q}$. The operator $L_{\varepsilon}$ satisfies the following minimum principle.

Lemma 2.1** (Minimum Principle).**

Let $w\in C^{2,1}(\overline{Q})$. If $w(x,t)\geq 0$, $\forall~{}(x,t)\in\Gamma$ and $L_{\varepsilon}w(x,t)\leq 0$, $\forall~{}(x,t)\in Q$ then $w(x,t)\geq 0$, $\forall(x,t)\in\overline{Q}$.

Proof.

The proof follows easily from [10]. ∎

Lemma 2.2**.**

Let $u(x,t)$ be the solution of the problem (1.2)-(1) then for all $\varepsilon>0$ the following bound holds

[TABLE]

Proof.

Using the barrier function

[TABLE]

the desired estimate can be obtained using the minimum principle. ∎

The problem data are assumed to be sufficiently smooth that guarantee the required smoothness of the solution on the set $\overline{Q}$. We assume that the data of the problem (1.2)-(1) satisfy the following conditions

[TABLE]

Also, the data of the problem (1.2)-(1) satisfy on $\Gamma^{c}=(\overline{\Gamma}_{0}\cup\overline{\Gamma}_{1})\cap\Gamma_{b}$ ( i.e. the corner points $(0,0)$, $(1,0)$, $(0,-\tau)$ and $(1,-\tau)$ ) the compatibility conditions for the derivatives in $t$ upto order $K_{0}=[l/2]+1$. In the case when the initial function $s(x,t)$ together with its derivatives vanish on the set $\Gamma^{c}$, the following conditions,

[TABLE]

guarantee the compatibility conditions for the derivatives in $t$ upto order $K_{0}=[l/2]+1$. These compatibility conditions ensure the existence of the unique solution $u(x,t)\in C^{K,K/2}_{\lambda}(\overline{Q})$, where $K=l+2$, for the problem (1.2)-(1) [9].

Lemma 2.3**.**

Let the solution $u(x,t)$ of the problem (1.2)-(1) satisfying the assumptions (2)-(2.2) for $K=6$, then

[TABLE]

where C is independent of $\varepsilon$.

Proof.

Using the method of steps we derive the bounds on the derivatives of the solution $u(x,t)$. First, we consider the case for $t\leq\tau$. Since, $u(x,t-\tau)$ is a known function in $[0,1]\times[0,\tau]$, the problem (1.2)-(1) becomes

[TABLE]

The bounds in the interval $Q_{1}$ are obtained as follows. The variable $x$ is transformed to the stretched variable $\widetilde{x}=x/\sqrt{\varepsilon}$, we write the problem (1.2)-(1) as

[TABLE]

where $\widetilde{Q}_{1}=(0,1/\sqrt{\varepsilon})\times(0,\tau]$ and $\widetilde{\Gamma}$ is the boundary analogous to $\Gamma$. Since, for $p>1$ the term $\varepsilon^{\frac{p-1}{2}}$ is very small it can be neglected and for $p=1$ its value is one, the differential equation (2.3) can be made independent of $\varepsilon$. Using [9, estimate (10.5)] we have, for all non negative integers $i,j$ such that $0\leq i+2j\leq 4$, and all $\widetilde{N}_{\delta}$ in $\widetilde{D}_{\varepsilon}$,

[TABLE]

Here, $C$ is a constant independent of ${\widetilde{N}_{\delta}}$ and for any $\delta>0$, ${\widetilde{N}_{\delta}}$ is a neighbourhood of diameter $\delta>0$ in $\widetilde{Q}_{1}$. Transforming back to the original variable $x$, we get

[TABLE]

Using the bounds of $u(x,t)$ given in Lemma 2.2 we get the desired estimates.

Next, we consider the case for $t\in[\tau,2\tau]$. In this case $u(x,t)$ is the solution of the following initial boundary value problem (IBVP):

[TABLE]

Again the RHS is a known function so the proof follows on the similar lines as discussed in the case for $[0,\tau]$. We proceed similarly to prove the result for $t\in[0,T]$. ∎

The bounds obtained in the Lemma 2 are not sufficient for proving $\varepsilon$-uniform error estimates. Therefore, stronger bounds on these derivatives are obtained by decomposing the solution $u(x,t)$ into the regular part $y(x,t)$ and the singular part $z(x,t)$. We define

[TABLE]

The regular component $y(x,t)$ is further decomposed into the sum

[TABLE]

where $y_{0}$, $y_{1}$, $y_{2}$ and $y_{3}$ are defined as

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Therefore, the regular component $y(x,t)$ satisfies

[TABLE]

The singular component $z(x,t)$ satisfies the following IBVP:

[TABLE]

Theorem 2.1**.**

For all non-negative integers $i,j$ such that $0\leq i+2j\leq 6,$ the regular component $y(x,t)$ satisfies

[TABLE]

and the singular component $z(x,t)$ satisfies

[TABLE]

where $m=\sqrt{\gamma}$.

Proof.

We first consider the interval $[0,1]\times[0,\tau]$. The data of the problems (2)-(2) are assumed to be sufficiently smooth and satisfy the appropriate compatibility conditions to ensure the existence of the unique solution $y_{0}$, $y_{1}$, $y_{2}$, $y$ $\in C_{\lambda}^{K,K/2}$, for $K=6$. Since $y_{0}$, $y_{1}$ and $y_{2}$ are solutions of first order hyperbolic equations (2), (2) and (2) over $Q_{1}$ as well we have the following estimates

[TABLE]

[TABLE]

and

[TABLE]

As $y_{3}$ is the solution of a problem similar to the initial boundary value problem (1.2)-(1) therefore, for all non-negative integer $i,j$ such that $0\leq i+2j\leq 4$, we have

[TABLE]

Using inequalities (2.11)-(2.14) we obtain the required estimates for the regular component $y(x,t)$ for $(x,t)\in\overline{Q}_{1}$.

To obtain the bound on the singular component, we define two barrier functions

[TABLE]

where $C$ is chosen sufficiently large such that we have

[TABLE]

Now,

[TABLE]

By Minimum principle, we have

[TABLE]

To obtain the bound on the derivatives of $z$ we transform the variable $x$ to the stretched variable $\widetilde{x}=\dfrac{x}{\sqrt{\varepsilon}}$. The transformed differential equation becomes independent of $\varepsilon$. For each neighbourhood $\widetilde{N}_{\delta}$ in $(2,1/\sqrt{\varepsilon})\times(0,\tau)$ using [9, §4.10], we have

[TABLE]

The required bounds can be obtained by transforming the inequality (2.15) in terms of the original variable $x$ and using the bound just obtained on $z(x,t)$. Similarly, for each neighbourhood $\widetilde{N}_{\delta}$ in $(0,2]\times(0,\tau)$ using [9, §4.10], we have

[TABLE]

Again transforming the inequality (2.16) in terms of the original variable $x$, using the bound on $z(x,t)$ and noting that $e^{-x/\sqrt{\varepsilon}}>C$ for $\widetilde{x}>2$, we have the required bounds. Next, consider the second interval $[\tau,2\tau]$. In this case we have $y(x,t)$ satisfy the problem (2) for $(x,t)\in Q_{2}$. The argument for rest of the proof is same as in the first case. The proof for $t\geq 2\tau$ also follows on the same lines. ∎

3 Discrete Problem

In this section, we discretize the problem (1.2)-(1) in both space and time direction. Firstly, a modified Shishkin mesh $S(L)$ is constructed to discretize the spatial domain.

Let $\overline{\Omega}^{N}:=\{x_{i}\}_{i=0}^{N}$ be the partition of the spatial domain $\overline{\Omega}$. We define the transition parameter $\sigma$ by

[TABLE]

where $L$ satisfies $\ln(\ln N)<L\leq\ln(N)$ and $e^{-L}\leq L/N$. The fitted piecewise uniform mesh $S(L)$ is constructed by dividing the domain $\overline{\Omega}$ into two subdomains $\overline{\Omega}=\overline{\Omega}_{1}\cup\overline{\Omega}_{2}$, where ${\Omega}_{1}=(0,\sigma]$ and $\Omega_{2}=(\sigma,1)$. A piecewise uniform mesh $\Omega_{\sigma}^{N}$ on $\Omega$ with $N$ mesh points is obtained by placing a uniform mesh with $N/2$ mesh points in each subintervals. The spatial step size $h_{i}=x_{i}-x_{i-1}$, for $i=1,2,\ldots,N$ is defined as

[TABLE]

where $h$ and $H$ are the spatial step size in $[0,\tau]$ and $(\tau,1],$ respectively.

For temporal discretization a uniform mesh $\Omega^{M}$ and $\Omega^{m_{\tau}}$ with $M$ and $m_{\tau}$ mesh points is considered by placing a uniform mesh with $M$ and $m_{\tau}$ mesh points in $[0,T]$ and $[-\tau,0)$, respectively. The uniform step size $\Delta t$ in time direction satisfies $\tau=m_{\tau}\Delta t$, where $m_{\tau}$ is a positive integer, $t_{n}=n\Delta t$, $n\geq-m_{\tau}$.

Piecewise uniform tensor product meshes $Q_{\sigma}^{N,M}$ on $Q$ and $\Gamma_{b,\sigma}^{N,M}$ on $\Gamma_{b}$ are defined as

[TABLE]

and the boundary points $\Gamma_{\sigma}^{N,M}$ of $Q_{\sigma}^{N,M}$ are defined as $\Gamma_{\sigma}^{N,M}=\overline{Q}^{N,M}\cap\Gamma$. We put $\Gamma_{l,\sigma}^{N,M}=\overline{Q}^{N,M}\cap\Gamma_{l}$ and $\Gamma_{r,\sigma}^{N,M}=\overline{Q}^{N,M}\cap\Gamma_{r}$. For $\sigma=1/2$, the mesh is uniform and for $\sigma=\sigma_{0}\sqrt{\varepsilon}L$ the mesh points get condensed at the left side of the domain.

3.1 The finite difference scheme

For any mesh function $v_{i}^{n}=v(x_{i},t_{n})$, the forward, backward and central difference operators $D_{x}^{+}$, $D_{x}^{-}$, $D_{x}^{0}$ in space and $D_{t}^{-}$ in time are defined as

[TABLE]

where $\widehat{h}_{i}=h_{i}+h_{i+1}$ for $i=1,\ldots,N-1$. We also define the second-order finite difference operator $\delta_{x}^{2}v_{i}^{n}$ in space by

[TABLE]

and $v_{i\pm 1/2}^{n}=\dfrac{v_{i\pm 1}^{n}+v_{i}^{n}}{2}$. Applying the central difference scheme in the interval $I=\{i\in\{1,2,...,N-1\},~{}a_{i}^{n}h_{i}<2\varepsilon\}$ and the mid-point upwind scheme in the remaining region, we get the following discrete problem

[TABLE]

where

[TABLE]

On simplifying the terms in the system of Eqns. (3.1), we obtain the following system of equations on the mesh $\overline{Q}_{\sigma}^{N,M}$

[TABLE]

where

[TABLE]

[TABLE]

and various coefficients are given by

[TABLE]

and

[TABLE]

Remark: As $a(x_{N/2},t_{n})>0$ and $a_{0}(x,t)\geq\alpha>0$ on $\overline{Q}$, we can conclude that there exists a constant $\kappa>0$ such that $a(x_{i},t_{n})\geq\kappa>0$, for $N/2\leq i\leq N$, $n\Delta t\leq T$.

Lemma 3.1**.**

Let $N_{0}$ be the smallest positive integer satisfying

[TABLE]

Then, for all $N\geq N_{0}$, we have

[TABLE]

Proof.

We first consider the case when $i\in I$. From (3.5) we clearly have

[TABLE]

and

[TABLE]

Also, from (3.5) we have

[TABLE]

Next, we consider the case when $i\notin I$. For $N\geq N_{0}$, where $N_{0}$ satisfies

[TABLE]

we have $a_{i}^{n}h_{i}<2\varepsilon$, for $i=1,2,\ldots,N/2$. It can be clearly seen that $\{1,2,\ldots,N/2\}\subset I$, for $N\geq N_{0}$ and therefore, $L_{\varepsilon,mu}^{N}$ is applied for $i>N/2$ where $i\notin I$.

Clearly, $r_{mu,i}^{-}=\dfrac{2\varepsilon\Delta t}{\widehat{h_{i}}h_{i}}>0$ and

[TABLE]

Using $\dfrac{1}{H}=\dfrac{N}{2(1-\tau)}\geq\dfrac{N_{0}}{2}$, we get

[TABLE]

Applying inequality (3.7) and taking $N\geq N_{0}$ we have $r^{+}_{mu,i}>0$. Also, from (3.6) we have

[TABLE]

From (3.5) and (3.6) we can easily get $|r_{cen,1}^{+}|<|r_{cen,1}^{0}|$, $|r_{cen,N-1}^{-}|<|r_{cen,N-1}^{0}|$ and $|r_{mu,N-1}^{-}|<|r_{mu,N-1}^{0}|$. ∎

Above Lemma establishes that the operator $L_{\varepsilon}^{N,M}$ satisfies the following discrete minimum principle

Lemma 3.2** (Discrete Minimum Principle).**

Let $W^{N}$ be any mesh function defined on $\overline{Q}_{\sigma}^{N,M}$. If $W^{N}(x_{i},t_{n})\geq 0$,  $\forall~{}(x_{i},t_{n})\in\Gamma_{\sigma}^{N,M}$ and $L_{\varepsilon}^{N,M}W^{N}(x_{i},t_{n})\leq 0$,  $\forall~{}(x_{i},t_{n})\in Q_{\sigma}^{N,M}$, then $W^{N}(x_{i},t_{n})\geq 0$,  $\forall~{}(x_{i},t_{n})\in\overline{Q}_{\sigma}^{N,M}$.

Lemma 3.3**.**

Let $W^{N}$ be any mesh function defined on $\overline{Q}_{\sigma}^{N,M}$. If $W^{N}(x_{i},t_{n})\geq 0,\thinspace\forall\thinspace(x_{i},t_{n})\in\Gamma_{\sigma}^{N,M}$ then

[TABLE]

Proof.

Constructing the following barrier function

[TABLE]

and using the discrete minimum principle we get the desired estimate. ∎

3.2 Error Analysis

In this section, we provide error estimates for the regular and the singular component of the numerical solution separately. Finally, they are combined to provide parameter uniform error estimates for the proposed hybrid scheme. To prove $\varepsilon$-uniform convergence of the proposed scheme we consider the barrier function

[TABLE]

where $\mu$ is a constant. Also,

[TABLE]

Lemma 3.4**.**

For each $0\leq i\leq N$ and $0<\mu<\dfrac{m}{2}$, the barrier function $\phi_{i}^{n}(\mu)$ satisfies the following inequalities

[TABLE]

Proof.

We first consider the case when $i\in I$. On applying the operator $L_{\varepsilon,cen}^{N,M}$ on the barrier function $\phi_{i}^{n}(\mu)$, we get

[TABLE]

where

[TABLE]

On Simplifying the Eqn. (3.11) using (3.9), we get

[TABLE]

Next, we consider the case when $i\notin I$. On applying the operator $L_{\varepsilon,mu}^{N,\Delta t}$ on the barrier function $\phi_{i}^{n}(\mu)$, we get

[TABLE]

where

[TABLE]

On simplifying Eqn. (3.12) using (3.9), we get

[TABLE]

∎

Lemma 3.5**.**

For each $0\leq i\leq N$ and $0<\mu<\dfrac{m}{2}$, we have the following inequalities

1. (i)

$\exp\left(\dfrac{-mx_{i}}{\sqrt{\varepsilon}}\right)\leq\phi_{i}^{n}(\mu),\quad 0\leq i\leq N,\quad n\Delta t\leq T;$** 2. (ii)

$\phi_{i}^{n}\leq\begin{cases}CL^{\displaystyle 2\mu\sigma_{0}\left(\frac{i}{N}\right)}N^{\displaystyle-2\mu\sigma_{0}\left(\frac{i}{N}\right)},\quad 1\leq i\leq\frac{N}{2},\quad n\Delta t\leq T,\\ CL^{\displaystyle\mu\sigma_{0}}N^{\displaystyle-\mu\sigma_{0}},\quad\frac{N}{2}\leq i\leq N,\quad n\Delta t\leq T.\end{cases}$**

Proof.

1. (i)

Using $e^{-x}=(1+x)^{-1},~{}x\geq 0,$ we obtain the desired inequality. 2. (ii)

Considering the barrier function $\phi_{i}^{n}(\mu)$ for $i\in\{1,2,\ldots,\frac{N}{2}\}$, we have

[TABLE]

Using $e^{-L}\leq L/N$, we get

[TABLE]

Next, consider the barrier function $\phi_{i}^{n}(\mu)$ for $i\in\{\frac{N}{2},\ldots,N\}$, we have

[TABLE]

Using $e^{-L}\leq L/N$, we get

[TABLE]

∎

To obtain $\varepsilon$-uniform error estimate we decompose the numerical solution $U_{i}^{n}=U(x_{i},t_{n})$ of the discrete problem (3.1) into a regular part $Y_{i}^{n}$ and a singular part $Z_{i}^{n}$ analogously to the decomposition of the continuous solution $u(x,t)$ as:

[TABLE]

where $Y_{i}^{n}$ satisfies the following non-homogeneous problem

[TABLE]

and $Z_{i}^{n}$ satisfies the following homogeneous problem

[TABLE]

As a result the pointwise error at the node $(x_{i},t_{n})$ in the discrete solution can be decomposed as

[TABLE]

Lemma 3.6** (Error in the Regular Component ).**

Under the assumption (3.7) of Lemma 3.1, the regular component at each mesh points $(x_{i},t_{n})\in\overline{Q}^{N,M}_{\sigma}$, satisfies the following error estimate:

[TABLE]

Proof.

On the interval $[0,\tau]$, the right-hand side of (1.2) becomes $f(x,t)+e(x,t)s(x,t-\tau)$ which is known and is independent of $\varepsilon$. We will consider two cases depending upon the relation between $\varepsilon$ and $N$:

1. Case (i)

When $\varepsilon>\|a\|_{\overline{Q}}/N$. In this case we have $a_{i}h_{i}<2\varepsilon$ for all $i\in\{1,\ldots,N-1\}$ which implies the set $\{1,\ldots,N-1\}\subseteq I$. We get

[TABLE]

Using $h_{i+1}+h_{i}\leq 2N^{-1}$ and the bounds on the derivatives of $y$ given in Theorem 2.1, we get

[TABLE]

Using Lemma 3.3 we can obtain the desired result. 2. Case (ii)

When $\varepsilon\leq\|a\|_{\overline{Q}}/N$. For the smooth component the truncation error is defined as

[TABLE]

Using $h_{i+1}+h_{i}\leq 2N^{-1}$ and the bounds on the derivatives of $y(x,t)$ given in Theorem 2.1, we get

[TABLE]

Using $\sqrt{\varepsilon}\leq C/N$, we get

[TABLE]

where $(x_{i},t_{n})\in Q_{\sigma,~{}\tau}^{N,M}=\Omega_{\sigma}^{N}\times\Omega_{1}^{M}$ ($\Omega_{1}^{M}$ is the uniform mesh with $M=m_{\tau}$ mesh elements in the interval $[0,\tau]$). Applying Lemma 3.3 we can obtain the desired result for the considered case.

Combining both the cases, we get

[TABLE]

The regular part $Y_{i}^{n}=Y(x_{i},t_{n})$ of the numerical solution on $Q_{\sigma,~{}\tau}^{N,M}$ is denoted by $Y_{\tau,i}^{n}=Y_{\tau}(x_{i},t_{n})$.

Next, we consider the second interval $[\tau,2\tau]$. On the second interval $[\tau,2\tau]$, the delay term $u(x,t-\tau)$ is the numerical solution obtained in the first interval $[0,\tau]$. To do the error analysis in the interval $[\tau,2\tau]$ consider the problem:-

[TABLE]

The regular component $Y(x_{i},t_{n})$ of the numerical solution at mesh point $(x_{i},t_{n})\in Q_{\sigma,~{}2\tau}^{N,M}=\Omega_{\sigma}^{N}\times\Omega_{2}^{M}$ ($\Omega_{2}^{M}$ is the uniform mesh with $M=m_{\tau}$ mesh elements in the interval $[\tau,2\tau]$) is determined by

[TABLE]

We observe that

[TABLE]

On using inequality (3.14) and again considering two cases when $(i)\varepsilon>\|a\|_{\overline{Q}}/N$ and when $(ii)\varepsilon\leq\|a\|_{\overline{Q}}/N$, we get

[TABLE]

Application of Lemma 3.3 gives us the desired result in the interval $[\tau,2\tau]$. Similarly, we can prove the result for $t\in[2\tau,3\tau]$ and so on. ∎

Lemma 3.7** (Error in the Singular Component).**

Under the assumption (3.7) and $0<\mu<\frac{m}{2}$ the following error estimate is satisfied by the singular component $Z(x_{i},t_{n})$ at each mesh points $(x_{i},t_{n})\in\overline{Q}^{N,M}_{\sigma}$,

[TABLE]

Proof.

We proceed by method of steps. Firstly, we compute the error in the interval $[0,\tau]$ and then consider the interval $[\tau,2\tau]$.

Case (i): In this case, for $t\in[0,\tau]$ the right hand side of Eqn. (1.2) is known and independent of $\varepsilon$. We first consider the outer region $(\sigma,1]\times(0,\tau]$. We know that $Z$ and $z$ are small in the outer region irrespective of the fact that whether $i\in I$ or $i\notin I$.

Consider the barrier functions

[TABLE]

where $C=|z(x_{0},t_{n})|$. We observe that

[TABLE]

Using the discrete minimum principle, we get

[TABLE]

Using triangle inequality and Lemma 3.5, we get

[TABLE]

Using the bounds of $\phi_{i}^{n}(\mu)$ given in Lemma 3.5, we get

[TABLE]

On taking $\sigma_{0}=\dfrac{2}{\mu}$, we get

[TABLE]

Next, we consider the inner region $[0,\sigma]\times(0,\tau]$. For $N\geq N_{0}$, satisfying (3.9) we have $a_{i}^{n}h_{i}<2{\varepsilon}$, for all $i=1,\ldots,N/2$ which implies $\{1,\ldots,N/2\}\subseteq I$. Therefore, $L_{\varepsilon,cen}^{N,M}$ is applied in the region $[0,\sigma]\times[0,\tau]$. We have,

[TABLE]

Since for inner region $i\in I$ so, using (3.9) we get $mh<\sqrt{\varepsilon}$ and sin$h\xi\leq C\xi$, for $0\leq\xi\leq 2$. This implies sinh$\left(\frac{mh}{\sqrt{\varepsilon}}\right)\leq\frac{Cmh}{\sqrt{\varepsilon}}$. Therefore, (3.2) becomes

[TABLE]

Also, in the inner region, we have $h_{i}=h$, so

[TABLE]

We have $|(Z-z)(x_{0},t_{n})|=0$, $\forall~{}t_{n}\in\Omega_{1}^{M}$ and $|(Z-z)(x_{i},t_{n})|=0$, $\forall~{}(x_{i},t_{n})\in\Gamma_{b,\sigma}^{N}$.

Also, from (3.15) we have

[TABLE]

Considering the barrier functions,

[TABLE]

we observe that $\Psi^{\pm}(x_{N/2},t_{n})>0$, $\Psi^{\pm}(x_{N},t_{n})>0$, $\forall~{}t_{n}\in\Omega_{1}^{M}$ and $\Psi^{\pm}(x_{i},t_{n})>0$, $\forall~{}(x_{i},t_{n})\in\Gamma_{b,\sigma}^{N}$. Using Lemma 3.4 we have,

[TABLE]

From discrete minimum principle we get,

[TABLE]

The singular part $Z(x_{i},t_{n})$ of the numerical solution on $Q_{\sigma,~{}\tau}^{N,M}$ is denoted by $Z_{\tau}(x_{i},t_{n})$.

Case (ii): On the second interval $[\tau,2\tau]$, the delay term $u(x,t-\tau)$ is the numerical solution obtained in the first interval $[0,\tau]$. We will do the error analysis over the interval $[\tau,2\tau]$ in the following way. We will consider the singularly perturbed delay parabolic partial differential equation (2) on the second interval $[\tau,2\tau]$. The singular component $Z(x_{i},t_{n})$ of the numerical solution at mesh points $(x_{i},t_{n})\in Q_{\sigma,~{}2\tau}^{N,M}=\Omega_{\sigma}^{N}\times\Omega_{2}^{M}$ is determined by

[TABLE]

We observe that

[TABLE]

First term on the right hand side can be approximated using (3.2) and the second term can be approximated using the same approach as discussed in the first interval. The proof is completed by introducing the barrier functions and applying the discrete minimum principle as in the previous case for $t\in[0,\tau]$.

[TABLE]

Similarly, the case, for $t\geq 2\tau$ also follows on the same lines. ∎

Theorem 3.1**.**

Let $u(x_{i},t_{n})$ be the exact solution of the problem (1.2)-(1) and $U(x_{i},t_{n})$ be the discrete solution of the system of Eqns.(3.1) at each mesh point $(x_{i},t_{n})\in\overline{Q}_{\sigma}^{N,M}$. Then under the assumption of Lemma 3.4 and $0<\mu<m/2$ for $N\geq N_{0}$ we have

[TABLE]

4 Richardson extrapolation

In this section, the Richardson extrapolation technique is used to obtain higher accuracy and order of convergence in time direction. We consider two meshes $\overline{Q}_{\sigma}^{N,M}=\overline{\Omega}_{\sigma}^{N}\times\overline{\Omega}^{M}$ and $\overline{Q}_{\sigma}^{N,2M}=\overline{\Omega}_{\sigma}^{N}\times\overline{\Omega}^{2M}$ where $\overline{\Omega}^{M}$ and $\overline{\Omega}^{2M}$ are uniform meshes with $M$ and $2M$ mesh points, respectively in the temporal direction. Both the considered meshes have same number of mesh points in the spatial direction. Let

[TABLE]

Then, $\overline{Q}_{\sigma,0}^{N,M}=\overline{Q}_{\sigma}^{N,M}$, as $\overline{Q}_{\sigma}^{N,M}\subseteq\overline{Q}_{\sigma}^{N,2M}$. Let $U^{k}$ denote the numerical solution of the problem (3.1) on the mesh $\overline{Q}_{\sigma}^{N,kM}=\overline{\Omega}_{\sigma}^{N}\times\overline{\Omega}^{kM}$ where $k=1,2$. Then, we approximate $u(x,t)$ by $U_{ext}(x_{i},t_{n})$ where

[TABLE]

The numerical approximation $U_{ext}$ has improved order of convergence in time. To verify this we use a technique similar to [13]. We have

[TABLE]

where $R_{n}^{k}$, $k=1,2$ is the remainder term and $\xi^{k}$ is the solution of the following problem:

[TABLE]

We need to derive the estimates for the remainder term $R^{k}$ on $\overline{Q}_{\sigma}^{N,kM}$, $k=1,2$. Also $R^{k}(x_{i},t_{n})=0$, $\forall(x_{i},t_{n})\in\Gamma_{\sigma}^{N,kM}$, where $\Gamma_{\sigma}^{N,kM}$, $k=1,2$ is the boundary of $\overline{Q}^{N,kM}$. We have

[TABLE]

From discrete minimum principle, we get

[TABLE]

Theorem 4.1**.**

Let $u(x_{i},t_{n})$ be the exact solution of the problem. (1.2)-(1) and $U_{ext}(x_{i},t_{n})$ be the discrete solution obtained using the Richardson extrapolation at each mesh point $(x_{i},t_{n})\in\overline{Q}_{\sigma}^{N,M}$. Then, for $N\geq N_{0}$ where $N_{0}$ satisfies the assumption (3.7) and $0<\mu<m/2$, we have the following $\varepsilon$-uniform error estimate

[TABLE]

5 Numerical results

In this section, we present the numerical results for two test problems to validate the theoretical results. They also verify the high accuracy and convergence rate of the proposed schemes.

Problem 1**.**

Consider the following singularly perturbed parabolic IBVP :-

[TABLE]

Problem 2**.**

Consider the following singularly perturbed parabolic IBVP :-

[TABLE]

Since exact solutions of the given problems are not known, the performances of the proposed schemes are illustrated by using the double mesh principle to calculate the maximum pointwise error. The maximum pointwise error is defined as

[TABLE]

where $\widetilde{U}=U$ for hybrid scheme (3.1) and $\widetilde{U}=U_{ext}$ if Richardson extrapolation is applied on the scheme (3.1). The corresponding order of convergence $q_{\varepsilon}^{N,M}$ is computed as

[TABLE]

Also, the $\varepsilon$-uniform maximum pointwise error $E^{N,M}$ is computed as

[TABLE]

and the corresponding $\varepsilon$-uniform order of convergence $q^{N,M}$ is given by

[TABLE]

For various values of $\varepsilon$, $N$ and $M$ the computed maximum pointwise errors $E_{\varepsilon}^{N,M}$ and the corresponding order of convergence $q_{\varepsilon}^{N,M}$ for the considered problems are tabulated in Tables 1 to 6.

In Table 1 and 3 we have given the results for upwind scheme on uniform mesh, upwind scheme on piecewise uniform Shishkin mesh and hybrid scheme (3.1) on modified Shishkin mesh for Problem 1, 2, respectively. It can be seen that the uniform mesh do not work. The upwind scheme on Shishkin mesh has almost first order of convergence. The proposed hybrid scheme gives better result than the upwind scheme with Shishkin mesh. The numerical results computed using hybrid scheme show monotonically decreasing behaviour as $N$ increases which confirms the $\varepsilon$-uniform convergence of the hybrid scheme (3.1). The order of convergence of hybrid scheme is not depicting the theoretical order of convergence of order two upto a logarithmic factor as proved in Theorem 3.1 as the error consists of two parts due to spatial and temporal discretization. The hybrid scheme improves the accuracy space but temporal order of convergence remains one. As a result, the numerical results display almost first order of convergence due to the much influence of the temporal error. Therefore, Richardson extrapolation is used to increase the order of convergence in time.

Tables 2 and 4 show that the use of Richardson extrapolation on hybrid scheme further improves the error and the rate of convergence. The resulting scheme provides almost second order of convergence in space and second order of convergence in time variable.

Tables 5 and 6 display the numerical results computed using Richardson extrapolation for Problem 1 and 2, respectively, with different values of $p$.

6 Conclusions

In this article, we proposed and analysed a higher order numerical scheme for the solution of singularly perturbed parabolic problems with time delay and degenerate coefficient. The proposed scheme is comprised of implicit Euler scheme for time discretization on uniform mesh and a hybrid scheme for space discretization on modified Shishkin mesh. Parameter uniform convergence of order one in time direction and order two upto a logarithmic factor in space direction is established for the proposed scheme (3.1). Further, to improve the order of convergence in time direction Richardson extrapolation is employed in the time direction. The resulting scheme increase the order of convergence to two in time direction. The numerical experiments are presented to illustrate the theoretical results.

Acknowledgements

The first author wish to acknowledge UGC Non-NET Fellowships for financial support vide Ref.No. Sch./139/Non-NET/Ext-152/2018-19/67.