# On mixed steps-collocation schemes for nonlinear fractional delay
differential equations

M. Mousa-Abadian<EMAIL_ADDRESS><EMAIL_ADDRESS>and  S. H.
Momeni-Masuleh∗ Department of Mathematics, Shahed University, P.O. Box
18151-159, Tehran, Iran

###### Abstract.

This research deals with the numerical solution of non-linear fractional
differential equations with delay using the method of steps and shifted
Legendre (Chebyshev) collocation method. This article aims to present a new
formula for the fractional derivatives (in the Caputo sense) of shifted
Legendre polynomials. With the help of this tool and previous work of the
authors, efficient numerical schemes for solving non-linear continuous
fractional delay differential equations are proposed. The proposed schemes
transform the nonlinear fractional delay differential equations to a non-delay
one by employing the method of steps. Then, the approximate solution is
expanded in terms of Legendre (Chebyshev) basis. Furthermore, the convergence
analysis of the proposed schemes is provided. Several practical model examples
are considered to illustrate the efficiency and accuracy of the proposed
schemes.

###### Key words and phrases:

Nonlinear fractional differential equations, delay differential equations,
method of steps, shifted Legendre (Chebyshev) basis

###### 2010 Mathematics Subject Classification:

34K37, 34A34, 26A33, 65M70, 33C45

∗Corresponding author

##  1\. Introduction

Delay differential equations (DDEs) belong to a broader class of functional
differential equations in which the rate of change of the unknown function at
a specific time is represented due to the values of the function at previous
times. DDEs are also known as differential-difference equations. Laplace and
Condorcet [4] first studied these equations, and naturally, appear in various
fields of science and engineering [8]. Fractional delay differential equations
(FDDEs) are considered as a generalization of DDEs which contain derivatives
of arbitrary fractional order. The integer order differential operator is a
local operator, while the fractional order differential operator is a non-
local operator. More precisely, the next state of a system which is modeled by
fractional differential equations (FDEs), depends not only upon its present
situation but also upon all of its past situations. The fractional order
differential operator enables us to describe a real event more accurately than
the classical integer order differential operator. In recent years, FDEs and
FDDEs are frequently used to model many real phenomena in the fields of
control theory [28], biology [14, 7], economy [2], and so on.

Because of the computational complexities of fractional delay derivatives, for
most of the FDEs and FDDEs, the exact solution is available and therefore, it
is necessary to employ numerical methods for solving such equations. Shakeri
and Dehghan [27] employ the homotopy perturbation method to solve delay
differential equations with integer order derivatives. Variational iteration
method is considered by Saadatmandi and Dehghan [24] to obtain the numerical
solution of the generalized pantograph equation. Moghaddam and Mostaghim [15,
20] introduced numerical methods in the frame work of the finite difference
method for solving FDDEs. They also presented a matrix approach using
fractional finite difference method to solve nonlinear FDDEs [16]. Prakash et
al. [23] proposed a numerical algorithm based on a modified He-Laplace method
to solve nonlinear FDDEs. Wang [29] combined the Adams-Bashforth-Moulton
method with the linear interpolation method to find an approximate solution
for FDDEs. Mohammed and Khlaif [17] applied the Adomian decomposition method
to get the numerical solution of FDDEs. Mousa-Abadian and Momeni-Masuleh
proposed a numerical scheme to solve linear fractional delay differential
systems [19]. Their scheme employs the method of steps to handle the delay
term and Chebyshev-Tau method to construct the approximate solution. Sedaghat
et al. [26] introduced a numerical scheme using Chebyshev polynomials for
solving FDDEs of pantograph type. Saeed et al. [25] developed Chebyshev
wavelet methods for solving FDDEs. Khader [12] derived an approximate formula
of the Laguerre polynomials for the numerical treatment of FDDEs. Daftardar-
Gejji et al. [6] extended a new predictor-corrector method to solve FDDEs.
Moghaddam et al. [21] developed a numerical method based on the Adams-
Bashforth-Moulton method for solving variable-order FDDEs. Yaghoobi et al.
[30] devised a numerical scheme based on a cubic spline interpolation for
solving variable-order FDDEs. Khosravian-Arab et al. [13] devised new Lagrange
basis functions to approximate fractional derivatives in unbounded domains.
Their approach is based on the pseudo-spectral, Galerkin, and Petrov-Galerkin
methods.

A numerical approach to solve nonlinear FDDEs can be a generalization of the
method which introduced in Ref. [19]. The Chebyshev collocation method can be
considered solving nonlinear FDDEs. Of course, employing collocation methods
are not restricted to use Chebyshev basis, but also Legendre basis can be
applied. Therefore, a significant part of this article deals with solving
nonlinear FDDEs by using Legendre basis. In this paper, we derive a new
formula for the fractional derivatives of shifted Legendre polynomials and
then present the efficient numerical schemes for solving nonlinear FDDEs.

The remainder of this article proceeds as follows. In the next section, the
properties of shifted collocation-type bases are discussed. In Section 3, a
formula for the fractional derivatives of shifted Legendre basis is derived.
Section 4 describes mixed steps-collocation schemes for solving nonlinear
FDDEs. The convergence analysis of the proposed schemes has been done in
Section 5. Section 6 concerns with applying the proposed schemes to several
nonlinear FDDEs. The conclusion is given in Section 7.

##  2\. Shifted collocation-type bases

The most common collocation methods are those based on Chebyshev and Legendre
basis. Properties of shifted Chebyshev basis have been investigated in Ref.
[19]. Here, we deduce the properties of shifted Legendre basis. The Legendre
basis $L_{\scriptscriptstyle\scriptscriptstyle{k}}(t)$ for $k=0,1,\dots,$ and
$-1\leq t\leq 1$, can be defined as the solution of following ordinary
differential equation [3]

$\frac{d}{dt}\Big{(}(1-t^{2})\frac{dL_{\scriptscriptstyle\scriptscriptstyle{k}}}{dt}(t)\Big{)}+k(k+1)L_{\scriptscriptstyle\scriptscriptstyle{k}}(t)=0,$

which satisfy $L_{\scriptscriptstyle\scriptscriptstyle{k}}(\pm 1)=(\pm
1)^{k}$. For $k\geq 1$, we have the following recurrence formula

(2.1)
$L_{\scriptscriptstyle\scriptscriptstyle{k+1}}(t)=\frac{2k+1}{k+1}tL_{\scriptscriptstyle\scriptscriptstyle{k}}(t)-\frac{k}{k+1}L_{\scriptscriptstyle\scriptscriptstyle{k-1}}(t),$

where $L_{\scriptscriptstyle\scriptscriptstyle{0}}(t)=1$ and
$L_{\scriptscriptstyle\scriptscriptstyle{1}}(t)=t$. The shifted Legendre basis
are defined on the interval $[\alpha,\beta]$ using the change of variable
$t=\frac{2}{\beta-\alpha}(x-\beta)+1$. For the simplicity, let us to denote
$L_{\scriptscriptstyle\scriptscriptstyle{k}}(\frac{2}{\beta-\alpha}(x-\beta)+1)$
by $L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)$. Therefore,
similar to (2.1), the following recurrence relation can be obtained

$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k+1}}(x)=\frac{2k+1}{k+1}\big{(}\frac{2}{\beta-\alpha}(x-\beta)+1\big{)}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)-\frac{k}{k+1}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k-1}}(x).$

The shifted Legendre basis
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)$ can be presented
in the following form

(2.2)
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=\frac{1}{2^{n}}\sum_{k=0}^{[\frac{n}{2}]}(-1)^{k}\binom{n}{k}\binom{2n-2k}{n}\big{(}\frac{2}{\beta-\alpha}(x-\beta)+1\big{)}^{n-2k}.$

By using the identity

(2.3)
$\big{(}\frac{2}{\beta-\alpha}(x-\beta)+1\big{)}^{n-2k}=\sum_{l=0}^{n-2k}\sum_{j=0}^{l}\frac{(-1)^{l-j}\binom{l}{j}\binom{n-2k}{l}2^{l}\beta^{l-j}}{(\beta-\alpha)^{l}}x^{j},$

the shifted Legendre basis
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)$ can be written
in terms of a power series in $x$ as

(2.4)
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=\sum_{k=0}^{[\frac{n}{2}]}\sum_{l=0}^{n-2k}\sum_{j=0}^{l}\frac{(-1)^{k-j+l}2^{l-n}\beta^{l-j}(2n-2k)!}{(\beta-\alpha)^{l}(n-k)!k!(n-2k-l)!j!(l-j)!}x^{j},$

which satisfy
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(\beta)=1$ and
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(\alpha)=(-1)^{k}$.

The next lemma describes integer order derivatives of the shifted Legendre
basis.

######  Lemma 2.1.

For $m\leq j\leq n$, we have

(2.5)
$D^{m}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=\sum_{k=0}^{[\frac{n}{2}]}\sum_{l=m}^{n-2k}\sum_{j=m}^{l}\frac{(-1)^{k-j+l}2^{l-n}(\beta-\alpha)^{-l}\beta^{l-j}(2n-2k)!}{(n-k)!k!(n-2k-l)!(l-j)!(j-m)!}x^{j-m}.$

###### Proof.

The proof is easily obtained from Eq. (2.4). ∎

The shifted Legendre polynomials satisfy the following relation

(2.6)
$\int_{\alpha}^{\beta}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,i}}(x)L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)\mathrm{d}x=\delta_{ik}(k+\frac{1}{2})^{-1},$

i.e., they are orthogonal with each other concerning the unit weight function.

The shifted Legendre basis form an orthogonal system of polynomials, which is
complete in the space of square integrable functions, i.e.,
${\mathcal{L}}^{2}(\alpha,\beta)$. Therefore, any
$u\in{\mathcal{L}}^{2}(\alpha,\beta)$ can be written as

