# The Analysis of Fractional-Order System Delay Differential Equations Using a Numerical Method

Authors: Pongsakorn Sunthrayuth; Hina M. Dutt; Fazal Ghani; Mohammad Asif Arefin
Journal: Complexity
 (2022)

---

## Abstract

To solve fractional delay differential equation systems, the Laguerre Wavelets Method (LWM) is presented and coupled with the steps method in this article. Caputo fractional derivative is used in the proposed technique. The results show that the current procedure is accurate and reliable. Different nonlinear systems have been solved, and the results have been compared to the exact solution and different methods. Furthermore, it is clear from the figures that the LWM error converges quickly when compared to other approaches. When compared with the exact solution to other approaches, it is clear that LWM is more accurate and gets closer to the exact solution faster. Moreover, on the basis of the novelty and scientific importance, the present method can be extended to solve other nonlinear fractional-order delay differential equations.


---

## Body

## 1. Introduction

In 1965, a mathematician named L’Hopital asked Leibniz what would be the solution to the problem if the derivatives and integrals were fractional order. This L’Hopital question has resulted in the creation of new mathematical knowledge, but no one has been able to deal with it for a long time [1]. Mathematicians began to conduct study in the field of fractional derivatives, integration, and the development of a new field of fractional calculus after a period of time. In mathematics, this domain is known as fractional calculus, and it is a significant branch of mathematics that deals with the study of fractional derivatives and integration. Mathematicians have recently started working on fractional calculus because of its wide applications in all fields of research such as economics [2], viscoelastic materials [3], dynamics of interfaces between soft nanoparticles and rough substrates [4], continuum and statistical mechanics [5], solid mechanics [6], and many other topics.Many natural problems can be solved using mathematical formulations by transforming physical facts into equation form. Differential equations (DEs) are a type of equation that is used to model a variety of phenomena. However, certain cases are too complicated to be solved using a differential equation. In this case, the researchers used fractional differential equations (FDEs), which are more accurate than differential equations with order integers in modelling the phenomenon. FDEs have realised the importance of real-world modelling challenges in recent years. Such as electrochemistry of corrosion [7], electrode-electrolyte polarization [8], heat conduction [9], optics and signal processing [10], diffusion wave [11], circuit systems [12], control theory of dynamical systems [6], probability and statistics [14, 15], fluid flow [16], and so on.Equations with delayed arguments are known as fractional delay differential equations (FDDEs). Time delay, spatial delay, step size delay, constant delay, and so on are examples of delayed arguments. Due to various delay arguments found in nature, FDDEs are classified into distinct types. FDDEs are time delay DDEs, which are equations in which the current time derivatives are dependent on the solution and possibly its derivatives at a previous time. In the last few decades, mathematicians have paid more attention to FDDEs for modelling than simple ODEs, because a small delay has a big impact. FDDEs are employed in a variety of domains of mathematics, including infection diseases, navigation control, population dynamics, circulating blood, and the body’s reaction to carbon dioxide [17–19], as well as some additional applications in advanced research studies.It is necessary to develop accurate, time-efficient, and computationally efficient numerical algorithms for solving FDDEs. Xu and Ma [20] investigated the SEIRS epidemic model with a saturation incidence rate and a time delay that defined the latent period. Rihan et al. [21] investigated a delay differential model, numerically analysed it, and established an effective method of combining chemotherapy with therapeutic immunotherapy in 2014. The global stability of the Lotka–Volterra autonomous model with diffusion and time delay was studied by Beretta and Takeuchi [22]. Lv and Gao [23] used the well-known reproducing kernel Hilbert space approach to solve neutral functional proportional delay differential equations (RKHSM). Galach [24] investigated the time delay in the model presented by Kuznetsov and Taylor, where the time delay was included to gain better compatibility with reality. Furthermore, some researchers discussed the behaviour of delay fractional differential equations or a system of delay fractional differential equations, as well as their stability and analysis. Some works, such as in [25, 26], demonstrate this style of research.In a number of situations, exact FDDEs solutions are difficult to get. As a result, the researchers’ key goal is to develop a numerical or analytical solution to FDDEs. As a result, many strategies have been employed such as the New Predictor Corrector Method (NPCM) [27], New Iterative Method (NIM) [28], Adomian Decomposition Method (ADM) [29], Backward Differentiation Formula (BDF) [30], Chebyshev Pseudospectral Method (CPM) [31], Legendre–Gauss Collocation Method (LGCM) [32], Adams–Bashforth–Moulton Algorithm (ABMA) [33], operational matrix based on poly-Bernoulli polynomials (OMM) [34], and Runge Kutta-type Method (RKM) [35]. Overall, some of the approaches used to obtain numerical or analytical solutions to FDDEs have low accuracy of convergence, while others have great accuracy. Among all of these approaches, the wavelet approximation family is one of the more recent methods for locating FDDE solutions. For the approximate solution of FDDEs systems in the current study, we implement Laguerre Wavelets Method (LWM) in combination with the steps method. The proposed solution is shown to be entirely compatible with the complexity of such problems and to be extremely user-friendly. The error comparison shows that the suggested technique has a very high level of accuracy.The structure of remaining paper is summarized as follows. Section2 defines some basic definitions related to our present work. The general methodology for solving FDDEs is provided in Section 3. Section 4 presents the main results, numerical simulations, and graphical representations. The conclusion along with future research directions is drawn in Section 5.


