Solving a System of Fractional-Order Volterra-Fredholm Integro-Differential Equations with Weakly Singular Kernels via the Second Chebyshev Wavelets Method

: In this paper, by means of the second Chebyshev wavelet and its operational matrix, we solve a system of fractional-order Volterra–Fredholm integro-differential equations with weakly singular kernels. We estimate the functions by using the wavelet basis and then obtain the approximate solutions from the algebraic system corresponding to the main system. Moreover, the implementation of our scheme is presented, and the error bounds of approximations are analyzed. Finally, we evaluate the efﬁciency of the method through a numerical example.


Introduction
Fractional calculus has been the interest of many scientists and engineers [1][2][3].Many engineering and science phenomena, such as the heat conduction problem, radiative equilibrium, elasticity and fracture mechanics [4], viscoelastic deformation, viscoelasticity, viscous fluid [5], continuous population [6] and so forth, are modeled using the fractional integro-differential equations with a weakly singular kernel, fractional differential equations, fractional integral equations and system of nonlinear Volterra integro-differential equations.Many applied problems are transformed into the system of fractional differential and integral equations by mathematical modeling [7][8][9][10].Consequently, it is essential to obtain the approximate solution of a system of integro-differential equations by numerical methods.
In recent years, researchers have proposed different methods for solving the system of differential equations.In 2015, Sahu and Ray [11] developed a numerical method based on the Legendre hybrid block pulse function to approximate the solution of nonlinear systems of Fredholm-Manhattan integral equations.In the same year, they presented another scheme to solve a system of nonlinear Volterra integro-differential equations using Legendre wavelets [6].Yüzbasi [12] solved the system of linear Fredholm-Volterra integro-differential equations, which includes the derivatives of unknown functions in integral parts using the Bessel collocation method.In 2016, Deif and Grace [13] developed a new iterative method to approximate the solution of a system of linear fractional differential integral equations.In 2019, Xie and Yi [14] developed a numerical method for solving a nonlinear system of fractional-order Volterra-Fredholm integro-differential equations based on block-pulse functions.In 2020, Saemi et al. [15] developed a solution for the system of fractional-order Volterra-Fredholm integro-differential equations based on Müntz-Legendre wavelets.
Wavelets are one of the most important tools used in various fields such as quantum mechanics, signal processing, image processing, time-frequency analysis and data compression [16].One of the methods that have been considered in recent years to solve various ordinary and fractional-order equations is the use of wavelets [17][18][19].Wavelets provide a detailed accurate representation of different types of functions and operators, and their relationship with numerical algorithms [20,21].The second Chebyshev wavelet is one of the wavelets, which has gained attention in solving many problems and is applicable for solving various types of Volterra integral equations with a weakly singular kernel [17], fractional-order nonlinear Fredholm integro-differential equations [20], fractional-order differential equations [22], a system of linear differential equations [23], fractional-order integro-differential equations with a weakly singular kernel [5], and Abel's integral equations [16].Approximation of equations using the second Chebyshev wavelets has been considered by many researchers, such as Zhu and Wang [17,24], Zhou and Xu [21], Wang and Fan [22], Tavassoli Kajani et al. [25], Zhou et al. [26], Yi et al. [27], Lal and Sharma [16], and Manchanda and Rani [23].
In this paper, we apply the second Chebyshev wavelets method to solve the system of fractional-order Volterra-Fredholm integro-differential equations with a weakly singular kernel.In fact, the main purpose of this study is solving the system of equations with singularity.The second Chebyshev wavelets method converts the system of fractionalorder Volterra-Fredholm integro-differential equations with a weakly singular kernel to a system of algebraic equations, which can be solved using the conventional linear methods.

Preliminaries
In this section, we introduce fractional-order operators, block-pulse functions and explain their features.Definition 1.The Riemann-Liouville fractional-order integral operator of order α is given by [17,27]: The properties of this operator are as follows: 1.
Definition 2. The Caputo fractional derivative operator of order α for a function f is given by [17,18,27]: where n ∈ N, and Γ(.) is the Gamma function.
The properties between the Caputo fractional-order derivative operator and the Riemann-Liouville fractional-order integral operator is given by the following expressions [17,18,24]: 2.
Definition 3. The set of block pulse function on [0, 1) is defined as: where i = 1, 2, . . ., m.Furthermore, the vector of block pulse functions is obtained as follows: and the important properties of these functions are as follows: Lemma 1.The block pulse function operational matrix of fractional-order integration F α is obtained by: where , For example with α = 0.5 and m = 6: .

The Second Chebyshev Wavelets and Function Approximation
In this section, we introduce the second Chebyshev wavelet and then use this basis to provide an approximation of functions.