$u(x)=\sum_{k=0}^{\infty}\hat{u}_{\scriptscriptstyle\scriptscriptstyle{k}}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x),$

where

$\hat{u}_{\scriptscriptstyle\scriptscriptstyle{k}}=(k+\frac{1}{2})\int_{\alpha}^{\beta}u(x)L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,k}}(x)\mathrm{d}x,\quad
k\geq 0.$

The associated inner product and norm are given by

$(f,g)=\int_{\alpha}^{\beta}f(x)g(x)\mathrm{d}x$

and

$||f||^{2}_{{\mathcal{L}}^{2}(\alpha,\beta)}=(f,f)=\int_{\alpha}^{\beta}|f(t)|^{2}\mathrm{d}t.$

We define $H^{m}(\alpha,\beta)$ to be the vector space of the functions
$g\in{\mathcal{L}}^{2}(\alpha,\beta)$ such that all the distributional
derivatives of $f$ of order up to $m$ can be represented by functions in
${\mathcal{L}}^{2}(\alpha,\beta)$. $H^{m}(\alpha,\beta)$ is endowed with the
norm

$||f||^{2}_{H^{m}(\alpha,\beta)}=\sum_{k=0}^{m}\Big{|}\Big{|}\frac{\partial^{k}}{\partial
x^{k}}f(x)\Big{|}\Big{|}^{2}_{{\mathcal{L}}^{2}(\alpha,\beta)}.$

Furthermore, the associated semi-norm is defined as follows

$|f|^{2}_{H^{m:N}(\alpha,\beta)}=\sum_{j=\min(m,N)}^{m}||f^{(j)}||^{2}_{{\mathcal{L}}^{2}(\alpha,\beta)},$

where $N$ is the number of nodal bases.

In the later sections of this paper, we will use the Gaussian integration
formula to approximate integrals such as

(2.7) $\int_{\alpha}^{\beta}f(x)\mathrm{d}x.$

Explicit formulas for the quadrature nodes and weights for discrete shifted
Chebyshev and Legendre basis are [6]

  * •

Chebyshev Gauss-Lobatto
For $j=0,1,\cdots,N$,

(2.8)
$x_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N,j}}=\frac{\beta-\alpha}{2}(x_{\scriptscriptstyle\scriptscriptstyle{N,j}}-1)+\beta,w_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N,j}}=\begin{cases}\frac{\pi}{2N},&j=0,N,\\\
\frac{\pi}{N},&j=1,\cdots,N-1,\end{cases}$

where

$x_{\scriptscriptstyle\scriptscriptstyle{N,j}}=\cos\frac{j\pi}{N}.$

  * •

Legendre Gauss-Lobatto

(2.9)
$x_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,0}}=\alpha,\,x_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N}}=\beta,\,x_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}(j=1,2,\cdots,N-1){\mbox{
roots of
}}L^{{}^{\prime}}_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N}}(x),$

and

$w_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N,j}}=\frac{2}{N(N+1)[L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N}}(x_{\scriptscriptstyle\scriptscriptstyle{j}})]^{2}},\quad
j=0,1,\cdots,N.$

For any $p(x)\in\mathbb{P}_{2N+1}$, where $\mathbb{P}_{2N+1}$ is the space of
polynomials of degree at most $2N+1$, we have

$\displaystyle\int_{\alpha}^{\beta}\frac{p(x)}{\sqrt{1-x^{2}}}\mathrm{d}x$
$\displaystyle=\sum_{j=0}^{N}p(x_{\scriptscriptstyle\scriptscriptstyle{j}})w_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N,j}},$
$\displaystyle(\text{Chebyshev Gauss-Lobatto}),$
$\displaystyle{\mbox{and}}\quad\int_{\alpha}^{\beta}p(x)\mathrm{d}x$
$\displaystyle=\sum_{j=0}^{N}p(x_{\scriptscriptstyle\scriptscriptstyle{j}})w_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,N,j}},$
$\displaystyle(\text{Legendre Gauss-Lobatto}).$

##  3\. Fractional derivatives of collocation bases

The shifted Chebyshev basis’ fractional derivatives have been discussed in
Ref. [1]. This section continues with obtaining a novel formula for fractional
derivatives of shifted Legendre basis
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)$ in the Caputo
sense [22]. Similar to the shifted Chebyshev basis [19] one can find the
following lemma and theorem.

######  Lemma 3.1.

Let $\nu$ be a positive real number. Then fractional derivative of order $\nu$
of shifted Legendre polynomials
$L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)$ can be given by

(3.1)
$D^{\nu}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=0,\quad
n=0,1,\cdots,\lceil\nu\rceil-1.$

######  Theorem 3.2.

The fractional derivative of order $\nu$ of the shifted Legendre basis is

(3.2)
$D^{\nu}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=\sum_{i=0}^{\infty}S_{\nu}(n,i)L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,i}}(x),\quad
n=\lceil\nu\rceil,\lceil\nu\rceil+1,\cdots,$

where the ceiling function $\lceil\nu\rceil$ stands for the smallest integer
greater than or equal to $\nu$ and

(3.3)
$S_{\nu}(n,i)=\displaystyle\sum_{k=0}^{[\frac{n}{2}]}\sum_{l=\lceil\nu\rceil}^{n-2k}\sum_{j=\lceil\nu\rceil}^{l}\frac{(-1)^{k-j+l}2^{l-n}\beta^{l-j}(\beta-\alpha)^{-l}(2n-2k)!j!}{(n-k)!k!(n-2k-l)!j!(l-j)!\Gamma(j+1-\nu)}c_{ij},$

in which

(3.4)
$c_{\scriptscriptstyle\scriptscriptstyle{ij}}=(i+\frac{1}{2})\int_{\alpha}^{\beta}x^{j-\nu}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,i}}(x)\mathrm{d}x.$

###### Proof.

As we know, the Caputo fractional derivative of $x^{k}$ of order $\nu$ is

$D^{\nu}x^{k}=\left\\{\begin{array}[]{ll}0,&k=0,1,\cdots\text{ and}\
k<\lceil\nu\rceil,\\\
\frac{\Gamma(k+1)}{\Gamma(k+1-\nu)}x^{k-\nu},&k=0,1,\cdots\text{ and}\
k\geq\lceil\nu\rceil.\end{array}\right.$

Considering (2.4), for $n=\lceil\nu\rceil,\lceil\nu\rceil+1,\cdots,$ we obtain

$\displaystyle
D^{\nu}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,n}}(x)=$
$\displaystyle\sum_{k=0}^{[\frac{n}{2}]}\sum_{l=0}^{n-2k}\sum_{j=0}^{l}\frac{(-1)^{k-j+l}2^{l-n}\beta^{l-j}(\beta-\alpha)^{-l}(2n-2k)!}{(n-k)!k!(n-2k-l)!j!(l-j)!}D^{\nu}x^{j}$
$\displaystyle=$
$\displaystyle\sum_{k=0}^{[\frac{n}{2}]}\sum_{l=\lceil\nu\rceil}^{n-2k}\sum_{j=\lceil\nu\rceil}^{l}\frac{(-1)^{k-j+l}2^{l-n}\beta^{l-j}(\beta-\alpha)^{-l}(2n-2k)!}{(n-k)!k!(n-2k-l)!(l-j)!\Gamma(j+1-\nu)}x^{j-\nu}.$

By expanding $x^{j-\nu}$ in terms of the shifted Legendre basis, we arrive at
the following form

$x^{j-\nu}=\sum_{i=0}^{\infty}c_{\scriptscriptstyle\scriptscriptstyle{ij}}L_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}(x),$

where $c_{ij}$ is given in (3.4), which completes the proof. ∎

##  4\. Mixed steps-collocation schemes

In this section, we propose new numerical schemes based on the method of steps
and Legendre (Chebyshev) collocation method to solve the following nonlinear
FDDE

(4.1)
$\left\\{\begin{array}[]{lc}\sum\limits_{j=0}^{n}A_{j}u^{(j)}(t)+D^{\nu}u(t)+\sum_{i=1}^{l}\lambda_{i}D^{\nu_{i}}u(t)=f(t,u(t),u(t-\tau)),&t\geq
0,\\\ u(t)=\phi(t),&-\tau\leq t\leq 0,\\\
u^{\prime}(0)=\phi^{\prime}(0)=\phi_{1},&\\\ \hskip 19.91684pt\vdots&\\\
u^{(\mu-1)}(0)=\phi^{(\mu-1)}(0)=\phi_{\mu-1},&\end{array}\right.$

where $\lambda_{i}$, $A_{j}\in\mathbb{R}$ are constants and $A_{n}\neq 0$,
$0<\nu_{1}<\nu_{2}<\cdots<\nu_{l}<\nu$, $m-1<\nu\leq m$ are real constants,
$D^{\nu}u(t)\equiv u^{\nu}(t)$ represents the Caputo fractional derivative of
order $\nu$ of function $u(t)$, $\mu=\max\\{m,n\\}$ and the function $f$ is
given nonlinear continuous function in $u$ which satisfies the following
Lipschitz conditions

(4.2) $\displaystyle|f(t,y_{1},u)-f(t,y_{2},u)|$ $\displaystyle\leq
L_{1}|y_{1}-y_{2}|,$ (4.3) $\displaystyle|f(t,y,u_{1})-f(t,y,u_{2})|$
$\displaystyle\leq L_{2}|u_{1}-u_{2}|.$

Also, throughout this paper, we shall assume the initial function $\phi(t)$ to
be continuous on $[-\tau,0]$. These conditions guarantee the existence and
uniqueness of the solution of problem (4.1) (see, e.g., [5, 31]).

Clearly, for $t\in[0,\tau]$, the nonlinear FDDE (4.1) is equivalent to the
following nonlinear non-delay FDEs

(4.4)
$\left\\{\begin{array}[]{lc}\sum\limits_{j=0}^{n}A_{j}u^{(j)}(t)+D^{\nu}u(t)+\sum_{i=1}^{l}\lambda_{i}D^{\nu_{i}}u(t)=f(t,u(t),\phi(t-\tau)),&t\in(0,\tau],\\\
u(0)=\phi(0)=\phi_{0},&\\\ u^{\prime}(0)=\phi^{\prime}(0)=\phi_{1},&\\\ \hskip
19.91684pt\vdots&\\\
u^{(\mu-1)}(0)=\phi^{(\mu-1)}(0)=\phi_{\mu-1}.&\end{array}\right.$

One way of solving the nonlinear FDE (4.4) is to use shifted Chebyshev basis,
which is an extension of the scheme presented in Ref. [19]. Another way is to
expand the approximate solution
$u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)$ in terms of truncated shifted
Legendre basis. The latter idea leads to

(4.5)
$u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)=\sum_{k=0}^{N}a_{\scriptscriptstyle\scriptscriptstyle{k}}L_{\scriptscriptstyle\scriptscriptstyle{0,\tau,k}}(t),\qquad
t\in[0,\tau],$

