Skip to main content

An accurate Chebyshev pseudospectral scheme for multi-dimensional parabolic problems with time delays

Abstract

In this paper, the Chebyshev Gauss-Lobatto pseudospectral scheme is investigated in spatial directions for solving one-dimensional, coupled, and two-dimensional parabolic partial differential equations with time delays. For the one-dimensional problem, the spatial integration is discretized by the Chebyshev pseudospectral scheme with Gauss-Lobatto quadrature nodes to provide a delay system of ordinary differential equations. The time integration of the reduced system in temporal direction is implemented by the continuous Runge-Kutta scheme. In addition, the present algorithm is extended to solve the coupled time delay parabolic equations. We also develop an efficient numerical algorithm based on the Chebyshev pseudospectral algorithm to obtain the two spatial variables in solving the two-dimensional time delay parabolic equations. This algorithm possesses spectral accuracy in the spatial directions. The obtained numerical results show the effectiveness and highly accuracy of the present algorithms for solving one-dimensional and two-dimensional partial differential equations.

1 Introduction

In recent years there has been a high level of interest in employing spectral methods for numerically solving many types of integral and differential equations, due to their high accuracy and ease of applying them for finite and infinite domains [1–9]. The spectral collocation method is a specific type of spectral methods that is more applicable and widely used to solve most types of differential equations [10–13].

In biology, physics, engineering, and computer design, many models can be attributed to time-delay partial differential equations (PDEs) (see, e.g., [14–17]). In the last few years, various analytical and numerical methods have been proposed for solving delay integral and PDEs (see, e.g., [18–20]). The Legendre pseudospectral algorithm [21] was implemented successfully to provide very accurate solutions of parabolic integro-differential equations in bounded and unbounded intervals. Kaushik et al. [22] proposed the uniform difference scheme for time-delayed PDEs. In the same direction, the singularly perturbed delay PDEs have been solved by an operator finite difference scheme [23]. The Chebyshev wavelet method was improved in [24] to obtain numerical solutions of time-varying delay systems.

Recently, Rashidinia et al. [25] proposed an efficient numerical technique for solving parabolic convection-diffusion problems, in which a finite difference scheme was used in the temporal direction to reduce the problem to a system of ordinary differential equations (ODEs) in a spatial direction, and the Sinc-Galerkin scheme was then used to solve this system. The Bessel collocation approximation was investigated in [26] to reduce the parabolic convection-diffusion equation to a system of algebraic equations. Moreover, Bhrawy et al. [27] proposed the Jacobi pseudospectral scheme for numerically solving the time-delayed Burgers’ equations and the same method was applied successfully for the generalized Fitzhugh-Nagumo equation in [28]; additionally, the Jacobi Gauss-Lobatto pseudospectral scheme has been proposed and extended to solve the complex generalized Zakharov system [29]. Ashyralyev and Agirseven [30] presented a difference scheme for the delay parabolic differential equation and studied the convergence of this scheme. Furthermore, Tian [31] discussed the asymptotic stability and convergence of three difference schemes which were applied to solve delay parabolic PDEs. More recently, the authors of [32] and [33] proposed a linearized compact multi-splitting method and compact difference method combined with an extrapolation scheme to solve convection-reaction-diffusion and neutral parabolic PDEs with delay, respectively.

As pointed out above, most researchers used finite difference schemes for the spatial discretizations of time-delay partial differential equations to get a system of ODEs with time delays. However, the accuracy of such methods is poor in the spatial directions. This motivated our interest to investigate the Chebyshev pseudospectral method for spatial discretizations for the initial-boundary value problems with time delays in one and two dimensions. This method is known for its ease in implementation along with the high accuracy and exponential convergence that can be achieved.

The main aim of this article is to investigate the Chebyshev pseudospectral scheme to solve the one-dimensional parabolic PDEs. We focus primarily on implementing this scheme in spatial independent variable. The spatial direction is collocated at \((N-1)\) collocation points of the Chebyshev Gauss-Lobatto quadrature nodes. This scheme has the advantage of reducing the one-dimensional parabolic PDEs into a system of \((N-1)\) ODEs with time delays in the time direction, that can be solved by continuous Runge-Kutta (RK) method.

We also extend the application of this algorithm to solve the coupled time delay parabolic equations at \((N-1)\) collocation nodes, which provides a system of \((2N-2)\) ODEs with time delays. Moreover, this algorithm is developed to solve the two-dimensional time delay parabolic equations, in which the two spatial variables are collocated at \((N-1)\times (M-1)\) Chebyshev Gauss-Lobatto quadrature nodes. This provides a system of \((N-1)\times(M-1)\) ODEs with time delays. Finally, several numerical simulations are given to confirm the high accuracy of the proposed algorithm.

This paper is outlined as follows: In the next section, we present some properties of Chebyshev polynomials. The Chebyshev pseudospectral approximation in spatial directions is obtained for solving the one-dimensional time-dependent PDEs, coupled parabolic PDEs, and two-dimensional time-dependent parabolic PDEs with time delays. In Section 4, some illustrative numerical experiments are given and some comparisons are made between our method and other methods. The paper ends with some conclusions and observations in Section 5.

2 Chebyshev polynomials interpolation

In this section, we present some properties of Chebyshev polynomials that will be used to construct the method. The Chebyshev polynomials are determined from the following recurrence formula:

$$T_{i+1}(x)=2 x T_{i}(x)-T_{i-1}(x), \quad i=1,2, \ldots,x \in [-1,1], $$

where \(T_{0}(x)=1\) and \(T_{1}(x)=x\). Also

$$ T_{k}(-1)=(-1)^{k}, \qquad T_{k}(1)=1. $$
(1)

Let \(w(x) = (1-x^{2})^{-\frac{1}{2}}\), and \(\mathrm{I}\equiv[-1,1]\), then we define the space \(L^{2}_{w}(\mathrm{I})\)-orthogonal system, i.e.,

$$ \int_{\mathrm{I}} T_{k}(x)T_{l}(x)w(x) \,dx=h_{k} \delta_{k,l}, $$
(2)

where

$$ h_{k}=\frac{c_{k}}{2}\pi, \qquad c_{0}=2, \qquad c_{k}=1, \quad k\geq1. $$
(3)

In the forthcoming discussions, we denote the norm and semi-norm of the weighted Sobolev space \(H_{w}^{r}(\mathrm{I})\) by \(\|\upsilon\|_{r,w,\mathrm{I}}\) and \(|\upsilon|_{r,w,\mathrm{I}}\), respectively. In particular, \(L_{w}^{2}(\mathrm{I})=H_{w}^{0}(\mathrm{I})\) and \(\|\upsilon\|_{w,\mathrm{I}}=\|\upsilon\|_{0,w,\mathrm{I}}\).

Let \(\zeta_{N,j}\), \(0\leqslant j\leqslant N\), be the nodes of the Chebyshev-Gauss interpolation over I. Their corresponding Christoffel numbers are \(\varpi_{N,j}\), \(0\leqslant j\leqslant N\), and

$$\varpi_{N,j}=\frac{\pi}{{\gamma}_{j} N}, \quad j=0,1,\dots,N, $$

where \({\gamma}_{0}={\gamma}_{N}=2\) and \({\gamma}_{1}\) for \(j=1,2,\dots,N-1\). Let \(S_{N}\) be the space of all polynomials of degree ≤N, then for any \(\phi\in S_{2N-1}(0,L)\),

$$ \int_{-1}^{1}\omega(x) \phi(x)\,dx= \sum_{j=0}^{N}\varpi_{N,j} \phi ( \zeta_{N,j}). $$
(4)

Let us define the following discrete inner product and norm:

$$ (u,v)_{w,N}=\sum_{j=0}^{N}u( \zeta_{N,j} ) v(\zeta_{N,j}) \varpi_{N,j},\qquad \| u \|_{w,N}=\sqrt{(u,u)_{w,N}}. $$
(5)

