# Numerical, Approximate Solutions, and Optimal Control on the Deathly Lassa Hemorrhagic Fever Disease in Pregnant Women

Authors: M. Higazy; A. El-Mesady; A. M. S. Mahdy; Sami Ullah; A. Al-Ghamdi
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

This paper is devoted to the model of Lassa hemorrhagic fever (LHF) disease in pregnant women. This disease is a biocidal fever and epidemic. LHF disease in pregnant women has negative impacts that were initially appeared in Africa. In the present study, we find an approximate solution to the fractional-order model that describes the fatal LHF disease. Laplace transforms coupled with the Adomian decomposition method (ADM) are applied. In addition, the fractional-order LHF model is numerically simulated in terms of a varied fractional order. Furthermore, a fractional order optimal control for the LHF model is studied.


---

## Body

## 1. Introduction

Infectious diseases remain a major threat to human health and welfare throughout history. Their spread is affected by several factors such as mode of transmission, infectious agent, infectious periods, incubation period, susceptibility, and resistance [1]. Some devastating infectious diseases like Lassa fever are endemic in many parts of the world and continue to emerge.Lassa fever, also known as Lassa hemorrhagic fever (LHF), is an animal-borne, or zoonotic, acute viral hemorrhagic illness caused by Lassa virus, a member of the familyArenaviridae. The reservoir is the rats Mastomys natalensis, M. erythroleucus, and Hylomyscus pamfi. Food contaminated with rodent urine, saliva, or feces is the main cause of human cases. Human-to-human transmission can occur via exposure to the blood, feces, urine, saliva, or vomitus of infected patients [2, 3]. The incubation period is ca. 2-21 days. LHF is basically endemic in some countries of West Africa, and outbreaks have occurred in Benin, Ghana, Guinea, Liberia, Nigeria, Sierra Leone, and Togo. Cases have been imported to Germany, United Kingdom, United States, and Sweden. This disease with high mortality rates (e.g., 80% among pregnant women) is likely affecting between 100,000 and 300,000 people every year, and it kills around 5,000 people, almost all in West Africa alone [4, 5].Fractional Calculus (FC) is a fruitful field of mathematical research with numerous applications in engineering [6], nanotechnology [7], optics [8], human diseases [9], and chaos soliton theory [10]. Application of FC is also adopted in biology, heat transfer, system identification, genetic algorithms, traffic systems, telecommunications, physics as well as finance, and economics [11–20].Herein, we are interested in the paradigm showing the nature of LHF, since it is lethal and transmissible [1, 2, 21, 22]. Using mathematical modeling to probe this fatal illness may provide insight, answers, and guidance useful for controlling its spread, and the originating applications. The given model is studied through Caputo fractional operator. In this model, the Caputo fractional operator is used since it permits both initial and boundary conditions to involve in modeling problems. The Caputo fractional operator can be applied on sufficiently differentiable functions only. One of the advantages of the fractional derivative of Caputo compared with the fractional derivative 3 of Riemann-Liouville is that the fractional derivative of Caputo for the constant is equal to zero. We develop an approximate solution to the fractional-order model describing the LHF affecting pregnant women using the Laplace transform coupled with the ADM. Numerical simulations are also presented. A fractional optimal control for the SIRD model is ultimately offered. For a good survey on the approximating of the differential equations, see [23–26]. Some recent research fractional optimal controls can be found in [27–29]. Atangana [21] studied the LHF using beta differential operator, where the model was proposed as follows in (1). Table 1 shows the meaning of the parameters in the model (1).
(1)D0AtβSt=−u1StIt+μ−lN+ωRt+u2−lSt,D0AtβIt=u1StIt−ω+u2It−u2St,D0AtβRt=u2It−ωRt,D0AtβDt=ωIt+l−μN+lSt,where D0Atβgt=t+1/Γβ1−βlimh⟶0gt+h−gt/h.Table 1
The meaning of the parameters in the model (1).
ParametersMeaningStSusceptible pregnant women populationItInfected pregnant women populationRtRecovery of pregnant women populationDtThe population of pregnant women dying in that countryu1The rate at which the women are infectedNThe total number of adult women in the given countryμThe rate at which the women are pregnantωThe rate of death of pregnant women with LHFlThe rate of death of pregnant women with another disease or with natural deathu2The rate of recovery of pregnant women from LHFAlso, the model has been studied by Goyal et al. [22]. The q-homotopy analysis transform method (q-HATM) has been applied for solving the LHF model [30]. The LHF model via Atangana-Baleanu fractional derivative has been studied in [31]. The SIR model has been studied through Caputo-Fabrizio fractional operator in [32].In this paper, we reconsider the system (1) by the following fractional sense of the system (1)
(2)Dα1St=−u1SI+μ−lN+ωR+u2−lS,(3)Dα2It=u1SI−ω+u2I−u2S,(4)Dα3Rt=u2I−ωR,(5)Dα4Dt=ωI+l−μN+lS,where Dr is the Caputo fractional derivative of order r.We organize the rest of our paper as follows. In Section2, we provide some basic mathematical concepts that will be necessary for our work. Also, the description of the proposed model is presented. We then follow by presenting the series solution for the proposed model in Section 3. We provide the approximate solution construction in Section 4. The numerical solution of the proposed model, having arbitrary order, is given in Section 5. The optimal control case is studied in Section 6. In Section 7, a concise conclusion of our paper results is presented.


## 2. Preliminaries

For studying the proposed model (2), we present in this section the basic definitions concerned with fractional calculus. For a good survey on the basic definitions and several applications on fractional calculus, see [6, 7, 33].


### 2.1. Some Important Definitions

Definition 1.
The Riemann–Liouville integral is represented by(6)Iαzt=1Γα∫0tt−xα−1zxdx,where Γ is the known gamma function and α∈0,1 is the order of this fractional integral.Definition 2.
The Caputo-type fractional derivative of orderα is defined as
(7)Dαzt=1Γk−α∫0tt−xk−α−1zkxdx,Definition 3.
The Laplace transform of Caputo fractional derivative of orderα [33] is given by
(8)LDαzt=sαZs−∑i=0m−1sα−i−1zi0;m−1<α<m,m∈N;Zs=Lzt.