## 2. Preliminaries Concept

This section introduces the basic concept and several important definitions from fractional calculus, which we will apply in our current research.


### 2.1. Definition

The following mathematical statement demonstrates Caputo’s definition for fractional derivatives of orderδ [36, 37].(1)Dδξψ=1Γm−δ∫0ψψ−τm−δ−1ξmτdτ,for n−1<δ≤m, m∈ℕ, ψ>0, ξ∈ℂ−1n.


### 2.2. Definition

The Riemann–Liouville integral operator for orderδ is given as [36, 37].(2)Iδξψ=1Γδ∫0ψψ−τδ−1ξτdτ.The following are the properties of the Caputo derivative and Riemann–Liouville integral operators.(3)DδIδξψ=ξψ,IδDδξψ=ξψ−∑k=0n−1ξk0+k!ψk,ψ≥0n−1<δ<n.


## 2.1. Definition

The following mathematical statement demonstrates Caputo’s definition for fractional derivatives of orderδ [36, 37].(1)Dδξψ=1Γm−δ∫0ψψ−τm−δ−1ξmτdτ,for n−1<δ≤m, m∈ℕ, ψ>0, ξ∈ℂ−1n.


## 2.2. Definition

The Riemann–Liouville integral operator for orderδ is given as [36, 37].(2)Iδξψ=1Γδ∫0ψψ−τδ−1ξτdτ.The following are the properties of the Caputo derivative and Riemann–Liouville integral operators.(3)DδIδξψ=ξψ,IδDδξψ=ξψ−∑k=0n−1ξk0+k!ψk,ψ≥0n−1<δ<n.


## 3. Laguerre Wavelets

