Numerical solutions to systems of fractional Voltera Integro differential equations, using Chebyshev wavelet method

ABSTRACT

The Chebyshev Wavelet Method (CWM) is applied to evaluate the numerical solutions of some systems of linear fractional Voltera integro differential equations (FVIDEs). The applicability and validity of the proposed method is ensured by discussing some illustrative examples. The numerical results obtained by this technique are compared with the exact solutions of the problems. The error analysis reveals that the accuracy of the present method is higher than any existing numerical method.

KEYWORDS:

• Fractional Voltera Integro differential equations
• Chebyshev wavelet method
• iterative method
• numerical solution

AMS subject classifications:

• 65L05
• 65L06

View correction statement:

Correction

1. Introduction

The theory and applications of fractional calculus can be observed in many fields of science and engineering such as nonlinear oscillation of earth quakes [1], fluid dynamic traffic [2] and signal processing [3].

Due to precious contribution of fractional calculus in various fields of science and engineering, the researchers have shown great interest to study fractional calculus. In this regard as in many cases, it is very difficult to find the exact or analytical solutions of fractional differential and integral equations. The numerical methods have gained importance to avoid this difficulty. Initially, the authors have used different numerical techniques to find the approximate solution of fractional differential and integral equations such as Adomian Decomposition Method (ADM) [4], Spline Collocation Method (SCM) [5], Fractional Transform Method (FTM) [6], Homotopy Perturbation Method (HPM) [7], Operational Tau Method (OTM) [8], Shifted Chebyshev Polynomial Method (SCPM) [9], Rationalized Haar Functions Method (RHFM) [10] and Reproducing Kernel Hilbert Space (RKHSM) [11,12].

Besides these methods, most of the authors have applied a comparatively new numerical techniques based on wavelets [13–16]. The methods based on Chebyshev wavelets have gained much importance during the last decade because of its simple, effective and straightforward implementation. Therefore, the mathematician paid great attention to these methods for solving problems in different fields of science and engineering. Examples are Chebyshev Wavelet Operational Matrix (CWOM) [17], Chebyshev Finite Difference Method (CFDM) [18], Shifted Chebyshev Polynomial Method (SCPM) [9] and Chebyshev Wavelet Method (CWM) [19,20].

The solution to fractional system of differential equations is also a point of investigation for many researchers. Therefore, they introduced and extended many numerical techniques for solving these fractional systems of differential equations. Examples are Adomian Decomposition Method (ADM) [21] and B-Spline method [22].

In this paper, we approximate the solution of fractional systems of Volterra Integro Differential Equations using an efficient Chebyshev Wavelet method (CWM). The simulations are done by the present method and a very useful Chebyshev Wavelet algorithm is developed. The numerical results found by the present method are compared with exact solution of the problem, showing the greatest degree accuracy.

2. Preliminaries and definitions

In this section, we present some definitions and other mathematical preliminaries for the completion of the current work.

Definition 2.1:

The Riemann fractional integral operator $I^{\mu }$ of order $\gamma$ on the usual Lebesgue space $L_1\left[a,b\right]$ is given by (1) $\left(I^{\mu }g\right)\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\mu \right)}{\int }_0^t(t-\xi {)}^{\mu -1}g\left(\xi \right)\text{d}\xi ,\gamma >0,$(1) $\left(I^{0}g\right)\left(t\right)=g\left(t\right).$ This integral operator has the following properties

1. $I^{\mu }{I}^{\eta }=I^{\mu +\eta }$,

2. $I^{\mu }{I}^{\eta }=I^{\eta }{I}^{\mu }$,

3. $I^{\mu }(t-a{)}^{\upsilon }=\frac{\mathrm{\Gamma }\left(\upsilon +1\right)}{\mathrm{\Gamma }\left(\mu +\upsilon +1\right)}(t-a{)}^{\mu +\upsilon }$

where $L_1\left[a,b\right]$, $\alpha$, $\mu \ge 0,$ and $\upsilon >-1$,

Definition 2.2:

The Caputo definition of fractional differential operator is given by $\begin{array}{rl}\left(D^{\mu }g\right)\left(t\right)& =\frac{1}{\mathrm{\Gamma }\left(n-\mu \right)}{\int }_0^t(t-\xi {)}^{n-\mu -1}{g}^{\left(n\right)}\left(\xi \right)\text{d}\xi ,\\ & n-1⟨\mu ⟩n,\end{array}$ where $t>0,$ $n$ is an integer

It has the following two basic properties

1. $\left(D^{\mu }{I}^{\mu }g\right)\left(t\right)=g\left(t\right)$ ,

2. $\left(I^{\mu }{D}^{\mu }g\right)\left(t\right)=\left(g\right)\left(t\right)-(x+a{)}^n=\sum _{k=0}^n{f}^{\left(k\right)}\left(0^{+}\right)\frac{{\left(t-a\right)}^k}{k!},t>0$