### 2.2. Mathematical Formulation and Description for the Proposed Model

The meaning of the parameters is shown in Table1. Now, the description and the mathematical formulation of the model can be shown as follows.With timet, we can express the rate of change in susceptible population as follows
(9)Dα1St=−u1SI+μ−lN+ωR+u2−lS.With timet, we can express the rate of change in the population of the infected pregnant women as follows
(10)Dα2It=u1SI−ω+u2I−u2S.Here,u1SI is the number of pregnant women excluded from the group of susceptible.With timet, the rate of change in recovery population is written as
(11)Dα3Rt=u2I−ωR.Finally, the rate of change of dying population with timet is given as
(12)Dα4Dt=ωI+l−μN+lS.Hence, the mathematical model for the LHF is given by model (2), with initial conditions
(13)S0=S0,I0=I0,R0=R0,D0=D0.Also0<αi<1,i∈1,2,3,4. Here, St,It,Rt, and Dt are sufficiently differentiable functions. Model (2) has not yet been studied by using the Laplace transform coupled with the ADM. Also, we introduce numerical simulations and fractional optimal controls for the SIRD model.The interaction between different compartments of LHF model (2) can be viewed via the following signal flow graph shown in Figure 1.Figure 1
Signal flow graph of the studied LHF model.


### 2.3. Calculating the Basic Reproduction NumberR0

Basic reproduction numberR0 is the expected production of newly infected females, in a totally susceptible women, by meeting a typical infective woman. The local stability and instability of the disease-free equilibrium (DFE) depend on the value of R0. Also, it defines the local stability threshold for the DFE. Moreover, it is a very important tool for controlling the disease, and it is an essential epidemiological criteria of disease. When R0<1, the DFE is locally asymptotically stable; a small number of infections into the population may cause it to evolve into an endemic prevalence. From another point of view, when R0>1, the DFE is locally unstable; a sufficiently small number of infected women will produce an outbreak. Here, R0 is estimated from the paradigm of nonlinear FDE’s (2) via the next-generation matrix method [34].Since the Jacobian matrix is(14)J=−u1I+u2−l−u1Sω0u1I−u2u1S−ω−u2000u2−ω0lω00.At the disease-free equilibrium pointμ−l/lN,0,0,0 and u2=0 then by putting the model in the F−V form [34] where
(15)F=u1μ−lNl000,V=ω00ω,V−1=1ω001ω,then(16)FV−1=u1μ−lNlω000.Then, the basic reproduction number of the given paradigm isR0=u1μ−l/ωlN. If R0<1, then, infection cannot develop because an infected woman infects less than one new woman during their infectious period on average. Conversely, if R0>1, the disease can enter the population because on average each infected woman will infect more than one new woman.


## 2.1. Some Important Definitions

Definition 1.
The Riemann–Liouville integral is represented by(6)Iαzt=1Γα∫0tt−xα−1zxdx,where Γ is the known gamma function and α∈0,1 is the order of this fractional integral.Definition 2.
The Caputo-type fractional derivative of orderα is defined as
(7)Dαzt=1Γk−α∫0tt−xk−α−1zkxdx,Definition 3.
The Laplace transform of Caputo fractional derivative of orderα [33] is given by
(8)LDαzt=sαZs−∑i=0m−1sα−i−1zi0;m−1<α<m,m∈N;Zs=Lzt.


## 2.2. Mathematical Formulation and Description for the Proposed Model

The meaning of the parameters is shown in Table1. Now, the description and the mathematical formulation of the model can be shown as follows.With timet, we can express the rate of change in susceptible population as follows
(9)Dα1St=−u1SI+μ−lN+ωR+u2−lS.With timet, we can express the rate of change in the population of the infected pregnant women as follows
(10)Dα2It=u1SI−ω+u2I−u2S.Here,u1SI is the number of pregnant women excluded from the group of susceptible.With timet, the rate of change in recovery population is written as
(11)Dα3Rt=u2I−ωR.Finally, the rate of change of dying population with timet is given as
(12)Dα4Dt=ωI+l−μN+lS.Hence, the mathematical model for the LHF is given by model (2), with initial conditions
(13)S0=S0,I0=I0,R0=R0,D0=D0.Also0<αi<1,i∈1,2,3,4. Here, St,It,Rt, and Dt are sufficiently differentiable functions. Model (2) has not yet been studied by using the Laplace transform coupled with the ADM. Also, we introduce numerical simulations and fractional optimal controls for the SIRD model.The interaction between different compartments of LHF model (2) can be viewed via the following signal flow graph shown in Figure 1.Figure 1
Signal flow graph of the studied LHF model.


## 2.3. Calculating the Basic Reproduction NumberR0

Basic reproduction numberR0 is the expected production of newly infected females, in a totally susceptible women, by meeting a typical infective woman. The local stability and instability of the disease-free equilibrium (DFE) depend on the value of R0. Also, it defines the local stability threshold for the DFE. Moreover, it is a very important tool for controlling the disease, and it is an essential epidemiological criteria of disease. When R0<1, the DFE is locally asymptotically stable; a small number of infections into the population may cause it to evolve into an endemic prevalence. From another point of view, when R0>1, the DFE is locally unstable; a sufficiently small number of infected women will produce an outbreak. Here, R0 is estimated from the paradigm of nonlinear FDE’s (2) via the next-generation matrix method [34].Since the Jacobian matrix is(14)J=−u1I+u2−l−u1Sω0u1I−u2u1S−ω−u2000u2−ω0lω00.At the disease-free equilibrium pointμ−l/lN,0,0,0 and u2=0 then by putting the model in the F−V form [34] where
(15)F=u1μ−lNl000,V=ω00ω,V−1=1ω001ω,then(16)FV−1=u1μ−lNlω000.Then, the basic reproduction number of the given paradigm isR0=u1μ−l/ωlN. If R0<1, then, infection cannot develop because an infected woman infects less than one new woman during their infectious period on average. Conversely, if R0>1, the disease can enter the population because on average each infected woman will infect more than one new woman.