where $a_{\scriptscriptstyle\scriptscriptstyle{k}}\in\mathbb{R}$ are the
unknown coefficients that we aim to find them. Thanks to Theorem 3.2, we can
express the derivatives $u^{(1)}(t)$, $u^{(2)}(t)$, $\cdots,u^{(n)}(t)$,
$D^{\nu}u(t)$, $D^{\nu_{1}}u(t)$, $D^{\nu_{2}}u(t),\cdots,D^{\nu_{l}}u(t)$ in
terms of unknowns coefficients $a_{\scriptscriptstyle\scriptscriptstyle{k}}$.

Now, we employ the Legendre (Chebyshev) collocation method to solve (4.4)
numerically. To do this, the following equation

(4.6)
$\sum\limits_{j=0}^{n}A_{j}u{\scriptscriptstyle\scriptscriptstyle{N}}^{(j)}(t)+D^{\nu}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)+\sum_{i=1}^{l}\lambda_{i}D^{\nu_{i}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)=f(t,u_{\scriptscriptstyle\scriptscriptstyle{N}}(t),\phi(t-\tau)),$

must be satisfied at the shifted Legendre collocation nodes (2.9) (shifted
Chebyshev collocation nodes (2.8)) exactly. In fact, by using (3.2) and (4.5),
for $j=0,1,\cdots,N-\mu$, we get the following equations

$\displaystyle\sum\limits_{j=0}^{n}\sum_{k=0}^{N}A_{j}a_{\scriptscriptstyle\scriptscriptstyle{k}}L^{(j)}_{\scriptscriptstyle\scriptscriptstyle{0,\tau,k}}(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}})+\sum_{k=0}^{N}a_{\scriptscriptstyle\scriptscriptstyle{k}}D^{\nu}$
$\displaystyle
L_{\scriptscriptstyle\scriptscriptstyle{0,\tau,k}}(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}})+\sum_{i=1}^{l}\sum_{k=0}^{N}\lambda_{i}a_{\scriptscriptstyle\scriptscriptstyle{k}}D^{\nu_{i}}L_{\scriptscriptstyle\scriptscriptstyle{0,\tau,k}}(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}})$
(4.7)
$\displaystyle=f(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}},\sum_{k=0}^{N}a_{\scriptscriptstyle\scriptscriptstyle{k}}L_{\scriptscriptstyle\scriptscriptstyle{0,\tau,k}}$
$\displaystyle(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}),\phi(t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}-\tau)),\
t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}\geq 0,$

where $t_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}$ are the same
as $x_{\scriptscriptstyle\scriptscriptstyle{\alpha,\beta,j}}$ which are
defined by (2.9). After imposing the initial conditions

(4.8) $u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(i)}(0)=\phi_{i},\quad
i=0,1,\cdots,\mu-1,$

we arrive at a nonlinear system of algebraic equations. Similarly, using
Theorem 1 in Ref. [19] and related shifted Chebyshev expansion, we obtain an
algebraic system of nonlinear equations. The nonlinear resulted systems can be
solved, for example, by Newton’s method. Therefore, the approximate solution
$u_{\scriptscriptstyle\scriptscriptstyle{N}}$ in the interval $[0,\tau]$ is
now available. To obtain the approximate solution of Eq. (4.1) in
$[\tau,2\tau]$, the presented procedure is used. Generally, if we want to
solve Eq. (4.1) in the interval $[(i-1)\tau,i\tau]$, $i\geq 1$, we need to
solve the following equation

(4.9)
$\sum\limits_{j=0}^{n}A_{j}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u^{(j)}(t)+D^{\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)+\sum_{p=1}^{l}\lambda_{p}D^{\nu_{p}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)=f(t,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t),\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i-1)}}}u(t-\tau)),$

with the initial conditions

$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u^{(j)}(0)=\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i-1)}}}u^{(j)}(\tau),\quad
j=0,1,\cdots,\mu,$

where

$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)=u((i-1)\tau+t),\quad
i\geq 1.$

Now, using the proposed procedure, we get the approximate solution
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}$
of Eq. (4.9).

##  5\. Convergence analysis

In this section, by a similar manner presented in Ref. [9], we show that the
obtained approximate solutions in the previous section are convergent to the
exact solutions. In order to investigate the exponential rate of convergence
of the proposed schemes, we consider the nonlinear FDDE (4.9) on the interval
$\Omega_{i}=[(i-1)\tau,i\tau]$.

Let us define
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)=\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)-\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)$
to be the error function of the proposed scheme, where
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)$ and
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)$
are the exact and Legendre (Chebyshev) collocation solution of (4.9) at $i-th$
step, respectively.

Hereafter, we use the subscript $w$ for the Legendre and Chebyshev weight
functions.

The orthogonal projection operator
$P_{\scriptscriptstyle\scriptscriptstyle{N}}$ from
${\mathcal{L}}_{\scriptscriptstyle\scriptscriptstyle{w}}^{2}(\Omega)$ onto
$\mathbb{P}_{\scriptscriptstyle\scriptscriptstyle{N}}$, where
$\Omega=[\alpha,\beta]$, satisfies

$\forall\psi_{N}\in\mathbb{P}_{\scriptscriptstyle\scriptscriptstyle{N}},\quad\int_{0}^{\tau}(f-P_{\scriptscriptstyle\scriptscriptstyle{N}}f)(r)dr=0,$

for any function $f$ in
${\mathcal{L}}_{\scriptscriptstyle\scriptscriptstyle{w}}^{2}(\Omega)$. It is
clear that $P_{\scriptscriptstyle\scriptscriptstyle{N}}$ belongs to
$\mathbb{P}_{\scriptscriptstyle\scriptscriptstyle{N}}$.

The following inequalities for the shifted Legendre (Chebyshev) polynomials
and shifted Legendre (Chebyshev)-Gauss-Lobatto nodes for $k\geq 1$ can be
obtained by a similar argument provided in Ref. [3]

(5.1)
$\displaystyle||y-P_{\scriptscriptstyle\scriptscriptstyle{N}}(y)||_{H^{l}_{w}(\Omega)}\leq$
$\displaystyle CN^{2l-1/2-k}|y|_{H^{k:N}_{w}(\Omega)},$

where $y\in H_{w}^{k}(\Omega)$.

######  Lemma 5.1 ([11]).

Let $F$ be a continuous function on the real line. Then the $n$-fold multiple
integral of $F$ based at $a$ is given by

$\int_{a}^{t}\int_{a}^{t_{n}}\cdots\int_{a}^{t_{2}}F(t_{1})\mathrm{d}T=\frac{1}{(n-1)!}\int_{a}^{t}(t-s)^{n-1}F(s)\mathrm{d}s,$

where $\mathrm{d}T=\mathrm{d}t_{1}\mathrm{d}t_{2}\cdots\mathrm{d}t_{n}$.

Now, we present the convergence theorem of the proposed schemes.

######  Theorem 5.2.

Suppose that the exact solution
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)$ at $i-th$ step
of Eq. (4.9) is smooth enough, i.e.
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)\in
H^{k}_{w}(\Omega)$ for $i,k\geq 1$, and the corresponding mixed steps-
collocation solution
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)$
is given by shifted Legendre or Chebyshev basis. Then for sufficiently large
$N$ we have

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}\leq$
$\displaystyle
C_{1}N^{-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}(\Omega)}+C_{2}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{1:N}(\Omega)}$
(5.2)
$\displaystyle+C_{3}N^{2(\eta_{1}-1)-1/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}}$
$\displaystyle+C_{4}N^{2(\eta_{2}-1)-1/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}}+C_{5}N^{-3/2},$

where

$\eta_{1}=\begin{cases}m,&m\leq 3,\\\
3,&m>3,\end{cases}\qquad\eta_{2}=\begin{cases}m_{p},&m_{p}\leq 3,\\\
3,&m_{p}>3,\end{cases}$

and the constants $C_{i}$ are independent of $N$ and only depend on $n,m$ and
$\nu$.

###### Proof.

As
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)$
is the mixed steps-collocation solution of Eq. (4.9) on the interval
$\Omega_{i}$, it satisfies the following equation

$\displaystyle\sum_{j=0}^{n}A_{j}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(j)}(t)+P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\frac{1}{\Gamma(m-\nu)}$
$\displaystyle\int_{(i-1)\tau}^{t}(t-s)^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s)\mathrm{d}s\big{)}$
$\displaystyle+P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\sum_{p=1}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}$
$\displaystyle\int_{(i-1)\tau}^{t}(t-s)^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s)\mathrm{d}s\big{)}$
$\displaystyle=P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}f(t,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t),$
$\displaystyle\leftidx{{}_{(i-1)}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t))\big{)}.$

By $n$-times integration of the above expression, we have