From [34], the following relation holds:

$$ \| u\|_{w}\leq\| u \|_{N,w}\leq\sqrt{2} \| u \|_{w},\quad \forall u\in S_{N}. $$
(6)

Theorem 2.1

Let \(u(x)\in H^{k}(\mathrm{I})\) Sobolev space and \(u_{N}(x)=\sum_{i=0}^{N}a_{i}T_{i}(x)\) be the best approximation polynomial of \(u(x)\) in \(L_{w}^{2}\)-norm. Thus, the truncation error is

$$ \bigl\| u(x)-u_{N}(x)\bigr\| _{L_{w}^{2}(\mathrm{I})}\leq C_{0}N^{-k} \bigl\| u(x)\bigr\| _{H^{k}(\mathrm{I})}, $$
(7)

where \(C_{0}\) is a positive constant, which is dependent on the selected norm and independent of \(u(x)\) and N (for a proof see [1]).

3 Chebyshev pseudospectral scheme

We propose a pseudospectral algorithm based on Chebyshev polynomials to integrate the spatial variable for the \((1+1)\) parabolic PDEs with discrete time delay. The problem is then transformed into a system of ODEs with time delay. The algorithm is extended to treat the coupled parabolic PDEs and \((2+1)\) parabolic PDEs with time delay. The spatial variables are integrated based on the Chebyshev pseudospectral scheme, in which we use the Gauss-Lobatto interpolation nodes, to convert the problems to those of ODEs with time delay. The continuous RK scheme is employed to solve the resulting systems of ODEs.

3.1 \((1+1)\) Parabolic PDEs with time delay

Consider the \((1+1)\) parabolic PDEs with time delay of the form

$$ \frac{\partial u(x,t)}{\partial{t}}= \eta \frac{\partial^{2} u(x,t)}{\partial{x^{2}}}+ \lambda_{1} u(x,t-\tau)+\lambda_{2} u^{2}(x,t-\tau) +f(x,t),\quad(x, t)\in\mathrm{I} \times\mathrm{T}, $$
(8)

with the boundary conditions

$$ u(-1,t)=c_{1}(t),\qquad u(1,t)=c_{2}(t),\quad t \in\mathrm{T}, $$
(9)

and the initial state

$$ u(x,t)=c_{3}(x,t),\quad (x, t)\in \mathrm{I} \times[-\tau,0], $$
(10)

where \(\mathrm{I}\equiv[-1,1]\), \(\mathrm{T}\equiv(0,T]\), and η and \(\tau>0\) are the diffusion coefficient and delay parameter.

The polynomial approximation of degree N to \(u(x,t)\) may be expressed in terms of the orthogonal family \(\{T_{i}(x)\} \) in the form

$$ u(x,t)=\sum_{i=0}^{N}a_{i}(t) T_{i}(x). $$
(11)

It follows from (4) and (5) that

$$ a_{i}(t)=\frac{1}{h_{i}} \int_{-1}^{1}u(x,t) w(x) T_{i}(x)\,dx. $$
(12)

It can be seen directly from (4) that \(a_{j}(t)\) may be expanded in the form

$$a_{i}(t)=\frac{1}{h_{i}} \sum _{j=0}^{N}T_{i}( \zeta_{N,j}) \varpi_{N,j} u(\zeta_{N,j},t), $$
(13)

where \(\zeta_{N,j}\) (\(0\leq j\leq N\)) are the zeros and \(\varpi_{N,j}\) (\(0\leq j\leq N\)) are the corresponding quadrature weights.

We can further rewrite (11) as

$$ u(x,t)=\sum_{j=0}^{N} \Biggl(\sum_{i=0}^{N}\frac{1}{h_{i}}T_{i} (\zeta_{N,j})T_{i}(x)\varpi_{N,j} \Biggr) u( \zeta_{N,j},t). $$
(14)

The first-order spatial partial derivative at a specific collocation node \(\zeta_{N,n}\) can be obtained from (14) as

$$ u_{x}(\zeta_{N,n},t)= \sum _{i=0}^{N}\mu_{ni}u( \zeta_{N,i},t),\quad 0\leq n \leq N, $$
(15)

where

$$ \mu_{ni}=\sum_{j=0}^{N} \frac{1}{h_{j}}T_{j} ( \zeta_{N,i})\partial_{x} \bigl(T_{j}( \zeta_{N,n}) \bigr) \varpi_{N,i}. $$
(16)

This result can be extended to compute the second-order spatial partial derivative at a specific collocation node \(\zeta_{N,n}\) as

$$ u_{xx}(\zeta_{N,n},t)= \sum _{i=0}^{N}\gamma_{ni}u( \zeta_{N,i},t),\quad 0\leq n \leq N, $$
(17)

where

$$ \gamma_{ni}=\sum _{j=0}^{N}\frac{1}{h_{j}}T_{j} ( \zeta_{N,i})\partial_{xx} \bigl(T_{j}( \zeta_{N,n}) \bigr) \varpi_{N,i}. $$
(18)

In the context of Chebyshev pseudospectral approximation, putting (14) and (17) in (8) gives

$$\begin{aligned} &\dot{u}_{n}(t)=\eta \sum _{i=0}^{N}\gamma_{ni}u_{i}(t)+ \lambda_{1} u_{n}(t-\tau)+\lambda_{2} u^{2}_{n}(t-\tau)+f_{n}(t), \\ &\quad 1\leq n \leq N-1, \end{aligned}$$
(19)

where

$$\begin{aligned}& u_{k}(t)=u(\zeta_{N,k},t),\qquad u_{k}(t- \tau)=u(\zeta _{N,k},-\tau), \\& f_{k}(t)=f(\zeta_{N,k},t),\quad 1\leq k \leq N-1. \end{aligned}$$

The boundary conditions (9) are satisfied exactly at the two collocation points \(\zeta_{N,0}=-1\) and \(\zeta_{N,N}=1\). Furthermore, the preceding equation provides a system of \((N-1)\) ODEs with discrete time delay Ï„, namely

$$ \dot{u}_{n}(t)=\eta \sum _{i=1}^{N-1}\gamma_{ni}u_{i}(t)+ \rho_{n}(t)+\lambda_{1} u_{n}(t-\tau)+ \lambda_{2} u^{2}_{n}(t-\tau)+f_{n}(t), $$
(20)

subject to

$$ u_{n}(t)=c_{3}( \zeta_{N,n},t),\quad n=1,\ldots, N-1, t\in[-\tau,0], $$
(21)

where

$$\rho_{n}(t)=\gamma_{n0} c_{1}(t)+ \gamma_{nN} c_{2}(t). $$

Finally, the system (20)-(21) can be written in a matrix form as

$$ \begin{aligned} &\dot{\mathbf{U}}(t)=\mathbf{G} \bigl(t,u(t), u(t-\tau) \bigr), \\ &\mathbf{U}(0)=\mathbf{c},\quad t\in[-\tau,0], \end{aligned} $$
(22)

where

$$\begin{aligned}& \dot{\mathbf{U}}(t)= \bigl[\dot{u}_{1}(t),\dot{u}_{2}(t), \ldots,\dot{u}_{N-1}(t) \bigr]^{T}, \\& \mathbf{c}= \bigl[c_{3}(\zeta_{N,1}),c_{3}( \zeta_{N,2}),\ldots,c_{3}(\zeta _{N,N-1}) \bigr]^{T}, \\& \mathbf{G} \bigl(t,u(t),u(t-\tau) \bigr)= \bigl[g_{1} \bigl(t,u(t), u(t- \tau) \bigr),\ldots ,g_{N-1} \bigl(t,u(t), u(t-\tau) \bigr) \bigr]^{T}, \end{aligned}$$

and

$$g_{i} \bigl(t,u(t), u(t-\tau) \bigr)=\eta \sum _{j=1}^{N-1} \gamma _{ij}u_{j}(t)+ \rho_{i}(t)+ \lambda_{1} u_{i}(t-\tau)+ \lambda_{2} u^{2}_{i}(t-\tau)+f_{i}(t). $$