## 3. Series Solution for the Proposed Model

Applying Laplace transformation of (2) gives
(17)LDα1St=L−u1SI+μ−lN+ωR+u2−lS,LDα2It=Lu1SI−ω+u2I−u2S,LDα3Rt=Lu2I−ωR,LDα4Dt=LωI+l−μN+lS,which gives
(18)sα1S∗s−sα1−1S0=L−u1SI+μ−lN+ωR+u2−lS,sα2I∗s−sα2−1I0=Lu1SI−ω+u2I−u2S,sα3R∗s−sα3−1R0=Lu2I−ωR,sα4D∗s−sα4−1D0=LωI+l−μN+lS,where S∗s=LSt,I∗s=LIt,R∗s=LRt,D∗s=LDt.By dividing theith equation of (18) by sαi, we get
(19)S∗s=1sS0+1sα1L−u1SI+μ−lN+ωR+u2−lS,I∗s=1sI0+1sα2Lu1SI−ω+u2I−u2S,R∗s=1sR0+1sα3Lu2I−ωR,D∗s=1sD0+1sα4LωI+l−μN+lS.Taking the inverse of Laplace transformation on (19) gives
(20)St=S0+L−11sα1L−u1SI+μ−lN+ωR+u2−lS,It=I0+L−11sα2Lu1SI−ω+u2I−u2S,Rt=R0+L−11sα3Lu2I−ωR,Dt=D0+L−11sα4LωI+l−μN+lS.Now, we will find the series solution by putting(21)St=∑k=0∞Skt,It=∑k=0∞Ikt,Rt=∑k=0∞Rkt,Dt=∑k=0∞Dkt.Decomposing the termsStIt gives
(22)StIt=∑k=0∞Gkt.The polynomialGkt is known as “Adomian polynomial” with the form
(23)Gkt=1k!ddεk∑r=0kεrSrt∑r=0kεrIrtε=0.By substituting from (21) and (22) into (20), we get
(24)LS0=S0s,LI0=I0s,LR0=R0s,LD0=D0s,(25)LS1=1sα1L−u1G0+μ−lN+ωR0+u2−lS0,LI1=1sα2Lu1G0−ω+u2I0−u2S0,LR1=1sα3Lu2I0−ωR0,LD1=1sα4LωI0+l−μN+lS0,(26)LS2=1sα1L−u1G1+ωR1+u2−lS1,LI2=1sα2Lu1G1−ω+u2I1−u2S1,LR2=1sα3Lu2I1−ωR1,LD2=1sα4LωI1+lS1,(27)LS3=1sα1L−u1G2+ωR2+u2−lS2,LI3=1sα2Lu1G2−ω+u2I2−u2S2,LR3=1sα3Lu2I2−ωR2,LD3=1sα4LωI2+lS2,(28)LS4=1sα1L−u1G3+ωR3+u2−lS3,LI4=1sα2Lu1G3−ω+u2I3−u2S3,LR4=1sα3Lu2I3−ωR3,LD4=1sα4LωI3+lS3,⋮(29)LSk=1sα1L−u1Gk−1+ωRk−1+u2−lSk−1,LIk=1sα2Lu1Gk−1−ω+u2Ik−1−u2Sk−1,LRk=1sα3Lu2Ik−1−ωRk−1,LDk=1sα4LωIk−1+lSk−1.Note thatk∈2,3,4,..⋯ in (29).


## 4. Approximate Solution

### 4.1. Construction of the Approximate Solution