3. Properties of the Chebyshev wavelets

Wavelets consist of family of functions generated from the dilation $m$ and translation $l$ of a single function $\psi \left(x\right)$ called the mother wavelet. When the dilation $a$ and translation $b$ change continuously, then we get the following continuous family of Wavelet [19] ${\psi }_{l,m}\left(x\right)=|a{|}^{\frac{1}{2}}\psi \left(\frac{x-l}{m}\right),l,m\in R,m\ne 0,$ If we restrict the parameters $l$ and $m$ to discrete values as $m={a_0}^{-k},l=n{b}_0{a_0}^{-k},a_0>1,b_0>1.$ We have the following family of discrete wavelets ${\psi }_{k,n}\left(x\right)=|a{|}^{\frac{k}{2}}\psi \left({a_0}^{k}x-n{b}_0\right),k,n\in Z,$ where ${\psi }_{k,n}$ form a wavelet basis for $L^2\left(R\right)$. Especially when $a_0=2$ and $b_0=1$, then ${\psi }_{k,n}\left(x\right)$ forms an orthogonal basis.

The second kind of Chebyshev wavelets is constituted of four parameters, ${\psi }_{n,m}\left(x\right)=\psi \left(k,n,m,x\right)$, where $n=1,2,\dots ,2^{k-1},$ $k$ is any positive integer, $m$ is the degree of the second Chebyshev polynomial. The Chebyshev wavelets are defined on the interval $0\le x<1$ as ${\psi }_{n,m}\left(x\right)=\left\{\begin{array}{c}{2}^{k/2}\tilde{T}\left(2^k-2n+1\right),\frac{n-1}{2}\le x\le \frac{n+1}{2}\\ 0,\text{otherwise}\end{array}\right.$ where (2) $\widetilde{T_m}\left(x\right)=\sqrt{\frac{2}{\pi }}{T}_m\left(x\right),m=0,1,2,\dots .,M-1$(2)

Here, $T_m\left(x\right)$ are the second Chebyshev polynomials of degree $m$ with respect to the weight function $w\left(x\right)=\sqrt{1-x^2}$ on the interval [−1,1], and satisfying the following recursive formula $T_0\left(x\right)=1,T_1\left(x\right)=2x,$ $T_{m+1}\left(x\right)=2x{T}_m\left(x\right)-T_{m-1}\left(x\right),m=1,2,3,\dots ,$

Lemma 3.1:

If the Chebyshev Wavelet expansion of a continuous function $f\left(x\right)$ converges uniformly, then the Chebyshev Wavelet expansion converges to the function $f\left(x\right)$.

Proof:

See [23].

Theorem 3.2:

A function $f\left(x\right)\in L_2\left[0,1\right]$, with bounded second derivative, say $|f^{\mathrm{\prime \prime }}\left(x\right)|\le N$, can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to $f\left(x\right)$, that is, $f\left(x\right)=\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{\mathrm{\infty }}{c}_{nm}{\psi }_{n,m}\left(x\right)$

Proof:

See [23].

4. Chebyshev wavelet method (CWM)

In this paper, we consider the fractional systems of Volterra integro differential equations $y^{\left(m\right)}\left(x\right)=f\left(x\right)+{\int }_0^{x}k(x,t,{y_1}^{\left(m\right)},{y_2}^{\left(m\right)},{y_3}^{\left(m\right)},\dots ,{y_n}^{\left(m\right)})\text{d}t,$ where $f\left(x\right)=[f_1\left(x\right),f_2\left(x\right),\dots ,f_n\left(x\right){]}^{\mathrm{T}}$ And (3) $\begin{array}{rl}{y}^{\left(m\right)}\left(x\right)& =[{y_1}^{\left(m\right)},{y_2}^{\left(m\right)},{y_3}^{\left(m\right)},\dots ,{y_n}^{\left(m\right)}{]}^{\mathrm{T}},\\ k& =[k_1,k_2,\dots ,k_n{]}^{\mathrm{T}},m=1,2,\dots ,\end{array}$(3) with the initial conditions $y^{\left(m\right)}\left(0\right)=b_m.$ The solution to system (3) can be expended by Chebyshev wavelets series as (4) $\begin{array}{rl}{y}_1\left(x\right)& =\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{\mathrm{\infty }}{c}_{n,m}{\psi }_{n,m}\left(x\right),y_2\left(x\right)\\ & =\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{\mathrm{\infty }}{d}_{n,m}{\psi }_{n,m}\left(x\right),y_n\left(x\right)\\ & =\sum _{n=1}^{\mathrm{\infty }}\sum _{m=0}^{\mathrm{\infty }}{e}_{n,m}{\psi }_{n,m}\left(x\right),\left(4\right)\end{array}$(4) where ${\psi }_{n,m}\left(x\right)$ is given by Equation (1)(1) $\left(I^{\mu }g\right)\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(\mu \right)}{\int }_0^t(t-\xi {)}^{\mu -1}g\left(\xi \right)\text{d}\xi ,\gamma >0,$(1) . The series in Equation (4)(5) $\begin{array}{rl}{y_1}_{k,M}\left(x\right)& =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(x\right),{y_2}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{d}_{n,m}{\psi }_{n,m}\left(x\right),{y_n}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{2^{k-1}}{e}_{n,m}{\psi }_{n,m}\left(x\right),\end{array}$(5) are truncated as (5) $\begin{array}{rl}{y_1}_{k,M}\left(x\right)& =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(x\right),{y_2}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{d}_{n,m}{\psi }_{n,m}\left(x\right),{y_n}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{2^{k-1}}{e}_{n,m}{\psi }_{n,m}\left(x\right),\end{array}$(5) This implies that there are $2^{k-1}M\times 2^{k-1}M\dots n-\text{time}\dots \times 2^{k-1}M$ conditions to determine $2^{k-1}M\times 2^{k-1}M\dots n-\text{time}\dots \times 2^{k-1}M$ coefficients $c_{i,j}$, $d_{i,j}$, … , $e_{i,j}$. We put these coefficients in Equation (5)(6) ${y_1}_{k,M}\left(0\right)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(0\right)={\alpha }_0,$(6) to obtain the approximate solution by Chebyshev Wavelet method (CWM).

