# Linear B-spline finite element method for the generalized diffusion equation with delay

## Abstract

### Objectives

The main aim of this paper is to develop a linear B-spline finite element method for solving generalized diffusion equations with delay. The linear B-spline basis function is used to discretize the space variable. The time discretization process is based on Crank-Nicolson. The benefit of the scheme is that the numerical solution is obtained as a smooth piecewise continuous function which empowers one to find an approximate solution at any desired position in the domain.

### Result

Sufficient and necessary conditions for the numerical method to be asymptotically stable are derived. The convergence of the numerical method is studied. Some numerical experiments are performed to verify the applicability of the numerical method.

## Introduction

In this paper, we consider a class of the generalized delay diffusion equation of the form

\begin{aligned} \left\{ \begin{array}{lll} \frac{\partial u(x,t)}{\partial t} = a_{1}\frac{\partial ^{2} u(x,t)}{\partial x^{2}}+a_{2}\frac{\partial ^{2} u(x,t-\tau )}{\partial x^{2}} ,t>0,0< x <\pi , \\ u(x,t) =\psi \ (x,t), ~~-\tau \le t\le 0,0 \le x \le \pi , \\ u(0,t)=u(\pi ,t)=0, t>0, \end{array} \right. \end{aligned}
(1)

with $$a_{1},a_{2} \in {\mathbb {R}}$$ are real numbers and $$\tau >0$$ is a delay constant. The delay diffusion equation has several applications in science and engineering [1,2,3,4,5]. The generalized delay diffusion equation has intrinsic complex nature because its exact solutions are difficult to obtain. Therefore, one has to mostly rely on numerical treatments. Jackiewicz and Zubik-Kowal  used spectral collocation and waveform relaxation methods to investigate nonlinear partial differential equations with delay. Chen and Wang  used the variational iteration method to study a neutral functional differential equation with delays. The numerical treatments of the generalized delay diffusion equations were studied by many authors(see for instance [8,9,10,11]). Test equation of the type Eq. (1) is also considered in [12, 13]. In these works, the authors applied the separation of the variables to solve analytically.

The finite element method (FEM) is a well-established numerical method for solving partial differential equations (PDEs). The method approximates the exact solution by using piecewise polynomials or B-spline basis functions. B-splines as finite element basis functions provide the required continuity and smoothness. The use of various degrees of B-spline functions to obtain the numerical solutions of some PDEs has been shown to provide easy and simple algorithms. For instance, B-spline finite elements have been widely applied to solve elliptic equations [14, 15], Korteweg-De Vries equation [16,17,18], Burgers’ equation [19,20,21,22], regularized long-wave equation [23, 24], Fokker-Planck equation , advection-diffusion equation , and generalized equal width wave equation , etc., successfully. However, to the best knowledge of the authors, the B-spline FEM method is not considered for finding the approximate solution of the diffusion equation with delay. In this paper, we have applied a linear B-spline FEM to find numerical solutions to the problem under consideration.

### Notations

Let $$H^{r}=H^{r}(\Omega )= W_2^r(\Omega )$$ denotes the Sobolev spaces of order r with respective to norm $$\left\| . \right\| _{r}$$ defined as

\begin{aligned} \quad \left\| \nu \right\| =\quad \left\| \nu \right\| _{L_{2}}:=\Bigg (\int _{\Omega } \nu (x)^{2}dx\Bigg )^{\frac{1}{2}} \end{aligned}

and

\begin{aligned} \quad \left\| \nu \right\| _{r} =\quad \left\| \nu \right\| _{H^{r}}:=\Bigg (\sum _{i \le r}\quad \left\| \frac{\partial ^{i} \nu (x)}{\partial x^{i}}\right\| ^{2}\Bigg )^{\frac{1}{2}}. \end{aligned}

Let $$\nu (x), w(x)(x\in \Omega )$$ be real valued functions.

\begin{aligned} (\nu (x),w(x)):=\int _{\Omega } \nu (x) w(x)dx,~~~~~~ (\nabla \nu (x),\nabla w(x)) := \int _{\Omega }\frac{\partial \nu (x)}{\partial x}\frac{\partial w(x)}{\partial x}dx. \end{aligned}

### Assumption

Assume $$u(t):=u(.,t), u_{t}(t):=u_{t}(.,t),u_{tt}(t):=u_{tt}(.,t), u_{ttt}(t):=u_{ttt}(.,t),\psi (t):=\psi \ (.,t)$$, and $$\psi _{t} (t):=\psi _{t}(.,t)$$.

## Description of the method

Let $$\Delta t = \tau /m$$ be a step size with $$m \ge 1$$, the grid points $$t_{n} = n\Delta t (n=0,1,\dots )$$ and be the approximation in $$S_h$$ of u(t) at $$t=t_{n} = n\Delta t$$. We partition the x -axis into N finite element by choosing a set of equally-spaced knots $$\{x_{k}\}_{k=0}^{N}$$ such at $$0 = x_{0}<x_{1}<\dots x_{N-1}< x_{N}=\pi$$ and $$x_{i+1}-x_{i}=h, i =0,1,2, \dots ,N -1$$.

The linear B-spline basis functions is chosen as follows:

\begin{aligned} Q_{j}(x) = \left\{ \begin{array}{lll} \frac{x-x_{j-1}}{x_{j}-x_{j-1}}, x \in [x_{j-1},x_{j}]\\ \frac{x_{j+1}-x}{x_{j+1}-x_{j}}, x \in [x_{j},x_{j+1}] \\ 0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,x \notin [x_{j-1},x_{j+1}]\\ \end{array}\right. j =1,2,\dots ,N-1. \end{aligned}
(2)

Applying Green’s formula to the second and third terms of equation (1) we obtain

\begin{aligned} (u_{t}(x,t),\nu )+a_{1}(\nabla u(x,t),\nabla \nu )+a_{2}( \nabla u(x,t-\tau ),\nabla \nu )=0,\forall \nu \in H_0^1(\Omega ),t >0. \end{aligned}
(3)

Define the space

\begin{aligned} S_{h} =\{\zeta :\zeta \in C^{2}([0,\pi ]),\zeta |_{[x_{n-1},x_{n}]} \in P^{1} ,1\le n\le N,\zeta (0) = \zeta (\pi ) = 0 \}, \end{aligned}

where $$P^{1}$$ is the space of all polynomials degree less or equal to 1.

We can find the approximate solution $$u_{h}(t): =u_{h}(.,t)$$ belonging to $$S_{h}$$ for each t, so that

\begin{aligned} \left\{ \begin{array}{lll} (u_{h,t}(t,\zeta ))+ a_{1}(\nabla u_{h}(t),\nabla \zeta )+a_{2}( \nabla u_{h}(t-\tau ),\nabla \zeta )=0,\forall \zeta \in S_{h},t>0,\\ u_{h}(x,t)=\psi _{h}(x,t)=0, t>0, \end{array} \right. \end{aligned}
(4)

where $$\psi _{h}(.,t)$$ is an approximation of $$\psi (.,t)$$ in $$S_{h}$$.

Let $$\Delta t =\tau / m$$ be a given step size with $$m \ge 1$$, the grid points $$t_{n} = n\Delta t ( n=0,1,\dots )$$ and $$U^{n}$$ be the approximation in $$S_{h}$$ of u(t) at $$t =t_{n} =n\Delta t$$.

Application of Galerkin Crank-Nicloson method to Eq. (4) gives a numerical scheme of the following type

\begin{aligned} \Bigg (\frac{U^{n}-U^{n-1}}{\Delta t},\zeta \Bigg ) +a_{1} \Bigg (\frac{\nabla U^{n}+\nabla U^{n-1}}{2},\nabla \zeta \Bigg ) +a_{2}\Bigg (\frac{\nabla U^{n-m}+\nabla U^{n-m-1}}{2 },\nabla \zeta \Bigg )=0, \end{aligned}
(5)

where $$U^{n}(.) = \psi (.,t_{n})$$ for $$-m \le n \le 0$$.

Let

\begin{aligned} U^{n}(x) :=\sum _{j=1}^{N-1}Q_{j}(x) \alpha _j^{n}. \end{aligned}
(6)

Substituting Eq. (6) into Eq. (5) and choosing $$\zeta =Q_{i},i=0,\dots ,N-1$$, we get

\begin{aligned} \frac{1}{\Delta t}\sum _{j=1}^{N-1}(\alpha _j^{n}-\alpha _j^{n-1})(Q_{i}(x),Q_{j}(x))=-\frac{a_{1}}{2}\sum _{j=1}^{N-1}(\alpha _j^{n}+\alpha _j^{n-1})(\nabla Q_{i}(x),\nabla Q_{j}(x)) \end{aligned}
\begin{aligned} -\frac{a_{2}}{2}\sum _{j=1}^{N-1}(\alpha _j^{n-m}+\alpha _j^{n-m-1})(\nabla Q_{i}(x),\nabla Q_{j}(x)), \end{aligned}
(7)

which can be rewritten as

\begin{aligned} \frac{1}{\Delta t}\sum _{j=1}^{N-1}(\alpha _j^{n}-\alpha _j^{n-1})\int _{0}^{\pi }Q_{i}(x)Q_{j}(x)dx =-\frac{a_{1}}{2}\sum _{j=1}^{N-1}(\alpha _j^{n}+\alpha _j^{n-1})\int _{0}^{\pi }Q_{i}'(x)Q_{j}'(x)dx \\ \quad -\frac{a_{2}}{2}\sum _{j=1}^{N-1}(\alpha _j^{n-m}+\alpha _j^{n-m-1})\int _{0}^{\pi }Q_{i}'(x)Q_{j}'(x)dx. \end{aligned}
(8)

Define the following matrices:

$$A= (a_{i,j})_{i,j=1}^{N-1}=\int _{0}^{\pi }Q_{i}'(x)Q_{j}'(x)dx,$$
(9)
$$B= (b_{i,j})_{i,j=1}^{N-1}=\int _{0}^{\pi }Q_{i}(x)Q_{j}(x)dx.$$
(10)

The $$(N -1)\times (N-1)$$ matrices A and B are given as follows

\begin{aligned} A= & {} \frac{1}{h}\left( \begin{array}{cccccc} 2 &{} -1 &{} 0 &{}\ldots &{} 0 &{} 0 \\ -1 &{} 2 &{} -1 &{}\ldots &{} 0 &{} 0 \\ 0 &{} -1 &{} 2 &{}\ldots &{} 0 &{} 0 \\ \vdots &{}\vdots &{}\vdots &{}\ddots &{}\vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{}\ldots &{} 2 &{} -1 \\ 0 &{} 0 &{} 0 &{}\ldots &{} -1 &{} 2 \\ \end{array} \right) \nonumber \\ B= & {} \frac{h}{6}\left( \begin{array}{cccccc} 4 &{} 1 &{} 0 &{}\ldots &{} 0 &{} 0 \\ 1 &{} 4 &{} 1 &{}\ldots &{} 0 &{} 0 \\ 0 &{} 1 &{} 4 &{}\ldots &{} 0 &{} 0 \\ \vdots &{}\vdots &{}\vdots &{}\ddots &{}\vdots &{} \vdots \\ 0 &{} 0 &{} 0 &{}\ldots &{} 4 &{} 1 \\ 0 &{} 0 &{} 0 &{}\ldots &{} 1 &{} 4 \\ \end{array} \right) \nonumber \\&\left\{ \begin{array}{lll} (B +\frac{1}{2}a_{1}\Delta t A)\alpha ^{n} =(B -\frac{1}{2}a_{1}\Delta t A )\alpha ^{n-1} -\frac{1}{2}a_{2}\Delta t A(\alpha ^{n-m} + \alpha ^{n-m-1}),\\ \alpha ^{n} = \gamma ^{n} ,for -m \le n \le 0. \end{array}\right. \end{aligned}
(11)

with $$\gamma ^{n}=\psi (t_{n})$$ an initial approximation and $$\alpha ^{n}:=(\alpha _1,\dots ,\alpha _{N} )^{T}$$, and $$B +\frac{1}{2}a_{1}\Delta t A$$ is positive definite and hence, in particular, invertible. Therefore, it has a unique solution.

## Stability analysis

### Definition 1

If the solution $$U^{n}$$ of Eq. (5) corresponding to any sufficiently differentiable function $$\psi _{h}(x,t)$$ with $$\psi _{h}(0,t)$$ =$$\psi _{h}(\pi ,t)$$ satisfies

\begin{aligned} \lim _{n\rightarrow \infty }{U}^{n} = 0,x \in [0,1], \end{aligned}
(12)

then the zero solution of Eq. (5) is called asymptotically stable.

Let $$K :=[x_{i},x_{i+1}]$$ be an element the finite element, and $${\tilde{K}}:=[-1,1]$$ be the reference element in $$\eta$$ -plane. Then

\begin{aligned} \int _{K}{\tilde{Q}}_{i}{\tilde{Q}}_{j}dx=\frac{h}{2}\int _{{\tilde{K}}}\tilde{\tilde{Q_{i}}}\tilde{\tilde{Q_{j}}}d\eta , \int _{K}\nabla \tilde{ Q_{i}}\nabla \tilde{Q_{j}}dx=\frac{2}{h}\int _{{\tilde{K}}}\nabla \tilde{\tilde{Q_{i}}}\nabla \tilde{\tilde{Q_{j}}}d\eta , \end{aligned}

where $${\tilde{B}}=\int _{{\tilde{K}}}\tilde{\tilde{Q_{i}}}\tilde{\tilde{Q_{j}}}d\eta$$ and $${\tilde{A}} = \int _{{\tilde{K}}}\nabla \tilde{\tilde{Q_{i}}}\nabla \tilde{\tilde{Q_{j}}}d\eta$$.

From Eq. (8),

\begin{aligned} \alpha ^{n}= & {} \Bigg (\frac{h}{2}{\tilde{B}} +\frac{a_{1} \Delta t}{h}{\tilde{A}} \Bigg )^{-1}\Bigg (\frac{h}{2}{\tilde{B}} -\frac{a_{1}\Delta t}{h}{\tilde{A}} \Bigg )\alpha ^{n-1} \nonumber \\&-\frac{a_{2}\Delta t}{h}\Bigg (\frac{h}{2}{\tilde{B}} +\frac{a_{1}\Delta t}{h}{\tilde{A}} \Bigg )^{-1}{\tilde{A}}(\alpha ^{n-m}+\alpha ^{n-m-1}) \end{aligned}
(13)
\begin{aligned} \alpha ^{n}= & {} \Bigg (I +\frac{2a_{1}\Delta t}{h^{2}}{\tilde{B}}^{-1}{\tilde{A}}\Bigg )^{-1}\Bigg (I -\frac{2a_{1}\Delta t}{h^{2}}{\tilde{B}}^{-1}{\tilde{A}}\Bigg )\alpha ^{n-1} \nonumber \\&-\frac{2a_{2}\Delta t}{h^{2}}\Bigg (I +\frac{2a_{1}\Delta t}{h^{2}}{\tilde{B}}^{-1}{\tilde{A}} \Bigg )^{-1}{\tilde{B}}^{-1}{\tilde{A}}(\alpha ^{n-m}+\alpha ^{n-m-1}). \end{aligned}
(14)

Let $$\alpha ^{n} =\gamma ^{n}C_{1}$$, where $$C_{1}$$ is a constant vector. The characteristic of Eq. (14) is:

\begin{aligned} \gamma ^{m} -\Bigg (\frac{1 -\frac{2a_{1}\Delta t}{h^{2}}\lambda _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1}\Delta t}{h^{2}}\lambda _{{\tilde{B}}^{-1}{\tilde{A}}}}\Bigg )\gamma ^{m-1}- \Bigg (\frac{\frac{2a_{2}\Delta t}{h^{2}}\lambda _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1}\Delta t}{h^{2}}\lambda _{{\tilde{B}}^{-1}{\tilde{A}}}}\Bigg ) (\gamma +1)=0, \end{aligned}
(15)

