Legendre wavelet collocation method combined with the Gauss--Jacobi quadrature for solving fractional delay-type integro-differential equations

TL;DR

This paper introduces a novel collocation method combining Legendre wavelets and Gauss--Jacobi quadrature to efficiently solve fractional delay-type integro-differential equations with high accuracy.

Contribution

The work develops a new numerical approach that integrates Legendre wavelets with Gauss--Jacobi quadrature for solving fractional integro-differential equations, including error analysis.

Findings

Demonstrates high accuracy in example problems

Provides error bounds for Legendre wavelet approximation

Shows efficiency of the method through numerical results

Abstract

In this work, we present a collocation method based on the Legendre wavelet combined with the Gauss--Jacobi quadrature formula for solving a class of fractional delay-type integro-differential equations. The problem is considered with either initial or boundary conditions and the fractional derivative is described in the Caputo sense. First, an approximation of the unknown solution is considered in terms of the Legendre wavelet basis functions. Then, we substitute this approximation and its derivatives into the considered equation. The Caputo derivative of the unknown function is approximated using the Gauss--Jacobi quadrature formula. By collocating the obtained residual at the well-known shifted Chebyshev points, we get a system of nonlinear algebraic equations. In order to obtain a continuous solution, some conditions are added to the resulting system. Some error bounds are given for the Legendre wavelet approximation of an arbitrary function. Finally, some examples are included to show the efficiency and accuracy of this new technique.

Tables

Table 1. Table 1: (Example 5.2 .) Numerical results with different values of M 𝑀 M and k 𝑘 k .

	$M=3$			$M=4$			$M=5$
$k$	$L^2$-Error	Ratio		$L^2$-Error	Ratio		$L^2$-Error	Ratio
$1$	$3.23e-01$	—		$5.53e-02$	—		$7.07e-03$	—
$2$	$1.99e-01$	$1.62$		$4.73e-02$	$1.17$		$8.06e-03$	$0.88$
$3$	$1.06e-01$	$1.88$		$2.08e-02$	$2.27$		$1.91e-03$	$4.22$
$4$	$5.23e-02$	$2.03$		$7.04e-03$	$2.95$		$4.01e-04$	$4.76$
$5$	$2.50e-02$	$2.09$		$2.04e-03$	$3.45$		$7.60e-05$	$5.28$
$6$	$1.18e-02$	$2.19$		$5.40e-04$	$3.78$		$1.36e-05$	$5.59$

Table 2. Table 2: (Example 5.3 .) Numerical results with M = 3 𝑀 3 M=3 , α = 2 𝛼 2 \alpha=2 and different values of k 𝑘 k .

$k$	$L^2$-Error	Ratio
$1$	$1.76e-01$	—
$2$	$6.89e-02$	$2.55$
$3$	$3.43e-02$	$2.01$
$4$	$1.68e-02$	$2.04$
$5$	$8.34e-03$	$2.01$
$6$	$4.15e-03$	$2.01$

Table 3. Table 3: (Example 5.4 .) Absolute error at some selected points with α = 1 𝛼 1 \alpha=1 and k = 1 𝑘 1 k=1 .

	Present method				Method of [38]
$x$	$M=10$	$M=20$			$M=10$	$M=20$
$0.0$	$0.00e+00$	$0.00e+00$			$0.00e+00$	$0.00e+00$
$0.6$	$1.21e-08$	$3.33e-21$			$1.46e-05$	$2.13e-13$
$1.2$	$2.45e-08$	$2.70e-21$			$2.25e-05$	$3.73e-14$
$1.8$	$1.79e-07$	$4.93e-21$			$4.21e-05$	$2.93e-13$
$2.4$	$2.18e-07$	$3.71e-20$			$7.66e-05$	$4.23e-13$
$3.0$	$4.35e-05$	$1.67e-17$			$1.45e-04$	$8.17e-13$

Table 4. Table 4: (Example 5.4 .) Numerical results with α = 1 𝛼 1 \alpha=1 and different values of M 𝑀 M and k 𝑘 k .

	$M=3$			$M=4$			$M=5$