In this article, we considered fractional volterra systems of order one or less than one, two or less than two and three or less than three.

First, we consider system of order one or less than one, which consists of two equations and two unknowns. $\begin{array}{rl}\frac{\text{d}{y}_1}{\text{d}x}& =1+x-\frac{1}{3}{x}^3+{\int }_0^x(\left(x-t\right)y_1\left(t\right)\\ & +\left(x-t\right)y_2\left(t\right))\text{d}t,\end{array}$ $\begin{array}{rl}\frac{\text{d}{y}_2}{\text{d}x}& =1-x-\frac{1}{12}{x}^4+{\int }_0^x(\left(x-t\right)y_1\left(t\right)\\ & -\left(x-t\right)y_2\left(t\right))\text{d}t\end{array}$ with initial conditions $y_1\left(0\right)=0,\frac{d{y}_1}{dx}=1,y_2\left(0\right)=0,\frac{\text{d}{y}_2}{\text{d}x}=1.$ The procedure is as follows, the initial conditions for both $y_1$ and $y_2$ are approximated as (6) ${y_1}_{k,M}\left(0\right)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(0\right)={\alpha }_0,$(6) (7) ${y_2}_{k,M}\left(0\right)=\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(0\right)={\beta }_0,$(7) (8) $\frac{\text{d}{y_1}_{k,M}\left(0\right)}{\text{d}x}=\frac{\text{d}}{\text{d}x}\left(\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(0\right)\right)={\alpha }_1,$(8) (9) $\frac{\text{d}{y_2}_{k,M}\left(0\right)}{\text{d}x}=\frac{\text{d}}{\text{d}x}\left(\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(0\right)\right)={\beta }_1.$(9)

The remaining $2^{k-1}M\times 2^{k-1}M-4$ conditions can be obtain by substituting Equation (4)(5) $\begin{array}{rl}{y_1}_{k,M}\left(x\right)& =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{c}_{n,m}{\psi }_{n,m}\left(x\right),{y_2}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-1}{d}_{n,m}{\psi }_{n,m}\left(x\right),{y_n}_{k,M}\left(x\right)\\ & =\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{2^{k-1}}{e}_{n,m}{\psi }_{n,m}\left(x\right),\end{array}$(5) in Equation (3)(3) $\begin{array}{rl}{y}^{\left(m\right)}\left(x\right)& =[{y_1}^{\left(m\right)},{y_2}^{\left(m\right)},{y_3}^{\left(m\right)},\dots ,{y_n}^{\left(m\right)}{]}^{\mathrm{T}},\\ k& =[k_1,k_2,\dots ,k_n{]}^{\mathrm{T}},m=1,2,\dots ,\end{array}$(3) , by taking $n=2$,$m=1$, and $0<\alpha \le 1$, we get (10) $\begin{array}{rl}& \frac{{\text{d}}^{\alpha }}{\text{d}{x}^{\alpha }}\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{c}_{n,m}{\psi }_{n,m}\left(x_i\right)\\ & =f\left(x\right)+{\int }_0^{x}k\left(x,t,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{c}_{n,m}{\psi }_{n,m}\left(x_i\right),\right.\\ & \times \left.\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{d}_{n,m}{\psi }_{n,m}\left(x_i\right)\right)\text{d}t,\end{array}$(10) (11) $\begin{array}{rl}& \frac{{\text{d}}^{\alpha }}{\text{d}{x}^{\alpha }}\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{\text{d}}_{n,m}{\psi }_{n,m}\left(x_i\right)\\ & =g\left(x\right)+{\int }_0^{x}k\left(x,t,\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{c}_{n,m}{\psi }_{n,m}\left(x_i\right),\right.\\ & \times \left.\sum _{n=1}^{2^{k-1}}\sum _{m=0}^{M-5}{d}_{n,m}{\psi }_{n,m}\left(x_i\right)\right)\text{d}t.\end{array}$(11)