The Second Chebyshev Wavelets and Their Properties
In this part, we introduce the second Chebyshev wavelet and its features.Definition 4. The second Chebyshev wavelet on the interval [0, 1) is defined as [16,23,27]: where n = 1, 2, . . ., 2 k−1 , m = 0, 1, . . ., M − 1, and k, M ∈ N. The coefficient 2 π is for orthonormality, and U m (t) is the Chebyshev polynomial of the second kind with degree m, which is as follows: Furthermore, the weight function of the second kind Chebyshev polynomials is , and with transmission and dilation, first we obtain ω(t) = ω(2t − 1), and then we get ω n (t) = ω(2 k t − 2n + 1) as the weight function of the second Chebyshev wavelets basis.
For example, with k = 2 and M = 3, we have n = 1, 2, m = 0, 1, 2, and for 0 ≤ t < 0.5, and also for 0.5 ≤ t < 1, The second Chebyshev wavelets have an orthonormal basis of L 2 [0, 1), i.e., where < ., .> ω n denotes the inner product.The second Chebyshev wavelet has compact support The Chebyshev wavelet charts for k = 3 and M = 4 are shown in Figure 1.According to the second Chebyshev wavelet, the vector of this wavelet is given by [16,27]: In other words, for 0 ≤ t < 0.5: , and for 0.5 ≤ t < 1: Moreover, for the collocation points [24] with m = 2 k−1 M, the second Chebyshev wavelets matrix is obtained as follows [24,25]: For k = 2 and M = 3, i.e., There is a relationship between the second Chebyshev wavelet and the block pulse function: If I α is the fractional-order integration operator of the second Chebyshev wavelets, one can achieve [17,24]: where P α is named as the operational matrix of fractional-order integration of the second Chebyshev wavelet.For example,

Function Approximation
Using the second Chebyshev wavelet, each function in L 2 can be approximated using the following lemma.Lemma 2. Any function f ∈ L 2 ([0, 1]) can be expanded into the second Chebyshev wavelet as [16,27,28]: where If Equation ( 6) is truncated, then with Equation (3)

Method Analysis
For solving the system (1), without reducing the generality of the equations under consideration, we assume that the initial conditions are zero, and we approximate D α y i (t), f i (t), and k i (t, s) for i = 1, 2 in terms of the second Chebyshev wavelet as follows: From Equations ( 2), ( 7) and ( 8), we obtain Thus, and from Equation ( 8) and 1 0 Ψ(s)Ψ(s) T ds = D, we have By substituting Equations ( 8)-( 11) into (1), we get and we obtain By solving system (13), one can get C 1 and C 2 .Then substituting them into (9), the unknown solutions can be obtained.

Error Analysis
In this section, we present an error estimation for the system of Equation ( 1).Theorem 1.Let ŷ1 (t) and ŷ2 (t) be the approximations of y 1 (t) and y 2 (t), obtained by the second Chebyshev wavelets basis (6), as the solutions of system (1) with 0 < β 1 , β 2 < 1/2.Assume also that there is a pair of constants k 1 and k 2 such that , then the approximate solutions of system (1) converge to the exact solutions with respect to L 2 norm.
Proof.Compute, by the assumptions, and Using the triangular inequality, we get and as a result, we obtain and Now, we denote e 2 (t), (14) and consequently by the assumption We have similar relations for the second approximation, i.e., and attention to assumption Substituting ( 17) into (14), and ( 15) into ( 16), we have where and .
On the other hand, according to the assumptions of this theorem, we get 1 − ε 1 > 0 and 1 − ε 2 > 0; also, we know that e 1 (t) ≥ 0 and e 2 (t) ≥ 0. Therefore, the relations in (18) are satisfied only for e 1 (t), e 2 (t) = 0.The conditions of this theorem are sufficient to find the appropriate approximate solution of system (1), but they can be met in certain cases rarely.Furthermore, according to the relation (6) and the expressed discretization, the approximation of the solution is done in some node points on [0, 1].In other words, according to Lemma 2, it is possible to get a suitable approximation so that e 1 (t) → 0 and e 2 (t) → 0 by the constant consideration of m and when k → ∞, even if the conditions of Theorem 1 do not hold.

Numerical Example
To demonstrate the efficiency of our proposed method, we consider the following example.
Example 1.Consider the system of fractional-order Volterra-Fredholm integro-differential equations with a weakly singular kernel: The exact solutions of this system are y 1 (t) = t and y 2 (t) = −t.The absolute errors of y 1 (t) and y 2 (t) for different values of t are listed in Tables 1 and 2.  3 shows the execution time.

Conclusions
In this paper, a numerical algorithm using the second Chebyshev wavelet was proposed for a system of fractional-order Volterra-Fredholm integro-differential equations with weakly singular kernels in the Volterra part.Applying the properties of the second Chebyshev wavelet, we transformed the main system into an algebraic system of equations with a sparse coefficient matrix.By solving this system, an approximate solution was obtained for the system of fractional integral equations.Furthermore, the error analysis of the proposed approach was presented.

Figures 2 -
Figures 2-5 display the results of comparing the errors of different approximations by the second Chebyshev wavelets method for various values of k and M. Furthermore, Table3shows the execution time.
t M =

Table 3 .
Run times of Example 1.