The previous system of first-order ODEs with time delay (22) is integrated by the continuous RK scheme [35, 36].

3.2 Coupled \((1+1)\) parabolic PDEs with time delay

The main objective of this section is to develop the Chebyshev pseudospectral method to numerically solve the coupled parabolic PDEs with discrete time delay. This converts the following coupled parabolic PDEs with time delay into system of ODEs with time delay:

$$ \left . \textstyle\begin{array}{l@{}} \frac{\partial u(x,t)}{\partial{t}}= \eta_{1} \frac{\partial^{2} u(x,t)}{\partial{x^{2}}} +u(x,t) (1+\lambda_{1} u(x,t- \tau)+\kappa_{1} v(x,t-\tau) )+r(x,t),\\ \frac{\partial v(x,t)}{\partial{t}}=\eta_{2} \frac{\partial^{2} v(x,t)}{\partial{x^{2}}} +v(x,t) (1+ \lambda_{2} u(x,t-\tau)+\kappa_{2} v(x,t-\tau) )+s(x,t),\\ \quad (x,t)\in \mathrm{I} \times\mathrm{T}, \end{array}\displaystyle \right \} $$
(23)

subject to

$$ \begin{aligned} &u(-1,t)=c_{1}(t),\qquad u(1,t)=c_{2}(t), \\ &v(-1,t)=c_{3}(t),\qquad v(1,t)=c_{4}(t),\quad t\in \mathrm{T}, \end{aligned} $$
(24)

and initial data

$$ u(x,t)=c_{5}(x,t),\qquad v(x,t)=c_{6}(x,t), \quad (x,t) \in\mathrm{I} \times[-\tau,0]. $$
(25)

Now, we outline the main step of the Chebyshev pseudospectral algorithm for the coupled parabolic PDEs with time delay. Let us expand the approximate solutions as

$$ \left . \textstyle\begin{array}{l@{}} u(x,t)=\sum_{i=0}^{N}a_{i}(t) T_{i}(x),\\ v(x,t)=\sum_{i=0}^{N}b_{i}(t) T_{i}(x). \end{array}\displaystyle \right \} $$
(26)

The expansion coefficients \(\{a_{j}(t)\mbox{ and }b_{j}(t)\}\), can be described in terms of the solution at the Chebyshev Gauss-Lobatto quadrature points \(\{\zeta_{N,j}\}\), \(j=0,1,\ldots,N\), as

$$ \left . \textstyle\begin{array}{l@{}} a_{i}(t)= \frac{1}{h_{i}}\sum_{j=0}^{N}T_{i}( \zeta_{N,j}) \varpi_{N,j} u(\zeta_{N,j},t),\\ b_{i}(t)=\frac{1}{h_{i}}\sum_{j=0}^{N}T_{i} (\zeta_{N,j})\varpi_{N,j} v(\zeta_{N,j},t). \end{array}\displaystyle \right \} $$
(27)

Therefore, the approximate solutions (26) can be expressed in terms of Chebyshev polynomials, at the Chebyshev Gauss-Lobatto quadrature points, in the form

$$ \left . \textstyle\begin{array}{l@{}} u(x,t)=\sum_{j=0}^{N} (\sum_{i=0}^{N}\frac{1}{h_{i}}T_{i} (\zeta_{N,j})T_{i}(x)\varpi_{N,j} ) u( \zeta_{N,j},t),\\ v(x,t)=\sum_{j=0}^{N} (\sum_{i=0}^{N}\frac {1}{h_{i}}T_{i} ( \zeta_{N,j})T_{i}(x)\varpi_{N,j} ) v( \zeta_{N,j},t). \end{array}\displaystyle \right \} $$
(28)

The first-order partial derivative with respect to x for the solutions (28), at a specific collocation node \(\zeta_{N,n}\) can be written as

$$ \left . \textstyle\begin{array}{l@{}} u_{x}( \zeta_{N,i},t)=\sum_{j=0}^{N} \mu_{ij}u(\zeta _{N,j},t),\\ v_{x}(\zeta_{N,i},t)=\sum_{j=0}^{N} \mu_{ij}v(\zeta_{N,j},t), \quad i=0,1,\ldots,N . \end{array}\displaystyle \right \} $$
(29)

According to the previous step, the second-order spatial derivatives of the approximate solutions (28) are

$$ \left . \textstyle\begin{array}{l@{}} u_{xx}( \zeta_{N,i},t)=\sum_{j=0}^{N} \gamma_{ij}u(\zeta_{N,j},t),\\ v_{xx}(\zeta_{N,i},t)=\sum_{j=0}^{N} \gamma_{ij}v(\zeta_{N,j},t),\quad i=0,1,\ldots,N. \end{array}\displaystyle \right \} $$
(30)

Applying the Chebyshev pseudospectral approximation [37, 38] for the coupled parabolic PDEs with discrete time delay (23), at the nodes of Chebyshev-Gauss-Lobatto quadrature. Thanks to (28)-(30), the coupled equations (23) may be written in the form

$$ \left . \textstyle\begin{array}{l@{}} \dot{u}_{n}(t)= \eta_{1} \sum_{i=0}^{N} \gamma_{ni}u_{i}(t)+u_{n}(t) (1+ \lambda_{1} u_{n}(t-\tau)+\kappa_{1} v_{n}(t-\tau) )+r_{n}(t),\\ \dot{v}_{n}(t)=\eta_{2} \sum_{i=0}^{N} \gamma _{ni}v_{i}(t)+v_{n}(t) (1+ \lambda_{2} u_{n}(t-\tau)+\kappa_{2} v_{n}(t-\tau) )+s_{n}(t), \end{array}\displaystyle \right \} $$
(31)

where

$$\begin{aligned}& u_{k}(t)=u(\zeta_{N,k},t),\qquad v_{k}(t)=v( \zeta_{N,k},t), \\& r_{k}(t)=r(\zeta_{N,k},t), \qquad s_{k}(t)=s( \zeta_{N,k},t), \quad k=1,\ldots,N-1, n=1,\ldots,N-1. \end{aligned}$$

The boundary conditions (24) have been satisfied exactly at the two collocation points \(\zeta_{N,0}=-1\) and \(\zeta_{N,N}=1\).

Let us denote

$$\rho_{n}(t)=\gamma_{n0} c_{1}(t)+ \gamma_{nN} c_{2}(t),\qquad \sigma _{n}(t)= \gamma_{n0} c_{3}(t)+\gamma_{nN} c_{4}(t), $$

then the coupled parabolic PDEs (23) with the set of boundary conditions (24) are transformed into the following ODEs with time delay:

$$ \left . \textstyle\begin{array}{l@{}} \dot{u}_{n}(t)= \eta_{1} \sum_{i=1}^{N-1} \gamma_{ni}u_{i}(t)+u_{n}(t) (1+ \lambda_{1} u_{n}(t-\tau)+\kappa_{1} v_{n}(t-\tau) )+r_{n}(t)+\eta_{1} \rho_{n}(t), \\ \dot{v}_{n}(t)=\eta_{2} \sum_{i=1}^{N-1} \gamma _{ni}v_{i}(t)+v_{n}(t) (1+ \lambda_{2} u_{n}(t-\tau)+\kappa_{2} v_{n}(t-\tau ) )+s_{n}(t)+\eta_{2} \sigma_{n}(t), \end{array}\displaystyle \right \} $$
(32)

and from (25) we get the initial values

$$ \left . \textstyle\begin{array}{l@{}} u_{n}(t)=c_{5}( \zeta_{N,n},t),\\ v_{n}(t)=c_{6}(\zeta_{N,n},t), \quad n=1,\ldots, N-1, t \in[-\tau,0]. \end{array}\displaystyle \right \} $$
(33)

