Skip to main content

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



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.


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.


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}$$

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 [6] used spectral collocation and waveform relaxation methods to investigate nonlinear partial differential equations with delay. Chen and Wang [7] 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 [25], advection-diffusion equation [26], and generalized equal width wave equation [27], 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.


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}$$


$$\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}$$


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)\).

Main text

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}$$

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}$$

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}$$

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}$$

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


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

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}$$

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}$$

Define the following matrices:

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

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}$$

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}$$

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}$$
$$\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}$$

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}$$

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

Lemma 1

[28] 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.


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}$$

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}$$


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

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\).



$$\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}$$


\(\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}$$


$$\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}$$


$$\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}$$


$$\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}$$


$$\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}$$


$$\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}$$


$$\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}$$

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}$$


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

[29] 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}$$

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\)).

Fig. 1
figure 1

Solution of (23) with parameter values a) \(N=10\) and \(m=40\). b) \(N=10\) and \(m=50\). c) \(N=10\) and \(m=200\). d) \(N=10\) and \(m=500\)

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

[30] 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}$$

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 Errors norms and the corresponding convergence orders (\(\Delta t\approx \Delta x^{2}\)) for example 2
Table 2 Comparison of the numerical solutions obtained with various values of m for \(N = 10 , T = 1\), and \(\tau =0.5\) with the exact solution for example 2

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 [30]. In Table  2, the comparison between the exact and approximation solution are given.


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.


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.



Finite element method


Partial differential equations


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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Article  Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Google Scholar 

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

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Google Scholar 

  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.

    Article  Google Scholar 

  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.

    Article  Google Scholar 

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

    Google Scholar 

  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.

    Google Scholar 

Download references


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


No funding organization for this research work.

Author information

Authors and Affiliations



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.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Lubo, G.T., Duressa, G.F. Linear B-spline finite element method for the generalized diffusion equation with delay. BMC Res Notes 15, 195 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


Mathematics Subject Classification