Wavelets [38–40] are a family of functions made up of dilation and translation of a single function called the mother wavelet, φψ. The family of continuous wavelets [41] is formed when the dilation parameter a and the translation parameter b vary continuously.(4)φa,bψ=a−1/2φψ−ba,a,b∈R,a≠0.The following family of discrete wavelets results from restricting the parametersa and b to discrete values as a=a0−p,a=nb0a0−p,a0>1,b0>0,(5)φp,nψ=a−p/2φa0pψ−nb0,p,n∈ℤ,where the wavelet basis for L2R is φp,n. When a0=2 and b0=1, for instance, φp,nψ forms an orthonormal basis. There are four arguments in the Laguerre wavelets Φn,mψ=φk,n,m,ψ, n=1,2,…,2k−1, where k is non-negative integer, m represents the Laguerre polynomials degree, and represents normalized time. Over the interval 0,1, they are defined as(6)φn,m=2p/2ℒ˜m2pψ−2n+1,n−12p−1≤ψ<n2p−1,0,Otherwise,where(7)ℒ˜m=1m!ℒmψm=0,1,2,…A−1.m=0,1,2,…,M−1. The coefficients are utilised in (10) to determine orthonormality. The Laguerre polynomials having degree m with regard to wψ=1 weight function on the interval 0,∞ are Lmψ and satisfy the recursive formula:(8)ℒ0ψ=1,ℒ1ψ=1−ψ,ℒm+2=2m+3−xℒm+1ψ−m+1ℒmm+2m=0,1,2,3,4,….whereModified Laguerre wavelets method (MLWM): Here, we consider the delay differential equation of the form:(9)yαψ=fψ+gψyψa−c,0<ψ<b,0<α≤1,yψ=pψ,−b≤ψ≤0,where fψ is a provided continuous linear or nonlinear function and gψ is a source term function. Using the proposed method, transform the delay differential equation (12) to an inhomogeneous ordinary differential equation by using the initial source, pψ, as shown in (12):(10)yαψ=fψ+gψpψa−c,0<ψ<b,1<α≤2.Equation (14) can be expanded as a Laguerre wavelets series as follows:(11)yψ=∑n=0∞∑m=0∞dn,mφn,mψ,where φn,mψ is determined by (9). The truncated series is used to approximate yψ.(12)yp,A=∑n=0∞∑m=0∞dn,mφn,mψ,Then, there should be a total of2p−1A conditions for determining the 2p−1A coefficient:(13)c10,c11…cA−1…c20,c2A−1…c2p−11…c2p−1A−1.Since the initial and boundary conditions, respectively, provide the conditions.(14)yp,A0=∑n=12p−1∑m=0A−1dn,mφn,m0=q0.(15)ddψyp,A1=ddψ∑n=12p−1∑m=0A−1dn,mφn,m1=q´1.We see that there should be2p−1A−2 extra condition to recover the unknown coefficient dn,m. These conditions can be obtained by substituting (14) in (12):(16)dαdψα∑n=12p−1∑m=0A−3dn,mφn,mψ=f∑n=12p−1∑m=0A−3dn,mφn,mψ+gψpψa−c.We, now assume equation (18) is exact at 2p−1A−3 points ψi as follows:(17)dαdψα∑n=12p−1∑m=0A−3dn,mφn,mψi=f∑n=12p−1∑m=0A−3dn,mφn,mψi+gψipψia−c.The best choice of theψi points are the zeros of the shifted Laguerre polynomials of degree 2p−1A−2 in the interval [0, 1] that is ψi=si−1/2, where si=cos2i−1π/2p−1A−1,i=1,2,3,…2p−1A−2. Since the initial and boundary conditions, respectively, provide the conditions. Combining equations (9) and (12) yields 2p−1A linear equations from which the unknown coefficients, dn,m, can be computed. The same technique is followed for first- and second-order delay differential equations.


## 4. Numerical Representation

### 4.1. Example

Consider the system of fractional ordinary delay differential equations [42],(18)Dδξψ=−ζψ−2e−34ψcos12ψsin14ψξ0.25ψ−e−ψcos12ψζ0.5ψ,Dδζψ=eψξ20.5ψ−ζ20.5ψ,with the initial sources ξ0=1, ζ0=0, and having exact solution at δ=1 as ξψ=e−ψcosψ, ζψ=sinψ.Table1 shows the exact solution and numerical results achieved using the proposed method. Table 2 shows the comparison on the basis of absolute error between our technique and those derived from FBPs. When δ=1, the behaviour of the exact solution and proposed method solution of this problem is shown in Figures 1 and 2, respectively, whereas the error comparison of CPM and FBPs is shown in Figures 3 and 4. Figures 5 and 6 show graphical representations for different fractional orders of δ, confirming that the proposed method solution converges to the exact solution as the value of δ approaches from fractional-order towards integer-order.Table 1
Comparison of the exact and MLWM solution for example 1 atm = 9.
ψExactξψExactζψMLWM solutionξψMLWM solutionζψ01.0000000000000000.0000000000000001.0000000000000000.0000000000000000.10.9003169998451940.9983341664682810.9003169998451940.9983341664682810.20.8024106473425200.1986693307950610.8024106473425270.1986693307950610.30.7077306780263510.2955202066613390.7077306780256620.2955202066613390.40.6174056479016460.3894183423086500.6174056479016530.3894183423086370.50.5322807302156710.4794255386042030.5322807302152730.4794255386042030.60.452953789145250.56464247339503530.452953789145497.56464247339470350.70.3798093899251540.6442176872376910.379809389925360.6442176872368760.80.3130505040044800.7173560908995220.3130505040043460.7173560908985010.90.2527277532911690.7833269096274830.2527277532918680.7833269096282491.00.1987661103464130.8414709848078960.1987661103464800.841470984813237Table 2
Error estimation of proposed method with FBPs for example 1 atm = 9.
ψErrorξMLWMErrorζMLWMErrorξFBPsErrorζFBPs0.27.1240361137E-152.1190342156276E-171.22E-113.55E-120.43.1090504913E-131.2920241091118E-149.91E-121.04E-110.64.0077803050E-123.3177676221433E-137.20E-121.59E-110.82.5602591191E-121.0209754593019E-126.56E-122.06E-111.04.9538867166E-115.3406078485367E-127.58E-101.47E-11Figure 1
Behaviour of the exact solution and proposed method solution forξψ of problem 1.Figure 2
Behaviour of the exact solution and proposed method solution forζψ of problem 1.Figure 3
FBPs and proposed method error analysis forξψ of example 1.Figure 4
FBPs and proposed method error analysis forζψ of example 1.Figure 5
The error comparison at various fractional-orders of forξψ example 1.Figure 6
The error comparison at various fractional-orders of forζψ example 1.