Hereafter, we will discuss an approximate solution for the proposed model. Taking the inverse of Laplace transformation for (24), (25), (26), and (27) gives
(30)L−1S0s=S0,L−1I0s=I0,L−1R0s=R0,L−1D0s=D0,(31)S1=atα1,a=−u1G0+μ−lN+ωR0+u2−lS0Γα1+1,(32)I1=btα2,b=u1G0−ω+u2I0−u2S0Γα2+1,(33)R1=ctα3,c=u2I0−ωR0Γα3+1,(34)D1=d1tα4,d1=ωI0+l−μN+lS0Γα4+1,(35)S2=dt2α1+etα1+α2+ftα1+α3,(36)d=au2−l−u1I0Γα1+1Γ2α1+1,e=−u1S0bΓα2+1Γα1+α2+1,f=cωΓα3+1Γα1+α3+1,(37)I2=gtα1+α2+ht2α2,g=au1I0−u2Γα1+1Γα1+α2+1,h=bu1S0−ω+u2Γα2+1Γ2α2+1,(38)R2=itα2+α3+jt2α3,i=u2bΓα2+1Γα2+α3+1,j=−ωcΓα3+1Γ2α3+1,(39)D2=ptα2+α4+qtα1+α4,p=ωbΓα2+1Γα2+α4+1,q=laΓα1+1Γα1+α4+1,(40)S3=rt2α1+α2+st3α1+r1tα1+2α2+ut2α1+α3+vtα1+α2+α3+wtα1+2α3,(41)r=−u1gS0−u1I0e−u1ab+eu2−lΓα1+α2+1Γ2α1+α2+1,s=−u1I0d+du2−lΓ2α1+1Γ3α1+1,r1=−u1hS0Γ2α2+1Γα1+2α2+1,u=−u1I0f+fu2−lΓα1+α3+1Γ2α1+α3+1,v=iωΓα2+α3+1Γα1+α2+α3+1,w=jωΓ2α3+1Γα1+2α3+1,(42)I3=xtα1+2α2+yt3α2+zt2α1+α2+a1tα1+α2+α3,(43)x=u1S0g+eu1I0−gω+u2+u1ab−u2eΓα1+α2+1Γα1+2α2+1,y=u1S0h−hω+u2Γ2α2+1Γ3α2+1,z=du1I0−u2dΓ2α1+1Γ2α1+α2+1,a1=fu1I0−u2Γα1+α3+1Γα1+α2+α3+1,(44)R3=b1tα1+α2+α3+c1t2α2+α3+d1tα2+2α3+e1t3α3,(45)b1=gu2Γα1+α2+1Γα1+α2+α3+1,c1=hu2Γ2α2+1Γ2α2+α3+1,d1=−ωiΓα2+α3+1Γα2+2α3+1,e1=−ωjΓ2α3+1Γ3α3+1,(46)D3=f1tα1+α2+α4+g1t2α2+α4+h1t2α1+α4+i1tα1+α3+α4,(47)f1=gω+leΓα1+α2+1Γα1+α2+α4+1,g1=ωhΓ2α2+1Γ2α2+α4+1,h1=ldΓ2α1+1Γ2α1+α4+1,i1=lfΓα1+α3+1Γα1+α3+α4+1.SinceSt=∑k=0∞Skt,It=∑k=0∞Ikt,Rt=∑k=0∞Rkt,Dt=∑k=0∞Dkt, then, if we take the first four terms from the above series and take into consideration the equations from (30) to (46), we get
(48)St=S0+S1+S2+S3=S0+atα1+dt2α1+etα1+α2+ftα1+α3+rt2α1+α2+st3α1+r1tα1+2α2+ut2α1+α3+vtα1+α2+α3+wtα1+2α3,(49)It=I0+I1+I2+I3=I0+btα2+gtα1+α2+ht2α2+xtα1+2α2+yt3α2+zt2α1+α2+a1tα1+α2+α3,(50)Rt=R0+R1+R2+R3=R0+ctα3+itα2+α3+jt2α3+b1tα1+α2+α3+c1t2α2+α3+d1tα2+2α3+e1t3α3,(51)Dt=D0+D1+D2+D3=D0+d1tα4+ptα2+α4+qtα1+α4+f1tα1+α2+α4+g1t2α2+α4+h1t2α1+α4+i1tα1+α3+α4.After the compensation by the values ofa,b,c,d1,d,e,f,g,h,i,j,p,q,r,r1,s,u,v,w,x,y,z,a1,b1,c1,d1,e1,f1,g1,h1,i1, we obtain St=S0+S1+S2+S3=(52)S0+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1tα1+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1t2α1−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1tα1+α2+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1tα1+α3+−u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1S0−u1I0−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1−u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1G0−ω+u2I0−u2S0Γα2+1+−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1u2−lΓα1+α2+1Γ2α1+α2+1t2α1+α2+−u1I0−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1u2−lΓ2α1+1Γ3α1+1t3α1+−u1I0u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1u2−lΓα1+α3+1Γ2α1+α3+1t2α1+α3+u2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1ωΓα2+α3+1Γα1+α2+α3+1tα1+α2+α3+−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1ωΓ2α3+1Γα1+2α3+1tα1+2α3−u1u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1S0Γ2α2+1Γα1+2α2+1tα1+2α2,(53)It=I0+I1+I2+I3=I0+u1G0−ω+u2I0−u2S0Γα2+1tα2+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1tα1+α2+u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1t2α2+u1S0−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1+−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1u1I0−−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1ω+u2+u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1G0−ω+u2I0−u2S0Γα2+1−u2−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1Γα1+α2+1Γα1+2α2+1tα1+2α2+u1S0u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1−u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1ω+u2Γ2α2+1Γ3α2+1t3α2+u1I0−u2−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1Γ2α1+1Γ2α1+α2+1t2α1+α2+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1u1I0−u2Γα1+α3+1Γα1+α2+α3+1tα1+α2+α3,(54)Rt=R0+R1+R2+R3=R0+u2I0−ωR0Γα3+1tα3+u2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1tα2+α3−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1t2α3+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1u2Γα1+α2+1Γα1+α2+α3+1tα1+α2+α3+u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1u2Γ2α2+1Γ2α2+α3+1t2α2+α3−ωu2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1Γα2+α3+1Γα2+2α3+1tα2+2α3−ω−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1Γ2α3+1Γ3α3+1t3α3,(55)Dt=D0+D1+D2+D3=D0+ωI0+l−μN+lS0Γα4+1tα4+ωu1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α4+1tα2+α4+l−u1G0+μ−lN+ωR0+u2−lS0Γα1+1Γα1+1Γα1+α4+1tα1+α4+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1ω+l−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1Γα1+α2+1Γα1+α2+α4+1tα1+α2+α4+ωu1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1Γ2α2+1Γ2α2+α4+1t2α2+α4+l−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1Γ2α1+1Γ2α1+α4+1t2α1+α4+lu2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1Γα1+α3+1Γα1+α3+α4+1tα1+α3+α4.Now, by taking the first three terms ofSt, It, Rt, and Dt and putting α1=α2=α3=α4=α,u1=0.4, μ=0.3, N=1000, l=0.2, ω=0.8, u2=0.2, S0=900, I0=10, R0= D0=0, in the proposed model yield (see [30]).
(56)St=S0+S1+S2=900−3500Γα+1tα−1213598.4Γ2α+1t2α,(57)It=I0+I1+I2=10+3410Γα+1tα+1210890Γ2α+1t2α,(58)Rt=R0+R1+R2=0+2Γα+1tα+680.4Γ2α+1t2α,(59)Dt=D0+D1+D2=0+88Γα+1tα+2028Γ2α+1t2α.Figure2 exhibits the dynamical behavior of St,It,Rt, and Dt against various fractional orders (α=0,7;0.8;0.9;1).Figure 2
The behavior ofSt,It,Rt, and Dt populations versus time with different fractional derivative orders.


### 4.2. Convergence and Error Analysis