Assume that Equations (10) and (11) are exact at $2^{k-1}M-4$ points. Then $x_i$ points are calculated by the following formula $x_i=\frac{i-0.5}{2^{k-1}M}fori=1,2,3,\dots 2^{k-1}M-4.$ The combination of Equations (6), (7), (8), (9), (10) and (11) forms the linear system of $2^{k-1}M\times 2^{k-1}M$ linear equations. The unknown $c_{i,j}$ and $d_{i,j}$are calculated through the solution of this system.

The same procedure can be repeated for other initial value fractional systems of Volterra Integro differential equations.

5. Numerical examples

Example 5.1 :

Consider the following fractional system of Volterra Integro differential equation of order, $0<\alpha \le 1$ $\begin{array}{rl}& \frac{{\text{d}}^{\alpha }{y}_1}{\text{d}{x}^{\alpha }}-1-x+\frac{1}{3}{x}^3-{\int }_0^x(\left(x-t\right)y_1\left(t\right)\\ & +\left(x-t\right)y_2\left(t\right))\text{d}t,\end{array}$ $\begin{array}{rl}& \frac{{\text{d}}^{\alpha }{y}_2}{\text{d}{x}^{\alpha }}-1+x+\frac{1}{12}{x}^4-{\int }_0^x(\left(x-t\right)y_1\left(t\right)\\ & -\left(x-t\right)y_2\left(t\right))\text{d}t,\end{array}$ with initial conditions $y_1\left(0\right)=0,\frac{d{y}_1}{dx}\left(0\right)=1,y_2\left(0\right)=0,\frac{d{y}_2}{dx}\left(0\right)=1.$ The exact solution of the system is $y_1\left(x\right)=x+\frac{x^2}{2}$ and $y_2\left(x\right)=x-\frac{x^2}{2}.$

In Table 1, the exact solution and Chebyshev wavelet method approximations are shown for the system solution $y_1\left(x\right)$ and$y_2\left(x\right)$. The approximate solutions are represented by $y_1$ (CWM) and $y_2$ (CWM). The corresponding errors are represented by Error $y_1$ and Error $y_2$ associated with $y_1$ and $y_2$, respectively. The method is applied for $M=19,K=1$. The absolute errors between the exact and approximate solutions are measured. The error investigation from the table shows that the proposed method has higher accuracy.

Table 2, analysed the approximate solutions for Example 5.1 for different values of $\alpha$ such that $0<\alpha \le 1$. This investigation shows that as the value of $\alpha$ increases from 0.75 to 1, the accuracy is increasing and attained its maximum accuracy at $\alpha =1$.

Table 1. The numerical results of Example 5.1.

$x_i$	$y_1$	$y_1$ (CWM)	$y_2$	$y_2$ (CWM)	Error $y_1$	Error $y_2$
0.0	0.000000	0.00000000000000000001	0.000	0.00000000000000000000	8.86E−21	5.22E−22
0.1	0.105000	0.10500000000000017529	0.095	0.09500000000000009781	1.75E−16	9.78E−17
0.2	0.220000	0.22000000000000017686	0.180	0.18000000000000009876	1.76E−16	9.87E−17
0.3	0.345000	0.34500000000000017641	0.255	0.25500000000000010000	1.76E−16	1.00E−16
0.4	0.480000	0.48000000000000017144	0.320	0.32000000000000010673	1.71E−16	1.06E−16
0.5	0.625000	0.62500000000000014947	0.375	0.37500000000000013167	1.49E−16	1.31E−16
0.6	0.780000	0.78000000000000008278	0.420	0.42000000000000020303	8.28E−17	2.03E−16
0.7	0.945000	0.94499999999999991864	0.455	0.45500000000000037359	8.12E−17	3.73E−16
0.8	1.120000	1.11999999999999956790	0.480	0.48000000000000073240	4.32E−16	7.32E−16
0.9	1.305000	1.30499999999999889160	0.495	0.49500000000000141871	1.10E−15	1.41E−15
1.0	1.500000	1.49999999999999769060	0.500	0.50000000000000264463	2.30E−15	2.64E−15

