Next Article in Journal
The Dynamics of a Fractional-Order Mathematical Model of Cancer Tumor Disease
Previous Article in Journal
A Symmetric Extensible Protocol for Quantum Secret Sharing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Computation of Mixed Volterra–Fredholm Integro-Fractional Differential Equations by Using Newton-Cotes Methods

by
Shazad Shawki Ahmed
1,* and
Hiwa Abdullah Rasol
2
1
Department of Mathematics, College of Science, University of Sulaimani, Sulaymaniyah 46001, Iraq
2
Department of Mathematics, College of Basic Education, University of Raparin, Raniyah 46012, Iraq
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(8), 1693; https://doi.org/10.3390/sym14081693
Submission received: 5 July 2022 / Revised: 10 August 2022 / Accepted: 11 August 2022 / Published: 15 August 2022
(This article belongs to the Section Mathematics)

Abstract

:
In this article, the numerical solution of the mixed Volterra–Fredholm integro-differential equations of multi-fractional order less than or equal to one in the Caputo sense (V-FIFDEs) under the initial conditions is presented with powerful algorithms. The method is based upon the quadrature rule with the aid of finite difference approximation to Caputo derivative usage collocation points. For treatments, our technique converts the V-FIFDEs into algebraic equations with operational matrices, some of which have the symmetry property, which is simple for evaluating. Furthermore, numerical examples are presented to show the technique’s validity and usefulness as well comparisons with previous results. The majority of programs are performed using MATLAB v. 9.7.

1. Introduction

One of the most significant areas of mathematics that deals with arbitrary order integrals and derivatives is fractional calculus (FC). In recent years, the idea of FC has been successfully researched to focus on a variety of issues in physics, signal processing, engineering, bio-science, and other domains. A few examples of practical uses for fractional calculus are fluid waft, electrical networks, control theory, electromagnetic, facts, optics, capacity, diffusion, and viscoelasticity (see, for instance, [1,2,3,4,5,6,7]). Many practical issues are converted into a series of fractional differential and integral equations via mathematical modeling [8]. Numerical approaches are used to approximate the solution of integro-fractional differential equations since it is challenging to obtain the analytical solution of fractional differential equations.
Then again, integro-fractional differential equations have numerous applications in the mentioned fields, and they pique the interests of many scientists. Volterra and Fredholm integro-differential equations and mixed Volterra–Fredholm integro-differential equations are three prominent types of such equations. In this regard, as in many cases, it is very difficult to find the correct analytical solutions of fractional differential and integral equations. The approximate methods have gained importance to prevent this difficulty. Many methods evaluate the approximate solution of integro-fractional differential equations; see references [9,10,11,12,13,14,15,16,17,18,19]. In [11], the author utilized modified Navot–Simpson’s quadrature for solving second-kind Volterra integral equations of singular type. Hamasalih et al. [14] proposed a technique based on quadrature methods for numerical treatment of the most general VI-FD equations of linear type. Additionally, in [17], the author applied the finite difference method together with homotopy perturbation method for numerically solving the linear VI-D equations, and Smadpour et al. [19] used the iterative method with the aid of the midpoint quadrature formula for solving fuzzy Fredholm integral equations.
The class of equations studied in this paper is the multi-order mixed Volterra–Fredholm integro-differential equations of the fractional orders that lies in the interval ( 0 , 1 ) :
= 0 n 1 p ( x ) a c D x α n u ( x ) + p n ( x ) u ( x ) = f ( x ) + i = 0 m λ i a x a b K i ( t , s ) a c D s β i u ( s ) d s d t ,
with the initial condition u ( a ) = u 0 .
Where x [ a , b ] , p 0 ( x ) = 1 , the fractional orders 0 = α 0 < α 1 < < α n < 1 , and 0 = β 0 < β 1 < < β m < 1 and put μ = max { α n , β m } = 1 . Furthermore, the functions f ( x ) , p ( x ) C [ a , b ] , R , K i C ( S , R ) , the  S = { ( t , s ) : a s t b } denotes the continuous function for all i = 0 , 1 , , m , and  u ( x ) is the unknown function to be determined of Equation (1). The  λ i are scalar parameters. The  a c D x γ indicates the γ -Caputo fractional derivative of u ( x ) on [ a , b ] and γ = { α n , β i } [ 0 , 1 ) for all = 0 , 1 , , n and i = 0 , 1 , , m .
The paper is organized as follows: Section 2 presents the necessary definitions and basic preliminaries of the FC and Lagrange interpolation with basic three techniques of quadrature rules. Section 3 is devoted to the formulation of Newton–Cots formulas: Trapezoidal, Simpson’s 1/3 and Midterm techniques for solving multi-order linear mixed IFDEs of V-F type with variable coefficients. Our results are illustrated throughout the examples in Section 4. Finally, Section 5 includes a discussion for this method.
This paper aims to apply three closed Newton–Cotes techniques with the aid of the finite difference approximation for Caputo derivative terms depending on collocation points converting to an algebraic system. It also intends to evaluate the approximate solution for the multi-order linear mixed V-FIFDEs.

2. Preliminaries

In this section, definitions of the fractional Riemann–Liouville (R-L) integrals and Caputo derivative are presented as well as their relationship and some properties. In addition, we stated some basic definitions and lemmas that are used in this paper.

2.1. Fractional Operators