The matrix formulation of the previous systems is

$$\begin{aligned} &\begin{pmatrix} \dot{u}_{1}(t)\\ \dot{u}_{2}(t)\\ \vdots\\ \dot{u}_{N-1}(t)\\ \dot{v}_{1}(t)\\ \dot{v}_{2}(t)\\ \vdots\\ \dot{v}_{N-1}(t) \end{pmatrix} \\ &\quad= \begin{pmatrix} \eta_{1} \sum_{i=0}^{N}\gamma_{1i}u_{i}(t)+u_{1}(t)(1+\lambda_{1} u_{1}(t-\tau)+\kappa_{1} v_{1}(t-\tau))+r_{1}(t)+\eta_{1} \rho_{1}(t)\\ \eta_{1} \sum_{i=0}^{N}\gamma_{2i}u_{i}(t)+u_{2}(t)(1+\lambda_{1} u_{2}(t-\tau)+\kappa_{1} v_{2}(t-\tau))+r_{2}(t)+\eta_{1} \rho_{2}(t)\\ \cdots\\ \cdots\\ \cdots\\ \eta_{1} \sum_{i=0}^{N}\gamma _{N-1i}u_{i}(t)+u_{N-1}(t)(1+\lambda_{1} u_{N-1}(t-\tau)+\kappa_{1} v_{N-1}(t-\tau))+r_{N-1}(t)+\eta_{1} \rho_{N-1}(t)\\ \eta_{2} \sum_{i=0}^{N}\gamma_{1i}v_{i}(t)+v_{1}(t)(1+\lambda_{2} u_{1}(t-\tau)+\kappa_{2} v_{1}(t-\tau))+s_{1}(t)+\eta_{2}\sigma_{1}(t) \\\eta_{2} \sum_{i=0}^{N}\gamma_{2i}v_{i}(t)+v_{2}(t)(1+\lambda_{2} u_{2}(t-\tau)+\kappa_{2} v_{2}(t-\tau))+s_{2}(t)\eta_{2}\sigma_{2}(t)\\ \cdots\\ \cdots\\ \cdots\\ \eta_{2} \sum_{i=0}^{N}\gamma_{N-1i}v_{i}(t)+v_{N-1}(t)(1+\lambda_{2} u_{N-1}(t-\tau)+\kappa_{2} v_{N-1}(t-\tau))+s_{N-1}(t)\eta_{2}\sigma_{N-1}(t) \end{pmatrix}, \end{aligned}$$
(34)

subject to the vector of initials

$$ \begin{pmatrix} u_{1}(t)\\ u_{2}(t)\\ \vdots\\ u_{N-1}(t)\\ v_{1}(t)\\ v_{2}(t)\\ \vdots\\ v_{N-1}(t) \end{pmatrix}= \begin{pmatrix} c_{5}(\zeta_{N,1},t)\\ c_{5}(\zeta_{N,2},t)\\ \vdots\\ c_{5}(\zeta_{N,N-1},t)\\ c_{6}(\zeta_{N,1},t)\\ c_{6}(\zeta_{N,2},t)\\ \vdots\\ c_{6}(\zeta_{N,N-1},t) \end{pmatrix}, \quad t\in[-\tau,0]. $$
(35)

3.3 \((2+1)\) Parabolic PDEs with time delay

In the present subsection, we extend the application of the Chebyshev pseudospectral approximation for a numerical treatment of the \((2+1)\) parabolic PDEs with time delay, namely

$$\begin{aligned} &\frac{\partial v(p,q,t)}{\partial t} =\mu_{1} \frac{\partial^{2} v(p,q,t)}{\partial p^{2}}+\mu_{2} \frac{\partial^{2} v(p,q,t)}{\partial q^{2}}+\mu_{3} v(p,q,t- \tau)+f_{1}(p,q,t), \\ &\quad (p,q, t)\in \Omega_{1}\times\Omega_{2}\times \Omega_{3} , \end{aligned}$$
(36)

where

$$\Omega_{1}\equiv [l_{0},l_{1}],\qquad \Omega_{2}\equiv[l_{2},l_{3}],\quad\mbox{and}\quad \Omega_{3}\equiv (0,T], $$

subject to the initial data

$$ v(p,q,t)=c_{1}(p,q,t),\quad (p,q,t) \in \Omega_{1} \times \Omega_{2}\times[-\tau,0], $$
(37)

and the boundary conditions

$$ \begin{aligned} &v(l_{0},q,t)=c_{2}(q,t), \qquad v(l_{1},q,t)=c_{3}(q,t),\quad(q,t) \in \Omega_{2}\times\Omega_{3}, \\ & v(p,l_{2},t)=c_{4}(p,t),\qquad v(p,l_{3},t)=c_{5}(p,t), \quad (p,t)\in\Omega_{1}\times\Omega_{3}. \end{aligned} $$
(38)

Firstly, we suppose we have the change of spatial variables

$$x=\frac{2}{l_{1}-l_{0}}p+\frac{l_{0}+l_{1}}{l_{0}-l_{1}},\qquad y=\frac {2}{l_{3}-l_{2}}q+ \frac{l_{2}+l_{3}}{l_{2}-l_{3}}, $$

and also \(v(p,q,t)=u(x,y,t)\). For the space variables, the above-mentioned problem transforms into another two-dimensional problem in \([-1,1]\times[-1,1]\),

$$\begin{aligned} \frac{\partial u(x,y,t)}{\partial t} ={}&\mu_{1} \biggl( \frac{2}{l_{1}-l_{0}} \biggr)^{2} \frac{\partial^{2} u(x,y,t)}{\partial x^{2}} \\ &{}+\mu_{2} \biggl(\frac{2}{l_{3}-l_{2}} \biggr)^{2} \frac{\partial^{2} u(x,y,t)}{\partial y^{2}}+\mu_{3} u(x,y,t-\tau)+f_{2}(x,y,t), \\ &{} (x,y, t)\in\mathrm{I}^{2}\times(0,T], \end{aligned}$$
(39)

subject to the initial data

$$ u(x,y,t)=c_{6}(x,y,t),\quad (x,y,t) \in \mathrm{I}^{2}\times[-\tau,0], $$
(40)

and the boundary conditions

$$ \begin{aligned} &u(-1,y,t)=c_{7}(y,t),\qquad u(1,y,t)=c_{8}(y,t),\quad (y,t)\in \mathrm{I}\times\mathrm{T}, \\ &u(x,-1,t)=c_{9}(x,t),\qquad u(x,1,t)=c_{10}(x,t),\quad (x,t)\in \mathrm{I}\times \mathrm{T}. \end{aligned} $$
(41)

Secondly, we construct an algorithm based on the Chebyshev pseudospectral approximation to convert the previous two-dimensional PPDE into a system of ODEs with time delay. We expand the approximate solution in the form

$$ u(x,y,t)=\sum_{i=0}^{N} \sum_{j=0}^{M}a_{i,j}(t) T_{i}(x) T_{j}(y). $$
(42)

Let \(\{\zeta_{N,i}; 0\leq i\leq N\}\) and \(\{\lambda _{M,j}; 0\leq j\leq M\}\) be the zeros of the Chebyshev polynomials \(T_{N}(x)\) and \(T_{M}(y)\), respectively. Making use of the Chebyshev Gauss-Lobatto quadrature and due to the orthogonality relation (4), we have the inversion formulas for \(a_{i,j}(t)\), thus

$$ a_{i,j}(t)=\frac{1}{h_{i} h_{j}} \sum _{l=0}^{N}\sum_{k=0}^{M} \bigl(T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} \bigr) u( \zeta_{N,l},\lambda _{M,k},t). $$
(43)

The approximate solution (42) may be expressed in the form

$$ u(x,y,t)=\sum_{i=0}^{N} \sum_{j=0}^{M}\sum _{l=0}^{N}\sum_{k=0}^{M} \frac{ (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i}h_{j}} T_{i}(x)T_{j}(y) u_{l,k}(t), $$
(44)