It is clear that the previous solution is a series. This solution converges rapidly and uniformly to the exact solution. We will use classical techniques to prove the convergence of the series (21), and hence, we can get the sufficient conditions for the convergence of the used method.Theorem 4.
Suppose℘ is the Banach space and σ:℘⟶℘ is a contractive nonlinear operator with for all x,x′∈℘,σx−σx′≤λx−x′,0<λ<1. Then, based on Banach contraction principle, there is a unique point x such that σx=x, where x=s,i,r,d. By applying the ADM, the series in (21) can be written as follows:
(60)xn=σxn‐1,xn‐1=∑k=1n‐1xk,n=1,2,3,⋯,and suppose that x0=x0∈Ajx where Ajx=x′∈℘:x−x′<j, then, we obtain
(61)xn∈Ajx;(62)limn⟶∞xn=x.Proof.
For the first partaxn∈Ajx; by the mathematical induction if n=1, we obtain
(63)x0−x=σx0−σx≤λx0−x.
Suppose the result is true atn−1, hence,
(64)x0−x≤λn−1x0−x.
We have(65)xn−x=σxn−1−σx≤λxn−1−x≤λnx0−x,(66)xn−x≤λnx0−x≤λnj<j.This shows that xn∈Ajx..☐For the second part (b)limn⟶∞xn=x.Sincexn−x≤λnx0−x and limn⟶∞λn=0. So, we obtain limn⟶∞xn−x=0⇒limn⟶∞xn=x.


## 4.1. Construction of the Approximate Solution

Hereafter, we will discuss an approximate solution for the proposed model. Taking the inverse of Laplace transformation for (24), (25), (26), and (27) gives
(30)L−1S0s=S0,L−1I0s=I0,L−1R0s=R0,L−1D0s=D0,(31)S1=atα1,a=−u1G0+μ−lN+ωR0+u2−lS0Γα1+1,(32)I1=btα2,b=u1G0−ω+u2I0−u2S0Γα2+1,(33)R1=ctα3,c=u2I0−ωR0Γα3+1,(34)D1=d1tα4,d1=ωI0+l−μN+lS0Γα4+1,(35)S2=dt2α1+etα1+α2+ftα1+α3,(36)d=au2−l−u1I0Γα1+1Γ2α1+1,e=−u1S0bΓα2+1Γα1+α2+1,f=cωΓα3+1Γα1+α3+1,(37)I2=gtα1+α2+ht2α2,g=au1I0−u2Γα1+1Γα1+α2+1,h=bu1S0−ω+u2Γα2+1Γ2α2+1,(38)R2=itα2+α3+jt2α3,i=u2bΓα2+1Γα2+α3+1,j=−ωcΓα3+1Γ2α3+1,(39)D2=ptα2+α4+qtα1+α4,p=ωbΓα2+1Γα2+α4+1,q=laΓα1+1Γα1+α4+1,(40)S3=rt2α1+α2+st3α1+r1tα1+2α2+ut2α1+α3+vtα1+α2+α3+wtα1+2α3,(41)r=−u1gS0−u1I0e−u1ab+eu2−lΓα1+α2+1Γ2α1+α2+1,s=−u1I0d+du2−lΓ2α1+1Γ3α1+1,r1=−u1hS0Γ2α2+1Γα1+2α2+1,u=−u1I0f+fu2−lΓα1+α3+1Γ2α1+α3+1,v=iωΓα2+α3+1Γα1+α2+α3+1,w=jωΓ2α3+1Γα1+2α3+1,(42)I3=xtα1+2α2+yt3α2+zt2α1+α2+a1tα1+α2+α3,(43)x=u1S0g+eu1I0−gω+u2+u1ab−u2eΓα1+α2+1Γα1+2α2+1,y=u1S0h−hω+u2Γ2α2+1Γ3α2+1,z=du1I0−u2dΓ2α1+1Γ2α1+α2+1,a1=fu1I0−u2Γα1+α3+1Γα1+α2+α3+1,(44)R3=b1tα1+α2+α3+c1t2α2+α3+d1tα2+2α3+e1t3α3,(45)b1=gu2Γα1+α2+1Γα1+α2+α3+1,c1=hu2Γ2α2+1Γ2α2+α3+1,d1=−ωiΓα2+α3+1Γα2+2α3+1,e1=−ωjΓ2α3+1Γ3α3+1,(46)D3=f1tα1+α2+α4+g1t2α2+α4+h1t2α1+α4+i1tα1+α3+α4,(47)f1=gω+leΓα1+α2+1Γα1+α2+α4+1,g1=ωhΓ2α2+1Γ2α2+α4+1,h1=ldΓ2α1+1Γ2α1+α4+1,i1=lfΓα1+α3+1Γα1+α3+α4+1.SinceSt=∑k=0∞Skt,It=∑k=0∞Ikt,Rt=∑k=0∞Rkt,Dt=∑k=0∞Dkt, then, if we take the first four terms from the above series and take into consideration the equations from (30) to (46), we get
(48)St=S0+S1+S2+S3=S0+atα1+dt2α1+etα1+α2+ftα1+α3+rt2α1+α2+st3α1+r1tα1+2α2+ut2α1+α3+vtα1+α2+α3+wtα1+2α3,(49)It=I0+I1+I2+I3=I0+btα2+gtα1+α2+ht2α2+xtα1+2α2+yt3α2+zt2α1+α2+a1tα1+α2+α3,(50)Rt=R0+R1+R2+R3=R0+ctα3+itα2+α3+jt2α3+b1tα1+α2+α3+c1t2α2+α3+d1tα2+2α3+e1t3α3,(51)Dt=D0+D1+D2+D3=D0+d1tα4+ptα2+α4+qtα1+α4+f1tα1+α2+α4+g1t2α2+α4+h1t2α1+α4+i1tα1+α3+α4.After the compensation by the values ofa,b,c,d1,d,e,f,g,h,i,j,p,q,r,r1,s,u,v,w,x,y,z,a1,b1,c1,d1,e1,f1,g1,h1,i1, we obtain St=S0+S1+S2+S3=(52)S0+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1tα1+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1t2α1−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1tα1+α2+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1tα1+α3+−u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1S0−u1I0−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1−u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1G0−ω+u2I0−u2S0Γα2+1+−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1u2−lΓα1+α2+1Γ2α1+α2+1t2α1+α2+−u1I0−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1u2−lΓ2α1+1Γ3α1+1t3α1+−u1I0u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1u2−lΓα1+α3+1Γ2α1+α3+1t2α1+α3+u2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1ωΓα2+α3+1Γα1+α2+α3+1tα1+α2+α3+−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1ωΓ2α3+1Γα1+2α3+1tα1+2α3−u1u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1S0Γ2α2+1Γα1+2α2+1tα1+2α2,(53)It=I0+I1+I2+I3=I0+u1G0−ω+u2I0−u2S0Γα2+1tα2+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1tα1+α2+u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1t2α2+u1S0−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1+−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1u1I0−−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1ω+u2+u1−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1G0−ω+u2I0−u2S0Γα2+1−u2−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1Γα1+α2+1Γα1+2α2+1tα1+2α2+u1S0u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1−u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1ω+u2Γ2α2+1Γ3α2+1t3α2+u1I0−u2−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1Γ2α1+1Γ2α1+α2+1t2α1+α2+u2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1u1I0−u2Γα1+α3+1Γα1+α2+α3+1tα1+α2+α3,(54)Rt=R0+R1+R2+R3=R0+u2I0−ωR0Γα3+1tα3+u2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1tα2+α3−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1t2α3+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1u2Γα1+α2+1Γα1+α2+α3+1tα1+α2+α3+u1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1u2Γ2α2+1Γ2α2+α3+1t2α2+α3−ωu2u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α3+1Γα2+α3+1Γα2+2α3+1tα2+2α3−ω−ωu2I0−ωR0Γα3+1Γα3+1Γ2α3+1Γ2α3+1Γ3α3+1t3α3,(55)Dt=D0+D1+D2+D3=D0+ωI0+l−μN+lS0Γα4+1tα4+ωu1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα2+α4+1tα2+α4+l−u1G0+μ−lN+ωR0+u2−lS0Γα1+1Γα1+1Γα1+α4+1tα1+α4+−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u1I0−u2Γα1+1Γα1+α2+1ω+l−u1S0u1G0−ω+u2I0−u2S0Γα2+1Γα2+1Γα1+α2+1Γα1+α2+1Γα1+α2+α4+1tα1+α2+α4+ωu1G0−ω+u2I0−u2S0Γα2+1u1S0−ω+u2Γα2+1Γ2α2+1Γ2α2+1Γ2α2+α4+1t2α2+α4+l−u1G0+μ−lN+ωR0+u2−lS0Γα1+1u2−l−u1I0Γα1+1Γ2α1+1Γ2α1+1Γ2α1+α4+1t2α1+α4+lu2I0−ωR0Γα3+1ωΓα3+1Γα1+α3+1Γα1+α3+1Γα1+α3+α4+1tα1+α3+α4.Now, by taking the first three terms ofSt, It, Rt, and Dt and putting α1=α2=α3=α4=α,u1=0.4, μ=0.3, N=1000, l=0.2, ω=0.8, u2=0.2, S0=900, I0=10, R0= D0=0, in the proposed model yield (see [30]).
(56)St=S0+S1+S2=900−3500Γα+1tα−1213598.4Γ2α+1t2α,(57)It=I0+I1+I2=10+3410Γα+1tα+1210890Γ2α+1t2α,(58)Rt=R0+R1+R2=0+2Γα+1tα+680.4Γ2α+1t2α,(59)Dt=D0+D1+D2=0+88Γα+1tα+2028Γ2α+1t2α.Figure2 exhibits the dynamical behavior of St,It,Rt, and Dt against various fractional orders (α=0,7;0.8;0.9;1).Figure 2
The behavior ofSt,It,Rt, and Dt populations versus time with different fractional derivative orders.


## 4.2. Convergence and Error Analysis

It is clear that the previous solution is a series. This solution converges rapidly and uniformly to the exact solution. We will use classical techniques to prove the convergence of the series (21), and hence, we can get the sufficient conditions for the convergence of the used method.Theorem 4.
Suppose℘ is the Banach space and σ:℘⟶℘ is a contractive nonlinear operator with for all x,x′∈℘,σx−σx′≤λx−x′,0<λ<1. Then, based on Banach contraction principle, there is a unique point x such that σx=x, where x=s,i,r,d. By applying the ADM, the series in (21) can be written as follows:
(60)xn=σxn‐1,xn‐1=∑k=1n‐1xk,n=1,2,3,⋯,and suppose that x0=x0∈Ajx where Ajx=x′∈℘:x−x′<j, then, we obtain
(61)xn∈Ajx;(62)limn⟶∞xn=x.Proof.
For the first partaxn∈Ajx; by the mathematical induction if n=1, we obtain
(63)x0−x=σx0−σx≤λx0−x.
Suppose the result is true atn−1, hence,
(64)x0−x≤λn−1x0−x.
We have(65)xn−x=σxn−1−σx≤λxn−1−x≤λnx0−x,(66)xn−x≤λnx0−x≤λnj<j.This shows that xn∈Ajx..☐For the second part (b)limn⟶∞xn=x.Sincexn−x≤λnx0−x and limn⟶∞λn=0. So, we obtain limn⟶∞xn−x=0⇒limn⟶∞xn=x.


## 5. Fractional-Order Numerical Simulations of the SIRD System

First, we recall the basics of the applied numerical technique that have been given for the numerical simulation of fractional IVPs with Caputo derivatives (2). The technique is an extension of the well-known Adams-Bashforth-Moulton (ABM) integrator which is common for the numerical simulation of 1st order differential equations [35]. This technique comes from the idea that the IVP is equal to the equation of Volterra integral. The fractional ABM technique is fully presented by the following formulas (all other states can be found same as S). Assuming that 0,T is the domain of the solution and n=0,1,2,⋯,N, h=T/N, tn=nh:
(67)Stn+1=∑i=0α−1S0itn+1i+1i!+hαΓα+2ψ1tn+1,Sptn+1+hαΓα+2∑j=0nΛj,n+1ψ1tj,Stj,Sptn+1=∑i=0α−1S0itn+1i+1i!+1Γα∑k=0nΞk,n+1ψ1tk,Stk,where
(68)Λj,n+1=nα+1−n−αn+1α,ifj=0,n−j+2α+1+n−jα+1−2n−j+1α+1,if1≤j≤n,1ifj=n+1,(69)Ξj,n+1=hααn+1−jα−n−jα.A complete study of the stability of the method has been produced in [36].In this part, we have found the numerical simulation of the proposed paradigm via the predictor-corrector PECE method of ABM [35] coded in MATLAB. The behavior of the obtained solution is shown by taking suitable values of the parameters as in [30]: u1=0.4, μ=0.3, N=1000, l=0.2, ω=0.8, u2=0.2, S0=900, I0=10, R0= D0=0, and α1=α2=α3=α4=α. In Figure 3, the behavior of St and It populations versus time with different fractional derivative orders is shown. In Figure 4, the behavior of Rt and Dt populations versus time with different fractional derivative orders is shown. In Figure 5, the dynamics of the populations in It,Rt and Dt 3D space with different fractional derivative orders are shown. In Figure 6, the dynamics of the populations in St,It and Dt 3D-space with different fractional derivative orders are shown. In Figure 7, the dynamics of the populations in St,It, and Rt 3D space with different fractional derivative orders are shown. In Figure 8, the dynamics of the populations in St,Rt, and Dt 3D space with different fractional derivative orders are shown. In Figure 9, the Lyapunov exponents concerning the fraction derivative order α are displayed. From Figure 9, we can conclude that the system is unstable in the transient period since LE1 is positive. But as shown LE1 loses its positivity after a short period which indicates that the system settles to its fixed point. In Figure 10, we show the 3D dynamics of the infected population concerning both time and infection rate at fractional derivative order 0.95. In Figure 11, we show the 3D dynamics of the susceptible population concerning both time and infection rate at fractional derivative order 0.95. There is a quick convergence between the obtained and the exact solution. In complex models, it is well known that a small variation in physical behavior stimulates several new results to understand and analyze nature in a systematic manner.Figure 3
The behavior ofSt and It populations versus time with different fractional derivative orders.Figure 4
The behavior ofRt and Dt populations versus time with different fractional derivative orders.Figure 5
The dynamics of the populations inIt,Rt, and Dt 3D space with different fractional derivative orders.Figure 6
The dynamics of the populations inSt,It, and Dt 3D space with different fractional derivative orders.Figure 7
The dynamics of the populations inSt,It, and Rt 3D space with different fractional derivative orders.Figure 8
The dynamics of the populations inSt,Rt, and Dt 3D space with different fractional derivative orders.Figure 9
The exponents of Lyapunov as a function of the fractional derivative orderα.Figure 10
3D dynamics of the infected population versus both time and infection rate at fractional derivative order 0.95.Figure 11
3D dynamics of the susceptible population versus both time and infection rate at fractional derivative order 0.95.


## 6. Fractional Optimal Control (FOCP) for SIRD Paradigm