Table 2. Numerical results of Example 5.1, for different fractional order $0<\alpha \le 1$.

$x_i$	$y_1\left(0.75\right)$	$y_2\left(0.75\right)$	$y_1\left(0.99\right)$	$y_2$ (0.99)	$y_1$ (0.99999)	$y_2$ (0.99999)
0.0	0.000000000	0.0000000	0.000000	0.000000	0.00000000	0.00000000
0.1	0.176199660	0.15340263	0.106888	0.102015	0.10500185	0.09500781
0.2	0.339855138	0.25283204	0.223632	0.258342	0.22000359	0.18008987
0.3	0.496712043	0.28935057	0.350566	0.592042	0.34500559	0.25538802
0.4	0.653952604	0.24973061	0.490544	1.245713	0.48001105	0.32106665
0.5	0.812194283	0.11645073	0.651549	2.355534	0.62502893	0.37728292
0.6	0.968110842	−0.1264186	0.851227	4.042822	0.78007915	0.42417707
0.7	1.113946717	−0.4917918	1.122739	6.413404	0.94519920	0.46187177
0.8	1.236233568	−0.9885288	1.522087	9.555548	1.12045253	0.49046985
0.9	1.314090673	−1.6188345	2.136990	13.52723	1.30593864	0.51003968
1.0	1.317517557	−2.3725543	3.096760	18.32412	1.50180445	0.52057835

Figure 1. The solution graph of Example 5.1, for $y_1\left(x\right)$ at different fractional order.

Figure 2. The solution graph of Example 5.1, for $y_2\left(x\right)$ at different fractional order.

In Table 3, Error$y_1\left(0.75\right)$, Error $y_1\left(0.99\right)$ and Error $y_1\left(0.99999\right)$ show the respective errors of $y_1$ at $\alpha =0.75$, 0.99 and 0.99999. Similarly, Error $y_2\left(0.75\right)$, Error $y_2\left(0.99\right)$ and Error $y_2\left(0.99999\right)$ are the errors generated with $y_2$ for $\alpha =0.75$, 0.99 and 0.99999. The error analysis again investigates that solutions at different values of $\alpha$ converges to integer order solutions (Figures 1 and 2).

Table 3. Errors of Example 5.1, for different fractional order $0<\alpha \le 1$.

$x_i$	Error $y_1\left(0.75\right)$	Error $y_2\left(0.75\right)$	Error $y_1\left(0.99\right)$	Error $y_2$(0.99)	Error $y_1$(0.99999)	Error $y_2$(0.99999)
0.0	1.34E−20	6.79E−21	3.30E−20	4.10E−19	2.14E−21	5.34E−21
0.1	7.11E−02	5.84E−02	1.88E−03	7.01E−03	1.85E−06	7.81E−06
0.2	1.19E−01	7.28E−02	3.63E−03	7.83E−02	3.59E−06	8.98E−05
0.3	1.51E−01	3.43E−02	5.56E−03	3.37E−01	5.59E−06	3.88E−04
0.4	1.73E−01	7.02E−02	1.05E−02	9.25E−01	1.10E−05	1.06E−03
0.5	1.87E−01	2.58E−01	2.65E−02	1.9805	2.89E−05	2.28E−03
0.6	1.88E−01	5.46E−01	7.12E−02	3.6228	7.91E−05	4.17E−03
0.7	1.68E−01	9.46E−01	1.77E−01	5.9584	1.99E−04	6.87E−03
0.8	1.16E−01	1.4685	4.02E−01	9.0755	4.52E−04	1.04E−02
0.9	9.09E−03	2.1138	8.31E−01	13.032	9.38E−04	1.50E−02
1.0	1.82E−01	2.8725	1.596760	17.824	1.80E−03	2.05E−02

Example 5.2:

Consider the following fractional system of Volterra integro differential equation of order, $1<\alpha \le 2$ $\frac{{\text{d}}^{\alpha }{y}_1}{\text{d}{x}^{\alpha }}+x^3+x^4-{\int }_0^x\left(3y_2\left(t\right)+4y_3\left(t\right)\right)\text{d}t=0,$ $\frac{{\text{d}}^{\alpha }{y}_2}{\text{d}{x}^{\alpha }}-2-x^2+x^4-{\int }_0^x\left(4y_3\left(t\right)-2y_1\left(t\right)\right)\text{d}t=0,$ $\frac{{\text{d}}^{\alpha }{y}_3}{\text{d}{x}^{\alpha }}-6x+x^2-x^3-{\int }_0^x\left(2y_1\left(t\right)-3y_2\left(t\right)\right)\text{d}t=0.$ With the initial conditions are $\begin{array}{rl}{y}_1\left(0\right)& =0,y_2\left(0\right)=0,y_3\left(0\right)=0,\frac{\text{d}{y}_1}{\text{d}x}\left(0\right)=1,\\ \frac{\text{d}{y}_2}{\text{d}x}\left(0\right)& =0,\frac{\text{d}{y}_3}{\text{d}x}\left(0\right)=0.\end{array}$