where

$$u(\zeta_{N,n},\lambda _{M,m},t)=u_{n,m}(t). $$

In what follows, the first-order partial derivative with respect to x for the solutions (44), at the specific collocation nodes \(\zeta_{N,n}\) and \(\lambda _{M,m}\) can be written as

$$\begin{aligned} &\partial_{x} u_{n,m}(t)=\sum _{l=0}^{N}\sum_{k=0}^{M} \sum_{i=0}^{N}\sum _{j=0}^{M}\frac{ (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} \partial_{x} \bigl(T_{i}(\zeta_{N,n}) \bigr)T_{j}( \lambda _{M,m}) u_{l,k}(t), \\ &\quad n=0,1,\ldots, N,m=0,1,\ldots, M, \end{aligned}$$
(45)

or

$$ \partial_{x} u_{n,m}(t)=\sum _{l=0}^{N}\sum_{k=0}^{M} \gamma^{n,m}_{l,k} u_{l,k}(t), $$
(46)

where \(\gamma^{n,m}_{i,j}\) are the expansion coefficients of the derivative

$$ \gamma^{n,m}_{l,k}=\sum _{i=0}^{N}\sum_{j=0}^{M} \frac{ (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} \partial_{x} \bigl(T_{i}( \zeta_{N,n}) \bigr)T_{j}(\lambda _{M,m}). $$
(47)

Accordingly, the first-order derivative with respect to y for \(u(x,y,t)\), at the nodes \(\zeta_{N,n}\) and \(\lambda _{M,m}\), is

$$\begin{aligned} &\partial_{y} u_{n,m}(t)=\sum _{l=0}^{N}\sum_{k=0}^{M} \sum_{i=0}^{N}\sum _{j=0}^{M}\frac{ (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} T_{i}( \zeta_{N,n})\partial_{y} \bigl(T_{j}( \lambda _{M,m}) \bigr) u_{l,k}(t), \\ &\quad n=0,1,\ldots, N,m=0,1,\ldots, M, \end{aligned}$$
(48)

or

$$ \partial_{y} u_{n,m}(t)=\sum _{l=0}^{N}\sum_{k=0}^{M} \delta^{n,m}_{l,k} u_{l,k}(t), $$
(49)

where

$$ \delta^{n,m}_{l,k}=\sum _{i=0}^{N}\sum_{j=0}^{M} \frac{ (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} T_{i}(\zeta_{N,n})\partial_{y} \bigl(T_{j}(\lambda _{M,m}) \bigr). $$
(50)

Also, the second-order partial derivatives with respect to x and y are, respectively, given by

$$ \partial^{2}_{x} u_{n,m}(t)=\sum_{l=0}^{N}\sum _{k=0}^{M} \xi^{n,m}_{l,k} u_{l,k}(t), $$
(51)

where

$$ \xi^{n,m}_{l,k}=\sum _{i=0}^{N}\sum_{j=0}^{M} \frac { (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} \partial^{2}_{x} \bigl(T_{i}( \zeta_{N,n}) \bigr) T_{j}(\lambda _{M,m}), $$
(52)

and

$$ \partial^{2}_{y} u_{n,m}(t)=\sum_{l=0}^{N}\sum _{k=0}^{M}\eta^{n,m}_{l,k} u_{l,k}(t), $$
(53)

where

$$ \eta^{n,m}_{l,k}=\sum _{i=0}^{N}\sum_{j=0}^{M} \frac { (T_{j}(\lambda _{M,k})\varpi_{M,k} T_{i}(\zeta_{N,l})\varpi_{N,l} )}{h_{i} h_{j}} T_{i}(\zeta_{N,n})\partial^{2}_{y} \bigl(T_{j}(\lambda _{M,m}) \bigr). $$
(54)

Let us denote

$$\begin{aligned}& \epsilon_{1}=\mu_{1} \biggl(\frac{2}{l_{1}-l_{0}} \biggr)^{2},\qquad \epsilon_{2}=\mu _{2} \biggl( \frac{2}{l_{3}-l_{2}} \biggr)^{2}, \\& f_{2}(\zeta_{N,n},\lambda _{M,m},t)=f_{n,m}(t), \qquad c_{6}( \zeta_{N,n},\lambda _{M,m},t)=c_{n,m}(t). \end{aligned}$$

Moreover, the values of \(u_{0,k}(t)\), \(u_{N,k}(t)\), \(u_{l,0}(t)\), and \(u_{l,N}(t)\) can be given by

$$ \begin{aligned} &u_{0,k}(t)=c_{7}( \lambda _{M,k},t),\qquad u_{N,k}(t)=c_{8}( \lambda _{M,k},t),\quad k=0,\ldots,M, \\ & u_{l,0}(t)=c_{9}(\zeta_{N,l},t),\qquad u_{l,N}(t)=c_{10}( \zeta _{N,l},t),\quad l=0, \ldots,N. \end{aligned} $$
(55)

In the Chebyshev pseudospectral approximation for the two-dimensional version of parabolic PDEs with time delay, the residual of (39) is set to zero at \((N-1)\times(M-1)\) of the nodes of Chebyshev Gauss-Lobatto quadrature. Therefore, adopting (44)-(55) enables one to write (39)-(41) as a \((N-1)\times(M-1)\) system of ODEs with time delay,

$$\begin{aligned} &\dot{u}_{n,m}(t)= \epsilon_{1} \sum_{l=0}^{N}\sum _{k=0}^{M} \xi ^{n,m}_{l,k} u_{l,k}(t)+\epsilon_{2} \sum_{l=0}^{N} \sum_{k=0}^{M}\eta^{n,m}_{l,k} u_{l,k}(t)+\mu_{3} u_{n,m}(t-\tau)+f_{n,m}(t), \\ &\quad n=1,\ldots,N-1, m=1,\ldots,M-1, \end{aligned}$$
(56)

subject to the initial data

$$ u_{n,m}(t)=c_{n,m}(t),\quad n=1,\ldots,N-1,m=1, \ldots,M-1. $$
(57)

Finally, the system (56)-(57) may be arranged in the matrix form, which can be solved by the continuous RK scheme,

$$\begin{aligned} & \begin{pmatrix} \dot{u}_{1,1}(t)&\cdots&\dot{u}_{1,M-1}(t)\\ \dot{u}_{2,1}(t)&\vdots&\dot{u}_{2,M-1}(t)\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ \dot{u}_{N-2,1}(t)&\vdots&\dot{u}_{N-2,M-1}(t)\\ \dot{u}_{N-1,1}(t)&\cdots&\dot{u}_{N-1,M-1}(t) \end{pmatrix} \\ &\quad= \begin{pmatrix} {\L}_{1,1}(t, u_{1},\ldots,u_{N-1})&\cdots&{\L}_{1,M-1}(t, u_{1},\ldots ,u_{N-1})\\ {\L}_{2,1}(t, u_{1},\ldots,u_{N-1})&\vdots&{\L}_{2,M-1}(t, u_{1},\ldots ,u_{N-1})\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ {\L}_{N-2,1}(t,u_{1},\ldots,u_{N-1})&\vdots&{\L }_{N-2,M-1}(t,u_{1},\ldots,u_{N-1})\\ {\L}_{N-1,1}(t,u_{1},\ldots,u_{N-1})&\cdots&{\L }_{N-1,M-1}(t,u_{1},\ldots,u_{N-1}) \end{pmatrix}, \end{aligned}$$
(58)
$$\begin{aligned} & \begin{pmatrix} u_{1,1}(t)&\cdots&u_{1,M-1}(t)\\ u_{2,1}(t)&\vdots&u_{2,M-1}(t)\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ u_{N-2,1}(t)&\vdots&u_{N-2,M-1}(t)\\ u_{N-1,1}(t)&\cdots&u_{N-1,M-1}(t) \end{pmatrix}= \begin{pmatrix} c_{1,1}(t)&\cdots&c_{1,M-1}(t)\\ c_{2,1}(t)&\vdots&c_{2,M-1}(t)\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ \cdots& \ddots& \cdots\\ c_{N-2,1}(t)&\vdots&c_{N-2,M-1}(t)\\ c_{N-1,1}(t)&\cdots&c_{N-1,M-1}(t) \end{pmatrix}, \quad t\in[- \tau,0], \end{aligned}$$
(59)

where

$$\begin{aligned} {\L}_{n,m}(t,u_{1}, \ldots,u_{N-1})={}&\epsilon_{1} \sum _{l=0}^{N} \sum_{k=0}^{M} \xi^{n,m}_{l,k} u_{l,k}(t) \\ &{}+\epsilon_{2} \sum_{l=0}^{N} \sum_{k=0}^{M}\eta ^{n,m}_{l,k} u_{l,k}(t)+\mu_{3} u_{n,m}(t-\tau)+f_{n,m}(t), \\ &{}n=1,\ldots,N-1,m=1,\ldots,M-1. \end{aligned}$$
(60)

This system of first-order ODEs can be solved by the continuous RK scheme.

4 Numerical results

To demonstrate the accuracy and applicability of the proposed method in the present paper, four numerical examples are carried out by the algorithms presented in the previous section with a comparison with a linearized compact difference scheme [39] in the one-dimensional case. In each of these examples, we shall highlight the accuracy and efficiency of the method.

Example 4.1

Consider the delay parabolic equation

$$\begin{aligned} u_{t}={}&2 u_{yy}+ \bigl(u(y,t-0.1) \bigr)^{2}-\frac{15}{8}y^{4} \bigl(t^{3}-2t-1 \bigr)- \biggl(\frac {y^{6}}{32}-y \biggr) \bigl(3t^{2}-2 \bigr) \\ &{}- \biggl(\frac{y^{6}}{32}-y \biggr)^{2} \bigl((t-0.1)^{3}-2(t-0.1)-1 \bigr)^{2},\quad (y, t)\in [1,2]\times(0,1] , \end{aligned}$$
(61)

subject to the boundary conditions

$$ u(1,t)=-\frac{31}{32} \bigl(t^{3}-2t-1 \bigr),\qquad u(1,t)=0,\quad t\in(0,1], $$
(62)

and the related initial data

$$ u(y,t)= \biggl(\frac{y^{6}}{32}-y \biggr)^{2} \bigl(t^{3}-2t-1 \bigr)^{2},\quad(y, t) \in [1,2]\times[-0.1,0]. $$
(63)

The exact solution of (61) is

$$ u(y,t)= \biggl(\frac{y^{6}}{32}-y \biggr) \bigl(t^{3}-2t-1 \bigr),\quad(y, t)\in [1,2]\times(0,1]. $$
(64)

For several points, a comparison between the absolute errors of problem (61) which were obtained using the linearized compact difference method (LCDM [39]) and the results obtained by our method is presented in Table 1 at \(N=10\). Table 2 lists the maximum absolute errors (MAEs) for different choices of N, while the best result for the maximum absolute error achieved by using the linearized compact difference method (LCDM [39]) with 80 space steps and \(6{,}400\) time steps was \(1.23\times10^{-8}\).

Table 1 Comparison between the absolute errors of problem ( 61 )
Table 2 MAEs using the Chebyshev collocation method for problem ( 61 )

Figure 1 shows the plot of the numerical solution \(\widetilde{u}(x,t)\), where \(\varepsilon=2^{-4}\) and \(N=16\). Figures 2 and 3 clearly demonstrate the coincidence of the curves of numerical and exact solutions of problem (61), with values of parameters listed in their captions. The absolute error of problem (61) with \(\varepsilon=2^{-4}\) and \(N=16\) is sketched in Figure 4. For the two values of \(t=1\) and \(y=1\), the curves of absolute errors of problem (61) are displayed in Figures 5 and 6, respectively.

Figure 1
figure 1

The numerical solution \(\pmb{\widetilde{u}(y,t)}\) at \(\pmb{N=10}\) .

Figure 2
figure 2

The curves of numerical and exact solutions of problem ( 61 ) at \(\pmb{N=10}\) .

Figure 3
figure 3

The curves of numerical and exact solutions of problem ( 61 ) at \(\pmb{N=10}\) .

Figure 4
figure 4

The absolute error of problem ( 61 ) at \(\pmb{N=8}\) .

Figure 5
figure 5

The absolute error of problem ( 61 ) at \(\pmb{t=1}\) and \(\pmb{N=10}\) .

Figure 6
figure 6

The absolute error of problem ( 61 ) at \(\pmb{y=1.5}\) at \(\pmb{N=10}\) .

Example 4.2

Consider the coupled nonlinear parabolic equation with time delay

$$\begin{aligned} \begin{aligned}u_{t}(x,t)={}& u_{xx}(x,t)+4u(x,t) \bigl(1-u(x,t-\tau) \bigr)- \bigl(v(x,t-\tau) \bigr)^{2} \\ &{}+\frac{e^{-1-\frac{2 t}{\tau}} (e \tau \cos^{2}(\pi x)+\sin(\pi x) (e^{t/\tau} (-1+ (-4+\pi^{2} ) \tau )+4 \tau \sin(\pi x) ) )}{\tau}, \\ v_{t}(x,t)={}& v_{xx}(x,t)+8v(x,t) \bigl(1-v(x,t-\tau) \bigr)- \bigl(u(x,t-\tau) \bigr)^{2} \\ &{}+\frac{e^{-1-\frac{2 t}{\tau}} (e^{t/\tau} (-1+ (-8+\pi^{2} ) \tau ) \cos(\pi x)+8 \tau \cos^{2}(\pi x)+e \tau \sin^{2}(\pi x) )}{\tau}, \\ &{} (x, t)\in[0,\pi]\times(0,1] , \end{aligned} \end{aligned}$$
(65)

subject to the boundary conditions

$$ u(0,t)=u(1,t)=0,\qquad v(0,t)=e^{- (1+\frac{t}{\tau} )}, \qquad v(\pi,t)=-e^{- (1+\frac{t}{\tau} )},\quad t\in(0,5], $$
(66)

and the related initial data

$$ u(x,t)=e^{- (1+\frac{t}{\tau} )} \sin(\pi x),\qquad v(x,t)=e^{- (1+\frac{t}{\tau} )}\cos(\pi x), \quad(x, t)\in [0,\pi]\times(-\tau,1]. $$
(67)

We implemented the algorithm given in Section 3.2 and in Table 3, and we list the MAEs between the exact and approximate solutions for \(u(x,t)\) and \(v(x,t)\) (denoted by \(M^{1}_{E}\) and \(M^{2}_{E}\)) at various choices of N. The errors are calculated through a comparison with the exact solution

$$ u(x,t)=e^{- (1+\frac{t}{\tau} )} \sin(\pi x),\qquad v(x,t)=e^{- (1+\frac{t}{\tau} )} \cos(\pi x),\quad(x, t)\in [0,\pi]\times(0,1]. $$
(68)

Moreover, Table 4 displays the absolute errors for the two solutions (denoted by \(E_{1}(x,t)\) and \(E_{2}(x,t)\)) for different points in the solution interval at \(N=12\). We see in these tables that the results are very accurate for even small choices of the number of nodes, N.

Table 3 MAEs of problem ( 65 )
Table 4 Absolute errors for problem ( 65 ) at \(\pmb{N=12}\)

Example 4.3

Consider the two-dimensional nonlinear parabolic equation with time delay

$$\begin{aligned} u_{t}(x,y,t)={}& u_{xx}(x,y,t)+ u_{yy}(x,y,t)+u(x,y,t-\tau) \\ &{}-e^{-\pi t} \bigl(e^{\pi \tau}+\pi-2 \pi^{2} \bigr) \cos( \pi y) \cos(\pi x), \\ &{} (x,y, t)\in[0,1]\times[0,1]\times(0,1] , \end{aligned}$$
(69)

subject to

$$ \begin{aligned} &u(0,y,t)=e^{-\pi t}\cos(\pi y), \qquad u(1,y,t)=-e^{-\pi t}\cos(\pi y), \\ &u(x,0,t)=e^{-\pi t}\cos(\pi x), \qquad u(x,1,t)=-e^{-\pi t}\cos(\pi x), \end{aligned} $$
(70)

and

$$ u(x,y,t)=e^{-\pi t}\cos(\pi x)\cos( \pi y), \quad(x,y, t)\in [0,1]\times[0,1]\times(-1,1]. $$
(71)

Assume that the absolute error is given by

$$ E_{3}(x,y,t)=\bigl|u(x,y,t)- \widetilde{u}(x,y,t)\bigr|, $$
(72)

where \(u(x,y,t)\) and \(\widetilde{u}(x,y,t)\) are, respectively, the exact and approximate solutions at a specific point \((x,y,t)\) in the solution domain. Also, the MAE is defined by

$$ M^{3}_{E}= \operatorname{Max} \bigl\{ E_{3}(x,y,t):\forall(x,y,t)\in [0,1]\times[0,1] \times(0,1] \bigr\} . $$
(73)

The norm infinity is given by

$$ M^{4}_{E}= \operatorname{Max} \bigl\{ E_{3}(x,y,t_{i}):\forall(x,y)\in [0,1] \times[0,1] \times(0,1] \bigr\} . $$
(74)

Absolute errors of \(u(x,y,t)\) related to (69)-(71) are introduced in Table 5 using the Chebyshev pseudospectral method for different points \((x,y,t)\) at \(N=M=8\), compared with the exact solution

$$ u(x,y,t)=e^{-\pi t}\cos(\pi x)\cos( \pi y). $$
(75)

Table 6 lists the MAEs and the norm infinities for \(u(x,y,t)\) for the two-dimensional problem (69)-(71) at various choices of N and M. Again and from these tables, we clearly observe that the Chebyshev pseudospectral method is accurate for even small values of N.

Table 5 Absolute errors for problem ( 69 )-( 71 ) at \(\pmb{N=M=8}\)
Table 6 The MAEs and the norm infinities of \(\pmb{u(x,y,t)}\) for problem ( 69 )-( 71 )

Example 4.4

Consider the two-dimensional nonlinear parabolic equation with time delay

$$\begin{aligned} u_{t}(x,y,t)={}& u_{xx}(x,y,t)+ u_{yy}(x,y,t)+u(x,y,t-\tau )u_{x}(x,y,t) \\ &{}-e^{\pi -2 \pi t} \cosh(x+y) \bigl(e^{\pi (-1+t)} (2+\pi )+\sinh(x+y) \bigr), \\ &{} (x,y, t)\in[0,1]\times[0,1]\times(0,1] , \end{aligned}$$
(76)

subject to

$$ \begin{aligned} &u(0,y,t)=e^{-\pi t} \cosh(y), \qquad u(1,y,t)=e^{-\pi t}\cosh(1+y), \quad(y,t)\in[0,1]\times(0,1], \\ &v(x,0,t)=e^{-\pi t}\cosh(x), \qquad v(x,1,t)=e^{-\pi t}\cosh(x+1), \quad(x,t) \in[0,1]\times(0,1], \end{aligned} $$
(77)

and

$$ u(x,y,t)=e^{-\pi t}\cosh(x+y),\quad(x,y, t)\in [0,1]\times[0,1]\times(-1,1]. $$
(78)

The exact solution is

$$ u(x,y,t)=e^{-\pi t}\cosh(x+y),\quad(x,y, t)\in [0,1] \times[0,1]\times(0,1]. $$
(79)

Table 7 reports the MAEs and the norm infinities of \(u(x,y,t)\) for problem (76)-(78) using the present method for various choices of N and M. The numerical results presented in this table show that the results are very accurate for small value of N and M.

Table 7 The MAEs and the norm infinities of \(\pmb{u(x,y,t)}\) for problem ( 76 )-( 78 )

5 Concluding remarks and future work

The Chebyshev Gauss-Lobatto pseudospectral method was investigated successfully in spatial discretizations to get accurate approximate solutions of one-dimensional, coupled, and two-dimensional parabolic PDEs with time delay. All of these problems were transformed to systems of ODEs with time delays, greatly simplifying the problems. The continuous RK scheme is then applied to the resulting semidiscrete delay systems. From the numerical experiments, by the obtained results the effectiveness and highly accuracy were demonstrated of the Chebyshev Gauss-Lobatto pseudospectral method for solving the mentioned problems. The present algorithm is also available for approximating the solution for singularly perturbed delay parabolic equations.

We have outlined the implementation of a Chebyshev pseudospectral approximation based on the Chebyshev Gauss-Lobatto quadrature nodes for solving the two-dimensional time delay parabolic PDEs. The technique can be extended to more sophisticated problems with variable and distributed time delays. In principle, this method may be extended to related problems in mathematical physics. It is possible to use other orthogonal polynomials, say Legendre polynomials, or Jacobi polynomials to solve the mentioned problems in this article. Furthermore, the proposed spectral method might be developed by considering the Chebyshev pseudospectral approximation in both temporal and spatial discretizations. We should note that, as a numerical method, we are restricted to solving problems over a finite domain. Also, the pseudospectral approximation might be employed based on generalized Laguerre or modified generalized Laguerre polynomials to solve similar problems in semi-infinite spatial intervals.

References

  1. Canuto, C, Hussaini, MY, Quarteroni, A, Zang, TA: Spectral Methods: Fundamentals in Single Domains. Springer, New York (2006)

    Google Scholar 

  2. Schamel, H, Elsässer, K: The application of the spectral method to nonlinear wave propagation. J. Comput. Phys. 22, 501-516 (1976)

    Article  MATH  Google Scholar 

  3. Bhrawy, AH, Abdelkawy, MA: A fully spectral collocation approximation for multi-dimensional fractional Schrödinger equations. J. Comput. Phys. 294, 462-483 (2015)

    Article  MathSciNet  Google Scholar 

  4. Abd-Elhameed, WM: On solving linear and nonlinear sixth-order two point boundary value problems via an elegant harmonic numbers operational matrix of derivatives. Comput. Model. Eng. Sci. 101(3), 159-185 (2014)

    MathSciNet  Google Scholar 

  5. Abd-Elhameed, WM, Doha, EH, Youssri, YH: Efficient spectral-Petrov-Galerkin methods for third- and fifth-order differential equations using general parameters generalized Jacobi polynomials. Quaest. Math. 36, 15-38 (2013)

    Article  MATH  MathSciNet  Google Scholar 

  6. Kayedi-Bardeh, A, Eslahchi, MR, Dehghan, M: A method for obtaining the operational matrix of the fractional Jacobi functions and applications. J. Vib. Control 20, 736-748 (2014)

    Article  MathSciNet  Google Scholar 

  7. Bhrawy, AH, Doha, EH, Ezz-Eldien, SS, Abdelkawy, MA: A numerical technique based on the shifted Legendre polynomials for solving the time-fractional coupled KdV equation. Calcolo (2015). doi:10.1007/s10092-014-0132-x

    Google Scholar 

  8. Bhrawy, AH, Zaky, MA: A method based on the Jacobi tau approximation for solving multi-term time-space fractional partial differential equations. J. Comput. Phys. 281, 876-895 (2015)

    Article  MathSciNet  Google Scholar 

  9. Eslahchi, MR, Dehghan, M, Parvizi, M: Application of the collocation method for solving nonlinear fractional integro-differential equations. J. Comput. Appl. Math. 257, 105-128 (2014)

    Article  MATH  MathSciNet  Google Scholar 

  10. Golbabai, A, Javidi, M: A numerical solution for non-classical parabolic problem based on Chebyshev spectral collocation method. Appl. Math. Comput. 190, 179-185 (2007)

    Article  MATH  MathSciNet  Google Scholar 

  11. Bhrawy, AH, Zaky, MA: Numerical simulation for two-dimensional variable-order fractional nonlinear cable equation. Nonlinear Dyn. 80(1), 101-116 (2015)

    Article  MathSciNet  Google Scholar 

  12. Bhrawy, AH, Doha, EH, Abdelkawy, MA, Hafez, RM: An efficient collocation algorithm for multidimensional wave type equations with nonlocal conservation conditions. Appl. Math. Model. (2015). doi:10.1016/j.apm.2015.01.029

    Google Scholar 

  13. Doha, EH, Bhrawy, AH, Hafez, RM: On shifted Jacobi spectral method for high-order multi-point boundary value problems. Commun. Nonlinear Sci. Numer. Simul. 17, 3802-3810 (2012)

    Article  MATH  MathSciNet  Google Scholar 

  14. Rezounenko, AV, Wu, JH: A non-local PDE model for population dynamics with state-selective delay: local theory and global attractors. J. Comput. Appl. Math. 190, 99-113 (2006)

    Article  MATH  MathSciNet  Google Scholar 

  15. Cooke, K, van den Driessche, P, Zou, X: Interaction of maturation delay and nonlinear birth in population and epidemic models. J. Math. Biol. 39, 332-352 (1999)

    Article  MATH  MathSciNet  Google Scholar 

  16. Reese, MG: Application of a time-delay neural network to promoter annotation in the Drosophila melanogaster genome. Comput. Chem. 26, 51-56 (2001)

    Article  Google Scholar 

  17. Adimy, M, Crauste, F: Global stability of a partial differential equation with distributed delay due to cellular replication. Nonlinear Anal. 54, 1469-1491 (2003)

    Article  MATH  MathSciNet  Google Scholar 

  18. Dehghan, M, Salehi, R: Solution of a nonlinear time-delay model in biology via semi-analytical approaches. Comput. Phys. Commun. 181, 1255-1265 (2010)

    Article  MATH  MathSciNet  Google Scholar 

  19. Kadalbajoo, MK, Kumar, D: Fitted mesh B-spline collocation method for singularly perturbed differential-difference equations with small delay. Appl. Math. Comput. 204, 90-98 (2008)

    Article  MATH  MathSciNet  Google Scholar 

  20. Wang, WS, Li, SF: Convergence of waveform relaxation methods for neutral delay differential equations. Math. Comput. Model. 48, 1875-1887 (2008)

    Article  MATH  Google Scholar 

  21. Fakhar-Izadi, F, Dehghan, M: The spectral methods for parabolic Volterra integro-differential equations. J. Comput. Appl. Math. 235, 4032-4046 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  22. Kaushik, A, Sharma, KK, Sharma, M: A parameter uniform difference scheme for parabolic partial differential equation with a retarded argument. Appl. Math. Model. 34, 4232-4242 (2010)

    Article  MATH  MathSciNet  Google Scholar 

  23. Bashier, EBM, Patidar, KC: A novel fitted operator finite difference method for a singularly perturbed delay parabolic partial differential equation. Appl. Math. Comput. 217, 4728-4739 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  24. Ghasemi, M, Kajani, MT: Numerical solution of time-varying delay systems by Chebyshev wavelets. Appl. Math. Model. 35, 5235-5244 (2011)

    Article  MATH  MathSciNet  Google Scholar 

  25. Rashidinia, J, Barati, A, Nabati, M: Application of Sinc-Galerkin method to singularly perturbed parabolic convection-diffusion problems. Numer. Algorithms (2014). doi:10.1007/s11075-013-9752-5

    MathSciNet  Google Scholar 

  26. Yuzbasi, S, Sahin, N: Numerical solutions of singularly perturbed one-dimensional parabolic convection.diffusion problems by the Bessel collocation method. Appl. Math. Comput. 220, 305-315 (2013)

    Article  MathSciNet  Google Scholar 

  27. Bhrawy, AH, Assas, LM, Alghamdi, MA: Fast spectral collocation method for solving nonlinear time-delayed Burgers-type equations with positive power terms. Abstr. Appl. Anal. 2013, Article ID 741278 (2013)

    MathSciNet  Google Scholar 

  28. Bhrawy, AH: A Jacobi-Gauss-Lobatto collocation method for solving generalized Fitzhugh-Nagumo equation with time-dependent coefficients. Appl. Math. Comput. 222, 255-264 (2013)

    Article  MathSciNet  Google Scholar 

  29. Bhrawy, AH: An efficient Jacobi pseudospectral approximation for nonlinear complex generalized Zakharov system. Appl. Math. Comput. 247, 30-46 (2014)

    Article  MathSciNet  Google Scholar 

  30. Ashyralyev, A, Agirseven, D: On convergence of difference schemes for delay parabolic equations. Comput. Math. Appl. 66, 1232-1244 (2013)

    Article  MathSciNet  Google Scholar 

  31. Tian, H: Asymptotic stability of numerical methods for linear delay parabolic differential equations. Comput. Math. Appl. 56, 1758-1765 (2008)

    Article  MATH  MathSciNet  Google Scholar 

  32. Zhang, Q, Zhang, C: A new linearized compact multisplitting scheme for the nonlinear convection-reaction-diffusion equations with delay. Commun. Nonlinear Sci. Numer. Simul. 18, 3278-3288 (2013)

    Article  MathSciNet  Google Scholar 

  33. Zhang, Q, Zhang, C: A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations. Appl. Math. Lett. 26, 306-312 (2013)

    Article  MATH  MathSciNet  Google Scholar 

  34. Shen, J, Tang, T, Wang, LL: Spectral Methods: Algorithms, Analysis and Applications. Springer Series in Computational Mathematics. Springer, Berlin (2011)

    Book  Google Scholar 

  35. Enright, WH, Hayashi, H: A delay differential equation solver based on a continuous Runge-Kutta method with defect control. Numer. Algorithms 16, 349-364 (1997)

    Article  MATH  MathSciNet  Google Scholar 

  36. Bellen, A, Vermiglio, R: Some applications of continuous Runge-Kutta methods. Appl. Numer. Math. 22, 63-80 (1996)

    Article  MATH  MathSciNet  Google Scholar 

  37. Odeyemi, T, Mohammadian, A, Seidou, O: Application of the Chebyshev pseudospectral method to van der Waals fluids. Commun. Nonlinear Sci. Numer. Simul. 17, 3499-3507 (2012)

    Article  MATH  MathSciNet  Google Scholar 

  38. Chen, S-S, Li, B-W, Sun, Y-S: Chebyshev collocation spectral method for solving radiative transfer with the modified discrete ordinates formulations. Int. J. Heat Mass Transf. 88, 388-397 (2015)

    Article  Google Scholar 

  39. Sun, Z-Z, Zhang, Z-B: A linearized compact difference scheme for a class of nonlinear delay partial differential equations. Appl. Math. Model. 37, 742-752 (2013)

    Article  MathSciNet  Google Scholar 

Download references

Acknowledgements

This article was funded by the Deanship of Scientific Research DSR, King Abdulaziz University, Jeddah. The authors, therefore, acknowledge with thanks DSR technical and financial support. The authors are very grateful to the referees for carefully reading the paper and for their comments and suggestions, which have improved the paper.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Ali H Bhrawy.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

The authors have made equal contributions to each part of this paper. All the authors read and approved the final manuscript.

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Bhrawy, A.H., Abdelkawy, M.A. & Mallawi, F. An accurate Chebyshev pseudospectral scheme for multi-dimensional parabolic problems with time delays. Bound Value Probl 2015, 103 (2015). https://doi.org/10.1186/s13661-015-0364-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13661-015-0364-y

Keywords