$\displaystyle\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{t_{n}}\cdots\int_{(i-1)\tau}^{t_{2}}A_{n}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(n)}(t_{1})\mathrm{d}T$
$\displaystyle+\sum_{j=0}^{n-1}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{t_{n}}\cdots\int_{(i-1)\tau}^{t_{2}}A_{j}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(j)}(t_{1})\mathrm{d}T$
$\displaystyle+\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{t_{n}}\cdots\int_{(i-1)\tau}^{t_{2}}P_{\scriptscriptstyle\scriptscriptstyle{N}}$
$\displaystyle\big{(}\int_{(i-1)\tau}^{t_{1}}(t_{1}-s)^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s)\mathrm{d}s\big{)}\mathrm{d}T$
$\displaystyle+\sum_{p=1}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{t_{n}}\cdots\int_{(i-1)\tau}^{t_{2}}$
$\displaystyle
P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\int_{(i-1)\tau}^{t_{1}}(t_{1}-s)^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s)\mathrm{d}s\big{)}\mathrm{d}T$
(5.3)
$\displaystyle=\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{t_{n}}\cdots\int_{(i-1)\tau}^{t_{2}}P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}f(t_{1},\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t_{1}),$
$\displaystyle\leftidx{{}_{(i-1)}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t_{1}))\big{)}\mathrm{d}T.$

Thanks to Lemma 5.1, we can rewrite each of multiple integrals appeared in
(5.3) as a single integral [9]

$\displaystyle
A_{n}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)+Q_{0}(t)+\sum_{j=0}^{n-1}\int_{(i-1)\tau}^{t}\frac{A_{j}}{(n-1-j)!}(t-s)^{n-1-j}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)\mathrm{d}s$
$\displaystyle+\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\big{)}\mathrm{d}s$
$\displaystyle+\sum_{p=1}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\big{)}\mathrm{d}s$
(5.4)
$\displaystyle=\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}f(s,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s),\leftidx{{}_{(i-1)}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s))\big{)}\mathrm{d}s,$

where $Q_{0}(t)$ contains initial conditions. Similarly, the exact solution
$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)$ satisfies the
following relation

$\displaystyle
A_{n}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(t)+Q_{0}(t)+\sum_{j=0}^{n-1}\int_{(i-1)\tau}^{t}\frac{A_{j}}{(n-1-j)!}(t-s)^{n-1-j}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)\mathrm{d}s$
$\displaystyle+\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle+\sum_{p=0}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.5)
$\displaystyle=\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}f(s,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s),\leftidx{{}_{(i-1)}}u(s))\mathrm{d}s.$

Subtracting (5.5) from (5.4) leads to

$\displaystyle
A_{n}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)$
$\displaystyle+\sum_{j=0}^{n-1}\frac{A_{j}}{(n-1-j)!}\int_{(i-1)\tau}^{t}(t-s)^{n-1-j}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)\mathrm{d}s$
$\displaystyle+\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P_{\scriptscriptstyle\scriptscriptstyle{N}}}\mathrm{d}s$
$\displaystyle+\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle+\sum_{p=1}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}\mathrm{d}s$
$\displaystyle+\sum_{p=1}^{l}\frac{\lambda_{p}}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.6)
$\displaystyle=\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)\mathrm{d}s,$

where

$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)$
$\displaystyle=P_{\scriptscriptstyle\scriptscriptstyle{N}}(f(s,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s),\leftidx{{}_{(i-1)}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)))-f(s,\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s),\leftidx{{}_{(i-1)}}u(s)),$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)$
$\displaystyle=P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\big{)}$
$\displaystyle\qquad-\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1},$

and

$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)$
$\displaystyle=P_{\scriptscriptstyle\scriptscriptstyle{N}}\big{(}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\big{)}$
$\displaystyle\qquad-\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}.$

After $(n-1)$-times integrating by parts of the fourth and sixth term of the
left-hand side of (5.6) we get (see Apendix A)

$\displaystyle\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}$
$\displaystyle\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.7) $\displaystyle=\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}$
$\displaystyle\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s,$

and

$\displaystyle\frac{1}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}$
$\displaystyle\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.8)
$\displaystyle=\frac{1}{\Gamma(m_{p}-\nu_{p})\prod_{j=1}^{n-1}(m_{p}-\nu_{p}+j-1)}$
$\displaystyle\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}+n-\nu_{p}-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s.$

Now we consider the three different cases:

(i):

$m\geq 4$ and $m_{p}\geq 4$,

(ii):

$m\geq 4$ and $m_{p}\leq 4$,

(iii):

$m\leq 4$.

Case $(i)$: $(m-3)$-times integrating by parts of right-hand side of equations
(5.7) and (5.8), for $n+1>\nu$, gives

$\displaystyle\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-2}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.9)
$\displaystyle=\frac{\prod_{j=1}^{m-3}(m+n-\nu-1-j)}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{n-\nu+1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s,$

and

$\displaystyle\frac{1}{\Gamma(m_{p}-\nu_{p})\prod_{j=1}^{n-1}(m_{p}-\nu_{p}+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}+n-\nu_{p}-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.10)
$\displaystyle=\frac{\prod_{j=1}^{m_{p}-3}(m_{p}+n-\nu_{p}-1-j)}{\Gamma(m_{p}-\nu_{p})\prod_{j=1}^{n-1}(m_{p}-\nu_{p}+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s.$

Substituting the right-hand side of (5.9) into the right-hand side of (5.7),
we obtain

$\displaystyle\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s$ (5.11)
$\displaystyle=\frac{\prod_{j=1}^{m-3}(m+n-\nu-1-j)}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{n-\nu+1}$
$\displaystyle\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s,$

and similarly, from equations (5.8) and (5.10) we have

$\displaystyle\frac{1}{\Gamma(m_{p}-\nu_{p})}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
(5.12)
$\displaystyle=\frac{\prod_{j=1}^{m_{p}-3}(m_{p}+n-\nu_{p}-1-j)}{\Gamma(m_{p}-\nu_{p})\prod_{j=1}^{n-1}(m_{p}-\nu_{p}+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s.$

Substituting (5.11) and (5.12) into (5.6), we arrive at

$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)|\leq$
$\displaystyle\sum_{j=0}^{n-1}\Big{|}\frac{A_{j}}{A_{n}(n-1-j)!}\Big{|}\int_{(i-1)\tau}^{t}\Big{|}(t-s)^{n-1-j}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\Big{|}\frac{1}{A_{n}(n-1)!\Gamma(m-\nu)}\Big{|}\int_{(i-1)\tau}^{t}\Big{|}(t-s)^{n-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\Big{|}\frac{\prod_{j=1}^{m-3}(m+n-\nu-1-j)}{A_{n}\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\Big{|}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}\Big{|}(s-s_{1})^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\Big{|}\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle+\sum_{p=1}^{l}\Big{|}\frac{\lambda_{p}}{A_{n}\Gamma(m_{p}-\nu_{p})}\Big{|}\int_{(i-1)\tau}^{t}\Big{|}(t-s)^{n-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\sum_{p=1}^{l}\Big{|}\frac{\prod_{j=1}^{m_{p}-3}(m_{p}+n-\nu_{p}-1-j)}{A_{n}\Gamma(m_{p}-\nu_{p})\prod_{j=1}^{n-1}(m_{p}-\nu_{p}+j-1)}\Big{|}\int_{(i-1)\tau}^{s}\Big{|}(s-s_{1})^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})\Big{|}\mathrm{d}s_{1}\mathrm{d}s$
(5.13)
$\displaystyle+\int_{(i-1)\tau}^{t}\Big{|}\frac{(t-s)^{n-1}}{A_{n}(n-1)!}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)\Big{|}\mathrm{d}s.$

We may rewrite Eq. (5.13) as

$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)|\leq$
$\displaystyle\gamma_{1}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|\mathrm{d}s+\gamma_{2}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)|\mathrm{d}s+\gamma_{3}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\gamma_{4}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)|\mathrm{d}s+\gamma_{5}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
(5.14)
$\displaystyle\qquad+\gamma_{6}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)|\mathrm{d}s,$

where $\gamma_{i}$ are independent of $N$ and only may depend on $n$, $m$,
$m_{p}$ and $\nu$.

By the Gronwall lemma [29] we get

$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)|\leq$
$\displaystyle
e^{(\int_{(i-1)\tau}^{t}(\gamma_{1})\mathrm{d}s}\Big{(}\gamma_{2}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)|\mathrm{d}s+\gamma_{3}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\gamma_{4}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)|\mathrm{d}s+\gamma_{5}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\gamma_{6}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)|\mathrm{d}s\Big{)}$
$\displaystyle\leq$
$\displaystyle\gamma_{7}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)|\mathrm{d}s+\gamma_{8}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
$\displaystyle+\gamma_{9}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)|\mathrm{d}s+\gamma_{10}\sum_{p=1}^{l}\int_{(i-1)\tau}^{t}\Big{|}s^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)\Big{|}\mathrm{d}s$
$\displaystyle\qquad+\gamma_{11}\int_{(i-1)\tau}^{t}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{f}(s)|\mathrm{d}s,$

where $\gamma_{7},\cdots,\gamma_{11}$ are some constants which related to
$\gamma_{1},\cdots,\gamma_{6}$.

From Lipschitz conditions (4.2) and (4.3), inequality (5.1) and generalized
Hardy’s inequality [10] we obtain

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}\leq$
$\displaystyle\gamma_{12}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}+\gamma_{13}||s^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle+\gamma_{14}\sum_{p=1}^{l}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
(5.15)
$\displaystyle+\gamma_{15}\sum_{p=1}^{l}||s^{n-\nu_{p}+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}+C_{5}N^{-3/2},$

where $\gamma_{12},\cdots,\gamma_{15}$ are some constants which related to
$\gamma_{7},\cdots,\gamma_{11}$ and are independent of $N$. From (5.1) we get