### 4.2. Example

Consider the system of fractional ordinary delay differential equations [43].(19)Dδξψ=−ξψ+ζψ2+34ψ2+2Γ3−δψ2−δ,Dδζψ=ζψ−ξψ2−34ψ2+2Γ3−δψ2−δ.The exact solution is given byξψ=ψ2 and ζψ=ψ2.The comparison among the exact solution and the Spline function polynomial technique solution are shown in Table3. In Table 4, the errors acquired by the current technique are compared to those obtained by the Spline function polynomial method. In Figures 7 and 8, we compare the exact and approximated solutions, which shows that they are very close to each other. In addition, Figures 9 and 10 show the MLWM and Spline function error comparisons, demonstrating that suggested approach is in best agreement with the exact solution.Table 3
Comparison at different fractional-order ofδ on the basis of error for example 2.
δψExactMLWMsolutionMLWMerrorSpline functions0.10.010.00010.000099999863751.3625E-108.2E-40.020.00040.00039999985531.447E-102.5E-30.030.00090.00089999986241.376E-104.7E-30.040.00160.0015999998831.17E-107.3E-30.050.00250.0024999998131.87E-101.0E-20.20.010.00010.00010000001111.11E-114.4E-40.020.00040.00040000010351.035E-101.4E-30.030.00090.00089999998451.55E-112.7E-30.040.00160.0016000000575.7E-114.4E-30.050.00250.0025000000222.2E-116.1E-30.30.010.00010.00010000017981.798E-102.1E-40.020.00040.00040000023032.303E-107.1E-40.030.00090.00090000017171.717E-101.4E-30.040.00160.0016000001061.06E-102.4E-30.050.00250.0025000001361.36E-103.5E-30.40.010.00010.00010000001891.89E-118.1E-50.020.00040.00039999995024.98E-112.9E-40.030.00090.00089999989161.084E-106.1E-40.040.00160.0015999999415.9E-111.0E-30.050.00250.0024999999973.000E-121.0E-30.50.010.00010.000099999880741.1926E-104.5E-60.020.00040.00039999977802.220E-102.6E-50.030.00090.00089999969063.094E-107.0E-50.040.00160.0015999997173.83E-101.4E-40.050.00250.0024999996573.43E-102.5E-4Table 4
Comparison at different fractional-orders ofδ on the basis of error for example 2.
δψExactMLWMsolutionMLWMerrorSpline functions0.010.00010.000099999863751.3625E-108.2E-40.020.00040.00039999985531.447E-102.5E-30.10.030.00090.00089999986241.376E-104.7E-30.040.00160.0015999998831.17E-107.3E-30.050.00250.0024999998131.87E-101.0E-20.010.00010.00010000001111.11E-114.4E-40.020.00040.00040000010351.035E-101.4E-30.20.030.00090.00089999998451.55E-112.7E-30.040.00160.0016000000575.7E-114.4E-30.050.00250.0025000000222.2E-116.1E-30.010.00010.00010000017981.798E-102.1E-40.020.00040.00040000023032.303E-107.1E-40.30.030.00090.00090000017171.717E-101.4E-30.040.00160.0016000001061.06E-102.4E-30.050.00250.0025000001361.36E-103.5E-30.010.00010.00010000001891.89E-118.1E-50.020.00040.00039999995024.98E-112.9E-40.40.030.00090.00089999989161.084E-106.1E-40.040.00160.0015999999415.9E-111.0E-30.050.00250.0024999999973.000E-121.0E-30.010.00010.000099999880741.1926E-104.5E-60.020.00040.00039999977802.220E-102.6E-50.50.030.00090.00089999969063.094E-107.0E-50.040.00160.0015999997172.83E-101.4E-40.050.00250.0024999996573.43E-102.5E-4Figure 7
Analysis of the exact and proposed method solution forξψ of problem 2.Figure 8
Analysis of the exact and proposed method solution forζψ of problem 2.Figure 9
Spline functions and proposed method error analysis forξψ of example 2.Figure 10
Spline functions and proposed method error analysis forζψ of example 2.


## 4.1. Example

Consider the system of fractional ordinary delay differential equations [42],(18)Dδξψ=−ζψ−2e−34ψcos12ψsin14ψξ0.25ψ−e−ψcos12ψζ0.5ψ,Dδζψ=eψξ20.5ψ−ζ20.5ψ,with the initial sources ξ0=1, ζ0=0, and having exact solution at δ=1 as ξψ=e−ψcosψ, ζψ=sinψ.Table1 shows the exact solution and numerical results achieved using the proposed method. Table 2 shows the comparison on the basis of absolute error between our technique and those derived from FBPs. When δ=1, the behaviour of the exact solution and proposed method solution of this problem is shown in Figures 1 and 2, respectively, whereas the error comparison of CPM and FBPs is shown in Figures 3 and 4. Figures 5 and 6 show graphical representations for different fractional orders of δ, confirming that the proposed method solution converges to the exact solution as the value of δ approaches from fractional-order towards integer-order.Table 1
Comparison of the exact and MLWM solution for example 1 atm = 9.
ψExactξψExactζψMLWM solutionξψMLWM solutionζψ01.0000000000000000.0000000000000001.0000000000000000.0000000000000000.10.9003169998451940.9983341664682810.9003169998451940.9983341664682810.20.8024106473425200.1986693307950610.8024106473425270.1986693307950610.30.7077306780263510.2955202066613390.7077306780256620.2955202066613390.40.6174056479016460.3894183423086500.6174056479016530.3894183423086370.50.5322807302156710.4794255386042030.5322807302152730.4794255386042030.60.452953789145250.56464247339503530.452953789145497.56464247339470350.70.3798093899251540.6442176872376910.379809389925360.6442176872368760.80.3130505040044800.7173560908995220.3130505040043460.7173560908985010.90.2527277532911690.7833269096274830.2527277532918680.7833269096282491.00.1987661103464130.8414709848078960.1987661103464800.841470984813237Table 2
Error estimation of proposed method with FBPs for example 1 atm = 9.
ψErrorξMLWMErrorζMLWMErrorξFBPsErrorζFBPs0.27.1240361137E-152.1190342156276E-171.22E-113.55E-120.43.1090504913E-131.2920241091118E-149.91E-121.04E-110.64.0077803050E-123.3177676221433E-137.20E-121.59E-110.82.5602591191E-121.0209754593019E-126.56E-122.06E-111.04.9538867166E-115.3406078485367E-127.58E-101.47E-11Figure 1
Behaviour of the exact solution and proposed method solution forξψ of problem 1.Figure 2
Behaviour of the exact solution and proposed method solution forζψ of problem 1.Figure 3
FBPs and proposed method error analysis forξψ of example 1.Figure 4
FBPs and proposed method error analysis forζψ of example 1.Figure 5
The error comparison at various fractional-orders of forξψ example 1.Figure 6
The error comparison at various fractional-orders of forζψ example 1.


## 4.2. Example