In this position, we discuss the optimal control of the fractional-order SIRD paradigm:(70)DαSt=−u1SI+μ−1N+ωR+u2−lS,DαIt=u1SI−ω+u2I−u2S,DαRt=u2I−ωR,DαDt=ωI+l−μN+lS.According to references [37–40], consider the model (70) inside R4. The proposed control functions can be formulated as follows:
(71)Ψ=u1.,u2.∈L∞×L∞0≤u1.,u2.≤1,∀t∈0,Tf=0,1,and Tf refers to the end of the time period, u1. provides the rate of susceptibility per year, and u2. provides the rate of curing from infection population to recovered population per year. u1. with u2. is the control functions. Here, the function space L∞ norm is defined as x∞=maxixi. This concept of the L∞ norm considered to be equal to taking the limit as p⟶∞ of the Lp norm.The following is the definition of the objective function (the control variables are squared)(72)Ju1,u2=∫0TfaIt+bu12t+cu22tdt.Such thata,b, and c clarify, respectively, the quantity of infectious, a per year rate of susceptible, and a per year rate of progress from infective to recovered.Finding the optimal controlsu1.andu2. is the basic goal in FOCPs, in order to minimize the predefined objective function:
(73)Ju1,u2=∫0TfψS,I,R,D,u1,u2,tdt,where ψS,I,R,D,u1,u2,t=aIt+bu12t+cu22t, subjected to the constraint
(74)DαS=η1,DαI=η2,DαR=η3,DαD=η4,where ηi=ηiS,I,R,D,u1,u2,t,i=1,2,3,4.The initial conditions are as follows:(75)S0=S1,I0=I1,R0=R1,D0=D1.To construct the FOCP, consider the selected objective function defined in [41–44]:
(76)J∗=∫0TfHS,I,R,D,u1,u2,t−∑i=14λiηiS,I,R,D,u1,u2,tdt,Such thati=1,2,3,4.In equation (77), (78), we write the Hamiltonian of the cost (objective) function (76) and the control of the fractional-order SIRD paradigm (70):
(77)HS,I,R,D,u1,u2,t=ψS,I,R,C,u1,u2,t+∑i=14λiηiS,I,R,D,u1,u2,t,then
(78)H=aI+bu12+cu22+λ1−u1SI+μ−1N+ωR+u2−lS+λ3u2I−ωR+λ2u1SI−ω+u2I−u2S+λ4ωI+l−μN+lS.From (76) and (78), following the method in [37, 39, 43–45], the FOPC’s necessary and sufficient conditions are formulated as
(79)Dαλ1=∂H∂S,Dαλ2=∂H∂I,Dαλ3=∂H∂R,Dαλ4=∂H∂D,(80)∂H∂uk=0,k=1,2⇒∂H∂u1=0,∂H∂u2=0,(81)DαS=∂H∂λ1,DαI=∂H∂λ2,DαR=∂H∂λ3,DαD=∂H∂λ4,additionally,
(82)λi,Tf=0,where λi,i=1,2,3,4 represent the Lagrange multipliers. In equations (79) and (80), the necessary conditions using FOPC Hamiltonian are represented.The main result is recorded in Theorem5.Theorem 5.
Ifu1andu2 are optimal controls related to the controlled states S∗,I∗,R∗andD∗, then, the adjoint variables λi∗,i=1,2,3,4 exist that satisfy the following conditions:
(i)
Adjoint conditions (costate conditions)Considering the related conditions represented in the recorded theorem and using equations (82), see [37, 39, 43, 44], we obtain the following four conditions that can be formulated as
(83)Dαλ1∗=λ1∗−u1I∗+u2−l+λ2∗u1I∗−u2+λ4∗l,(84)Dαλ2∗=a+λ1∗−u1S∗+λ2∗u1S∗−ω+u2+λ3∗u2+λ4∗ω,(85)Dαλ3∗=λ1∗ω+λ3∗−ω,(86)Dαλ4∗=0.(ii)
Transversality conditions(87)λi∗Tf=0,i=1,2,3,4.(iii)
Conditions of optimality(88)HS∗,I∗,R∗,D∗,u1∗,u2∗,λi=min0≤u1∗,u2∗≤1H∗S∗,I∗,R∗,D∗,u1∗,u2∗,λi,moreover, by applying equations (80), the control functions u1∗,u2∗ are given as follows:
(89)∂H∂u1=0⇒u1=λ1∗−λ2∗S∗I∗2b,(90)∂H∂u2=0⇒u2=λ2∗−λ3∗I∗+λ2∗−λ1∗S∗2c,(91)u1∗=min1,max0,λ1∗−λ2∗S∗I∗2b,(92)v2∗=min1,max0,λ2∗−λ3∗I∗+λ2∗−λ1∗S∗2c.Proof.
The system of costate conditions (83)–(86) is produced from equation (81) such that the Hamiltonian H∗ is presented as
(93)H∗=aI∗+bu1∗2+cu2∗2+λ1∗DαS∗+λ2∗DαI∗+λ3∗DαR∗+λ4∗DαD∗.
Furthermore, the conditions in equation (82) are guaranteed, and the optimal control given by equations (91)–(92) can be produced via equation (80).☐Substitutingvi∗,i=1,2 in (70), the state system is formulated as follows:
(94)DαS∗t=−u1∗S∗I∗+μ−1N+ωR∗+u2∗−lS∗,DαI∗t=u1∗S∗I∗−ω+u2∗I∗−u2∗S∗,DαR∗t=u2∗I∗−ωR∗,DαD∗t=ωI∗+l−μN+lS∗.For more details about FOC, see [37, 39, 40, 44].The existence of the control system (83)–(86) is to be proved here, the same can be found in [37–43] as follows:Let(95)g1λ1∗,λ2∗,λ3∗,λ4∗=λ1∗−u1∗I∗+u2∗−l+λ2∗u1∗I∗−u2∗+λ4∗l,(96)g2λ1∗,λ2∗,λ3∗,λ4∗=a+λ1∗−u1∗S∗+λ2∗u1∗S∗−ω+u2∗+λ3∗u2∗+λ4∗ω,(97)g3λ1∗,λ2∗,λ3∗,λ4∗=λ1∗ω+λ3∗−ω,g4λ1∗,λ2∗,λ3∗,λ4∗=0.LetΩ=λr∗∈R:λr∗≤k,r=1,2,3,4,t∈0,T.We have, atΩ:∂g1/∂λ1∗=−u1∗I∗+u2∗−l,∂g1/∂λ2∗=u1∗I∗−u2∗,∂g1/∂λ3∗=0,∂g1/∂λ4∗=l, ∂g2/∂λ1∗=−u1∗S∗, ∂g2/∂λ2∗=u1∗S∗−ω+u2∗, ∂g2/∂λ3∗=u2∗, ∂g2/∂λ4∗=ω,(98)∂g3∂λ1∗=ω,∂g3∂λ2∗=0,∂g3∂λ3∗=−ω,∂g3∂λ4∗=0,(99)∂g4∂λ1∗=∂g4∂λ2∗=∂g4∂λ3∗=∂g4∂λ4∗=0,(100)∂g1∂λ3∗≤k1⇒g1=g13,∂g1∂λ4∗≤g1⇒g1=k13,∂g1∂λ1∗≤k11,∂g1∂λ2∗≤k12,(101)∂g2∂λ1∗≤k21,⋯∂g2∂λ3∗≤k23,∂g2∂λ4∗≤k24,∂g2∂λ2∗≤k22,⋯(102)∂g3∂λ1∗≤k2,∂g3∂λ2∗≤k2⇒g3=k31,∂g3∂λ3∗≤k33,∂g3∂λ4∗≤k34,(103)∂g4∂λ1∗≤k41,∂g4∂λ3∗≤k43,∂g4∂λ4∗≤k44,∂g4∂λ2∗≤k42,where k1,k2,k3,k4,k5k22,k23,k24,k2,k31,k33,k34,k41,k42,k43, and k44 are constants greater than zero. Hence, each of the four functions f1,f2, and f4 agrees with the Lipschitz z conditions. Then, each of the four kernels f1,f2 is continuous respecting the four features, then. See for instance [37–43] for more information about the existence and uniqueness.


## 7. Conclusion

In this work, we have provided the approximate solution constructions for the fractional-order LHF model. In addition, the numerical simulations of the fractional-order LHF model have been experimented with many arbitrary orders. Also, an optimal control case has been investigated. The outcomes acquired by utilizing Laplace transform combined with the ADM are very compatible with results existing in the previous research. Further, the applied method for finding the numerical solution of the proposed model is introduced without considering any perturbation, discretization, or transformations. The current examination enlightens the considered nonlinear manner that rely upon the time history and the time moment by applying the idea of fractional order differentiation. Moreover, the obtained results show that Laplace transform combined with the ADM is amazingly deliberate, more powerful, and exceptionally precise, and which can be applied to break down the different classes of nonlinear marvels that emerge in science and innovation. For future work, we suggest generalizing the studied LHF model to show the effect of the infection from other sources as animals, men, and children.



---

*Source: 2444920-2021-12-15.xml*