$\displaystyle||s^{n-\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq||s^{n-\nu}||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq\gamma_{16}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)||_{H^{2}_{w}(\Omega)}$
(5.16)
$\displaystyle\leq\gamma_{17}N^{7/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{k:N}_{w}(\Omega)},$

and

$\displaystyle||s^{n-\nu_{p}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq||s^{n-\nu_{p}}||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{{}^{\prime\prime}}(s)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq\gamma_{18}||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)||_{H^{2}_{w}(\Omega)}$
(5.17)
$\displaystyle\leq\gamma_{19}N^{7/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{k:N}_{w}(\Omega)},$

where $\gamma_{16},\cdots,\gamma_{19}$ do not depend on $N$. From (5.1) we
have

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P_{\scriptscriptstyle\scriptscriptstyle{N}}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq
C_{6}N^{-3/2}\Big{|}\int_{(i-1)\tau}^{t}(t-s)^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s)\mathrm{d}s\Big{|}_{H^{1:N}_{w}(\Omega)}$
(5.18)
$\displaystyle=C_{7}N^{-3/2}|D^{\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{1:N}_{w}(\Omega)},$
$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq
C_{8}N^{-3/2}\Big{|}\int_{(i-1)\tau}^{t}(t-s)^{m_{p}-\nu_{p}-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m_{p})}(s)\mathrm{d}s\Big{|}_{H^{1:N}_{w}(\Omega)}$
(5.19)
$\displaystyle=C_{9}N^{-3/2}|D^{\nu_{p}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{1:N}_{w}(\Omega)}.$

Linear operators $D^{\nu}:\mathbb{P}_{N}\to\mathbb{P}_{N}^{\nu}$ and
$D^{\nu_{p}}:\mathbb{P}_{N}\to\mathbb{P}_{N}^{\nu_{p}}$ are bounded (see
Appendix C), so that the constants $C_{10}$ and $C_{11}^{p}$ can be found such
that

(5.20)
$|D^{\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}(\Omega)}\leq
C_{10}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}(\Omega)},$

and

(5.21)
$|D^{\nu_{p}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}(\Omega)}\leq
C_{11}^{p}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}(\Omega)}.$

Therefore, from (5.18) and (5.20) we have

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq
C_{12}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|_{H^{1:N}_{w}(\Omega)}$
$\displaystyle=C_{12}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)-\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)+\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)|_{H^{1:N}_{w}(\Omega)}$
(5.22) $\displaystyle\leq
C_{12}N^{-3/2}\big{(}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|_{H^{1:N}_{w}(\Omega)}+|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)|_{H^{1:N}_{w}(\Omega)}\big{)}.$

Since
$u_{\scriptscriptstyle\scriptscriptstyle{N}}=P_{\scriptscriptstyle\scriptscriptstyle{N}}u$,
we can write

$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|_{H^{1:N}_{w}(\Omega)}=$
$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)-\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)+P_{\scriptscriptstyle\scriptscriptstyle{N}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)-P_{\scriptscriptstyle\scriptscriptstyle{N}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)|_{H^{1:N}_{w}(\Omega)}$
$\displaystyle=$
$\displaystyle|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)-P_{\scriptscriptstyle\scriptscriptstyle{N}}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u(s)|_{H^{1:N}_{w}(\Omega)}.$

So that

(5.23)
$||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{P}_{\scriptscriptstyle\scriptscriptstyle{N}}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}{\overset{\scriptscriptstyle\scriptscriptstyle{\eqref{hef}}}{\leq}}C_{13}N^{-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{k:N}_{w}(\Omega)}+C_{12}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{1:N}_{w}(\Omega)}.$

From (5.19) and (5.21) we have

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{P^{p}_{N}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}$
$\displaystyle\leq
C_{14}^{p}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|_{H^{1:N}_{w}(\Omega)}$
$\displaystyle\leq
C_{14}^{p}N^{-3/2}\big{(}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(s)|_{H^{1:N}_{w}(\Omega)}+|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{1:N}_{w}(\Omega)}\big{)}$
(5.24)
$\displaystyle{\overset{\scriptscriptstyle\scriptscriptstyle{\eqref{hef}}}{\leq}}C_{15}^{p}N^{-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{k:N}_{w}(\Omega)}+C_{14}^{p}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H^{1:N}_{w}(\Omega)}.$

From Eqs. (5)-(5.24), for $m\geq 4$ and $m_{p}\geq 4$, we have

$\displaystyle||\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}(t)||_{{\mathcal{L}^{2}_{\scriptscriptstyle\scriptscriptstyle{w}}}(\Omega)}\leq$
$\displaystyle
C_{1}N^{-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}(\Omega)}+C_{2}N^{-3/2}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{1:N}(\Omega)}$
$\displaystyle+C_{3}N^{2(\eta_{1}-1)-1/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}}$
$\displaystyle+C_{4}N^{2(\eta_{2}-1)-1/2-k}|\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}u|_{H_{\scriptscriptstyle\scriptscriptstyle{w}}^{k:N}}+C_{5}N^{-3/2},$

where $C_{1},\cdots,C_{5}$ are constants which related to earlier
$\gamma_{i}$’s and $C_{i}$’s.

The same argument can be done for the cases $(ii)$ and $(iii)$ by assuming
that $n+1>\nu$, $m_{p}-\nu_{p}+n>2$, and $m-\nu+n>2$ respectively. ∎

##  6\. Numerical results

In this section, we consider several practical examples which, in general, do
not have an exact solution. Here, we made use of the $Matlab~{}2014a$ package.
The computational codes were conducted on an Intel(R) Core(TM) i7-6700K
processor, equipped with 8 GB of RAM. Also, We use the fix-point iteration
method for solving nonlinear systems, and stopping criterion is set to be
$10^{-15}$.

Example 1. Consider the following FDDE

(6.1) $D^{\nu}u(t)=h(t)-u(t)-u(t-\tau),\quad t\in(0,T],\\\ $

with the boundary condition

$u(t)=0,\quad t\in[-\tau,0],$

where $\tau$ is taken as a fraction of the length of time interval $[0,1]$.
Now, two different cases for the forcing term $h(t)$ may be considered:
Case $(i)$:

(6.2)
$h(t)=\begin{cases}\frac{\Gamma(11)}{\Gamma(11-\nu)}t^{10-\nu}+t^{10},&t\in[0,\tau],\\\
\frac{\Gamma(11)}{\Gamma(11-\nu)}t^{10-\nu}+t^{10}+(t-\tau)^{10},&t\in(\tau,1],\end{cases}$

which corresponding exact solution is

(6.3) $u(t)=\begin{cases}0,&t\in[-\tau,0],\\\ t^{10},&t\in(0,1].\end{cases}$

Case $(ii)$:

(6.4)
$h(t)=\begin{cases}\sum_{j=0}^{\infty}\gamma_{j}\frac{\Gamma(\beta_{j})}{\Gamma(\beta_{j}-\nu)}t^{\xi_{j}-\nu}+t^{\frac{13}{2}}\sin(\pi
t^{\frac{4}{3}}),&t\in[0,\tau],\\\
\sum_{j=0}^{\infty}\gamma_{j}\frac{\Gamma(\beta_{j})}{\Gamma(\beta_{j}-\nu)}t^{\xi_{j}-\nu}+t^{\frac{13}{2}}\sin(\pi
t^{\frac{4}{3}})+(t-\tau)^{\frac{13}{2}}\sin(\pi(t-\tau)^{\frac{4}{3}}),&t\in(\tau,1],\end{cases}$

where

$\displaystyle\gamma_{j}=$ $\displaystyle\frac{(-1)^{j}}{(2j+1)!}\pi^{2j+1},$
$\displaystyle\beta_{j}=$ $\displaystyle\frac{53+16j}{6},$
$\displaystyle\xi_{j}=$ $\displaystyle\frac{47+16j}{6}.$

The related exact solution is

(6.5) $u(t)=\left\\{\begin{array}[]{lr}0,\hskip 64.01869ptt\in[-\tau,0],\\\
t^{\frac{13}{2}}\sin(\pi t^{\frac{4}{3}}),\quad t\in(0,1].\end{array}\right.$

Zayernouri et al. [32] used Petrov-Galerkin spectral method to solve (6.1).
They employed Reimann-Liouville fractional derivatives while we use the
Caputo’s fractional derivatives. As we know, these are related together by the
following relation [18]

$\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{R}}}D^{\nu}u(t)=D^{\nu}u(t)+\sum_{k=0}^{m-1}\frac{t^{k-\nu}}{\Gamma(k+1-\nu)}u^{(k)}(0^{+}),$

where $\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{R}}}D^{\nu}$ stands
for Reimann-Liouville fractional derivative. Since $u(0)=0$, both fractional
derivatives are the same and consequently the resulted approximate solution
are comparable.

The $\mathcal{L}^{2}$-Error of the Case $(i)$ and Case $(ii)$ for $\tau=0.5$
and different values of $N$ and $\nu=0.1$ are reported in Table 1 and Table 2,
respectively. The time responses of the two cases are plotted in Fig. 1 and
Fig. 2. As one can observe, the results are in remarkable agreement with the
results of Zayernouri et al. [32]. The results are obtained by employing
shifted Chebyshev basis are more accurate than the others.

Table 1. $\mathcal{L}^{2}$-Error of Example 1 (Case $(i)$) for different values of $N$ with $\tau=0.5$ and $\nu=0.1$. $N$ | Shifted Legendre Basis | Shifted Chebyshev Basis
---|---|---
3 | $0.871781$ | $0.022220$
5 | $0.268627$ | $0.002870$
7 | $0.020670$ | $1.262008e{-4}$
9 | $1.863461e{-4}$ | $7.422971e{-7}$
11 | $2.232967e{-8}$ | $9.407349e{-9}$
13 | $1.954103e{-8}$ | $1.039491e{-12}$
15 | $1.020807e{-8}$ | $9.904553e{-14}$
17 | $6.124110e{-10}$ | $9.425039e{-15}$
19 | $2.681213e{-12}$ | $1.025904e{-15}$
Table 2. $\mathcal{L}^{2}$-Error of Example 1 (Case $(ii)$) for different values of $N$ with $\tau=0.5$ and $\nu=0.1$. $N$ | Shifted Legendre Basis | Shifted Chebyshev Basis
---|---|---
3 | $0.343192$ | $0.016382$
5 | $0.299374$ | $0.002711$
7 | $0.006547$ | $1.763871e{-4}$
9 | $0.003438$ | $1.300984e{-5}$
11 | $1.724315e{-4}$ | $7.690661e{-7}$
13 | $1.106626e{-6}$ | $8.797766e{-9}$
15 | $4.971592e{-7}$ | $2.055257e{-10}$
17 | $1.892239e{-8}$ | $5.953249e{-12}$
19 | $9.786602e-{11}$ | $1.270578e{-13}$
Figure 1. Time response of Example 1 (Case $(i)$) by the shifted Chebyshev
basis with $N=15$ and $\tau=0.5$. Figure 2. Time response of Example 1 (Case
$(ii)$) by the shifted Chebyshev basis with $N=15$ and $\tau=0.5$.

Example 2. Consider the following FDDE [32]

(6.6) $D^{\nu}u(t)=h(t)-A(t)u(t)-B(t)u(t-\tau),\quad t\in(0,1],\\\ $

with the initial condition

$u(t)=0,\quad t\in[-\tau,0].$

Now, two cases are taken into consideration. Case $(i)$: Take
$A(t)=B(t)=t^{2}-t^{3}$. The corresponding exact solution is given in (6.3).
Case $(ii)$: Put $A(t)=B(t)=\sin(\pi t)$, where the analytical solution is
given in (6.5). The $\mathcal{L}^{2}$-Error of the Case $(i)$ and Case $(ii)$
for different values of $N$, $\tau=0.5$ and $\nu=0.1$ are provided in Table 3
and Table 4, respectively. The $\mathcal{L}^{2}$-Error of the Case $(i)$ in
the current work is at least of order $10^{-13}$ for $N\geq 11$, while it
happened only when $N\geq 17$ in Ref. [32]. In the Case $(ii)$, the results
are in remarkable agreement with the results of Zayernouri et al. [32]. Again,
the results are obtained by employing shifted Chebyshev basis are more
accurate than the others.

Table 3. $\mathcal{L}^{2}$-Error of Example 2 (Case $(i)$) for different values of $N$ with $\tau=0.5$ and $\nu=0.1$. $N$ | Shifted Legendre Basis | Shifted Chebyshev Basis
---|---|---
3 | $0.708611$ | $0.655433$
5 | $0.113928$ | $0.082273$
7 | $0.005367$ | $0.003304$
9 | $3.324243e{-5}$ | $1.697280e{-5}$
11 | $6.758311e{-13}$ | $3.050396e{-14}$
13 | $5.424470e{-13}$ | $3.024794e{-14}$
15 | $3.985062e{-13}$ | $3.019384e{-14}$
17 | $2.606128e{-13}$ | $3.009684e{-14}$
19 | $1.112418e{-13}$ | $2.018803e{-14}$
Table 4. $\mathcal{L}^{2}$-Error of Example 2 (Case $(ii)$) for different values of $N$ with $\tau=0.5$ and $\nu=0.1$. $N$ | Shifted Legendre Basis | Shifted Chebyshev Basis
---|---|---
3 | $0.447450$ | $0.025522$
5 | $0.114932$ | $0.003248$
7 | $0.005888$ | $2.277593e-{4}$
9 | $6.020148e-{4}$ | $1.489728e{-5}$
11 | $3.318842e{-5}$ | $9.040363e{-7}$
13 | $3.284365e{-7}$ | $1.067213e{-8}$
15 | $9.812455e{-9}$ | $2.358068e{-10}$
17 | $3.035429e{-10}$ | $6.944012e{-12}$
19 | $2.880012e{-11}$ | $1.692653e{-13}$

Example 3. Consider Houseflies model as following [15]

(6.7) $D^{\nu}u(t)=-du(t)+cu(t-\tau)(k-czu(t-\tau)),\quad t>0,\\\ $

with the initial condition

$u(t)=160,\quad-\tau\leq t\leq 0.$

By taking $c=1.81$, $k=0.5107$, $d=0.147$ and $z=0.000226$, numerical results
of the shifted Chebyshev basis for different values of $\nu$, $\tau=3$ and
$\tau=5$ are presented in Table 5 and Table 6, respectively. Table 7 and Table
8 describe the numerical results of the shifted Legendre basis with the same
parameters. The approximate solutions are sketched in Fig. 3 and Fig. 4.
Comparison between the second and third columns of the Table 5 and Table 6
(Table 7 and Table 8) reveals that the maximum absolute error (MAE) is
$2\times 10^{-6}$ while the MAE reported in Ref. [15] which employed the
finite difference method was of order $10^{-5}$ using $N=100$. Additionally,
phase-space solution for $\nu=0.5$ and $\tau=5$ is plotted in Fig. 5.
Moreover, log plot of MAE for different values of $N$, $\tau=3$ and $\tau=5$
are plotted in Fig. 6 and Fig. 7, respectively.

Table 5. Numerical results of Example 3 by the shifted Chebyshev basis with $N=15$ and $\tau=3$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$
0.75 | $234.865602$ | $234.865602$ | $239.333361$ | $245.219089$ | $252.133912$
2.25 | $361.967021$ | $361.967021$ | $352.369942$ | $336.687494$ | $308.795525$
3.75 | $481.825305$ | $481.825304$ | $472.760850$ | $459.813354$ | $439.073221$
5.25 | $670.725206$ | $670.725204$ | $650.373599$ | $614.508045$ | $543.799282$
6.0 | $776.578086$ | $776.578085$ | $742.142667$ | $685.419046$ | $583.383909$
Table 6. Numerical results of Example 3 by the shifted Chebyshev basis with $N=15$ and $\tau=5$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$
1.25 | $280.382021$ | $280.382021$ | $281.116516$ | $280.559015$ | $275.428893$
3.75 | $463.917311$ | $463.917311$ | $437.587058$ | $399.856049$ | $342.930615$
6.25 | $636.682069$ | $636.682068$ | $609.967593$ | $571.739341$ | $510.393022$
8.75 | $950.311527$ | $950.311525$ | $890.874639$ | $798.685328$ | $645.755656$
10.0 | $1107.006007$ | $1107.006006$ | $1022.495184$ | $895.594663$ | $695.358908$
Table 7. Numerical results of Example 3 by the shifted Legendre basis with $N=15$ and $\tau=3$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$
0.75 | $234.865602$ | $234.865601$ | $239.333372$ | $245.219085$ | $252.133922$
2.25 | $361.967021$ | $361.967023$ | $352.369941$ | $336.687497$ | $308.795514$
3.75 | $481.825305$ | $481.825304$ | $472.760849$ | $459.813353$ | $439.073219$
5.25 | $670.725206$ | $670.725202$ | $650.373596$ | $614.508046$ | $543.799279$
6.0 | $776.578086$ | $776.578087$ | $742.142661$ | $685.419042$ | $583.383904$
Table 8. Numerical results of Example 3 by the shifted Legendre basis with $N=15$ and $\tau=5$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$ | $160.000000$
1.25 | $280.382021$ | $280.382023$ | $281.116518$ | $280.559021$ | $275.428895$
3.75 | $463.917311$ | $463.917314$ | $437.587055$ | $399.856060$ | $342.930617$
6.25 | $636.682069$ | $636.682069$ | $609.967581$ | $571.739351$ | $510.393025$
8.75 | $950.311527$ | $950.311526$ | $890.874641$ | $798.685332$ | $645.755659$
10.0 | $1107.006007$ | $1107.006008$ | $1022.495191$ | $895.594671$ | $695.358910$
Figure 3. Approximate solution of Example 3 using the shifted Chebyshev basis
with $N=15$ and $\tau=3$. Figure 4. Approximate solution of Example 3 using
the shifted Chebyshev basis with $N=15$ and $\tau=5$. Figure 5. Phase-space
solution of Example 3 using the shifted Chebyshev basis with $N=15$, $\nu=0.5$
and $\tau=5$. Figure 6. Log plot of MAE of Example 3 with $\nu=1$ and
$\tau=3$. Figure 7. Log plot of MAE of Example 3 with $\nu=1$ and $\tau=5$.

Example 4. The following model example concerns with the effect of noise on
light which is introduced by Pieroux [15]

(6.8)
$D^{\nu}u(t)=-\frac{1}{\epsilon}u(t)+\frac{1}{\epsilon}u(t)u(t-\tau),\quad
t>0,\\\ $

with the initial condition

$u(t)=0.9,\quad-\tau\leq t\leq 0.$

The obtained results of the shifted Chebyshev basis for various values of
$\nu$ and $\tau$ with $\epsilon=0.1$ are reported in Table 9 and Table 10.
Also, numerical results of the shifted Legendre basis are given in Table 11
and Table 12. In this model example, the MAEs related to the current works are
order of $10^{-6}$ using $N=15$, while the MAE reported in Ref. [15] which
employed the finite difference method achieved to this order of accuracy using
$N=100$ nodes. Fig. 8 and Fig. 9 show the approximate solutions, whereas Fig.
10 demonstrates the phase-space solution for $\nu=0.5$ and $\tau=5$. Also, log
plot of MAE for $\tau=1$ and $\tau=3$ are plotted in Fig. 11 and Fig. 12,
respectively.