$k$	$L^2$-Error	Ratio		$L^2$-Error	Ratio		$L^2$-Error	Ratio
$1$	$5.93e+00$	—		$1.85e+00$	—		$4.43e-01$	—
$2$	$4.63e+00$	$1.28$		$1.14e+00$	$1.62$		$1.88e-01$	$2.36$
$3$	$3.18e+00$	$1.46$		$4.86e-01$	$2.35$		$4.23e-02$	$4.44$
$4$	$1.94e+00$	$1.64$		$1.48e-01$	$3.28$		$6.85e-03$	$6.17$
$5$	$1.10e+00$	$1.76$		$4.28e-02$	$3.46$		$9.76e-04$	$7.02$
$6$	$5.89e-01$	$1.88$		$1.11e-02$	$3.86$		$1.29e-04$	$7.57$

Equations

$${}^c{D}^{\alpha }y(x)=f\left(x,y(x),y^{\prime }(x),y(px-\tau ),y^{\prime }(px-\tau ),{\int }_{px-\tau }^{x}g(x,s,y(s))ds\right),$$

$$y(x)=y_0(x),x\le 0,y^{\prime }(0)=y_0^{\prime },$$

$$y(x)=y_0(x),x\le 0,y(l)=y_1,$$

$$\rho =\underset{x\ge 0}{\min }(px-\tau )=-\tau ;\text{ moreover y0′(0)=y0′}.$$

