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

TL;DR

This paper introduces a new numerical approach combining the method of steps and shifted Legendre collocation for solving nonlinear fractional delay differential equations, with proven convergence and demonstrated accuracy.

Contribution

It presents a novel formula for fractional derivatives of shifted Legendre polynomials and develops efficient schemes for nonlinear fractional delay equations.

Findings

Schemes effectively solve nonlinear fractional delay differential equations.

Convergence analysis confirms the reliability of the methods.

Numerical examples demonstrate high accuracy and efficiency.

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 nonlinear 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.

Tables

Table 1. Table 1. ℒ 2 superscript ℒ 2 \mathcal{L}^{2} -Error of Example 1 (Case ( i ) 𝑖 (i) ) for different values of N 𝑁 N with τ = 0.5 𝜏 0.5 \tau=0.5 and ν = 0.1 𝜈 0.1 \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. Table 2. ℒ 2 superscript ℒ 2 \mathcal{L}^{2} -Error of Example 1 (Case ( i ​ i ) 𝑖 𝑖 (ii) ) for different values of N 𝑁 N with τ = 0.5 𝜏 0.5 \tau=0.5 and ν = 0.1 𝜈 0.1 \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$

Table 3. Table 3. ℒ 2 superscript ℒ 2 \mathcal{L}^{2} -Error of Example 2 (Case ( i ) 𝑖 (i) ) for different values of N 𝑁 N with τ = 0.5 𝜏 0.5 \tau=0.5 and ν = 0.1 𝜈 0.1 \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. Table 4. ℒ 2 superscript ℒ 2 \mathcal{L}^{2} -Error of Example 2 (Case ( i ​ i ) 𝑖 𝑖 (ii) ) for different values of N 𝑁 N with τ = 0.5 𝜏 0.5 \tau=0.5 and ν = 0.1 𝜈 0.1 \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$

Table 5. Table 5. Numerical results of Example 3 by the shifted Chebyshev basis with N = 15 𝑁 15 N=15 and τ = 3 𝜏 3 \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. Table 6. Numerical results of Example 3 by the shifted Chebyshev basis with N = 15 𝑁 15 N=15 and τ = 5 𝜏 5 \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. Table 7. Numerical results of Example 3 by the shifted Legendre basis with N = 15 𝑁 15 N=15 and τ = 3 𝜏 3 \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. Table 8. Numerical results of Example 3 by the shifted Legendre basis with N = 15 𝑁 15 N=15 and τ = 5 𝜏 5 \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$

Table 9. Table 9. Numerical results of Example 4 by the shifted Chebyshev basis with N = 15 𝑁 15 N=15 and τ = 1 𝜏 1 \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. Table 10. Numerical results of Example 4 by the shifted Chebyshev basis with N = 15 𝑁 15 N=15 and τ = 3 𝜏 3 \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. Table 11. Numerical results of Example 4 by the shifted Legendre basis with N = 15 𝑁 15 N=15 and τ = 1 𝜏 1 \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. Table 12. Numerical results of Example 4 by the shifted Legendre basis with N = 15 𝑁 15 N=15 and τ = 3 𝜏 3 \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$

Table 13. Table 13. Numerical results of Example 5 by the shifted Chebyshev basis with N = 15 𝑁 15 N=15 and τ = 5 𝜏 5 \tau=5 .

$t$
	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. Table 14. Numerical results of Example 5 by the shifted Legendre basis with N = 15 𝑁 15 N=15 and τ = 5 𝜏 5 \tau=5 .

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

Table 15. Table 15. CPU time (s) for the proposed schemes with N = 15 𝑁 15 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$

Equations

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

$$L_{k+1}(t)=\frac{2k+1}{k+1}t{L}_k(t)-\frac{k}{k+1}{L}_{k-1}(t),$$

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

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

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

$$L_{\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,$$

$$D^m{L}_{\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}.$$

$${\int }_{\alpha }^{\beta }{L}_{\alpha ,\beta ,i}(x)L_{\alpha ,\beta ,k}(x)\mathrm{d}x={\delta }_{ik}(k+\frac{1}{2}{)}^{-1},$$

$$u(x)=\sum _{k=0}^{\infty }{\hat{u}}_k{L}_{\alpha ,\beta ,k}(x),$$

$${\hat{u}}_k=(k+\frac{1}{2}){\int }_{\alpha }^{\beta }u(x)L_{\alpha ,\beta ,k}(x)\mathrm{d}x,k\ge 0.$$

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

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

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

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

$${\int }_{\alpha }^{\beta }f(x)\mathrm{d}x.$$

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

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

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

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

$${\int }_{\alpha }^{\beta }\frac{p(x)}{\sqrt{1-x^2}}\mathrm{d}x$$

$$\text{mbox}and{\int }_{\alpha }^{\beta }p(x)\mathrm{d}x$$

$$D^{\nu }{L}_{\alpha ,\beta ,n}(x)=0,n=0,1,\cdots ,\lceil \nu \rceil -1.$$

$$D^{\nu }{L}_{\alpha ,\beta ,n}(x)=\sum _{i=0}^{\infty }{S}_{\nu }(n,i)L_{\alpha ,\beta ,i}(x),n=\lceil \nu \rceil ,\lceil \nu \rceil +1,\cdots ,$$

$$S_{\nu }(n,i)=\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},$$

$$c_{ij}=(i+\frac{1}{2}){\int }_{\alpha }^{\beta }{x}^{j-\nu }{L}_{\alpha ,\beta ,i}(x)\mathrm{d}x.$$

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

$$D^{\nu }{L}_{\alpha ,\beta ,n}(x)=$$

$$=$$

$$x^{j-\nu }=\sum _{i=0}^{\infty }{c}_{ij}{L}_{\alpha ,\beta ,j}(x),$$

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

$$|f(t,y_1,u)-f(t,y_2,u)|$$

Full text

∗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]

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

By using the identity

[TABLE]

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

[TABLE]

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

[TABLE]

Proof.

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

The shifted Legendre polynomials satisfy the following relation

[TABLE]

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

[TABLE]

where

[TABLE]

The associated inner product and norm are given by

[TABLE]

and

[TABLE]

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

[TABLE]

Furthermore, the associated semi-norm is defined as follows

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

• •

Legendre Gauss-Lobatto

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

Theorem 3.2**.**

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

[TABLE]

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

[TABLE]

in which

[TABLE]

Proof.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with the initial conditions

[TABLE]

where

[TABLE]

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

[TABLE]

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]

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Subtracting (5.5) from (5.4) leads to

[TABLE]

where

[TABLE]

and

[TABLE]

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)

[TABLE]

and

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

We may rewrite Eq. (5.13) as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

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

[TABLE]

So that

[TABLE]

From (5.19) and (5.21) we have

[TABLE]

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

[TABLE]

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

[TABLE]

with the boundary condition

[TABLE]

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)$:

[TABLE]

which corresponding exact solution is

[TABLE]

Case $(ii)$:

[TABLE]

where

[TABLE]

The related exact solution is

[TABLE]

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]

[TABLE]

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.

Example 2. Consider the following FDDE [32]

[TABLE]

with the initial condition

[TABLE]

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.

Example 3. Consider Houseflies model as following [15]

[TABLE]

with the initial condition

[TABLE]

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.

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

[TABLE]

with the initial condition

[TABLE]

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.

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

[TABLE]

equipped with the conditions

[TABLE]

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.

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Therefore,

[TABLE]

Hence, for $k\geq 1$ we have

[TABLE]

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