The exact solution of the system is $y_1\left(x\right)=x$, $y_2\left(x\right)=x^2$, $y_3\left(x\right)=x^3$

Table 4 explained the numerical results of Example 5.2. The exact solutions for $y_1$,$y_2$ and $y_3$ are given by $y_1$ (exact), $y_2$ (exact) and $y_3$ (exact), respectively. The approximate solutions obtained by Chebyshev wavelet method, for $y_1$,$y_2$ and $y_3$, are $y_1$ (CWM), $y_2$ (CWM) and $y_3$ (CWM), respectively.

Table 4. The numerical solutions of Example 5.2.

$x_i$	$y_1$ (exact)	$y_1$ (CWM)	$y_2$ (exact)	$y_2$ (CWM)	$y_3$ (exact)	$y_3$ (CWM)
0.0	0.0000	0.0000000000000	0.0000	0.000000000000000	0.0000	0.00000000000000000000
0.1	0.1000	0.0999999999999	0.0100	0.010000000000013	0.0010	0.00100000000003054812
0.2	0.2000	0.1999999999990	0.0400	0.040000000000029	0.0080	0.00800000000092897704
0.3	0.3000	0.2999999999864	0.0900	0.090000000000049	0.0270	0.02700000001354235333
0.4	0.4000	0.3999999999091	0.1600	0.160000000000118	0.0640	0.06400000009091300146
0.5	0.5000	0.4999999996089	0.2500	0.250000000000588	0.1250	0.12500000039130475950
0.6	0.6000	0.5999999987301	0.3600	0.360000000003600	0.2160	0.21600000127162549054
0.7	0.7000	0.6999999966051	0.4900	0.490000000022462	0.3430	0.34300000340941206265
0.8	0.8000	0.7999999921595	0.6400	0.640000000141051	0.5120	0.51200000795371335706
0.9	0.9000	0.8999999840488	0.8100	0.810000000839143	0.7290	0.72900001670629419916
1.0	1.0000	0.9999999719332	1.0000	1.000000004487292	1.0000	1.00000003232985799840

Table 5 analysed the errors associated with the solutions $y_1$, $y_2$ and $y_3$ for$\alpha =2$. These errors are denoted by Error $y_1$, Error $y_2$ and Error $y_3$. The absolute error between the exact and approximate solution is obtained which shows the desired degree of accuracy. The numerical simulations are handled by using $k=1$ and $M=19$ in the current method.

Table 5. The errors of Example 5.2 for $\alpha =2$.

	Error $y_1$	Error $y_2$	Error $y_3$
$x_i$	$\alpha =2$	$\alpha =2$	$\alpha =2$
0.0	1.85E−21	3.34E−21	2.91E−21
0.1	2.22E−14	1.30E−14	3.05E−14
0.2	9.10E−13	2.95E−14	9.28E−14
0.3	1.35E−11	4.95E−14	1.35E−11
0.4	9.08E−11	1.18E−13	9.09E−11
0.5	3.91E−10	5.88E−13	3.91E−10
0.6	1.26E−09	3.60E−12	1.27E−09
0.7	3.39E−09	2.24E−11	3.40E−09
0.8	7.84E−09	1.41E−10	7.95E−09
0.9	1.59E−08	8.39E−10	1.67E−08
1.0	2.80E−08	4.48E−09	3.23E−08

In Table 6, the numerical solutions of Example 5.2 for $y_1$, $y_2$ and $y_3$ are given at different fractional orders $\alpha$ such that $1<\alpha \le 2$. Here, ${1.9}^4=1.9999$, ${1.9}^7=1.9999999$ and ${1.9}^{12}=1.999999999999$.

Table 6. Numerical results of Example 5.2 for different orders $1<\alpha \le 2$.

	$y_1$	$y_1$	$y_1$	$y_2$	$y_2$	$y_2$	$y_3$	$y_3$	$y_3$
$x_i$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$
0.0	0.000000	0.000000	0.0000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
0.1	0.099999	0.099999	0.0999999	0.008722	0.009998	0.009999	0.001000	0.001000	0.0010000
0.2	0.199989	0.199999	0.1999999	0.022193	0.039982	0.039999	0.008013	0.008000	0.0080000
0.3	0.299831	0.299999	0.2999999	0.011349	0.089921	0.089999	0.027175	0.027000	0.0270001
0.4	0.398859	0.399998	0.3999999	−0.05789	0.159782	0.159999	0.065154	0.064001	0.0640011
0.5	0.495066	0.499995	0.4999999	−0.21936	0.249532	0.249999	0.129958	0.125004	0.1250049
0.6	0.583904	0.599983	0.5999999	−0.50564	0.359137	0.359999	0.232116	0.216016	0.2160160
0.7	0.656706	0.699956	0.6999999	−0.94931	0.488565	0.489998	0.386155	0.343043	0.3430430
0.8	0.698732	0.799899	0.7999999	−1.58401	0.637783	0.639998	0.611985	0.512099	0.5120996
0.9	0.686822	0.899787	0.8999998	−2.44344	0.806757	0.809997	0.935064	0.729205	0.7292053
1.0	0.586575	0.999588	0.9999996	−3.55600	0.995458	0.999996	1.382024	1.000380	1.0003806