$$I^{\alpha}y(x)=\left\{\begin{array}[]{ll}\frac{1}{\Gamma(\alpha)}\int_{0}^{x}(x-s)^{\alpha-1}y(s)ds,&\alpha>0,\\ y(x),&\alpha=0,\end{array}\right.$$

$$\Gamma (\alpha )={\int }_0^{\infty }{x}^{\alpha -1}{e}^{-x}dx.$$

$${}^c{D}^{\alpha }y(x)=\frac{1}{\Gamma (n-\alpha )}{\int }_0^x(x-s{)}^{n-\alpha -1}{y}^{(n)}(s)ds,n-1<\alpha \le n,n\in \mathbb{N},$$

$$w^{\lambda ,\nu }(x)=(1-x{)}^{\lambda }(1+x{)}^{\nu }.$$

$${\int }_{-1}^1{w}^{\lambda ,\nu }(x)J_n^{(\lambda ,\nu )}(x)J_m^{(\lambda ,\nu )}(x)dx=h_n^{\lambda ,\nu }{\delta }_{mn},$$

$$h_n^{\lambda ,\nu }=\frac{2^{\lambda +\nu +1}\Gamma (\lambda +n+1)\Gamma (\nu +n+1)}{n!(\lambda +\nu +2n+1)\Gamma (\lambda +\nu +n+1)}.$$

$${\int }_{-1}^1(1-x{)}^{\lambda }(1+x{)}^{\nu }f(x)dx=\sum _{s=1}^N{w}_{s}f(x_s)+R_N(f),$$

$$w_s=\frac{2^{\lambda +\nu +1}\Gamma (\lambda +N+1)\Gamma (\nu +N+1)}{N!\Gamma (\lambda +\nu +N+1)(\frac{d}{dx}{J}_N^{(\lambda ,\nu )}(x_s){)}^2(1-x_s^2)},$$

$$\begin{array}{rl}{R}_n(f)=& \frac{2^{\lambda +\nu +2N+1}N!\Gamma (\lambda +N+1)\Gamma (\nu +N+1)\Gamma (\lambda +\nu +N+1)}{(\lambda +\nu +2N+1)(\Gamma (\lambda +\nu +2N+1){)}^2}\\ & \times \frac{f^{(2N)}(\xi )}{(2N)!},\xi \in (-1,1).\end{array}$$

$$\psi_{n,m}(x)=\left\{\begin{array}[]{ll}2^{\frac{k}{2}}\sqrt{\frac{2m+1}{2l}}P_{m}\left(\frac{2^{k}}{l}x-2n+1\right),&\frac{n-1}{2^{k-1}}l\leq x<\frac{n}{2^{k-1}}l,\\ 0,&\text{otherwise},\end{array}\right.$$

$$P_m(x)=2P_{m-1}(x)-P_{m-2}(x),m=2,3,\dots ,$$

$${\int }_{-1}^1{P}_i(x)P_j(x)dx=\frac{2}{2i+1}{\delta }_{ij}.$$

$$\psi_{n,m}^{{}^{\prime}}(x)=\left\{\begin{array}[]{ll}\frac{2^{\frac{3k}{2}}}{l}\sqrt{\frac{2m+1}{2l}}P_{m}^{{}^{\prime}}\left(\frac{2^{k}}{l}x-2n+1\right),&\frac{n-1}{2^{k-1}}l\leq x<\frac{n}{2^{k-1}}l,\\ 0,&\text{otherwise},\end{array}\right.$$

$$\psi_{n,m}^{{}^{\prime\prime}}(x)=\left\{\begin{array}[]{ll}\frac{2^{\frac{5k}{2}}}{l^{2}}\sqrt{\frac{2m+1}{2l}}P_{m}^{{}^{\prime\prime}}\left(\frac{2^{k}}{l}x-2n+1\right),&\frac{n-1}{2^{k-1}}l\leq x<\frac{n}{2^{k-1}}l,\\ 0,&\text{otherwise}.\end{array}\right.$$

$$f(x)\simeq {\Psi }_{k,M-1}(f)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(x)=A^T\psi (x),$$

$$a_{n,m}=⟨f(x),{\psi }_{n,m}(x)⟩={\int }_0^{l}f(x){\psi }_{n,m}(x)dx,$$

$$A=[a_{1,0},a_{1,1},\dots ,a_{1,M-1},\dots ,a_{2^{k-1},0},a_{2^{k-1},1},\dots ,a_{2^{k-1},M-1}{]}^T,$$

$$\psi (x)=[{\psi }_{1,0}(x),{\psi }_{1,1}(x),\dots ,{\psi }_{1,M-1}(x),\dots ,{\psi }_{2^{k-1},0}(x),{\psi }_{2^{k-1},1}(x),\dots ,{\psi }_{2^{k-1},M-1}(x){]}^T.$$

$$y(x)\simeq \sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(x),$$

$$y(px-\tau )\simeq \sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(px-\tau ).$$

$$\begin{array}{r}{R}_1(x)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{}^c{D}^{\alpha }{\psi }_{n,m}(x)-f(x,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(x),\\ \sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}^{\prime }(x),\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(px-\tau ),\\ \sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}^{\prime }(px-\tau ),{\int }_{px-\tau }^{x}g(x,s,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(s))ds).\end{array}$$

$$\begin{array}{r}{R}_2(x)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{}^c{D}^{\alpha }{\psi }_{n,m}(x)-f(x,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(x),\\ \sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}^{\prime }(x),y_0(px-\tau ),y_0^{\prime }(px-\tau ),\\ {\int }_{px-\tau }^{0}g(x,s,y_0(s))ds+{\int }_0^{x}g(x,s,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{a}_{n,m}{\psi }_{n,m}(s))ds)\end{array}$$

$${}^c{D}^{\alpha }{\psi }_{n,m}(x)=\frac{1}{\Gamma (2-\alpha )}{\int }_0^x(x-s{)}^{1-\alpha }{\psi }_{n,m}^{{}^{\prime \prime }}(s)ds.$$

$$t=2\left(\frac{s}{x}\right)-1,dt=\frac{2}{x}ds.$$

$${}^c{D}^{\alpha }{\psi }_{n,m}(x)=\frac{1}{\Gamma (2-\alpha )}{\left(\frac{x}{2}\right)}^{2-\alpha }{\int }_{-1}^1(1-t{)}^{1-\alpha }{\psi }_{n,m}^{{}^{\prime \prime }}(\frac{x}{2}(t+1))dt.$$

$${}^c{D}^{\alpha }{\psi }_{n,m}(x)=\frac{1}{\Gamma (2-\alpha )}{\left(\frac{x}{2}\right)}^{2-\alpha }\sum _{i=1}^N{w}_i{\psi }_{n,m}^{{}^{\prime \prime }}(\frac{x}{2}(t_i+1)),$$

Full text

1 Introduction

In recent decades, fractional calculus, which includes theories of derivatives and integrals of any arbitrary non-integer order, has been used for describing many real world phenomena, in such different domains as wave mechanics, viscoelasticity, physiology or epidemiology [1, 2, 3, 4, 5, 6]. Nevertheless, obtaining analytic solutions for fractional differential and integro-differential equations is very difficult. Thus, many researchers have introduced and developed numerical methods in order to get numerical solutions for these two classes of equations (see for example [7, 8, 9, 10, 11, 12, 13]). Moreover, partial differential equations of fractional order have also been investigated and various numerical methods for their approximation have been introduced recently [14, 15].

Delay differential and integro-differential equations are used in the modeling of a large number of problems in different fields of sciences such as ecology and epidemiology [16], immunology [17], physiology [18] and electrodynamics [19]. Fractional delay differential equations have been used to model certain processes and systems with memory and heritage properties. Discussions about the existence of the solution of such equations can be found in [20, 21, 22]. Because of the computational complexity of fractional derivatives in the delay case, the exact analytical solution of the fractional delay differential equations is hardly available. Therefore, in the last decades many researchers have been attracted to deal with the numerical solution of this class of problems (see for example [23, 24, 25, 26, 27, 28]).

Wavelets are a set of functions built by dilation and translation of a single function $\varphi(t)$ which is called the mother wavelet. These functions are known as a powerful tool to get the numerical solution of equations. The primary idea of wavelets (translation and dilation) returns to the early 1960’s [29, 30]. Some developments exist concerning the multiresolution analysis algorithm based on wavelets [31] and the construction of compactly supported orthonormal wavelet bases [32]. Wavelets form an unconditional (Riesz) basis for $L^{2}(\mathbb{R})$ which is the space of all square integrable functions on the real line, i.e. any function in $L^{2}(\mathbb{R})$ can be decomposed and reconstructed in terms of wavelets [33]. Some researchers have constructed and used a number of wavelets for solving fractional differential equations such as B-spline wavelets [34], Haar wavelets [35], Chebyshev wavelets [36], Legendre wavelets [37] and Bernoulli wavelets [23]. The advantage of wavelets in comparison with other basis functions is that when using them we can improve the accuracy of the approximation in two different ways: 1) by increasing the degree of the mother function (assuming that it is a polynomial); 2) by increasing the level of resolution, i.e. reducing the support of each basis function. Unfortunately no all the researchers who use wavelets for approximating fractional equations take this advantage; for example, in [38] the wavelets are used only with resolution level $k=1$, and only the degree of the polynomial is changed.

In this paper, we consider the following fractional order nonlinear delay-type integro-differential equation

[TABLE]

subject to the initial conditions

[TABLE]

or boundary conditions

[TABLE]

where $0\leq x\leq l$, $1<\alpha\leq 2$, ${}^{c}D^{\alpha}$ denotes the Caputo derivative, $f$ and $g$ are continuous linear or nonlinear functions and $y^{\prime}_{0}$, $y_{1}$ and $p$ are constants. The parameter $\tau$ is called the delay component and is always non-negative. The delay component may be just constant (the constant delay case) or a function of $x$, $\tau=\tau(x)$ (the variable dependent delay case). The term $px-\tau$ is called delayed argument, $y(x)$ is the unknown function to be determined and $y(px-\tau)$ is the value of this function corresponding to the delayed argument. The function $y_{0}(x)$ is defined and continuously differentiable in $[\rho,0]$, where

[TABLE]

Delay problems of this type include as particular case the pantograph equation, when $p\neq 1$ and $\tau=0$. The existence and uniqueness of solution for this kind of problems has been investigated by many authors. In the non-delay case (when there is not delayed argument) we have just a nonlinear fractional equation; the existence and uniqueness of solution for such problems has been discussed for example in [39]; among other restrictions the function $f$ is supposed to satisfy the Lipschitz condition. In the case of delayed argument, authors usually apply the method of steps to reduce the problem to a non-delay one; then existence and uniqueness is analyzed for the new problem. This was done, for example, in [27], for the case of linear fractional equations with delayed arguments.

We present a collocation spectral method based on the Legendre wavelet basis functions combined with the Gauss–Jacobi quadrature formula to solve the above problem. We underline that the application of wavelets as basis functions has special advantages when problems of this class are considered. Actually in many cases the solution is not regular, in the sense that its first or higher derivatives may have discontinuities. In such cases, increasing the degree of polynomial approximation is not a solution to improve the accuracy of the approximation; however such improvement can be obtained by increasing the level of resolution. This is why in the numerical examples in Section 5 we vary both the polynomial degree and the resolution level of wavelets, and we check the effect on the quality of the approximation. To this aim, the rest of the paper is organized as follows: in Section 2, we give some preliminaries of fractional calculus and introduce the Gauss–Jacobi quadrature formula and the Legendre wavelet basis functions. Section 3 is devoted to suggesting a numerical method based on the Legendre wavelet and using the Gauss–Jacobi quadrature for solving the fractional order delay-type integro-differential equations. In Section 4, some error bounds are given. Section 5 includes some numerical examples. Finally, conclusions are given in Section 6.

2 Basic concepts

In this section, we present some preliminaries which will be used further in the next section.

2.1 Preliminaries of fractional calculus

Here, we briefly give some preliminaries and notations of fractional calculus. Two of the most important definitions for fractional integral and derivative operators are Riemann-Liouville integral and Caputo derivative. The Riemann-Liouville fractional integral operator $I^{\alpha}$ of order $\alpha\geq 0$ is defined as follows [40]:

[TABLE]

where $\Gamma(\alpha)$ is the gamma function defined by

[TABLE]

Also, the Caputo fractional derivative operator ${}^{c}D^{\alpha}$ of order $\alpha$ is defined as follows [40]:

[TABLE]

where $n=\lceil\alpha\rceil$ is the smallest integer greater than or equal to $\alpha$.

2.2 Gauss–Jacobi quadrature formula

The Jacobi polynomials $\{J_{n}^{(\lambda,\nu)}\}_{n=0}^{\infty}$, $\lambda,\nu>-1$, $x\in[-1,1]$, built a set of functions which are orthogonal with respect to the weight function

[TABLE]

The orthogonality property is satisfied for these polynomials as follows:

[TABLE]

where $\delta_{mn}$ is the Kronecker delta and

[TABLE]

The Gauss–Jacobi quadrature formula is given by

[TABLE]

where $x_{s}$, $s=1,\ldots,N$, are the roots of $J_{N}^{(\lambda,\nu)}$. Moreover, $w_{s}$, $s=1,\ldots,N$, are the weights and $R_{N}(f)$ is the remainder term which are, respectively, given by (see [41]):

[TABLE]

[TABLE]

Taking into account the remainder term, formula (5) is exact for all polynomials of degree less than or equal to $2N-1$ and is valid if $f$ has no singularity in $(-1,1)$. We use this formula further in our numerical method.

2.3 Legendre wavelets

We define Legendre wavelet functions on the interval $[0,l)$ as

[TABLE]

where $k=1,2,3,\ldots$, is the level of resolution, $n=1,2,3,\ldots,2^{k-1}$, $x$ is the normalized time, and $m=0,1,2,\cdots$, is the degree of the Legendre polynomial which is defined on the interval $[-1,1]$ and is given by the aid of the following recursive formula:

[TABLE]

with $P_{0}(x)=1$ and $P_{1}(x)=x$. Furthermore, in the special case of Jacobi polynomials with $\lambda=\nu=0$, the Legendre polynomials are given. These polynomials are orthogonal with respect to the weight function $w(x)=1$, that is

[TABLE]

The first and second derivatives of the Legendre wavelet function $\psi_{n,m}(x)$ are given, respectively, by

[TABLE]

[TABLE]

An arbitrary function $f(x)\in L^{2}[0,l)$ may be approximated using Legendre wavelets as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

3 Numerical method

In order to introduce the numerical method for solving equation (1), we assume that its solution can be approximated by Legendre wavelets up to $k$-th level of resolution:

[TABLE]

where $M\geq 3$. Then, for all $x$ such that $px-\tau\geq 0$, we can write

[TABLE]

Therefore, if $px-\tau>0$, by substituting equations (11) and (12) in (1), the following residual is obtained:

[TABLE]

On the other hand, if $px-\tau\leq 0$, we define the residual as follows:

[TABLE]

The functions $\psi_{n,m}^{\prime}(\cdot)$ are obtained using (8) and are substituted in $R_{1}$ and $R_{2}$. In a particular case, if $\alpha=2$, we have ${{}^{c}D^{2}}\psi_{n,m}(x)=\psi_{n,m}^{{}^{\prime\prime}}(x)$ which are given by (9). Suppose that $1<\alpha<2$, then using (4), we have

[TABLE]

In order to compute the integral part of (13), first the interval $[0,x]$ is transformed into $[-1,1]$ using the following variable transformation

[TABLE]

Therefore, (13) can be rewritten as

[TABLE]

Since $1-\alpha>-1$, we introduce $\lambda=1-\alpha$ and $\nu=0$ as a special case of Jacobi polynomials and use the Gauss–Jacobi quadrature to compute the integral part of (14). We set $N:=\lceil\frac{M-2}{2}\rceil$, and get

[TABLE]

where $t_{i}$, $i=1,\ldots,N$, are the zeros of $J_{N}^{(1-\alpha,0)}$ and the weights $w_{i}$, $i=1,\ldots,N$, are given using (6), as follows:

[TABLE]

For computing the nodes and weights of the fractional Gauss–Jacobi quadrature numerical algorithms can be used such as Golub-Welsch algorithm introduced by Pang et al. [42].

Remark 3.1**.**

By setting $N:=\lceil\frac{M-2}{2}\rceil$, according to the error of the Gauss–Jacobi quadrature formula, (15) is exact for computing ${{}^{c}D^{\alpha}}\psi_{n,m}(x)$, $n=1$, $m=0,\ldots,M-1$.