where $$\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}$$ denotes the corresponding eigenvalue of $${\tilde{B}}^{-1}{\tilde{A}}$$.

### Lemma 1

 Let $$\kappa _{m}(z) =\alpha (z)z^{m} -\beta (z)$$ be a polynomial, with $$\alpha (z)$$ and $$\beta (z)$$ are polynomials of zero degree. Then $$\kappa _{m}(z)$$ is a Schur polynomial for $$m \ge 1$$ if and only if the following conditions hold

1. (i)

$$\alpha (z) =0 \Rightarrow \left| z\right| < 1,$$

2. (ii)

$$\left| \beta (z)\right| \le \left| \alpha (z)\right| ,\forall z \in {\mathbb {C}} , \left| z\right| = 1,$$ and

3. (iii)

$$\kappa _{m}(z) \ne 0 ,\forall z\in {\mathbb {C}}, \left| z\right| = 1.$$

### Theorem 1

Suppose that $$0 \le a_{2} <a_{1}$$. Then the zero solution of the B-spline finite element method is delay-independently asymptotically stable.

### Proof

Let $$\alpha {(\gamma )} =\gamma -\frac{1 -\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}$$ and $$\beta {(\gamma })=\frac{\frac{2a_{2} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}(\gamma +1)$$.

(i) If $$\alpha {(\gamma )}=0$$, then $$\left| \gamma \right| =\left| \frac{1 -\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}\right| < 1.$$

(ii) For $$\forall \gamma \in {\mathbb {C}}$$, $$\left| \gamma \right| =1$$, represent $$\gamma =\cos \varrho +i\sin \varrho$$, then we get

\begin{aligned} \frac{\gamma -1}{\gamma +1}=\frac{\cos \varrho -1 +i\sin \varrho }{\cos \varrho +1 +i\sin \varrho } =\frac{2i\sin \varrho }{2+2\cos \varrho }. \end{aligned}

We obtain

\begin{aligned}&\left| \frac{\alpha (\gamma ) }{\gamma +1}\right| =\left| \frac{\gamma -\frac{1 -\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1 +\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}}{\gamma +1} \right| =\left| \frac{(\gamma -1)}{(\gamma +1)(1+\frac{2a_{1}\Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}})}+\frac{\frac{2a_{1}\Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1+\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}\right| \\&\quad \ge \frac{\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1+\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}>\frac{\frac{2a_{2} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}{1+\frac{2a_{1} \Delta t}{h^{2}}\gamma _{{\tilde{B}}^{-1}{\tilde{A}}}}=\left| \frac{\beta (\gamma ) }{\gamma +1}\right| . \end{aligned}

(iii) By (ii), it is straightforward. $$\square$$

## Convergence Analysis

In this section, we present the convergence analysis for the proposed method.

The Ritz projection $$R_{h}:H_0^{1}(\Omega ) \rightarrow S_{h}$$ is a mapping for any $$\nu \in H_0^{1}(\Omega )$$ such that

\begin{aligned} (\nabla R_{h}\nu -\nu ,\nabla w)=0, \forall w \in S_{h}. \end{aligned}
(16)

### Lemma 2

Assume that for any $$v \in H^{s}(\Omega )\cap H_0^{1}(\Omega )$$,

\begin{aligned} \inf _{\zeta \in {S}_{h}}\{\left\| \nu -\zeta \right\| + h\left\| \nabla (\nu -\zeta )\right\| \} \le Ch^{s}\left\| \nu \right\| _{s},for ~~ 1\le s \le r . \end{aligned}

holds. Then, with $$R_{h}$$ defined by Eq. (16), we have

\begin{aligned} \left\| R_{h}\nu -\nu \right\| +h\left\| \nabla (R_{h}\nu -\nu )\right\| \le Ch^{s}\left\| \nu \right\| _{s}, for~~ any ~ \nu ~ \in H^{s}(\Omega )\cap H_0^{1}(\Omega ), 1\le s \le r. \end{aligned}

The number r is referred to as the order of accuracy of the family $$\{S_{h}\}$$. For the case of piecewise linear B-spline basis function, $$r =2$$.

Define $$u(t):= u(.,t)$$ and $$u:[0,+\infty ) \rightarrow H_0^{1}(\Omega )$$. Let $$D_{h}: H_0^{1}(\Omega )\rightarrow S_h$$ by

\begin{aligned} a_{1}(\nabla D_{h}u(t) -\nabla u(t), \nabla \zeta ) +a_{2} ( \nabla D_{h} u(t-\tau )-\nabla u(t-\tau ),\nabla \zeta ) =0, \forall \zeta \in S_h \end{aligned}
(17)

and

\begin{aligned} D_{h}u(t)=R_{h}u(t)=R_{h}\psi (t),for ~ -\tau \le t \le 0. \end{aligned}
(18)

### Theorem 2

Let u and $$U^{n}$$ be the solution of (3) and (5), respectively. Assume that $$\left\| u(t) -R_{h} u (t)\right\| \le Ch^{2} \left\| u (t)\right\| _{2}$$, $$\left\| u_{t}(t) -R_{h} u_{t} (t)\right\| \le Ch^{2} \left\| u _{t}(t)\right\| _{2}$$, $$-\tau \le t\le 0$$ and $$\left\| \psi _{h}(t)-\psi (t)\right\| \le Ch^{2}$$, then

\begin{aligned} \left\| U^{n} -u (t_{n})\right\| \le C(h^{2}+(\Delta t)^{2}) ,for~~ n=1,2,... \end{aligned}

where C is a positive constant independent of h and $$\Delta t$$.

### Proof

Define

\begin{aligned} e^{n}=U^{n} -u (t_{n})=(U^{n} -D_{h}u (t_{n}))+(D_{h}u (t_{n})- u (t_{n} )) =\mu ^{n}+\sigma ^{n}, \end{aligned}

where

$$\mu ^{n}=U^{n} -D_{h}u (t_{n})$$, $$\sigma ^{n} = D_{h}u (t_{n})- u (t_{n} )$$, so that

\begin{aligned} \left\| U^{n} -u (t_{n})\right\| \le \left\| \mu ^{n}\right\| +\left\| \sigma ^{n}\right\| . \end{aligned}

The term $$\sigma ^{n}(t) =\sigma (t_{n})$$ is easily bounded by lemma 2.

\begin{aligned}&\Bigg (\frac{\mu ^{n}-\mu ^{n-1}}{\Delta t},\zeta \Bigg )+a_{1}\Bigg (\frac{\nabla \mu ^{n}+\nabla \mu ^{n-1}}{2},\nabla \zeta \Bigg )+a_{2}\Bigg (\frac{\nabla \mu ^{n-m}+\nabla \mu ^{n-m-1}}{2},\nabla \zeta \Bigg ) \nonumber \\&\quad =-(W^{n},\zeta ), \forall \zeta \in S_h, \end{aligned}
(19)

where

\begin{aligned} W^{n}= & {} \frac{D_{h} u(t_{n})-D_{h} u(t_{n-1})}{\Delta t} -\frac{ u_{t}(t_{n})+ u_{t}(t_{n-1})}{2} \\= & {} (D_{h}-I){\bar{\partial }} u(t_{n}) +\Bigg ({\bar{\partial }} u(t_{n})-\frac{u_{t}(t_{n})+u_{t}(t_{n-1})}{2}\Bigg )=:W_1^{n}+W_2^{n}. \end{aligned}

Setting $$\zeta =\frac{\mu ^{n}+\mu ^{n-1}}{2}$$, gives

\begin{aligned}&\Bigg ( \frac{\mu ^{n}-\mu ^{n-1}}{\Delta t},\frac{\mu ^{n}+\mu ^{n-1}}{2}\Bigg )+a_{1} \left\| \frac{\mu ^{n}+\mu ^{n-1}}{2} \right\| _1^{2} +a_{2}\Bigg (\frac{\nabla \mu ^{n-m}+\nabla \mu ^{n-m-1}}{2},\frac{\nabla \mu ^{n}+\nabla \mu ^{n-1}}{2}\Bigg ) \\&\quad =-\Bigg (W^{n},\frac{\mu ^{n}+\mu ^{n-1}}{2}\Bigg ). \end{aligned}

By applying Schwartz inequality,

\begin{aligned}&\Bigg ( \frac{\mu ^{n}-\mu ^{n-1}}{\Delta t},\frac{\mu ^{n}+\mu ^{n-1}}{2}\Bigg )+ \left\| \frac{\mu ^{n}+\mu ^{n-1}}{2} \right\| _1^{2} \le C\Bigg ( \left\| \frac{\mu ^{n-m}+\mu ^{n-m-1}}{2} \right\| _1^{2} \nonumber \\&\qquad +\left\| W^{n}\right\| \left\| \frac{\mu ^{n}+\mu ^{n-1}}{2} \right\| \Bigg ). \end{aligned}

So

\begin{aligned} \left\| \mu ^{n} \right\| ^{2} +\Delta t \left\| \frac{\mu ^{n}+\mu ^{n-1}}{2} \right\| _1^{2}\le C\Bigg (\left\| \mu ^{n-1} \right\| ^{2} +\Delta t \left\| \frac{\mu ^{n-m}+\mu ^{n-m-1}}{2} \right\| _1^{2} +(\Delta t )^{2}\left\| W^{n}\right\| ^{2}\Bigg ). \end{aligned}

We can assume that $$n \in ((k-1)m,km],k \in N$$. Then

\begin{aligned}&\Delta t \left\| \frac{\mu ^{n}+\mu ^{n-1}}{2} \right\| _1^{2}\le C\Bigg (\left\| \mu ^{n-1} \right\| ^{2} +\Delta t \left\| \frac{\mu ^{n-m}+\mu ^{n-m-1}}{2} \right\| _1^{2} +(\Delta t )^{2}\left\| W^{n}\right\| ^{2}\Bigg ) \\&\quad \le C\Bigg (\left\| \mu ^{n-1} \right\| ^{2} + \left\| \mu ^{n-m-1} \right\| ^{2}+\Delta t \left\| \frac{\mu ^{n-2m}+\mu ^{n-2m-1}}{2} \right\| _1^{2}+(\Delta t )^{2}( \left\| W^{n}\right\| ^{2}+\left\| W^{n-m}\right\| ^{2})\Bigg ) \\&\quad \le ...\le C\Bigg (\sum _{i=0}^{k-1}\left\| \mu ^{n-im-1} \right\| ^{2}+\Delta t \left\| \frac{\mu ^{n-km}+\mu ^{n-km-1}}{2} \right\| _1^{2}+(\Delta t )^{2} \sum _{i=0}^{k-1}\left\| W^{n-im}\right\| ^{2}\Bigg ). \end{aligned}

Therefore

\begin{aligned} \left\| \mu ^{n} \right\| ^{2}\le C\Bigg ( \sum _{i=0}^{k-1}\left\| \mu ^{n-im-1} \right\| ^{2}+\Delta t \left\| \frac{\mu ^{n-km}+\mu ^{n-km-1}}{2}\right\| _1^{2}+(\Delta t )^{2} \sum _{i=0}^{k-1}\left\| W^{n-im}\right\| ^{2}\Bigg ). \end{aligned}

By applying Gronwall inequality,

\begin{aligned} \left\| \mu ^{n} \right\| ^{2}\le C\Bigg ( \left\| \mu ^{0} \right\| ^{2}+\Delta t \left\| \frac{\mu ^{n-km}+\mu ^{n-km-1}}{2} \right\| _1^{2}+(\Delta t )^{2} \sum _{i=0}^{k-1}\left\| W ^{n-im}\right\| ^{2}\Bigg ). \end{aligned}
(20)

Write

\begin{aligned} W_1^{n} =(D_{h}-I){\tilde{\partial }}u(t_{n}) = \Delta t ^{-1}\int _{t_{n-1}}^{t_{n}}(D_{h}-I)u_{t}(t)dt, \end{aligned}

so

\begin{aligned} (\Delta t )^{2} \sum _{i=1}^{k-1}\left\| W_1^{n-im}\right\| ^{2}\le \sum _{i=1}^{k-1}\Bigg (\int _{t_{n -im-1}}^{t_{n-im}} Ch^{2}\left\| u_{t}(t) \right\| _{2} dt\Bigg )^{2} \le Ch^{2(2)}. \end{aligned}
(21)

Further

\begin{aligned} \left\| \Delta t W_2^{i} \right\| =\left\| u(t_{i}) -u(t_{i-1})-\Delta t \frac{ u_{t}(t_{i})+ u_{t}(t_{i-1})}{2} \right\| \le C(\Delta t )^{2}\int _{t_{i-1}}^{t_{i}}\left\| u_{ttt}(t) \right\| dt , \end{aligned}

so that

\begin{aligned} (\Delta t )^{2} \sum _{i=1}^{k-1}\left\| W_2^{n-im}\right\| ^{2}\le C(\Delta t )^{4} \sum _{i=1}^{k-1}\Bigg (\int _{t_{n-im-1}}^{t_{n-im}}\left\| u_{ttt}(S) \right\| dt\Bigg )^{2} \le C(\Delta t )^{4}. \end{aligned}
(22)

From Eq. (21) and Eq. (22), we have

\begin{aligned} \left\| U^{n} -u (t_{n})\right\| \le C(h^{2}+(\Delta t)^{2}) ,for~~ n=1,2,\dots \end{aligned}

$$\square$$

## Numerical experiments

The performance of the proposed methods is tested by using numerical experiments. To evaluate errors, $$L_{\infty }$$ and $$L_{2}$$ error norms are applied as follows:

\begin{aligned} L_{\infty }= \mathop {\max }\limits _{1 \le n \le N } \left| {u(t_{n})-(U^{n})} \right| ,L_{2}=\sqrt{ h\sum _{i=1}^{N}\left| {u(t_{n})-(U^{n})} \right| ^{2}} \end{aligned}

Order of convergence is obtained by

\begin{aligned} Order =\frac{\log (E^{h_{1}}/E^{h_{2}})}{\log (h_{1}/h_{2})} \end{aligned}

where $$E^{h_{1}}$$ and $$E^{h_{2}}$$ represent the errors at step sizes $$h_{1}$$ and $$h_{2}$$, respectively.

### Example 1

 Consider

\begin{aligned} \left\{ \begin{array}{lll} \frac{\partial u(x,t)}{\partial t} = a_{1}\frac{\partial ^{2} u(x,t)}{\partial x^{2}}+a_{2}\frac{\partial ^{2} u(x,t-\tau )}{\partial x^{2}} ,t>0,0< x <\pi , \\ u(x,t) =\psi \ (x,t), ~~-\tau \le t\le 0,0 \le x \le \pi ,\\ u(0,t)=u(\pi ,t)=0, t>0. \end{array} \right. \end{aligned}
(23)

First, we take the initial function as $$\psi (x,t) =sin(x),\tau =1,a_{1}=1.5,a_{2}=1$$ such that the trivial solution of Eq.(1) is asymptotically stable. Numerical results are obtained and plotted at time $$T=5$$ using different $$(\Delta t =\tau /m$$,$$h=\pi /N$$).

We apply the proposed method with different step sizes to solve the problem. The graph of numerical results is shown in Fig. 1. This graph shows that the numerical solution is asymptotically stable. And these confirm the theoretical results in Theorem 1.

### Example 2

 Consider

\begin{aligned} \left\{ \begin{array}{lll} \frac{\partial u(x,t)}{\partial t} = a_{1}\frac{\partial ^{2} u(x,t)}{\partial x^{2}}+a_{2}\frac{\partial ^{2} u(x,t-\tau )}{\partial x^{2}}+h(x,t),t>0,0< x <\pi , \\ u(x,t) =\psi \ (x,t), ~~-\tau \le t\le 0,0 \le x \le \pi , \\ u(0,t)=u(\pi ,t)=0, t>0, \end{array}\right. \end{aligned}
(24)

with the initial condition we take the initial function as $$\psi (x,t) =\sin (x)$$, and the added term h(xt) where that is the exact solution is $$u(x,t)=\exp ({-t})sin(x)$$. Here, we take the parameters $$a_{1}=1,a_{2}=0.5,\tau =0.5$$ and compute the problem on $$[0,\pi ]\times [0,2]$$ for different space and temporal step sizes $$(\Delta x=\pi /N,\Delta t=\tau /m)$$.

Table  1 shows the numerical errors and the corresponding orders. When the grid size is reduced, both error norms are significantly reduced. These results show the convergence of the linear B-spline finite element method. The given results suggest that the proposed method has order 2 of accuracy. The calculated error norms are also compared with the result obtained using the central difference method . In Table  2, the comparison between the exact and approximation solution are given.

## Conclusion

In this paper, a finite element method is constructed based on linear B-spline basis functions for solving the generalized diffusion equations with delay. The detailed description of results through tables and graphs proves that the proposed numerical method is working efficiently. For all the test cases, simulations at a different set of data points are carried out to check the applicability of the numerical scheme. Based on these observations, our expectation that the given method is well suited to the generalized diffusion with the delay is confirmed.

## Limitations

The linear B-spline basis functions yields an order 2 of accuracy. One can use higher polynomial basis functions in order to increase the order of accuracy in space.

## Data availibility

No additional data is used for this research work.

## Abbreviations

FEM:

Finite element method

PDEs:

Partial differential equations

## References

1. Adomian G, Rach R. Nonlinear stochastic differential delay equations. J Math Anal Appl. 1983;91(1):94–101.

2. Baranowski J, Legendre polynomial approximations of time delay systems. In: Proceeding of 12th International PhD Workshop. 2010; 1: 15–20.

3. Asl F.M, Ulsoy A.G. Analysis of a system of linear delay differential equations. J Dyn Sys Meas Control. 2003;125(2):215–23.

4. Wang XT. Numerical solution of delay systems containing inverse time by hybrid functions. Appl Math Comput. 2006;173(1):535–46.

5. Hale JK, Retarded functional differential equations: basic theory. In: Theory of functional differential equations. 1977; 36–56.

6. Jackiewicz Z, Zubik-Kowal B. Spectral collocation and waveform relaxation methods for nonlinear delay partial differential equations. Appl Numer Math. 2006;56(3–4):433–43.

7. Chen X, Wang L. The variational iteration method for solving a neutral functional-differential equation with proportional delays. Comput Math Appl. 2010;59(8):2696–702.

8. Garcia P, Castro M, Martín JA, Sirvent A. Numerical solutions of diffusion mathematical models with delay. Math Comput Model. 2009;50(5–6):860–8.

9. Garcia P, Castro M, Martín JA, Sirvent A. Convergence of two implicit numerical schemes for diffusion mathematical models with delay. Math Comput Model. 2010;52(7–8):1279–87.

10. Adam A, Bashier E, Hashim M, Patidar K. Fitted galerkin spectral method to solve delay partial differential equations. Math Methods Appl Sci. 2016;39(11):3102–15.

11. Zhang Q, Chen M, Xu Y, Xu D. Compact $$\theta$$-method for the generalized delay diffusion equation. Appl Math Comput. 2018;316:357–69.

12. Martín JA, Rodríguez F, Company R. Analytic solution of mixed problems for the generalized diffusion equation with delay. Math Comput Model. 2004;40(3–4):361–9.

13. Reyes E, Rodríguez F, Martín JA. Analytic-numerical solutions of diffusion mathematical models with delays. Comput Math Appl. 2008;56(3):743–53.

14. Sana M, Mustahsan M. Finite element approximation of optimal control problem with weighted extended b-splines. Mathematics. 2019;7(5):452.

15. Wu J, Zhang X. Finite element method by using quartic b-splines. Numer Methods Partial Differ Equ. 2011;27(4):818–28.

16. Zaki S. A quintic b-spline finite elements scheme for the kdvb equation. Comput Methods Appl Mech Eng. 2000;188(1–3):121–34.

17. Aksan E, Özdeş A. Numerical solution of korteweg-de vries equation by galerkin b-spline finite element method. Appl Math Comput. 2006;175(2):1256–65.

18. Ak T, Dhawan S, Karakoc SBG, Bhowmik SK, Raslan KR. Numerical study of rosenau-kdv equation using finite element method based on collocation approach. Math Model Anal. 2017;22(3):373–88.

19. Özis T, Esen A, Kutluay S. Numerical solution of burgers’ equation by quadratic b-spline finite elements. Appl Math Comput. 2005;165(1):237–49.

20. Aksan EN. An application of cubic b-spline finite element method for the burgersequation. Therm Sci. 2018;22(Suppl. 1):195–202.

21. Kutluay S, Ucar Y. Numerical solutions of the coupled burgers’ equation by the galerkin quadratic b-spline finite element method. Math Methods Appl Sci. 2013;36(17):2403–15.

22. Soliman A, A galerkin solution for burgers’ equation using cubic b-spline finite elements. In: Abstract and applied analysis. London: Hindawi; 2012. 2012.

23. Dağ İ, Saka B, Irk D. Galerkin method for the numerical solution of the rlw equation using quintic b-splines. J Comput Appl Math. 2006;190(1–2):532–47.

24. Saka B, Dağ İ. Quartic b-spline collocation algorithms for numerical solution of the rlw equation. Numer Methods Partial Differ Equ. 2007;23(3):731–51.

25. Lakestani M, Dehghan M. Numerical solution of fokker-planck equation using the cubic b-spline scaling functions. Numer Methods Partial Differ Equ. 2009;25(2):418–29.

26. Dhawan S, Bhowmik SK, Kumar S. Galerkin-least square b-spline approach toward advection-diffusion equation. Appl Math Comput. 2015;261:128–40.

27. Bhowmik SK, Karakoc SBG. Numerical solutions of the generalized equal width wave equation using the petrov-galerkin method. Appl Anal. 2021;100(4):714–34.

28. Tian H. Asymptotic stability analysis of the linear $$\theta$$-method for linear parabolic differential equations with delay. J Differ Equ Appl. 2009;15(5):473–87.

29. Liang H. Convergence and asymptotic stability of galerkin methods for linear parabolic equations with delays. Appl Math Comput. 2015;264:160–78.

30. Wu F, Li D, Wen J, Duan J. Stability and convergence of compact finite difference method for parabolic problems with delay. Appl Math Comput. 2018;322:129–39.

## Acknowledgements

The authors would like to appreciate the anonymous referees for their constructive suggestions.

## Funding

No funding organization for this research work.

## Author information

Authors

### Contributions

GTL carried out scheme development, MATLAB coding, and numerical experimentation. GFD formulated the problem, designed, and drafted the manuscript. Both authors read and approved the final manuscript.

### Corresponding author

Correspondence to Gemeda Tolessa Lubo.

## Ethics declarations

### Competing interests

The authors declare that they have no competing interests. 