Table 7 emphasizes on error analysis. The errors for all $y_1$, $y_2$ and $y_3$ are computed at different fractional orders. This error analysis shows the error decreases as the fractional order approaches to integer order (Figures 3–5).

Figure 3. Solution graph of Example 5.2 for $y_1\left(x\right)$ at different fractional order.

Figure 4. Solution graph of Example 5.2 for $y_2\left(x\right)$ at different fractional order.

Figure 5. Solution graph of Example 5.2 for $y_3\left(x\right)$ at different fractional order.

Table 7. The error analysis of Example 5.2 for different order$\alpha ={1.9}^4,{1.9}^7,{1.9}^{12}$.

	Error $y_1$	Error $y_1$	Error $y_1$	Error $y_2$	Error $y_2$	Error $y_2$	Error $y_3$	Error $y_3$	Error $y_3$
$x_i$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$	$\alpha ={1.9}^4$	$\alpha ={1.9}^7$	$\alpha ={1.9}^{12}$
0.0	4.71E−21	5.84E−21	7.49E−21	7.41E−20	4.28E−21	1.84E−21	4.40E−21	3.11E−21	1.68E−21
0.1	9.47E−08	9.44E−11	8.37E−14	1.27E−03	1.27E−06	1.11E−09	4.65E−07	4.65E−11	8.44E−14
0.2	1.09E−05	1.08E−08	9.48E−12	1.78E−02	1.77E−05	1.54E−08	1.32E−05	1.32E−08	9.50E−12
0.3	1.68E−04	1.67E−07	1.46E−10	7.86E−02	7.83E−05	6.82E−08	1.75E−04	1.74E−07	9.89E−10
0.4	1.14E−03	1.13E−06	9.89E−10	2.17E−01	2.17E−04	1.89E−07	1.15E−03	1.15E−06	9.90E−10
0.5	4.93E−03	4.91E−06	4.28E−09	4.69E−01	4.67E−04	4.07E−07	4.95E−03	4.94E−06	4.28E−09
0.6	1.60E−02	1.60E−05	1.39E−08	8.65E−01	8.62E−04	7.51E−07	1.61E−02	1.60E−05	1.39E−08
0.7	4.32E−02	4.31E−05	3.75E−08	1.4393	1.43E−03	1.24E−06	4.31E−02	4.30E−05	3.73E−08
0.8	1.01E−01	1.00E−04	8.78E−08	2.2240	2.21E−03	1.93E−06	9.99E−02	9.96E−05	8.66E−08
0.9	2.31E−01	2.12E−04	1.84E−07	3.2534	3.24E−03	2.82E−06	2.06E−01	2.05E−04	1.78E−07
1.0	4.13E−01	4.11E−04	3.58E−07	4.5560	4.54E−03	3.95E−06	3.82E−01	3.80E−04	3.31E−07

Example 5.3:

Consider the following fractional system of Volterra Integro differential equation of order, $2<\alpha \le 3$ $\frac{{\text{d}}^{\alpha }}{\text{d}{x}^{\alpha }}{y}_1+2x+2x^3+\frac{2}{5}{x}^5-{\int }_0^x\left({y_1}^2\left(t\right)+{y_2}^2\left(t\right)\right)\text{d}t=0,$ $\frac{{\text{d}}^{\alpha }}{\text{d}{x}^{\alpha }}{y}_2+\frac{2}{3}{x}^3+\frac{1}{5}{x}^5-{\int }_0^x\left(x-t\right)\left({y_1}^2\left(t\right)-{y_2}^2\left(t\right)\right)\text{d}t=0,$ with the initial conditions $\begin{array}{rl}{y}_1\left(0\right)& =1,\frac{d{y}_1}{dx}\left(0\right)=1,\frac{{\text{d}}^2{y}_1}{\text{d}{x}^2}\left(0\right)=2,y_2\left(0\right)=1,\\ \frac{\text{d}{y}_2}{\text{d}x}\left(0\right)& =-1,\frac{{\text{d}}^2{y}_2}{\text{d}{x}^2}\left(0\right)=2.\end{array}$ The exact solution of the system is $y_1\left(x\right)=1+x+x^2,y_2\left(x\right)=1-x+x^2.$