Remark 3.2**.**

In the case of complicated functions $g(\cdot,\cdot,\cdot)$ in $R_{1}(\cdot)$ and $R_{2}(\cdot)$, first we transform the interval of integration into the interval $[-1,1]$ and then use the Gauss–Legendre quadrature formula.

We use the following shifted Chebyshev points as the collocation points

[TABLE]

where

[TABLE]

Then, set the residual, $R_{l}(x)$, $l=1,2$, equal to zero at these points. So that, for $i=1,2,\ldots,2^{k-1}$ and $j=2,3,\ldots,M-1$, if $px_{i,j}-\tau>0$, we set

[TABLE]

otherwise, if $px_{i,j}-\tau\leq 0$, we set

[TABLE]

which (16) and (17) yield a system of $2^{k-1}(M-2)$ nonlinear equations in terms of the $2^{k-1}M$ unknown parameters $a_{n,m}$. To find all the values $a_{n,m}$, it is necessary to complete the system in (16) with additional conditions in order to obtain a continuous solution. To do this, let us set

[TABLE]

[TABLE]

Now, two more equations are needed. These equations can be obtained from the initial conditions (2) or boundary conditions (3), respectively, as

[TABLE]

[TABLE]

Finally, by considering (16)–(19)–(20) or (16)–(19), (21) in the case of initial or boundary value problem, respectively, we have a system of $2^{k-1}M$ nonlinear equations which can be solved for computing the unknown parameters $a_{n,m}$.

4 Error bounds

Suppose $\mathbb{P}_{M-1}(J_{k,n})$ denotes the set of all functions whose restriction on each mesh interval $J_{k,n}=(\frac{n-1}{2^{k-1}}l,\frac{n}{2^{k-1}}l)$, $n=1,2,...,{2^{k-1}}$, are polynomials of degree at most $M-1$. In this section, we give some estimates for the Legendre wavelets truncated error $f(x)-\sum\limits_{n=1}^{{2^{k-1}}}{\sum\limits_{m=0}^{M-1}{{a_{n,m}}{\psi_{n,m}}(x)}}$ in terms of the Sobolev norms where $\psi_{n,m}(x)$ and $a_{n,m}$ are given by (7) and (10), respectively. The Sobolev norm of integer order $m\geq 0$ in the interval (a,b), is given by

[TABLE]

where $f^{(j)}$ denotes the derivative of $f$ of order $j$. Here, we recall the following results which have been given in [43].

Lemma 4.1**.**

[43]** Assume that $f(x)\in H^{m}(-1,1),$ ($H^{m}(-1,1)$ is a Sobolev space) with $m\geq 1$ and $L_{M-1}(f)\in\mathbb{P}_{M-1}$ denotes the truncated Legendre series of $f$. Then the following estimate holds:

[TABLE]

and for $1\leq r\leq m$, we have

[TABLE]

where

[TABLE]

and $C$ is a positive constant independent of the function $f$ and integer $M$.

Lemma 4.2**.**

Let $f_{n,k}:J_{k,n}\rightarrow\mathbb{R}$, $n=1,2,...,2^{k-1}$, be functions in $H^{m}(J_{k,n})$. Consider the function $\bar{f}_{n,k}:(-1,1)\rightarrow\mathbb{R}$ defined by $(\bar{f}_{n,k})(x)=f_{n,k}(\frac{l}{2^{k}}(x+2n-1))$ for all $x\in(-1,1)$, then for $0\leq j\leq m$, we have

[TABLE]

Proof. For $0\leq j\leq m$, using the definition of the $L^{2}$- norm and the change of variable by setting $x^{\prime}=\frac{l}{2^{k-1}}(x+2n-1)$, we have

[TABLE]

Now we introduce the following semi-norm of $f\in H^{m}(0,l)$:

[TABLE]

for $0\leq r\leq m$, $M\geq 1$ and $k\geq 1$. Note that for $M\geq m$, we have

[TABLE]

Theorem 4.1**.**

Suppose that $f\in H^{m}(0,l)$ with $m\geq 1$ and $\linebreak{\Psi_{k,M-1}(f)=\sum\limits_{n=1}^{{2^{k-1}}}{\sum\limits_{m^{\prime}=0}^{M-1}{{a_{n,m^{\prime}}}{\psi_{n,m^{\prime}}}(x)}}}$ is the best approximation of $f$ based on the Legendre wavelets. Then,