Consider the system of fractional ordinary delay differential equations [43].(19)Dδξψ=−ξψ+ζψ2+34ψ2+2Γ3−δψ2−δ,Dδζψ=ζψ−ξψ2−34ψ2+2Γ3−δψ2−δ.The exact solution is given byξψ=ψ2 and ζψ=ψ2.The comparison among the exact solution and the Spline function polynomial technique solution are shown in Table3. In Table 4, the errors acquired by the current technique are compared to those obtained by the Spline function polynomial method. In Figures 7 and 8, we compare the exact and approximated solutions, which shows that they are very close to each other. In addition, Figures 9 and 10 show the MLWM and Spline function error comparisons, demonstrating that suggested approach is in best agreement with the exact solution.Table 3
Comparison at different fractional-order ofδ on the basis of error for example 2.
δψExactMLWMsolutionMLWMerrorSpline functions0.10.010.00010.000099999863751.3625E-108.2E-40.020.00040.00039999985531.447E-102.5E-30.030.00090.00089999986241.376E-104.7E-30.040.00160.0015999998831.17E-107.3E-30.050.00250.0024999998131.87E-101.0E-20.20.010.00010.00010000001111.11E-114.4E-40.020.00040.00040000010351.035E-101.4E-30.030.00090.00089999998451.55E-112.7E-30.040.00160.0016000000575.7E-114.4E-30.050.00250.0025000000222.2E-116.1E-30.30.010.00010.00010000017981.798E-102.1E-40.020.00040.00040000023032.303E-107.1E-40.030.00090.00090000017171.717E-101.4E-30.040.00160.0016000001061.06E-102.4E-30.050.00250.0025000001361.36E-103.5E-30.40.010.00010.00010000001891.89E-118.1E-50.020.00040.00039999995024.98E-112.9E-40.030.00090.00089999989161.084E-106.1E-40.040.00160.0015999999415.9E-111.0E-30.050.00250.0024999999973.000E-121.0E-30.50.010.00010.000099999880741.1926E-104.5E-60.020.00040.00039999977802.220E-102.6E-50.030.00090.00089999969063.094E-107.0E-50.040.00160.0015999997173.83E-101.4E-40.050.00250.0024999996573.43E-102.5E-4Table 4
Comparison at different fractional-orders ofδ on the basis of error for example 2.
δψExactMLWMsolutionMLWMerrorSpline functions0.010.00010.000099999863751.3625E-108.2E-40.020.00040.00039999985531.447E-102.5E-30.10.030.00090.00089999986241.376E-104.7E-30.040.00160.0015999998831.17E-107.3E-30.050.00250.0024999998131.87E-101.0E-20.010.00010.00010000001111.11E-114.4E-40.020.00040.00040000010351.035E-101.4E-30.20.030.00090.00089999998451.55E-112.7E-30.040.00160.0016000000575.7E-114.4E-30.050.00250.0025000000222.2E-116.1E-30.010.00010.00010000017981.798E-102.1E-40.020.00040.00040000023032.303E-107.1E-40.30.030.00090.00090000017171.717E-101.4E-30.040.00160.0016000001061.06E-102.4E-30.050.00250.0025000001361.36E-103.5E-30.010.00010.00010000001891.89E-118.1E-50.020.00040.00039999995024.98E-112.9E-40.40.030.00090.00089999989161.084E-106.1E-40.040.00160.0015999999415.9E-111.0E-30.050.00250.0024999999973.000E-121.0E-30.010.00010.000099999880741.1926E-104.5E-60.020.00040.00039999977802.220E-102.6E-50.50.030.00090.00089999969063.094E-107.0E-50.040.00160.0015999997172.83E-101.4E-40.050.00250.0024999996573.43E-102.5E-4Figure 7
Analysis of the exact and proposed method solution forξψ of problem 2.Figure 8
Analysis of the exact and proposed method solution forζψ of problem 2.Figure 9
Spline functions and proposed method error analysis forξψ of example 2.Figure 10
Spline functions and proposed method error analysis forζψ of example 2.


## 5. Conclusion

We used the MLWM to solve fractional delay differential equations systems in this research. The proposed method’s convergence is given special consideration. As demonstrated in Figures1–12, the fractional-order delay differential equation solution approaches towards the solution of the integer-order delay differential equation. The results obtained by implementing the proposed method are in great agreement with the exact solution and are more accurate than those obtained by implementing other techniques. The proposed method (MLWM) is extremely user-friendly but extremely accurate, according to computational effort and numerical results. The computations work in this article are done using Maple.Figure 11
The error comparison at various fractional-orders forξψ of example 2.Figure 12
The error comparison at various fractional-orders forζψ of example 2.



---

*Source: 3570667-2022-05-28.xml*