Table 9. Numerical results of Example 4 by the shifted Chebyshev basis with $N=15$ and $\tau=1$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$
0.25 | $0.700921$ | $0.700921$ | $0.670849$ | $0.625076$ | $0.555064$
0.75 | $0.425130$ | $0.425130$ | $0.417343$ | $0.413280$ | $0.419605$
1.25 | $0.198977$ | $0.198979$ | $0.172599$ | $0.135910$ | $0.097493$
1.75 | $0.021138$ | $0.021138$ | $0.024889$ | $0.032106$ | $0.044666$
2.0 | $0.004444$ | $0.004444$ | $0.011581$ | $0.021070$ | $0.036340$
Table 10. Numerical results of Example 4 by the shifted Chebyshev basis with $N=15$ and $\tau=3$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$
0.75 | $0.425130$ | $0.425130$ | $0.417917$ | $0.414396$ | $0.421697$
2.25 | $0.094859$ | $0.094859$ | $0.137716$ | $0.197791$ | $0.288194$
3.75 | $0.002861$ | $0.002859$ | $0.006645$ | $0.014154$ | $0.030316$
5.25 | $0.000000$ | $0.000000$ | $0.000637$ | $0.003305$ | $0.013782$
6.0 | $0.000000$ | $0.000000$ | $0.000429$ | $0.002413$ | $0.011416$
Table 11. Numerical results of Example 4 by the shifted Legendre basis with $N=15$ and $\tau=1$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$
0.25 | $0.700921$ | $0.700923$ | $0.670842$ | $0.625061$ | $0.555051$
0.75 | $0.425130$ | $0.425132$ | $0.417353$ | $0.413275$ | $0.419613$
1.25 | $0.198977$ | $0.198976$ | $0.172538$ | $0.135923$ | $0.097481$
1.75 | $0.021138$ | $0.021134$ | $0.024892$ | $0.032112$ | $0.044652$
2.0 | $0.004444$ | $0.004446$ | $0.011572$ | $0.021083$ | $0.036361$
Table 12. Numerical results of Example 4 by the shifted Legendre basis with $N=15$ and $\tau=3$. $t$ | exact $\nu=1$ | $\nu=1$ | $\nu=0.9$ | $\nu=0.75$ | $\nu=0.5$
---|---|---|---|---|---
0.0 | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$ | $0.900000$
0.75 | $0.425130$ | $0.425130$ | $0.417921$ | $0.414410$ | $0.421730$
2.25 | $0.094859$ | $0.094855$ | $0.137725$ | $0.197820$ | $0.288230$
3.75 | $0.002861$ | $0.002854$ | $0.006653$ | $0.014168$ | $0.030334$
5.25 | $0.000000$ | $0.000001$ | $0.000648$ | $0.003319$ | $0.013797$
6.0 | $0.000000$ | $0.000003$ | $0.000437$ | $0.002443$ | $0.011443$
Figure 8. Approximate solution of Example 4 by the shifted Legendre basis with
$N=15$ and $\tau=1$. Figure 9. Approximate solution of Example 4 by the
shifted Legendre basis with $N=15$ and $\tau=3$. Figure 10. Phase-space
solution of Example 4 by the shifted Legendre basis with $N=15$, $\nu=0.5$ and
$\tau=3$. Figure 11. Log plot of MAE of Example 4 with $\nu=1$ and $\tau=1$.
Figure 12. Log plot of MAE of Example 4 with $\nu=1$ and $\tau=3$.

Example 5. As a final model example, consider the following FDDE which is
introduced by Huang and Williams [20]

(6.9) $D^{\nu}u(t)+\delta D^{\nu_{1}}u(t)=-u(t)+\frac{\mu
q^{2}}{u^{3}(t)}(u(t)-\gamma u(t-\tau)),\quad t\in(0,1),\\\ $

equipped with the conditions

$u(t)=1,\quad t\in[-\tau,0],\quad u(1)=3,$

where $\delta=0.3$, $\mu=1$, $q=0.4$ and $\gamma=0.2$. Computational results
of the shifted Legendre basis with $N=15$ and $\tau=5$ is reported in Table 14
and Table 13 demonstrates the results using the shifted Chebyshev basis. A
comparison between the second and third columns of Table 13 and Table 14 show
that the present work is in remarkable agreement with the function bvp4c of
the Matlab software. However, the MAE of the finite difference method at $t=1$
is of order $10^{-2}$ [20], but in both presented schemes we get the exact
results. The graph of the numerical solutions of (6.9) for different values of
$\nu$ and $\nu_{1}$ is sketched in Fig. 13. Furthermore, Fig. 14 confirms the
agreement between the present work and the function bvp4c of the Matlab.

Table 13. Numerical results of Example 5 by the shifted Chebyshev basis with $N=15$ and $\tau=5$. $t$ | $\nu=2\&\nu_{1}=1$ | $\nu=1.5\&\nu_{1}=0.5$ | $\nu=1.75\&\nu_{1}=0.75$ | $\nu=1.95\&\nu_{1}=0.95$
---|---|---|---|---
| Current work | bvp4c | Current work
0 | $1.000000$ | $1.000000$ | $1.000000$ | $1.000000$ | $1.000000$
0.25 | $1.758262$ | $1.758281$ | $1.830901$ | $1.800878$ | $1.767155$
0.75 | $2.777401$ | $2.777411$ | $2.802636$ | $2.800559$ | $2.783206$
1 | $3.000000$ | $3.000000$ | $3.000000$ | $3.000000$ | $3.000000$
Table 14. Numerical results of Example 5 by the shifted Legendre basis with $N=15$ and $\tau=5$. $t$ | $\nu=2\&\nu_{1}=1$ | $\nu=1.5\&\nu_{1}=0.5$ | $\nu=1.75\&\nu_{1}=0.75$ | $\nu=1.95\&\nu_{1}=0.95$
---|---|---|---|---
| Current work | bvp4c | Current work
0 | $1.000000$ | $1.000000$ | $1.000000$ | $1.000000$ | $1.000000$
0.25 | $1.758275$ | $1.758281$ | $1.830918$ | $1.800884$ | $1.767174$
0.75 | $2.777414$ | $2.777411$ | $2.802647$ | $2.800557$ | $2.783224$
1 | $3.000000$ | $3.000000$ | $3.000000$ | $3.000000$ | $3.000000$
Figure 13. Approximate solution of Example 5 by the shifted Legendre basis
with $N=15$ and $\tau=5$. Figure 14. Comparison between the function bvp4c and
shifted Legendre basis solutions of Example 5 with $N=15$ and $\tau=5$, where
$\nu=2$ and $\nu_{1}=1$.

The CPU time of the above examples is reported in Table 15. As we see from the
table, the CPU time of the shifted Chebyshev basis is less than Legendre one.

Table 15. CPU time (s) for the proposed schemes with $N=15$. Example | Shifted Legendre basis | Shifted Chebyshev basis | $\tau$ | $\nu$ | $\nu_{1}$
---|---|---|---|---|---
1(case i) | $0.765$ | $0.621$ | $0.5$ | $0.1$ | $-$
1(case ii) | $1.852$ | $1.631$ | $0.5$ | $0.1$ | $-$
2(case i) | $0.462$ | $0.283$ | $0.5$ | $0.1$ | $-$
2(case ii) | $0.542$ | $0.321$ | $0.5$ | $0.1$ | $-$
3 | $1.024$ | $0.826$ | $3$ | $0.9$ | $-$
4 | $4.657$ | $3.245$ | $1$ | $0.9$ | $-$
5 | $2.241$ | $2.798$ | $5$ | $1.5$ | $0.5$

##  7\. Conclusion

In this article, a new formula for fractional derivatives of shifted Legendre
polynomials is derived. All the fractional derivatives are considered in the
Caputo sense. By using the formula and the formula based on shifted Chebyshev
polynomials for fractional derivatives [19], numerical schemes for solving
nonlinear FDDEs are proposed. The proposed schemes exploit the method of steps
and shifted Legendre (Chebyshev) basis to generate an approximate solution. A
mathematical analysis shows that the proposed schemes have an exponential rate
of convergence. Moreover, practical examples are taken to demonstrate the
effectiveness of the obtained results. MAE reveals that the approximate
solution has acceptable conformity with the available literature. Further
development of the proposed schemes should be concentrated on solving
nonlinear fractional delay differential problems with more than one delay. It
would also be interesting to extend an approximate solution in which a
discontinuous nonlinear $f$ is considered.

##  Appendix A Proof of relation (31)

Integrating by parts of the left-hand side of (5.7) yields

$\displaystyle\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s\hskip
170.71652pt$
$\displaystyle=\frac{1}{\Gamma(m-\nu)(n-1)!(m-\nu)}\int_{(i-1)\tau}^{s}(t-s)^{n-1}(s-s_{1})^{m-\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})|_{s=(i-1)\tau}^{s=t}\mathrm{d}s_{1}\hskip
113.81102pt$
$\displaystyle+\frac{1}{\Gamma(m-\nu)(n-2)!(m-\nu)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(t-s)^{n-2}(s-s_{1})^{m-\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s.\hskip
113.81102pt$

As we know, one can replace
$\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}F\mathrm{d}s_{1}\mathrm{d}s$ by
$\int_{(i-1)\tau}^{t}\int_{s_{1}}^{t}F\mathrm{d}s\mathrm{d}s_{1}$. Therefore,
the first term of the right hand side of the above relation becomes zero when
we substitute the limits of integration. Thus, we have

$\displaystyle\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s\hskip
170.71652pt$
$\displaystyle=\frac{1}{\Gamma(m-\nu)(n-2)!(m-\nu)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(t-s)^{n-2}(s-s_{1})^{m-\nu}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s.\hskip
113.81102pt$

By repeating the above steps, i.e., after $k$ times integrating by parts we
have

$\displaystyle\frac{1}{\Gamma(m-\nu)}\int_{(i-1)\tau}^{t}\frac{(t-s)^{n-1}}{(n-1)!}\int_{(i-1)\tau}^{s}(s-s_{1})^{m-\nu-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle=M\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(t-s)^{n-k-1}(s-s_{1})^{m-\nu+k-1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})\mathrm{d}s_{1}\mathrm{d}s,$

where

$M=\frac{1}{\Gamma(m-\nu)(n-k-1)!\prod_{j=1}^{k}(m-\nu+j-1)}.$

Now, by putting $k=n-1$ we conclude the relation (5.7).

##  Appendix B Proof of relation (33)

Integrating by parts of the left-hand side of (5.9) yields

$\displaystyle\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle=\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}(s-s_{1})^{m+n-\nu-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m-1)}(s_{1})|_{s_{1}=(i-1)\tau}^{s_{1}=s}$
$\displaystyle\mathrm{d}s$
$\displaystyle+\frac{(m+n-\nu-2)}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-3}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m-1)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s.$

After substitution of the limits of integration, the first term of right-hand
side of the above relation becomes zero, considering the initial conditions of
Eq. (4.9). Therefore, we have

$\displaystyle\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle=\frac{(m+n-\nu-2)}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-3}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m-1)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s.$

By repeating the above steps, i.e., after $m-3$ times integrating by parts we
have

$\displaystyle\frac{1}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{m+n-\nu-2}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s$
$\displaystyle=\frac{\prod_{j=1}^{m-3}(m+n-\nu-1-j)}{\Gamma(m-\nu)\prod_{j=1}^{n-1}(m-\nu+j-1)}\int_{(i-1)\tau}^{t}\int_{(i-1)\tau}^{s}(s-s_{1})^{n-\nu+1}\leftidx{{}_{\scriptscriptstyle\scriptscriptstyle{(i)}}}e_{\scriptscriptstyle\scriptscriptstyle{N}}^{(3)}(s_{1})$
$\displaystyle\mathrm{d}s_{1}\mathrm{d}s.$

The obtained result completes the proof of relation (5.9).

##  Appendix C Boundedness of $D^{\nu}$

By the definition of Caputo fractional derivative, we have

(C.1)
$D^{\nu}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)=\frac{1}{\Gamma(m-\nu)}\int_{0}^{t}(t-s)^{m-\nu-1}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(s)ds.$

As we know, if
$u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}:[0,t]\rightarrow\mathbb{R}$
is continuous and $(t-s)^{m-\nu-1}$ is an integrable function that does not
change sign on $[0,t]$, then there exists $\xi$ in $(0,t)$ such that

$\int_{0}^{t}(t-s)^{m-\nu-1}u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(t)ds=u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(\xi)\int_{0}^{t}(t-s)^{m-\nu-1}ds=\frac{u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(\xi)}{m-\nu}t^{m-\nu}.$

Therefore,

$|D^{\nu}u_{\scriptscriptstyle\scriptscriptstyle{N}}(t)|\leq
M|u_{\scriptscriptstyle\scriptscriptstyle{N}}^{(m)}(\xi)|\leq
M|u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{1:N}_{w}}.$

Hence, for $k\geq 1$ we have

$|D^{\nu}u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}}\leq
C|u_{\scriptscriptstyle\scriptscriptstyle{N}}|_{H^{k:N}_{w}},$

which means $D^{\nu}$ is bounded.

## References

  * [1] A. H. Bhrawy, M. A. Zaky, and J. A. T. Machado, _Numerical solution of the two-sided space time fractional telegraph equation via Chebyshev tau approximation_ , J. Optim. Theory Appl. 174 (2017), 321–341.
  * [2] T. $\breve{\text{S}}$kovránek, I. Podlubny, and I. Petrá$\breve{\text{s}}$, _Modeling of the national economies in state-space: A fractional calculus approach_ , Econ. Model. 29 (2012), no. 4, 1322–1327.
  * [3] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, _Spectral methods: Fundamentals in single domains_ , Springer-Verlag, Berlin Heidelberg, 2006.
  * [4] Y. Q. Chen and K. L. Moore, _Analytical stability bound for delayed second-order systems with repeating poles using Lambert function W_ , Automatica 38 (2002), 891–895.
  * [5] S. Choudhary and V. Daftardar-Gejji, _Existence uniqueness theorems for multi-term fractional delay differential equations_ , Fract. Calc. Appl. Anal. 18 (2015), no. 5, 1113–1127.
  * [6] V. Daftardar-Gejji, Y. Sukale, and S. Bhalekar, _Solving fractional delay differential equations: a new approach_ , Fract. Calc. Appl. Anal. 18 (2015), no. 2, 400–418.
  * [7] M. Dehghan and R. Salehi, _Solution of a nonlinear time-delay model in biology via semi-analytical approaches_ , Comput. Phys. Commun. 181 (2010), 1255–1265.
  * [8] T. Erneux, _Applied delay differential equations_ , Springer-Verlag, New York, 2009.
  * [9] F. Ghoreishi and S. Yazdani, _An extension of the spectral Tau method for numerical solution of multi-order fractional differential equations with convergence analysis_, Comput. Math. Appl. 61 (2011), 30–43.
  * [10] A. Gogatishvill and J. Lang, _The generalized hardy operator with kernel and variable integral limits in Banach function spaces_ , J. Inequal. Appl. 4 (1999), no. 1, 1–16.
  * [11] R. P. Kanwal, _Linear integral equations_ , Birkhauser Boston, Boston, 1997.
  * [12] M. M. Khader, _The use of generalized Laguerre polynomials in spectral methods for solving fractional delay differential equations_ , J. Comput. Nonlinear Dyn. 8 (2013), no. 4, 041018.
  * [13] H. Khosravian-Arab, M. Dehghan, and M. R. Eslahchi, _Fractional spectral and pseudo-spectral methods in unbounded domains: Theory and applications_ , J. Comput. Phys. 338 (2017), 527–566.
  * [14] R. L. Magin, _Fractional calculus models of complex dynamics in biological tissues_ , Comput. Math. Appl. 59 (2010), no. 5, 1586–1593.
  * [15] B. P. Moghaddam and Z. S. Mostaghim, _A numerical method based on finite difference for solving fractional delay differential equations_ , J. Taibah Univ. Sci. 7 (2013), 120–127.
  * [16] by same author, _Novel matrix approach to fractional finite difference for solving models based on nonlinear fractional delay differential equations_ , Ain Shams Engineering Journal 5 (2014), no. 2, 585–594.
  * [17] O. H. Mohammed and A. I. Khlaif, _Adomian decomposition method for solving delay differential equations of fractional order_ , IOSR Journal of Mathematics 10 (2014), 1–5.
  * [18] C. A. Monje, Y. Chen, B. M. Vinagre, D. Xue, and V. Feliu, _Fractional-order systems and controls: Fundamentals and applications_ , 1 ed., Advances in Industrial Control, Springer-Verlag London, 2010.
  * [19] M. Mousa-Abadian and S. H. Momeni-Masuleh, _Steps Chebyshev-Tau scheme for solving linear fractional differential equations with time delay_ , Submitted for publication.
  * [20] B. Parsa Moghaddam and Z. Salamat Mostaghim, _Modified finite difference method for solving fractional delay differential equations_ , Bol. Soc. Paran. Mat. 35 (2017), 49–58.
  * [21] B. Parsa Moghaddam, Sh. Yaghoobi, and J. A. T. Machado, _An extended predictor-corrector algorithm for variable-order fractional delay differential equations_ , J. Comput. Nonlinear Dyn. 11 (2016), 061001\.
  * [22] I. Podlubny, _Fractional differential equations_ , Academic Press, San Diego, 1999.
  * [23] J. Prakash, M. Kothandapani, and V. Bharathi, _Numerical approximations of nonlinear fractional differential difference equations by using modified He-Laplace method_ , Alexandria Engineering Journal 55 (2016), no. 1, 645–651.
  * [24] A. Saadatmandi and M. Dehghan, _Variational iteration method for solving a generalized pantograph equation_ , Comput. Math. Appl. 58 (2009), 2190–2196.
  * [25] U. Saeed, M. ur Rehman, and M. A. Iqbal, _Modified Chebyshev wavelet methods for fractional delay-type equations_ , Appl. Math. Comput. 264 (2015), 431–442.
  * [26] S. Sedaghat, Y. Ordokhani, and M. Dehghan, _Numerical solution of the delay differential equations of pantograph type via Chebyshev polynomials_ , Commun. Nonlinear Sci. Numer. Simulat. 17 (2012), 4815–4830.
  * [27] F. Shakeri and M. Dehghan, _Solution of delay differential equations via a homotopy perturbation method_ , Math. Comput. Model. 48 (2008), 486–498.
  * [28] A. Si-Ammour, S. Djennoune, and M. Bettayeb, _A sliding mode control for linear fractional systems with input and state delays_ , Commun. Nonlinear Sci. Numer. Simulat. 14 (2009), no. 5, 2310–2318.
  * [29] Z. Wang, _A numerical method for delayed fractional-order differential equations_ , Journal of Applied Mathematics 2013 (2013), 1–7.
  * [30] Sh. Yaghoobi, B. Parsa Moghaddam, and K. Ivaz, _An efficient cubic spline approximation for variable-order fractional differential equations with time delay_ , Nonlinear Dyn. 87 (2017), 815–826.
  * [31] Z. Yang and J. Cao, _Initial value problems for arbitrary order fractional differential equations with delay_ , Commun. Nonlinear Sci. Numer. Simulat. 18 (2013), no. 11, 2993–3005.
  * [32] M. Zayernouri, W. Cao, Z. Zhang, and G. E. Karniadakis, _Spectral and discontinuous spectral element methods for fractional delay equations_ , SIAM J. Sci. Comput. 36 (2014), B904–B929.