[TABLE]

and for $1\leq r\leq m$,

[TABLE]

where $C$ denotes a positive constant that is independent of $M$ and $k$ but depends on the length $l$.

Proof. Consider the function $f_{n,k}:J_{k,n}\rightarrow\mathbb{R}$ such that $f_{n,k}(x)=f(x)$, for all $x\in J_{k,n}$. Then from equation (22) and Lemma 4.2 for $r\geq 0$ we have

[TABLE]

Setting $r=0$, we obtain

[TABLE]

here, we used equation (23) for the second inequality and Lemma 4.2 for the third inequality. This yields equation (29). Furthermore, for $1\leq r\leq m$ and $k\geq 1$, the following inequality is observed:

[TABLE]

where for the third inequality, equation (24) is used. Therefore, we have proved (30).

Theorem 4.2**.**

Suppose that $y(x)\in H^{m}(0,l)$ is the exact solution of equation (1), $\Psi_{k,M-1}(y)=\sum\limits_{n=1}^{2^{k-1}}\sum\limits_{m^{\prime}=0}^{M-1}a_{n,m^{\prime}}\psi_{n,m^{\prime}}(x)=A^{T}\psi(x)$ is its best approximation using the Legendre wavelets ($M\geq 3$) and $y_{k,M-1}(x)=\sum\limits_{n=1}^{2^{k-1}}\sum\limits_{m^{\prime}=0}^{M-1}\bar{a}_{n,m^{\prime}}\psi_{n,m^{\prime}}(x)=\bar{A}^{T}\psi(x)$ is the approximate solution obtained by the proposed method in this paper. Then, we have

[TABLE]

where the norm in $||A-\bar{A}||_{2}$ is the usual Euclidean norm for vectors.

Proof. Using the triangle property of $L^{2}$-norm and (29), we obtain

[TABLE]

We have

[TABLE]

Taking (32) and (33) into consideration, we obtain (34).

Remark 4.1**.**

Theorem 4.2 suggests that if the values of $M$ and $k$ are large simultaneously, the numerical solution obtained by the proposed method in this paper would be divergent, specially when the exact solution $y(x)$ belongs to the space $H^{m}(0,l)$ with small $m$.

5 Numerical examples

In this section, we consider some examples and solve them using the present method. In our implementation, the method was carried out using Mathematica 11.3. The error $||{y-y_{k,M-1}}||$ is computed in the $L^{2}$-norm, where $y$ is the exact solution and $y_{k,M-1}$ is the numerical solution obtained by the present method.

Example 5.1**.**

As the first test example, we consider the following nonlinear fractional delay integro-differential equation

[TABLE]

with conditions

[TABLE]

The exact solution is $y(x)=x^{2}$. This very simple example, that can be solved manually, illustrates how the proposed method can be applied to a boundary value problem. We employ the method with $k=1$, $M=3$, and $l=2$. So, we have

[TABLE]

where

[TABLE]

The first and second derivatives of the functions $\psi_{1,m}$, $i=0,1,2$, are given by

[TABLE]

[TABLE]

Moreover, we have

[TABLE]

We set $N=1$, and use the Jacobi polynomial $J_{1}^{(-0.5,0)}(x)=\frac{1}{4}(3x-1)$ in oreder to employ the Gauss-Jacobi quadrature formula for computing ${{}^{c}D^{1.5}}\psi_{1,2}$. The node and weight of the quadrature rule are easily given by

[TABLE]

Hence, using (15), we obtain

[TABLE]

Finally, using the collocation point $x_{1,2}=1-\frac{\sqrt{3}}{2}$ and the boundary conditions given by (21), the following system of algebraic equations are obtained

[TABLE]

By solving this system, we get

[TABLE]

By substituting these values into (34), we have $y(x)=x^{2}$ which is the exact solution.

Example 5.2**.**

Consider the following linear fractional delay pantograph equation

[TABLE]

with the initial conditions $y(0)=y^{\prime}(0)=0$. The exact solution of this problem is $y(x)=x^{7/2}$. We have solved this problem with different values of $k$ and $M$. The numerical results are shown in Table 1 and Figure 1. In Table 1 the error norm is displayed for different values of $k$ and $M$. In the case $M=3$ when $k$ is increased we see that the error norm tends to [math] like $2^{(-k)}$ (as $k\rightarrow\infty$), which agrees with the approximation error by a Legendre wavelet when the approximated function belongs to $H^{3}([0,1])$. For higher values of $M$ the convergence rate is even faster. In Figure 1, the exact solution and the approximate solutions given by $M=3$, $k=2,3,4,5$ (left) and $k=1$, $M=3,4,5,6$ (right) are displayed.

Example 5.3**.**

Consider the following nonlinear fractional pantograph equation [38]

[TABLE]

subject to the boundary conditions $y(0)=y(1)=1$. The exact solution, when $\alpha=2$, is $y(x)=1+x-x^{3}$.

We have solved this problem with different values of $k$ and $M$. We consider $M=3$, $k=1,2,\ldots,6$, $\alpha=2$ and display the numerical results for the errors in Table 2. A plot of the numerical solutions with $M=3$, $k=2,3,4,5$, $\alpha=2$ together with the exact solution is displayed in Figure 2 (left). We can see the convergence of the numerical solutions to the exact solution as the value of $k$ increases. The exact solution for the case $\alpha=2$ and the approximate solutions with $M=4$, $k=1$, and different values of $\alpha$ are seen in Figure 2 (right). As it could be expected, when $\alpha$ is close to $2$, the approximate solution gets close to the exact solution of the case $\alpha=2$. Note that in the case $\alpha=2$, since the exact solution is a third-degree polynomial, our method produces the exact solution with $M\geq 4$. In the case $M=3$, the error tends to zero approximately as $2^{(-k)}$ (as $k\rightarrow\infty$), as it happens in the previous example.

Example 5.4**.**

Consider the following fractional Volterra delay integro-differential equation [38]

[TABLE]

The exact solution of this problem, when $\alpha=1$, is $y(x)=e^{x}$.

Considering $\alpha=1$ and substituting $x=0$ in (35), we get $y^{\prime}(0)=1$. By taking the derivative of (35), we obtain the following equation

[TABLE]

We consider $l=3$ and solve this problem. A comparison between the absolute error given by our method and the first kind Chebyshev wavelet method introduced in [38] is displayed in Table 3 for $\alpha=1$ and with $k=1$, $M=10,20$. Furthermore, Table 4 demonstrates the relationship of the errors with the change of $M$ with considering $k=1,2,\ldots,6$ for the classical non-fractional equations corresponding to Example 5.4. The results of this table suggest that the error tends to zero very fast as $k$ and $M$ tend to infinity. This could be expected since according to Theorem 4.2 when a given function is sufficiently smooth the approximation by the Legendre wavelets has spectral convergence. The comparison between the exact solution and approximate solutions obtained in the case $M=3$ and different values of $k$ for $\alpha=1$ is shown in Figure 3 (left). From this figure, it can be seen that the numerical solutions converge to the exact solution as $k$ increases. Finally, the numerical solutions resulted by different values of $\alpha$ and $M=8$, $k=1$ together with the exact solution of the case $\alpha=1$ are plotted in Figure 3 (right). As it could be seen, the approximate solution gets close to the exact solution of the case $\alpha=1$ when $\alpha$ is close to $1$.

6 Concluding remarks

In this work, a collocation method based on the Legendre wavelets combined with the Gauss–Jacobi quadrature formula is introduced to solve a class of nonlinear delay-type fractional integro-differential equations. First, an approximation of the unknown function using the Legendre wavelets basis functions is considered. Then, by substituting this approximation in the problem and using the initial condition, we define suitable residual functions. We use the Gauss–Jacobi quadrature formula in order to compute the fractional derivative of the Legendre wavelets. By collocating the residual at the shifted Gauss–Chebyshev points, a system of nonlinear algebraic equations is obtained. In order to get a continuous solution at the subinterval division nodes, we consider some additional conditions. By using the considered approximation of the solution, the initial or boundary conditions are taken into account. Finally, by solving the resulting system, the unknown coefficients of the solution are computed. Some error bounds are given for the error of the best approximation of a function using the Legendre wavelets functions. To illustrate the efficiency and high accuracy of the suggested method in this paper, five numerical examples have been solved. The numerical examples show the great advantage of using the wavelets method for solving this type of equations. Actually if the solution of the equation is sufficiently smooth by increasing the degree of polynomial ($M$) we can achieve spectral accuracy. But even if the solution has discontinuous derivatives of high order, by choosing moderate value of $M$, fast convergence can be achieved by increasing the level of resolution ($k$).

Acknowledgements

P.M. Lima acknowledges support from Fundação para a Ciência e a Tecnologia (the Portuguese Foundation for Science and Technology) through the grant SFRH/BSAB/135130/2017 and UID/Multi/04621/2019.