Table 8 displayed the numerical results of Example 5.3 using Chebyshev wavelet method. The exact solutions for $y_1$ and $y_2$ are represented by $y_1$ (exact) and $y_2$ (exact). The Chebyshev approximate solutions are denoted by $y_1$ (CWM) and $y_2$ (CWM). The absolute of the errors between the exact and approximate solutions are measured which shows the highest degree of accuracy.

Table 8. Numerical results for Example 5.3.

$x_i$	$y_1$ (exact)	$y_1$ (CWM)	Error $y_1$	$y_2$ (exact)	$y_2$(CWM)	Error $y_2$
0.0	1.000	0.99999999999999999996	4.0E−20	1.000	0.99999999999999999997	3.0E−20
0.1	1.110	1.11000000000000000000	0.0000	0.910	0.91000000000000000001	1.0E−20
0.2	1.240	1.24000000000000000000	0.0000	0.840	0.84000000000000000023	2.3E−19
0.3	1.390	1.39000000000000000010	1.0E−19	0.790	0.79000000000000000063	6.3E−19
0.4	1.560	1.56000000000000000020	2.0E−19	0.760	0.76000000000000000126	1.2E−18
0.5	1.750	1.75000000000000000030	3.0E−19	0.750	0.75000000000000000212	2.1E−18
0.6	1.960	1.96000000000000000050	5.0E−19	0.760	0.76000000000000000316	3.1E−18
0.7	2.190	2.19000000000000000070	7.0E−19	0.790	0.79000000000000000446	4.4E−18
0.8	2.440	2.44000000000000000100	1.0E−18	0.840	0.84000000000000000597	5.9E−18
0.9	2.710	2.71000000000000000140	1.4E−18	0.910	0.91000000000000000767	7.6E−18
1.0	3.000	3.00000000000000000170	1.7E−18	1.000	1.00000000000000000960	9.6E−18

Table 9 expresses the error analysis of Example 5.3 for different fractional orders. The error analysis shows that there is a very small change between the numerical solution at fractional order as compared to integer order (Figures 6–9).

Figure 6. The graph of $y_1\left(exact\right)$ and $y_1\left(CWM\right)$ of Example 5.3.

Figure 7. The solution graph of $y_2\left(exact\right)$ and $y_2\left(CWM\right)$.

Figure 8. The error graph of $Error{y}_1\left(x\right)$ and for other fractional orders.

Figure 9. The error graph of $y_2\left(x\right)$ for different fractional order $\alpha$, where $0\le \alpha \le 2$.

Table 9. The numerical results of Example 5.3, for fractional orders,$2<\alpha \le 3$_.

	Error $y_1$	Error $y_1$	Error $y_1$	Error$y_2$	$y_2$	$y_2$
$x_i$	$\alpha =2.25$	$\alpha =2.50$	$\alpha =2.75$	$\alpha =2.25$	$\alpha =2.50$	$\alpha =2.75$
0.0	3.00E−20	1.00E−19	0.00000	0.000000	1.00E−20	0.000000
0.1	0.00000	1.00E−19	0.0000	0.000000	2.00E−20	3.00E−20
0.2	1.00E−19	2.00E−19	0.0000	0.000000	1.80E−19	2.00E−19
0.3	2.00E−19	2.00E−19	1.00E−19	2.00E−19	4.30E−19	5.10E−19
0.4	3.00E−19	3.00E−19	1.00E−19	6.70E−19	7.40E−19	9.50E−19
0.5	5.00E−19	5.00E−19	2.00E−19	9.70E−19	1.15E−18	1.52E−18
0.6	6.00E−19	6.00E−19	5.00E−19	1.28E−18	1.59E−18	2.20E−18
0.7	8.00E−19	8.00E−19	7.00E−19	1.59E−18	2.09E−18	2.96E−18
0.8	1.10E−18	9.00E−19	1.00E−18	1.91E−18	2.61E−18	3.85E−18
0.9	1.10E−18	1.20E−18	1.40E−18	2.26E−18	3.17E−18	4.80E−18
1.0	1.40E−18	1.40E−18	1.30E−18	2.60E−18	3.80E−18	5.90E−18

6. Conclusion

In this work, we have fully attempted to find the numerical solution of the fractional system of Volterra Integro differential equations by using Chebyshev wavelet method. The numerical procedure and methodology are done in a very straight forward and effective manner. The numerical accuracy is also a point of interest. During numerical simulations, we observed that the current method has the highest degree of accuracy. On the bases of current work, the researchers can extended this technique to some other fractional systems of ordinary and partial differential equations.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Hassan Khan http://orcid.org/0000-0001-6417-1181

Muhammad Arif http://orcid.org/0000-0003-1484-7643