Definition 1
([20]). A real valued function u defined on closed bounded interval [ a , b ] = I in the space C γ ( I ) , γ R , if there exists a real number k > γ , such that u ( x ) = ( x a ) k u 0 ( x ) , where u 0 ( x ) C ( I ) , and it is said to be in the space C γ n ( I ) if and only if u ( n ) ( x ) is in the space C γ ( I ) , where n Z + { 0 } .
Definition 2
([1,3]). The Riemann–Liouville (R-L) fractional integral operator, a J x α , of order α 0 of a function u C γ ( I ) , γ > 1 is defined as:
a J x α u ( x ) = { 1 Γ ( α ) a x ( x t ) α 1 u ( t ) d t , if α R + . u ( x ) if α = 0 .
Hence, for all α , β 0 , γ > 1 and μ C μ ( μ 1 ) on the closed interval [ a , b ] , we have
a J x α a J x β u ( x ) = a J x β a J x α u ( x ) = a J x α + β u ( x ) .
and
a J x α ( x a ) γ = Γ ( γ + 1 ) Γ ( γ + α + 1 ) ( x a ) γ + α .
Definition 3
([8,21]). The Caputo fractional derivative operator, a c D x α , of order α 0 of a function u C 1 μ ( I ) , μ = α is defined as:
a c D x α u ( x ) = { μ u ( x ) x μ if α = μ , μ N . 1 Γ ( μ α ) a x ( x t ) μ α 1 μ u ( t ) t μ d t if μ 1 < α < μ . .
For all p , q R and u , ψ C 1 μ [ a , b ] , μ = α , so the linear property of the Caputo derivative is formulated as:
a c D x α p u ( x ) + q ψ ( x ) = p a c D x α u ( x ) + q a c D x α ψ ( x ) , x [ a , b ] .
In addition, the  α -fractional derivative in the Caputo sense vanishes, i.e, a c D x α A = 0 , for any constant A R .
The relationship between the integral operator of R-L ( a J x α ) and the differential operator of Caputo ( a c D x α ) is expressed in the following terms:
a c D x α a J x α u ( x ) = u ( x )
is the left inverse, and not the right,
a J x α a c D x α u ( x ) = k = 0 μ 1 u ( k ) ( x ) k ! ( x a ) k , x > a .
Lemma 1
([22]). The finite difference approximation for the Caputo derivative of order 0 < α 1 at specific points x = x ι + 1 , ι = 0 , 1 , , N 1 , where h = ( b a ) / N , is of the form
a c D x α u ( x ι + 1 ) = h α Γ ( 2 α ) d = 0 ι [ u ( x ι d + 1 ) u ( x ι d ) ] b d α ,
where b d α = ( d + 1 ) 1 α d 1 α .
Lemma 2
([22]). Let α 0 , α N , μ = α and u C 1 μ [ a , b ] . Then, a c D x α u ( x ) is continuous on [ a , b ] , and  a c D x α u ( x ) x = x 0 = 0 , i.e.,
lim x a a c D x α u ( x ) = 0 .
Lemma 3
([8,23]). Let u ( x ) = ( x a ) β , β 0 , α 0 and μ = α . Then, the Caputo derivative of order α of u is formed by:
a c D x α u ( x ) = { Γ ( μ + 1 ) Γ ( μ α + 1 ) ( x a ) β α if β N and β > μ 1 , or β N and β μ 0 if β { 0 , 1 , , μ 1 }

2.2. Lagrange Interpolation

A Lagrange polynomial, that interpolates a function u ( x ) as an algebraic polynomial of degree n 1 , at some selected nodes x 0 , x 1 , , x n is denoted by (see [24])
L ( x ) = ι = 0 n u ( x ι ) ι ( x ) ,
where
ι ( x ) = τ = 0 τ ι n ( x x τ ) ( x ι x τ )
The Lagrange polynomial L ( x ) of degree n 1 and u ( x ) are approximately equivalent on the interval [ x 0 , x n ] , i.e.,
L ( x ) u ( x ) , x [ x 0 , x n ] .
If we fix x 0 = x ι 1 = 0 , x 1 = x ι = 1 2 , and x 2 = x ι + 1 = 1 for n = 2 on the interval [ 0 , 1 ] , the Lagrange polynomial will obtain the form
u ( x ) L ( x ) = ( 2 x 2 3 x + 1 ) u ( x ι 1 ) 4 ( x 2 x ) u ( x ι ) + ( 2 x 2 x ) u ( x ι + 1 ) .
then for x ι 1 2 = 1 4 , we have
u ( x ι 1 2 ) = 3 8 u ( x ι 1 ) + 3 4 u ( x ι ) 1 8 u ( x ι + 1 ) .

2.3. Closed Newton–Cots Formulas

2.3.1. Trapezoidal Rule

The trapezoidal rule is the most commonly used simple technique for the numerical integration of definite integrals. The formula is based on the linear interpolation of the positive valued function u ( x ) between two adjacent points x ι and x ι + 1 . So, if an interval [ a , b ] is divided into N equal sub-intervals with the length h = ( b a ) / N , it yields that there is N + 1 equally spaced nodes x ι = a + ι h for each ι = 0 , 1 , , N . Interpolating over all nodes x ι , the summation of the integrated interpolations can be collected in a formula as: [25]
a b u ( t ) d t h ι = 0 N ω ι u ( t ι ) ,
where ω 0 = ω N = 1 2 , and  ω ι = 1 for all ι = 1 , 2 , , N 1 .The order of this approximation is O ( h 2 ) .
The above formula can be expanded, for approximating definite double integrals, by applying Formula (13) to the inner and outer integral, respectively. It can be represented as
a b c d u ( t , s ) d s d t h 1 h 2 ι = 0 N 2 τ = 0 N 1 ω ι ω τ u ( t ι , s τ )
where h 1 = ( d c ) / N 1 , h 2 = ( b a ) / N 2 , s τ = c + τ h 1 for all τ = 1 , 2 , , N 1 , t ι = a + ι h 2 for all ι = 1 , 2 , , N 2 1 , ω 0 = ω N 1 = ω N 2 = 1 2 , and  ω τ = ω ι = 1 for all τ = 1 , 2 , , N 1 1 and ι = 1 , 2 , , N 2 1 . The order of this approximation is O ( h 1 2 + h 2 2 ) .

2.3.2. Simpson’s 1/3 Rule

Simpson’s rule is another unweighted Newton–Cotes method, which approximates the definite integrals by integrating the interpolated parabola of the function u ( t ) between points t ι 1 , t ι and t ι + 1 . The composed Simpson’s 1/3 rule for nodes t ι = a + ι h for all ι = 0 , 1 , , N with setting h = ( b a ) / N , and  N Z + such that N 2 , can be written as [25,26,27]
a b u ( t ) d t h ι = 0 N / 2 ω ι N u ( t ι ) ,
where
  • For odd N: ω 0 N = ω N N = 1 3 , ω 2 τ 1 = 4 3 for all τ = 1 , 2 , , N / 2 , and  ω 2 τ N = 2 3 for all τ = 2 , 4 , , ( N 2 ) / 2 .
  • For even N: ω 0 N = 1 3 , ω N 1 N = 5 3 , ω 2 τ 1 N = 4 3 for all τ = 1 , 2 , , ( N 2 ) / 2 , ω 2 τ N = 2 4 for all τ = 1 , 2 , , ( N 3 ) / 2 , and  ω N N = 1 2 .
Approximation (15) is of order O ( h 4 ) . The techniques discussed in (15) can be modified to be used for approximating multiple integrals. Let h 2 = ( b a ) / N 2 , and  h 1 = ( d c ) / N 1 with t ι = a + ι h 2 , s τ = c + τ h 1 , for all ι = 1 , 2 , , N 2 and τ = 1 , 2 , , N 1 , respectively. Then, we have
a b c d u ( t , s ) d s d t h 1 h 2 ι = 0 N 2 τ = 0 N 1 ω ι N 2 ω τ N 1 u ( t ι , s τ ) ,
where ω ι N 2 and ω τ N 1 are defined in (15). Furthermore, the order of this approximation is O ( h 1 4 + h 2 4 ) .

2.3.3. Midpoint Rule

The midpoint is another Newton–Cotes formula in which one or both integration endpoints are omitted from the interpolation and the integration node points. It is derived from interpolating the integrated u ( t ) by the constant u ( ( t ι + t ι + 1 ) / 2 ) where t ι and t ι + 1 are two adjacent nodes. However, for  the composite formula, we define h = ( b a ) / N , t ι 1 2 = a + ( ι 1 2 ) h for all ι = 1 , 2 , , N , so, it takes the form [25,28]
a b u ( t ) d t h ι = 1 N u ( t ι 1 2 )
The above techniques can be modified to use for approximating definite double integrals by setting h 2 = ( b a ) / N 2 , t ι 1 2 = a + ( ι 1 2 ) h 2 for all ι = 1 , 2 , , N 2 and h 1 = ( d c ) / N 1 , s τ 1 2 = c + ( τ 1 2 ) h 1 for all τ = 1 , 2 , , N 1 . Then
a b c d u ( t , s ) d s d t h 1 h 2 ι = 1 N 2 τ = 1 N 1 u ( t ι 1 2 , s τ 1 2 ) .

3. The Solution Method

By applying each proposed method on (1) for fractional order ( 0 , 1 ) , we can obtain a system of linear algebraic equations that can be represented in a matrix form
[ L R ] U = F
then, solving the matrices Equation (19) for U obtains the approximated solutions at each selected node. All matrices are defined at concerned methods.
Throughout this section, we will use the following notations:
K ρ , q i = K i ( t ρ , s q ) , u r = u ( x r ) , p r = p ( x r )

3.1. A Numerical Method Based on the Trapezoidal Rule

By setting x = x r for all r = 0 , 1 , , N , and then applying trapezoidal Formula (14) to the integral terms of (1) at nodes x ρ = t ρ = a + ρ h , and  s q = a + q h for all ρ = 0 , 1 , , r and q = 0 , 1 , , N , respectively, with setting h = ( x r a ) / r = ( b a ) / N , and then using finite difference approximation for fractional order ( 0 , 1 ) (Lemma 1) together with Lemma 2, we can write (1) as follows:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 4 λ 0 K 0 : r ¯ , 0 0 u 0 + 2 q = 1 N 1 K 0 : r ¯ , q 0 u q + K 0 : r ¯ , N 0 u N + h 2 4 i = 1 m λ i A i β 2 q = 1 N 1 K 0 : r ¯ , q i d = 0 q 1 [ u q d u q d 1 ] b d β i + K 0 : r ¯ , N i d = 0 N 1 [ u N d u N d 1 ] b d β i
where
K 0 : r ¯ , τ i = K 0 , τ i + 2 ρ = 1 r 1 K ρ , τ i + K r , τ i , for all τ = 1 , 2 , , N .
In addition, we assumed that
A n α ( ) = h α n Γ ( 2 α n )
and
A i β = h β i Γ ( 2 β i )
Simply manipulating (20) for all r = 1 , 2 , , N yields a system of N linear algebraic equations of u q s, where u q s are the approximated solutions of the system at x = x q for all q = 1 , 2 , , N . Now, we can write the results in matrix form (19),  where L = [ L r , j ] N × N , R = [ R r , j ] N × N , U = [ u 1 , u 2 , , u N ] T and F = [ F r ] N × 1 . The structure of each matrix used in this technique is as follows:
Firstly, we assume that
B j σ = b j σ b j 1 σ , for all j = 0 , 1 , , N , and σ = { α n , β i } .
with B 0 σ = b 0 σ = 1 and B j σ = b j σ = 0 .
Then, we have
L r , j = 0 if r < j = 0 n 1 p r A n α ( ) + p r n if r = j = 0 n 1 p r A n α ( ) B r j α n , if r > j
In addition, the elements of matrix R are defined as follows:
R r , j = h 2 4 λ 0 ω j N ρ = 0 r ω ρ r K ρ , j 0 + i = 1 m λ i A i β d = j N ρ = 0 r ω d N ω ρ r K ρ , d i B d j β i ,
where
ω τ σ = { 1 if τ = 0 , σ 2 o . w .
Moreover, U and F are column matrices of order N × 1 and are defined as follows:
U = [ u 1 , u 2 , , u N ] T , for all r = 1 , , N .
and
F r = f r + = 0 n 1 p r A n α ( ) b r 1 α n + h 2 4 λ 0 ρ = 0 r ω ρ r K ρ , 0 0 i = 1 m λ i A i β d = 1 N ρ = 0 r ω d N ω ρ K ρ , d i b d 1 β i u 0
where u 0 = u ( a ) and ω τ σ = 1 if τ = 0 , σ 2 o . w . .

3.2. The Algorithm (TV-FIFDE)

The approximate solution for the mixed Volterra–Fredholm integro-differential equations of the fractional order lies in ( 0 , 1 ) , which is based on the trapezoidal rule and can be summarized by the following Algorithm 1:
Algorithm 1: TV-FIFDE
Input: N Z + , N, a, b and the initial condition u 0 .
Set: h = b a N .
for r 1 to N do
    for  j 0 to N do
         x r = a + r h .
        for  0 to n do
            A n α ( ) h α n Γ ( 2 α n )
            B j α n b j b j 1
           if  r < j  then
                C r , j 0
           else if  r = j  then
                C r , j p r A n α ( ) .
           else
                C r , j p r A n α ( ) B r j α n
           end if
        end for
        if  r < j  then
            L r , j 0
        else
            L r , j sum ( C r , j ) + p r n .
        end if
        for  i 0 to m do
            A i β h β i Γ ( 2 β i )
            B j α n b j b j 1
           for  d j to N do
                S 0 0 , C 1 0
               for  ρ 0 to r do
                   if i = 0 then
                        C 0 C 0 + λ 0 ω j N ω ρ r K ρ , j 0 .
                   else
                        C 1 C 1 + λ i A i β ω d N ω ρ r K ρ , d i B d j β i
                   end if
               end for
                C d , 1 0 S 0 .
                C d , 1 i S 1 .
           end for
        end for
         R r , j = h 2 4 sum ( C d , 1 0 + C d , 1 i )
    end for
    for  0 to n 1  do
         E r p r A n α ( ) b r 1 α n
         J 0 0 , J 1 0
    end for
    for  d 1 to N do
        for  ρ 0 to r do
           if i = 0 then
                J 0 J 0 + λ 0 ω ρ r K ρ , 0 0 .
           else
                J 1 J 1 + λ i A i β ω d N ω ρ r K ρ , d i b d j β i
           end if
        end for
         D d , 1 0 J 0
         D d , 1 i J 0
    end for
     F r f r + ( E r ) + h 2 4 sum ( D d , 1 0 D d , 1 i )
end for
Output: U [ L R ] F

3.3. A Numerical Method Based on the Simpson’s 1/3 Rule

By setting x = x r and applying Simpson’s 1/3 Formula (16) to the integral terms of (1) at nodes x ρ = t ρ = a + ρ h , and  s q = a + q h for all ρ = 0 , 1 , , r and q = 0 , 1 , , N , respectively, with setting h = ( x r a ) / r = ( b a ) / N , then, using the finite difference (Lemma 1) together with Lemma 2 and depending on the value of N, we can classify the method into the following two cases:
Case I: N is even
If r is even, then (1) becomes:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 9 i = 1 m λ i A i β q = 1 N / 2 K 0 : r ¯ , 2 q i d = 0 2 q 1 [ u 2 q d u 2 q d 1 ] b d β i + 4 q = 1 N / 2 K 0 : r ¯ , 2 q 1 i d = 0 2 q 2 [ u 2 q d 1 u 2 q d 2 ] b d β i + q = 2 N / 2 K 0 : r ¯ , 2 q 2 i d = 0 2 q 3 [ u 2 q d 2 u 2 q d 3 ] b d β i + λ 0 h 2 9 q = 1 N / 2 ( K 0 : r ¯ , 2 q 0 u 2 q + 4 K 0 : r ¯ , 2 q 1 0 u 2 q 1 + K 0 : r ¯ , 2 q 2 0 u 2 q 2 ,
where
K 0 : r ¯ , q i = ρ = 0 r / 2 K 2 ρ , q i + 4 K 2 ρ 1 , q i + K 2 ρ 2 , q i , for all r = 2 , 4 , , N and i = 0 , 1 , , m .
If r is odd, then Equation (1) takes the following form:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 9 i = 1 m λ i A i β q = 1 N / 2 K 0 : r 1 ¯ , 2 q i + 3 2 K r 1 : r ¯ , 2 q i d = 0 2 q 1 [ u 2 q d u 2 q d 1 ] b d β i + 4 q = 1 N / 2 K 0 : r 1 ¯ , 2 q 1 i + 3 2 K r 1 : r ¯ , 2 q 1 i d = 0 2 q 2 [ u 2 q d 1 u 2 q d 2 ] b d β i + q = 2 N / 2 K 0 : r 1 ¯ , 2 q 2 i + 3 2 K r 1 : r ¯ , 2 q 2 i d = 0 2 q 3 [ u 2 q d 2 u 2 q d 3 ] b d β i + λ 0 h 2 9 q = 1 N / 2 K 0 : r 1 ¯ , 2 q 0 + 3 2 K r 1 : r ¯ , 2 q 0 u 2 q + 4 K 0 : r 1 ¯ , 2 q 1 0 + 3 2 K r 1 : r ¯ , 2 q 1 0 u 2 q 1 + K 0 : r 1 ¯ , 2 q 2 0 + 3 2 K r 1 : r ¯ , 2 q 2 0 u 2 q 2 ,
where
K 0 : r 1 ¯ , q i = ρ = 0 ( r 1 ) / 2 K 2 ρ , q i + 4 K 2 ρ 1 , q i + K 2 ρ 2 , q i , for all r = 3 , 5 , , N 1 , and i = 0 , 1 , , m . K r 1 : r ¯ , q i = K r 1 , q + K r , q , for all r = 1 , 3 , , N 1 and i = 0 , 1 , , m
We have to repeat the process in Equations (28) and (29) until r reaches N, and then, we try to combine both equations together to construct a system of linear equations. Then, it can be written in matrix form (19), where L = [ L r , j ] N × N and U = [ U r ] N × 1 are the same as (24) and (26) for all r , j = 1 , 2 , , N . Recall that N is even in this case. In addition, R = [ R r , j ] N × N is a square matrix and its construction was illustrated for all r , j = 1 , 2 , , N as follows:
R r , j = h 2 9 λ 0 ρ = 0 r ϖ j ω ρ r K ρ , j 0 + i = 1 m λ i A i β d = j N ρ = 0 r ϖ d ω ρ r K ρ , d i B d j β i .
F r = f r + ρ = 0 n 1 p r A n α ( ) b r 1 α n + h 2 9 λ 0 ρ = 0 r ω ρ r K ρ , 0 0 i = 1 m λ i A i β d = 1 N ρ = 0 r ϖ d ω ρ r K ρ , d i b d 1 β i u 0 .
where
u 0 = u ( a ) , ϖ d = 1 if d = N 2 if d N and d is even 4 if d N and d is odd
ω 0 1 = ω 1 1 = 3 2 , ω ρ 2 τ = 1 if ρ = 0 , 2 τ 2 if ρ 2 τ and ρ is even , τ = 1 , 2 , , N 2 , 4 if ρ 2 τ and ρ is odd ω ρ 2 τ + 1 = 1 if ρ = 0 2 if ρ 2 τ , 2 τ + 1 and ρ is even 4 if ρ 2 τ + 1 and ρ is odd , τ = 1 , 2 , , N 2 2 , 3 2 if ρ = 2 τ 3 2 if ρ = 2 τ + 1
Case II: N is odd
For even r, Equation (1) will obtain the following form:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 9 i = 1 m λ i A i β { q = 1 ( N 1 ) / 2 K 0 : r ¯ , 2 q i d = 0 2 q 1 [ u 2 q d u 2 q d 1 ] b d β i + 4 q = 1 ( N 1 ) / 2 K 0 : r ¯ , 2 q 1 i d = 0 2 q 2 [ u 2 q d 1 u 2 q d 2 ] b d β i + q = 2 ( N 1 ) / 2 K 0 : r ¯ , 2 q 2 i d = 0 2 q 3 [ u 2 q d 2 u 2 q d 3 ] b d β i + 3 2 K 0 : r ¯ , N 1 i d = 0 N 2 [ u N d 1 u N d 2 ] b d β i + K 0 : r ¯ , N i d = 0 N 1 [ u N d u N d 1 ] b d β i } + λ 0 h 2 9 { q = 1 ( N 1 ) / 2 ( K 0 : r ¯ , 2 q 0 u 2 q + 4 K 0 : r ¯ , 2 q 1 0 u 2 q 1 + K 0 : r ¯ , 2 q 2 0 u 2 q 2 + 3 2 K 0 : r ¯ , N 1 0 u N 1 + K 0 : r ¯ , N 0 u N } ,
where,
K 0 : r ¯ , q i = ρ = 0 r / 2 K 2 ρ , q i + 4 K 2 ρ 1 , q i + K 2 ρ 2 , q i , r = 2 , 4 , , N .
If r is odd, Equation (1) becomes:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 9 i = 1 m λ i A i β { q = 1 N / 2 K 0 : r 1 ¯ , 2 q i + 3 2 K r 1 : r ¯ , 2 q i d = 0 2 q 1 [ u 2 q d u 2 q d 1 ] b d β i + 4 q = 1 N / 2 K 0 : r 1 ¯ , 2 q 1 i + 3 2 K r 1 : r ¯ , 2 q 1 i d = 0 2 q 2 [ u 2 q d 1 u 2 q d 2 ] b d β i
+ q = 2 N / 2 K 0 : r 1 ¯ , 2 q 2 i + 3 2 K r 1 : r ¯ , 2 q 2 i d = 0 2 q 3 [ u 2 q d 2 u 2 q d 3 ] b d β i + 3 2 [ K 0 : r 1 ¯ , N 1 i + 3 2 K r 1 : r ¯ , N 1 i d = 0 N 2 [ u N d 1 u N d 2 ] b d β i + K 0 : r 1 ¯ , N i + 3 2 K r 1 : r ¯ , N i d = 0 N 1 [ u N d u N d 1 ] b d β i ] } + λ 0 h 2 9 { q = 1 N / 2 K 0 : r 1 ¯ , 2 q 0 + 3 2 K r 1 : r ¯ , 2 q 0 u 2 q + 4 K 0 : r 1 ¯ , 2 q 1 0 + 3 2 K r 1 : r ¯ , 2 q 1 0 u 2 q 1 + K 0 : r 1 ¯ , 2 q 2 0 + 3 2 K r 1 : r ¯ , 2 q 2 0 u 2 q 2 + 3 2 K 0 : r 1 ¯ , N 1 0 + 3 2 K r 1 : r ¯ , N 1 0 u N 1 + K 0 : r 1 ¯ , N 0 + 3 2 K r 1 : r ¯ , N 0 u N } ,
where
K 0 : r 1 ¯ , q i = ρ = 0 ( r 1 ) / 2 K 2 ρ , q i + 4 K 2 ρ 1 , q i + K 2 ρ 2 , q i ,
for all r = 3 , 5 , , N 1 , and  i = 0 , 1 , , m . and K r 1 : r ¯ , q i = K r 1 , q + K r , q , for all r = 1 , 3 , , N 1 , and  i = 0 , 1 , , m .
Again, by repeating the process and combining both Formulas (32) and (33), we will obtain a system of linear equations. Moreover, it can be written in matrix form (19), with the only difference being that in this case, N is odd. Furthermore, the key formulas for constructing R r , j and F r for all r , j = 1 , 2 , N are as follows:
R r , j = h 2 9 λ 0 ρ = 0 r ϖ j ω ρ r K ρ , j 0 + i = 1 m λ i A i β d = j N ρ = 0 r ϖ d ω ρ r K ρ , d i B d j β i .
F r = f r + ρ = 0 n 1 p r A n α ( ) b r 1 α n + h 2 9 λ 0 ρ = 0 r ω ρ r K ρ , 0 0 i = 1 m λ i A i β d = 1 N ρ = 0 r ϖ d ω ρ r K ρ , d i b d 1 β i u 0 .
where
u 0 = u ( a ) , ϖ d = 2 if d N , N 1 and d is even 4 if d N , N 1 and d is odd 3 2 if d = N 1 3 2 if d = N
ω 0 1 = ω 1 1 = 3 2 , ω ρ 2 τ = 1 if ρ = 0 , 2 τ 2 if ρ 2 τ and ρ is even , τ = 1 , 2 , , N 1 2 , 4 if ρ 2 τ and ρ is odd
ω ρ 2 τ + 1 = 1 if ρ = 0 2 if ρ 2 τ , 2 τ + 1 and ρ is even 4 if ρ 2 τ + 1 and ρ is odd , τ = 1 , 2 , , N 1 2 , 3 2 if ρ = 2 τ 3 2 if ρ = 2 τ + 1

3.4. The Algorithm (SV-FIFDE)

The following Algorithm 2 outline the approximate solution for a mixed Volterra–Fredholm integro-differential equations of fractional order in ( 0 , 1 ) based on the Simpson’s 1/3 rule:
Algorithm 2: SV-FIFDE
Input: N Z + , N, a, b and the initial condition u 0 .
Set: h = b a N .
for r 1 to N do
    for  j 0 to N do
         x r = a + r h .
        for  0 to n do
            A n α ( ) h α n Γ ( 2 α n )
            B j α n b j b j 1
           if  r < j  then
                C r , j 0
           else if r = j then
                C r , j p r A n α ( ) .
           else
                C r , j p r A n α ( ) B r j α n
           end if
        end for
        if  r < j  then
            L r , j 0
        else
            L r , j = sum ( C r , j ) + p r n .
        end if
        for  i 0 to m do
            A i β h β i Γ ( 2 β i )
            B j α n b j b j 1
           for  d j to N do
                S 0 0 , C 1 0
               for  ρ 1 to r do
                   if  i = 0  then
                        C 0 C 0 + λ 0 ϖ j ω ρ r K ρ , j 0 .▹ if N is even, use relations in (30), and use relations in (34) otherwise.
                   else
                        C 1 C 1 + λ i A i β ϖ d ω ρ r K ρ , d i B d j β i ▹ if N is even, use relations in (30), and use relations in (34) otherwise.
                   end if
               end for
                C d , 1 0 S 0 .
                C d , 1 i S 1 .
           end for
        end for
         R r , j = h 2 9 sum ( C d , 1 0 + C d , 1 i )
    end for
    for  0 to n 1  do
         E r p r A n α ( ) b r 1 α n
         J 0 0 , J 1 0
    end for
    for  d 1 to N do
        for  ρ 0 to r do
           if i = 0 then
                J 0 J 0 + λ 0 ω ρ r K ρ , 0 0 . ▹ if N is even, use relations in (31), and use relations in (35) otherwise.
           else
                J 1 J 1 + λ i A i β ϖ d ω ρ r K ρ , d i b d j β i ▹ if N is even, use relations in (31), and use relations in (35) otherwise.
           end if
        end for
         D d , 1 0 J 0
         D d , 1 i J 0
    end for
     F r f r + ( E r ) + h 2 9 sum ( D d , 1 0 D d , 1 i )
end for
Output: U [ L R ] F

3.5. A Numerical Method Based on the Midpoint Rule

Set x = x r , then h = ( x r a ) / r = ( b a ) / N , with  x ρ = t ρ = a + ρ h , and  s j = a + q h for all r , j = 1 , , N and ρ = 1 , , r . Then, we can write Equation (1) as
= 0 n 1 p ( x r ) a c D x α n u ( x ) x = x r + p n ( x r ) u ( x r ) = f ( x r ) + i m λ i a x r a s N 1 K i ( t , s ) a c D s β i u ( s ) d s d t + a x r s N 1 s N K i ( t , s ) a c D s β i u ( s ) d s d t ,
since
a x r a s N 1 K i ( t , s ) a c D s β i u ( s ) d s d t = a x r 1 a s N 1 K i ( t , s ) a c D s β i u ( s ) d s d t + x r 1 x r a s N 1 K i ( t , s ) a c D s β i u ( s ) d s d t ,
and
a x r s N 1 s N K i ( t , s ) a c D s β i u ( s ) d s d t = a x r 1 s N 1 s N K i ( t , s ) a c D s β i u ( s ) d s d t + x r 1 x r s N 1 s N K i ( t , s ) a c D s β i u ( s ) d s d t ,
The main idea is to use midpoint rule (18) for approximating the terms that lie in the intervals [ a , x r 1 ] and [ a , s N 1 ] and trapezoidal rule (14) for the terms in [ x r 1 , x r ] and [ s N 1 , s N ] . Furthermore, by using relation (12) on the interval [ s q 1 , s q + 1 ] , all the terms of u q 1 2 and [ a c D s β i u ( s ) ] s = s q 1 2 can be approximated as
a c D s β i u ( s ) s = s q 1 2 = 3 8 a c D s β i u ( s ) s = s q 1 + 1 4 a c D s β i u ( s ) s = s q 1 8 a c D s β i u ( s ) s = s q + 1
u q 1 2 = 3 8 u q 1 2 + 1 4 u q 1 8 u q + 1
It allows us to compute initial and final terms of u q for all q = 1 , 2 , , N . In addition, it is helpful for obtaining better results due to repeating the process in each step.
Finally, from relations (39)–(42) and using Lemma 1 together with Lemma 2, the main Equation (38) can be written as follows:
= 0 n 1 p r A n α ( ) d = 0 r 1 [ u r d u r d 1 ] b d α n + p r n u r = f r + h 2 i = 1 m λ i A i β { 3 8 q = 2 N 1 K 1 2 : r ¯ , q 1 2 i d = 0 q 2 [ u q d 1 u q d 2 ] b d β i + 3 4 q = 1 N 1 K 1 2 : r ¯ , q 1 2 i d = 0 q 1 [ u q d u q d 1 ] b d β i 1 8 q = 1 N 1 K 1 2 : r ¯ , q 1 2 i d = 0 q [ u q d + 1 u q d ] b d β i + 1 2 K 1 2 : r ¯ , N 1 i d = 0 N 2 [ u N d 1 u N d 2 ] b d β i + K 1 2 : r ¯ , N i d = 0 N 1 [ u N d u N d 1 ] b d β i } + h 2 λ 0 { q = 1 N 1 K 1 2 : r ¯ , q 1 2 0 3 8 u q 1 + 3 4 u q 1 8 u q + 1 + 1 2 K 1 2 : r ¯ , N 1 0 u N 1 + K 1 2 : r ¯ , N 0 u N }
where
K 1 2 : r ¯ , q 1 2 i = ρ = 1 r 1 K ρ 1 2 , q 1 2 i + 1 2 K r 1 , q 1 2 i + K r , q 1 2 i , i = 0 , 1 , , m .
After some simple manipulation and repeating (43) for all r from 1 to N, it allows obtaining a system of N linear algebraic equations of u q values for  q = 1 , 2 , . . , N . Now, we are ready to write the obtained in the matrix form (19).
Recall that L = [ L r , j ] N × N , R = [ R r , j ] N × N , U = [ U r ] N × 1 and F = [ F r ] N × 1 . Each matrix’s structure for this method is defined as follows:
R r , j = h 2 λ 0 3 8 K 1 2 : r ¯ , j + 1 2 0 + 3 4 K 1 2 : r ¯ , j 1 2 0 1 8 K 1 2 : r ¯ , j 3 2 0 + h 2 i = 1 m λ i A i β i [ d = 1 N 1 K 1 2 : r ¯ , d 1 2 i 3 8 B d j 1 β i + 3 4 B d j β i 1 8 B d j + 1 β i + 1 2 K 1 2 : r ¯ , N 1 i B N j 1 β i + K 1 2 : r ¯ , N i B N j β i ] ,
R r , N 1 = h 2 λ 0 3 4 K 1 2 : r ¯ , N 3 2 0 ω N 1 1 8 K 1 2 : r ¯ , N 5 2 0 + 1 2 K 1 2 : r ¯ , N 1 0 + h 2 i = 1 m λ i A i β i d = N 1 N K 1 2 : r ¯ d 3 2 i 3 4 B d N β i ϖ d 2 1 8 B d N + 1 β i + 1 2 K 1 2 : r ¯ , N 1 i + K 1 2 : r ¯ , N i B 1 β i ,
R r , N = h 2 λ 0 ω N 1 8 K 1 2 : r ¯ , N 3 2 0 + 1 2 K 1 2 : r ¯ , N 0 + h 2 i = 1 m λ i A i β i 1 8 K 1 2 : r ¯ , N 3 2 i + 1 2 K 1 2 : r ¯ , N i .
where ω d = 0 if d = 1 1 o . w , ϖ d = 0 if d < 0 1 o . w , B d β i = b d β i = 0 and B 0 β i = b 0 β i = 1 for all d 1 and β i R + . Furthermore, F r is defined as follows:
F r = f r + h 2 { λ 0 3 8 K 1 2 : r ¯ , 1 2 0 i = 1 m λ i A i β [ d = 1 N 1 K 1 2 : r ¯ , d 1 2 i 3 8 b d 2 β i + 3 4 b d 1 β i 1 8 b d β i + 1 2 K 1 2 : r ¯ , N 1 b N 2 β i + K 1 2 : r ¯ , N b N 1 β i ] } u 0 .
L and U will be obtained in the same way before being described in the trapezoidal and Simpson’s 1/3 methods.

3.6. The Algorithm (MV-FIFDE)

We can summarize the Midpoint rule technique for solving the mixed V-FIDEs of the α -fractional orders, for  α ( 0 , 1 ) , at the following Algorithm 3:
Algorithm 3: MV-FIFDE
Input: N Z + , N, a, b and the initial condition u 0 .
Set: h = b a N .
for r 1 to N do
    for  j 0 to N do
         x r = a + r h .
        for  0 to n do
            A n α ( ) h α n Γ ( 2 α n )
            B j α n b j b j 1
           if  r < j  then
                C r , j = 0
           else if  r = j  then
                C r , j p r A n α ( ) .
           else
                C r , j p r A n α ( ) B r j α n
           end if
        end for
        if  r < j  then
            L r , j 0
        else
            L r , j sum ( C r , j ) + p r n .
        end if
        for  i 0 to m do
            A i β h β i Γ ( 2 β i )
            B j α n b j b j 1
           if  i = 0 & j N , N 1  then
                C 0 λ 0 3 8 K 1 2 : r ¯ , j + 1 2 0 + 3 4 K 1 2 : r ¯ , j 1 2 0 1 8 K 1 2 : r ¯ , j 3 2 0
           else if  i = 0 & j = N 1  then
                C 0 λ 0 ω N 1 8 K 1 2 : r ¯ , N 3 2 0 + 1 2 K 1 2 : r ¯ , N 0
           else if  j = N  then
                C 0 0 ω N 1 8 K 1 2 : r ¯ , N 3 2 0 + 1 2 K 1 2 : r ¯ , N 0
           end if
           for  d j to N do
               if  i 0 & j N , N 1  then
                    C 1 C 1 + λ i A i β i [ K 1 2 : r ¯ , d 1 2 i 3 8 B d j 1 β i + 3 4 B d j β i 1 8 B d j + 1 β i
                    + 1 2 K 1 2 : r ¯ , N 1 i B N j 1 β i + K 1 2 : r ¯ , N i B N j β i ]
               else if  i 0 & j = N 1  then
                    C 1 λ i A i β i [ d = N 1 N K 1 2 : r ¯ d 3 2 i 3 4 B d N β i ϖ d 2 1 8 B d N + 1 β i
                    + 1 2 K 1 2 : r ¯ , N 1 i + K 1 2 : r ¯ , N i B 1 β i ]
               else if  j = N  then
                    C 1 A i β i 1 8 K 1 2 : r ¯ , N 3 2 i + 1 2 K 1 2 : r ¯ , N i
               end if
                J 0 J 0 + λ i A i β [ K 1 2 : r ¯ , d 1 2 i 3 8 b d 2 β i + 3 4 b d 1 β i 1 8 b d β i
                + 1 2 K 1 2 : r ¯ , N 1 b N 2 β i + K 1 2 : r ¯ , N b N 1 β i ] }
           end for
        end for
         R r , j = h 2 sum ( C 0 + C 1 )
    end for
     F r f r + h 2 λ 0 3 8 K 1 2 : r ¯ , 1 2 0 J 0
end for
Output: U [ L R ] F

4. Numerical Examples

To highlight the precision of our work, we provide various examples where the exact solution already exists. For this purpose, we have written different MATLAB codes for each method (Appendix A), which generally covers the whole procedure of solution methods. They all have been performed on a computer using MATLAB v.9.7.
Example 1.
Consider a mixed Volterra–Fredholm integro-fractional differential as:
0 c D x 0.5 u ( x ) + ( x 1 ) 2 cosh ( x ) 0 c D x 0.34 u ( x ) + u ( x ) = f ( x ) + 0 x 0 1 t 0 c D s 0.2 u ( s ) d s d t ,
where
f ( x ) = x 0.66 ( x 1 ) 2 cosh ( x ) Γ ( 1.66 ) x 2 2 Γ ( 2.8 ) + x + x Γ ( 1.5 ) + 1 ,
with initial condition u ( 0 ) = 1 , where u ( x ) = x + 1 is the exact solution.
For this example, we set h = 0.1 ( N = 10 ) . Applying Algorithms 1, 2 and 3 and then running MATLAB program (see Appendix A.1 and Appendix A.2: Listing A5) for each method gives the approximate solution to the above equation. The results are presented in Table 1 in terms of running time and least square error (L.S.Error). The L.S.Errors for N = 10 ,   20 ,   100 and N = 200 are presented in Table 2. Additionally, Figure 1 depicts a comparison between the least square error of approximate solutions for N = 10 on [ 0 , 1 ] .
Example 2.
Consider the MV-FIFDE equation on the closed bounded interval [ 1 , 2 ] as:
0 c D x 0.9 u ( x ) + e x 0 c D x 0.8 u ( x ) + ln ( x ) 0 c D x 0.7 u ( x ) x 2 0 c D x 0.6 u ( x ) + sinh ( x ) u ( x ) = f ( x ) + 1 x 1 2 sin ( t + s ) u ( s ) d s d t + 1 x 1 2 e t s 1 0 c D s 0.5 u ( s ) d s d t ,
where
f ( x ) = 66 x 1 23 10 ln x 5 2 e x e 7 + 64 x 1 11 5 e x 5 68 x 1 12 5 x 2 5 + 62 x 1 21 10 5 5 sin x + 2 3 cos x + 2 + 6 cos x + 1 + x 1 3 sinh x + 5 sin 3 + 3 cos 3 6 cos 2
with initial condition, u ( 1 ) = 0 , where u ( x ) = x 3 3 x 2 + 3 x 1 , is the exact solution.
For the above example, numerical results are obtained by following the Algorithms 1, 2 and 3 and running the MATLAB codes (Appendix A.1 and Appendix A.2: Listing A5) that implement the steps of the algorithms. Table 3 shows the results at some selected nodes for N = 200 in terms of running time and L.S.Error.
The L.S.Errors are compared in Table 4 for N = 10 , 20 , 100 and 200. Additionally, Figure 2 depicts a comparison between the least square error of approximate solutions.
Example 3.
Consider the Mixed V-FIFD equation on [ 0 , 1 ] as:
0 c D x 0.12 u ( x ) + cosh ( x ) 0 c D x 0.09 u ( x ) + x cos ( π x 4 ) u ( x ) = f ( x ) + 0 x 0 2 t 2 0 c D s 0.88 u ( s ) d s d t ,
where
f ( x ) = 2 x cos ( π x 4 ) x cos ( π x 4 ) e x ι = 0 η x ι + 0.88 Γ ( ι + 1.88 ) + cosh ( x ) x ι + 0.91 Γ ( ι + 1.91 ) + x 3 Γ ( ι + 2.12 ) .
with initial condition u ( 0 ) = 1 , where u ( x ) = 2 e x is the exact solution.
For the above example, numerical results are obtained by following the algorithms TV-FIFDE, SV-FIFDE and MV-FIFDE and running the MATLAB codes that implement the steps of the Algorithms 1, 2 and 3 (see Appendix A.1 and Appendix A.2). Table 5 shows the results at some selected nodes for N = 20 , in terms of running time and L.S.Error.
For N = 20 , 100 , 200 and 400 , the L.S.Errors are compared in Table 6. Additionally, Figure 3 depicts a comparison between the least square error of approximate solutions.

5. Conclusions

We presented three numerical approaches for solving Volterra–Fredholm integro-fractional equations, and each method was described in detail. In addition, a good algorithm was created for summarizing each technique’s process. In addition, three examples are provided to test the validity and applicability of the techniques, and they are compared in terms of total least square errors and running times. In practice, we conclude with the following remarks:
  • The Midpoint approach is almost more accurate than Simpson’s 1/3 and trapezoidal.
  • All three techniques usually provide similar outcomes.
  • The Midpoint and Simpson’s 1/3 approaches are similar, and in some cases, the Simpson’s 1/3 method gives accurate results.
  • In a problem that contains Mittag–Leffler terms (see Example 3), the accuracy of the results depends on the choice N being sufficiently large (small step size h) and the number of Mittag–Leffler terms ( η ). Due to saving time, we approximated the Mittag–Leffler into 21 terms instead of infinite terms.
  • As N and η are sufficiently large, the error rate decreases, and the answer approaches the exact solution.
Feature Works: The mixed Volterra–Fredholm integro-fractional differential equation studied in this article is a novel equation. We recommend using generalized Bernstein polynomials with discrete residual weighted techniques such as collocation, subdomain, moment, and least square approaches with certain adjustments. We believe they offer decent results. Additionally, it could be treated numerically by applying all quadrature techniques.

Author Contributions

Conceptualization, S.S.A. and H.A.R.; methodology, S.S.A.; software, H.A.R.; validation, S.S.A. and H.A.R.; formal analysis, S.S.A.; investigation, H.A.R.; resources, S.S.A.; data curation, H.A.R.; writing—original draft preparation, S.S.A.; writing—review and editing, H.A.R.; visualization, H.A.R.; supervision, S.S.A.; project administration, S.S.A.; funding acquisition, S.S.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Ethics Committees of Mathematics Department, College of Basic Educations, University of Raparin, Ranya 46012, Kurdistan Region, Iraq; and Mathematics Department, College of Science, University of Sulaimani, Sulaymaniyah 46001, Kurdistan Region, Iraq.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data used during the study are available fromthe corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
FCFractional Calculus
VI-DEVolterra Integro-Differential Equations
V-FIFDEMixed Volterra–Fredholm Integro-Fractional Differential Equations
TV-FIFDETrapezoidal Algorithm of Mixed V-FIFDEs
SV-FIFDESimpson’s 1/3 Algorithm of Mixed V-FIFDEs
MV-FIFDEMidpoint Algorithm of Mixed V-FIFDEs

Appendix A. The MATLAB Codes

Appendix A.1. The Main MATLAB Codes

Listing A1. Main MATLAB Code of Trapezoidal Rule.
Symmetry 14 01693 i001a Symmetry 14 01693 i001b Symmetry 14 01693 i001c Symmetry 14 01693 i001d
Listing A2. Main MATLAB Code of Simpson’s 1/3 Rule.
Symmetry 14 01693 i002a Symmetry 14 01693 i002b Symmetry 14 01693 i002c Symmetry 14 01693 i002d Symmetry 14 01693 i002e
Listing A3. Main MATLAB Code of Midpiont Rule.
Symmetry 14 01693 i003a Symmetry 14 01693 i003b Symmetry 14 01693 i003c Symmetry 14 01693 i003d

Appendix A.2. The MATLAB Sub-Code and Functions

Listing A4. The Mittag–Leffler Function.
Symmetry 14 01693 i004
Listing A5. The Kernels Substituter.
Symmetry 14 01693 i005

References

  1. Miller, K.S.; Ross, B. An Introduction to the Fractional Calculus and Fractional Differential Equations; Wiley: New York, NY, USA, 1993. [Google Scholar]
  2. Panwar, V.S.; Uduman, P.S.; Gomez-Aguilar, J.F. Mathematical modeling of coronavirus disease covid-19 dynamics using cf and abc non-singular fractional derivatives. Chaos Solitons Fractals 2011, 145, 110757. [Google Scholar] [CrossRef]
  3. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: San Diego, CA, USA, 1998. [Google Scholar]
  4. Ross, B. The development of fractional calculus 1695–1900. Hist. Math. 1977, 4, 75–89. [Google Scholar] [CrossRef]
  5. Sun, H.G.; Zhang, Y.; Baleanu, D.; Chen, W.; Chen, Y.Q. A new collection of real world applications of fractional calculus in science and engineering. Commun. Nonlinear Sci. Numer. Simul. 2018, 64, 213–231. [Google Scholar] [CrossRef]
  6. Tuan, N.H.; Mohammadi, H.; Rezapour, S. A mathematical model for covid-19 transmission by using the caputo fractional derivative. Chaos Solitons Fractals 2020, 140, 110107. [Google Scholar] [CrossRef]
  7. Valério, D.; Machado, T.J.; Kiryakova, V. Some pioneers of the applications of fractional calculus. Fract. Calc. Appl. Anal. 2014, 17, 552–578. [Google Scholar] [CrossRef]
  8. Roohollahi, A.; Ghazanfari, B.; Akhavan, S. Numerical solution of the mixed volterra-fredholm integro-differential multi-term equations of fractional order. J. Comput. Appl. Math. 2008, 376, 112828. [Google Scholar] [CrossRef]
  9. Ahmed, S.S.; MohammedFaeq, S.J. Bessel Collocation Method for Solving Fredholm–Volterra Integro-Fractional Differential Equations of Multi-High Order in the Caputo Sense. Symmetry 2021, 13, 2354. [Google Scholar] [CrossRef]
  10. Al-Saif, N.S.M.; Ameen, S.A. Numerical solution of mixed volterra—Fredholm integral equation using the collocation method. Baghdad Sci. J. 2020, 17, 0849. [Google Scholar]
  11. Araghi, M.A.F.; Hamed, D.K. Numerical solution of the second kind singular Volterra integral equations by modified Navot-Simpson’s quadrature. Int. J. Open Probl. Compt. Math. 2008, 1, 201–2013. [Google Scholar]
  12. Filiz, A. A fourth-order robust numerical method for integro-differential equations. Asian J. Fuzzy Appl. Math. 2013, 1, 21–33. [Google Scholar]
  13. Hamasalih, S.A.; Ahmed, M.R.; Ahmed, S.S. Solution Techniques Based on Adomian and Modified Adomian Decomposition for Nonlinear Integro-Fractional Differential Equations of the Volterra-Hammerstein Type. J. Univ. Babylon Pure Appl. Sci. 2020, 28, 194–216. [Google Scholar]
  14. Hamasalih, S.A.; Ahmed, S.S. Numerical treatment of the most general linear volterra integro-fractional differential equations with caputo derivatives by quadrature methods. J. Math. Comput. Sci. 2012, 2, 1293–1311. [Google Scholar]
  15. Hasan, P.M.A.; Sulaiman, N.A. Numerical treatment of mixed volterra-fredholm integral equations using trigonometric functions and laguerre polynomials. Zanco J. Pure Appl. Sci. 2018, 30, 97–106. [Google Scholar]
  16. Kamoh, N.M.; Gyemang, D.G.; Soomiyol, M.C. Comparing the efficiency of Simpson’s 1/3 and 3/8 rules for the numerical solution of first order Volterra integro-differential equations. Int. J. Math. Comput. Sci. 2019, 13, 136–140. [Google Scholar]
  17. Raftari, B. Numerical solutions of the linear volterra integro-differential equations: Homotopy perturbation method and finite difference method. World Appl. Sci. J. 2010, 9, 7–12. [Google Scholar]
  18. Rihan, F.A.; Doha, E.H.; Hassan, M.I.; Kamel, N.M. Numerical treatments for volterra delay integro-differential equations. Comput. Methods Appl. Math. 2009, 9, 292–318. [Google Scholar] [CrossRef]
  19. Samadpour, V.; Mahaleh, M.; Ezzati, R. Numerical solution of linear fuzzy fredholm integral equations of second kind using iterative method and midpoint quadrature formula. J. Intell. Fuzzy Syst. 2017, 33, 1293–1302. [Google Scholar] [CrossRef]
  20. Tunç, O.; Atan, Ö.; Tunç, C.; Yao, J.C. Qualitative analyses of integro-fractional differential equations with Caputo derivatives and retardations via the Lyapunov–Razumikhin method. Axioms 2021, 10, 58. [Google Scholar] [CrossRef]
  21. Kilbas, A.A.; Srivastava, H.M.; Juan, J.T. Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  22. Ahmed, S. On System of Linear Volterra Integro-Fractional Differential Equations. Ph.D. Thesis, University of Sulaimani, Sulaymaniyah, Iraq, July 2009. [Google Scholar]
  23. Ahmed, S.S.; Amin, M.B. Solving Linear Volterra Integro-Fractional Differential Equations in Caputo Sense with Constant Multi-Time Retarded Delay by Laplace Transform. Zanco J. Pure Appl. Sci. 2019, 31, 80–89. [Google Scholar]
  24. Dzyadyk, V.K.; Shevchuk, I.A. Theory of Uniform Approximation of Functions by Polynomials; Walter de Gruyter: Berlin, Germany, 2008. [Google Scholar]
  25. Burden, R.L.; Faires, J.D.; Burden, A.M. Numerical Analysis; Cengage Learning: Boston, MA, USA, 2015. [Google Scholar]
  26. Zahir, D.C. Numerical Solutions for the Most General Multi-Higher Fractional Order Linear Integro-Differential Equations of Fredholm Type in Caputo Sense. Master’s Thesis, University of Sulaimani, Sulaymaniyah, Iraq, March 2017. [Google Scholar]
  27. Hamasalih, S.A. Some Computational Methods for Solving Linear Volterra Integro-Fractional Differential Equations. Master’s Thesis, University of Sulaimani, Sulaymaniyah, Iraq, November 2011. [Google Scholar]
  28. Atkinson, K.E. An Introduction to Numerical Analysis; John Wiley & Sons: New York, NY, USA, 2008. [Google Scholar]
Figure 1. (a)Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 1 with N = 10 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error comparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Figure 1. (a)Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 1 with N = 10 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error comparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Symmetry 14 01693 g001
Figure 2. (a) Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 2 with N = 200 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error omparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Figure 2. (a) Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 2 with N = 200 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error omparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Symmetry 14 01693 g002
Figure 3. (a) Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 3 with N = 20 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error comparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Figure 3. (a) Least square error comparison of Trapezoidal, Simpson’s 1/3 and Midpoint methods for Example 3 with N = 20 . (b) L.S.Error comparison of Trapezoidal and Simpson’s 1/3 methods. (c) L.S.Error comparison of Trapezoidal and Midpoint methods. (d) L.S.Error comparison of Simpson’s 1/3 and Midpoint methods.
Symmetry 14 01693 g003
Table 1. Numerical results for Example 1, by using Trapezoidal, Simpson’s 1/3 and Midpoint methods, with setting h = 0.1 ( N = 10 ) .
Table 1. Numerical results for Example 1, by using Trapezoidal, Simpson’s 1/3 and Midpoint methods, with setting h = 0.1 ( N = 10 ) .
ExactApproximate Solution
x i SolutionTrapezoidalSimpson’s 1/3Midpoint
01111
0.11.11.0999988805105451.0999997053163491.099999738543476
0.21.21.1999947359588971.1999986143443831.199998770584320
0.31.31.2999861995212661.2999963672945371.299996776901128
0.41.41.3999720524239421.3999926433485271.399993472849557
0.51.51.4999511629608371.4999871446069111.499988594119893
0.61.61.5999225114778621.5999796026657421.599981902572140
0.71.71.6998852783756191.6999698017815461.699973206789032
0.81.81.7998390042126151.7999576210153601.799962399468106
0.91.91.8997838277314311.8999430968884461.899949513012665
121.9997207781720541.9999265003280531.999934787801486
L.S.Error1.731659846192182 × 10 7 1.199871274392243 × 10 8 9.44542575567171 × 10 9
Running Time/s1.5319001.0839980.573723
Table 2. Numerical results for Example 1, by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 10 ,   20 ,   100 and 200 .
Table 2. Numerical results for Example 1, by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 10 ,   20 ,   100 and 200 .
L.S.Error
NTrapezoidalSimpson’s 1/3Midpoint
101.731659846192182 × 10 7 1.199871274392243 × 10 8 9.44542575567171 × 10 9
201.59706048132134 × 10 8 9.71114646403669 × 10 10 5.31559621620084 × 10 10
1006.01739779418656 × 10 11 2.92811067355869 × 10 12 7.00371425335641 × 10 13
2005.322519055906852 × 10 12 2.413151704413131 × 10 13 3.98572790064252 × 10 14
Table 3. Numerical results for Example 2 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods, with setting h = 0.005   ( N = 200 ) .
Table 3. Numerical results for Example 2 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods, with setting h = 0.005   ( N = 200 ) .
ExactApproximate Solution
x i SolutionTrapezoidalSimpson’s 1/3Midpoint
10000
1.10.0010.0010571874752290.0010570805637760.001057808414434
1.20.0080.0010571874752290.0082033075639590.008204615780975
1.30.0270.0274248474663590.0274240820944010.027424404330519
1.40.0640.0647117298487930.0647104622397050.064707324828548
1.50.1250.1260572771077040.1260554135888370.126045784833891
1.60.2160.2174553457523410.2174528060446280.217433376969525
1.70.3430.3449002500346190.3448969662587180.344864386404779
1.80.5120.5143865395184670.5143824556351030.514333531797560
1.90.7290.7319088988248440.7319039702581170.731835830084353
211.0034621233276151.0034563166162221.003366536719375
L.S.Error3.372218807263587 × 10 5 3.360748127151794 × 10 5 3.217068140177817 × 10 5
Running Time/s33.72421544.85798573.434670
Table 4. Numerical results for Example 2 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 10 ,   20 ,   100 and 200 .
Table 4. Numerical results for Example 2 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 10 ,   20 ,   100 and 200 .
L.S.Error
NTrapezoidalSimpson’s 1/3Midpoint
102.973649696929000 × 10 2 2.828859549187800 × 10 2 2.781452375954820 × 10 2
206.255313580691000 × 10 3 6.094665228386000 × 10 3 5.921832956186230 × 10 3
1001.639219388468831 × 10 4 1.629077055504340 × 10 4 1.566048582826550 × 10 4
2003.372218807263587 × 10 5 3.360748127151794 × 10 5 3.217068140177817 × 10 5
Table 5. Numerical results for Example 3 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods with setting h = 0.05   ( N = 20 ) .
Table 5. Numerical results for Example 3 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods with setting h = 0.05   ( N = 20 ) .
ExactApproximate Solution
x i SolutionTrapezoidalSimpson’s 1/3Midpoint
01111
0.10.8948290819243520.8947153591955660.8947363246155770.894731092758491
0.20.7785972418398300.7784638657188390.7785043998008400.778512152261663
0.30.6501411924239970.6500534392976600.6501052285036000.650128318440008
0.40.5081753023587300.5082174048773390.5082683228180410.508309294002346
0.50.3512787292998720.3515528324680860.3515879589110750.351649915261491
0.60.1778811996094910.1785044495882640.1785067389726020.178593362632636
0.7−0.013752707470477−0.012650823016870−0.012700071994150−0.012584564326513
0.8−0.225540928492468−0.223820712965856−0.223941528402110−0.223792441192181
0.9−0.459603111156950−0.457116340563048−0.457329844884369−0.457142059569808
1−0.718281828459046−0.714873159147915−0.715201393706980−0.714969423867113
L.S.Error2.248011363124974 × 10 5 1.883703212204261 × 10 5 2.213048214188453 × 10 5
Running Time/s1.0154040.9356830.928798
Table 6. Numerical results for Example 3 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 20 ,   100 ,   200 and 400 .
Table 6. Numerical results for Example 3 by using Trapezoidal, Simpson’s 1/3 and Midpoint methods for N = 20 ,   100 ,   200 and 400 .
L.S.Error
NTrapezoidalSimpson’s 1/3Midpoint
202.248011363124974 × 10 5 1.883703212204261 × 10 5 2.213048214188453 × 10 5
1008.492239220911016 × 10 7 6.391725685307116 × 10 7 7.017017864489802 × 10 7
2001.895285546764007 × 10 7 1.407627874937719 × 10 7 1.528127122158101 × 10 7
4004.137313887688308 × 10 8 3.052288381134047 × 10 8 3.292994636418602 × 10 8
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ahmed, S.S.; Rasol, H.A. Numerical Computation of Mixed Volterra–Fredholm Integro-Fractional Differential Equations by Using Newton-Cotes Methods. Symmetry 2022, 14, 1693. https://doi.org/10.3390/sym14081693

AMA Style

Ahmed SS, Rasol HA. Numerical Computation of Mixed Volterra–Fredholm Integro-Fractional Differential Equations by Using Newton-Cotes Methods. Symmetry. 2022; 14(8):1693. https://doi.org/10.3390/sym14081693

Chicago/Turabian Style

Ahmed, Shazad Shawki, and Hiwa Abdullah Rasol. 2022. "Numerical Computation of Mixed Volterra–Fredholm Integro-Fractional Differential Equations by Using Newton-Cotes Methods" Symmetry 14, no. 8: 1693. https://doi.org/10.3390/sym14081693

APA Style

Ahmed, S. S., & Rasol, H. A. (2022). Numerical Computation of Mixed Volterra–Fredholm Integro-Fractional Differential Equations by Using Newton-Cotes Methods. Symmetry, 14(8), 1693. https://doi.org/10.3390/sym14081693

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop