Contents lists available atScienceDirect
Applied Numerical Mathematics
www.elsevier.com/locate/apnum
Operator splitting for space-dependent epidemic model
Petra Csomós
∗, Bálint Takács
InstituteofMathematics,EötvösLorándUniversityandMTA–ELTENumericalAnalysisandLargeNetworksResearchGroup,PázmányPéter sétány1/C,H-1117Budapest,Hungary
a rt i c l e i nf o a b s t ra c t
Articlehistory:
Received3December2019
Receivedinrevisedform15September 2020
Accepted16September2020 Availableonline21September2020
Keywords:
Space-dependentepidemicmodel Operatorsplittingprocedures Qualitativeproperties Non-negativitypreservation
Wepresentandanalysenumericalmethodswithoperatorsplittingprocedures,appliedto anepidemicmodel whichtakesintoaccountthespace-dependenceoftheinfection.We deriveconditionsonthetimestep,underwhichthenumericalmethodspreservethenon- negativityandmonotonicitypropertiesoftheexactsolution.Ourresultsareillustratedby numericalexperiments.
©2020TheAuthors.PublishedbyElsevierB.V.onbehalfofIMACS.Thisisanopenaccess articleundertheCCBY-NC-NDlicense (http://creativecommons.org/licenses/by-nc-nd/4.0/).
1. Introduction
Epidemicmodellingplaysanincreasinglyimportantrolenotonlyinappliedmathematicsbutalsoinmedicineandpublic health.Thereis,forinstance,ahighdemandonplanningtherightplaceandtimeofvaccination.Themorecomplexthese models are,the lesshope we haveinobtaining their analytical solution.Thus, thederivation andanalysisofbiologically adequatenumericalmethodsmeansavitalchallenge.
Epidemicmodels originatefromtheseminal workofKermackandMcKendrick[9] publishedin1927,whoconstructed a compartment model to study the process of epidemic propagation. The population is split into three classes: healthy butsusceptibleindividuals,infectedpeople whocaninfect otherindividuals, andalreadyrecovered orotherwiseimmune individuals.Thefirstattemptsdescribetwowaystheindividualscan“change”classes:(i)susceptibleindividualsgetinfected with some possibility, and (ii) infectedpeople recoverwith some other rateof change. There are several directions the original model can be generalised: by considering birth and death processes, by addingmore classes of individuals, by consideringalatentperiod,or,asinthepresentpaper,bytakingintoaccounttheeffectofvaccination.Inthepresentwork we analyseanepidemicmodelwhichalsotreats thespace-dependencyoftheeffectoftheinfection,that is,thedistance betweenthesusceptibleandinfectedindividuals.
Thenovelty ofourwork istoapplyoperatorsplittingprocedures whendiscretisingintime. Theyallowustosplitthe modelintotwosub-problems,andsolvethemoneaftertheother.Withthehelpofoperatorsplitting,thedifficultycaused by thespace-dependencycan behandledseparately.Moreover,theexactsolutionoftheremainingpartcanbe computed leadingtoamorestableandaccuratenumericalsolution.
Since epidemicmodels describe real-life phenomena, it isvital to study whetherthey reflect theproperties expected fromthebiologicalpointofview.Therefore,weputanefforttoinvestigateunderwhichconditionsthemodel’snumerical
*
Correspondingauthor.E-mailaddress:csomos@cs.elte.hu(P. Csomós).
https://doi.org/10.1016/j.apnum.2020.09.010
0168-9274/©2020TheAuthors.PublishedbyElsevierB.V.onbehalfofIMACS.ThisisanopenaccessarticleundertheCCBY-NC-NDlicense (http://creativecommons.org/licenses/by-nc-nd/4.0/).
solutionowns thenon-negativityandmonotonicityproperties.Wearealsoafterthecaseswhenourmethodgiveshigher boundonthetimestepthantheonealreadypresentedintheliterature.Weillustrateourtheoreticalresultsbynumerical experiments.
Section 2givesan overviewon basicepidemicmodelsandshowshow wetreat thespace-dependencyofinfection.In Section3weintroducethespaceandtimediscretisationmethodswhichareusedlater.InSection4wedefinethequalitative propertiestobeinvestigatedintherestofthepaper.Section5containssome necessarytechnicaltools.Sections6,7,and 8 aredevoted tothe analysisofthesequential,weighted, andStrangsplittings, respectively.In Section 9we presentour numericalexperimentsillustratingthetheoreticalresults.Section10brieflysummarisesourresults.
2. Space-dependentSIRmodel
Mostofthecurrentlyused andanalysed modelsare derived fromthe ideaofKermack andMcKendrick[9],who con- structed thecompartment model,introducedabove, tostudy theprocess ofepidemicpropagation.Let S
,
I,
R:R+0 →R+0 denotethedensityofsusceptible, infected,andrecoveredindividualsamongthetotalpopulation,respectively, andletthe constant parametersa,
b>
0 describetherateofinfectionandrecovery,respectively.Let S0,
I0,
R0≥0 be givennumbers.Thenthesusceptible–infected–recovered(SIR)epidemicmodelhastheform
⎧ ⎪
⎨
⎪ ⎩
S
(
t) = −
aS(
t)
I(
t),
I(
t) =
aS(
t)
I(
t) −
bI(
t),
R(
t) =
bI(
t)
(1)
forallt
>
0 withtheinitialconditionS
(
0) =
S0,
I(
0) =
I0,
R(
0) =
R0.
(2)The SIRmodel(1)–(2) is aninitial valueproblembeinga systemofthree ordinarydifferentialequations.Althoughmodel (1) alreadydescribesseveralimportantfeaturesrealepidemicsposses,itdoesnottakeintoaccountthespatialdistribution ofthedifferentspecies,andsupposesahomogeneousdistributiononthedomain.Amodelwhichconsiderstheaforemen- tionedpropertieswasintroducedbyKendall[8] inthefollowingway.
Foranarbitrarym∈N,weconsideraboundeddomain
⊂Rmandtheopenball Bδ
(
x)
aroundthepoint x∈with radius
δ >
0.Let|Bδ(
x)
|denoteitsLebesguemeasure(orvolume),and XBδ(x)(
t)
thenumberofindividualsinthisballfor each X∈ {S,
I,
R}attimet≥0.Thenthedensityofclass X atpointx∈andattimet
>
0 isdefinedas X(
t,
x) :=
limδ→0 1
|
Bδ(
x) |
XBδ(x)(
t).
Toeasethenotation,wewillomitthetilde,anddenotethedensitybyX
(
t,
x)
foreach X∈ {S,
I,
R}.Theconsiderationabove leadstoaspace-dependentSIRmodelwhichis,however,atthispointnotsobeneficialbecausethedensityfunctionsbehave independentlyateachpoint x∈.Sincetheinfectiontakesplacepointwise,itcannotspreadinspace,beinghoweverthe maingoalofthegeneralization.Thus,itismorenaturaltosupposethattheinfectedindividualshavean influenceonthe susceptibles inacertain distancearoundthemselves insuch awaythat theylesslikelyinfect healthyindividualsfurther away fromthemselves.Thatis,a susceptiblecan getinfectedonlyinapredefineddomain,e.g.,a circle.Wenote thatthe radius
δ >
0 oftheinfectiousdomaincanvary dependingonthediseaseconsidered.We furthersupposethatthedisease processisthesameateverypointx∈.
Sinceit isthemostcommonway,we alsoformulate ourmodelintwodimensionsandsuppose arectangulardomain
= [0
,
A]× [0,
B]⊂R2 withsome A,
B>
0 arbitrarynumbers(althoughtheresultsofthispapercanbeextendedtomore generaldomains).Aroundthepointx=(
x,
y)
∈,wedenotetheinfectiousdomainbyBδ
(
x,
y)
,beingthecirclewithorigin(
x,
y)
andradiusδ >
0.Tothisend,letr≥0 denoteanarbitrarypoint’sdistancefrom(
x,
y)
andϑ
∈ [0,
2π )
itsangle.The mainideaistoreplace I(
t)
intheterms±S(
t)
I(
t)
inthespace-dependentversionofmodel(1) bythefollowingweighted integralontheball Bδ(
x,
y)
:I
(
t,
x,
y) :=
∞ 02π
0
G
(
x,
y,
r, ϑ)
I(
t,
x+
rcos(ϑ),
y+
rsin(ϑ))
rdϑ
dr (3)where function G:
×R+0 × [0
,
2π )
→R+0 describesthe disease process in.Since we aim atimitating theeffect of an infectedindividual at the point
(
x,
y)
∈on its
δ
-radius neighbourhood Bδ(
x,
y)
, we want function G to represent the combinedeffect of(i) thenon-negative andmonotonically decreasingfunction g1: [0, δ]
→R+0 whichdescribes the dependenceontheradiusr,and(ii)thenon-negativefunction g2: [0,
2π )
→R+0 describingthedependenceontheangle.Forthesakeofsimplicity,webuildtheinfectiousratea
>
0 intofunctiong1.Weremarkthatthecaseofconstantfunction g2 iswidelystudiedin[5] and[6].Anon-constant g2 maybeusefulformodellingthespreadofdiseaseswhenthereisaconstant windblowinginonedirectionwhichwas describedin[14].Wealsosuppose thatthefunction g2 is periodicin thewayg2
(
0)
= limϑ→2πg2
(ϑ)
.Sinceitisanaturalassumption,wetakethefunctionG beingseparableinr andϑ
: G(
x,
y,
r, ϑ) =
g1(
r)
g2(ϑ),
if x+
rcos(ϑ),
y+
rsin(ϑ)
∈
Bδ(
x,
y),
0
,
otherwise. (4)ThenthetermIinrelation(3) hasthefollowingform:
I
(
t,
x,
y) =
δ 02π
0
g1
(
r)
g2(ϑ)
I(
t,
x+
rcos(ϑ),
y+
rsin(ϑ))
rdϑ
dr (5)forallt≥0,
(
x,
y)
∈.
Toconsideramorerealisticmodelthan(1),wetakeintoaccounttheeffectofvaccinationaswell.Letc
>
0 denotethe raterelatedtothevaccinatedpopulationgettingimmune.ThenforthenewunknownfunctionsS,
I,
R:R+0 ×→R+0,we getthefollowingsystemofintegro-differentialequations:
⎧ ⎪
⎨
⎪ ⎩
∂
tS(
t,
x,
y) = −
S(
t,
x,
y)
I(
t,
x,
y) −
c S(
t,
x,
y),
∂
tI(
t,
x,
y) =
S(
t,
x,
y)
I(
t,
x,
y) −
bI(
t,
x,
y),
∂
tR(
t,
x,
y) =
c S(
t,
x,
y) +
bI(
t,
x,
y)
(6)
forallt≥0,
(
x,
y)
∈,andwiththeinitialcondition
S
(
0,
x,
y) =
S0(
x,
y),
I(
0,
x,
y) =
I0(
x,
y),
R(
0,
x,
y) =
R0(
x,
y),
(7) where S0,
I0,
R0:→R+0 aregivencontinuousfunctionssuchthat
S0
(
x,
y) +
I0(
x,
y) +
R0(
x,
y) =
0 holds for all(
x,
y) ∈ .
(8) Asalreadymentioned,onecannothopetofindananalyticalsolutiontosystem(6),althoughitwasprovedin[13] thatsuch asolutionexists,whichisalsounique.Therefore,weare goingtousenumericalmethodstosolvetheseequations.Inthe nextsectionweintroducethespaceandtimediscretisationmethodstobeusedinthepresentstudy.3. Discretisationmethods
Thepresentsectionaimsatintroducingthespaceandtimediscretisationmethodsofthespace-dependentSIRmodel(6) aswellastheoperatorsplittingprocedures.Sinceitisthemostchallengingpartofthenumericalmethodbeingconstructed, firstweshowhowweapproximatetheintegralappearingin(5).
3.1. Approximatingtheintegral
Thekeypointofthenumericalsolutionofproblem(6) istheapproximationofthedoubleintegralin(5),whichcanbe doneindifferentways.OneapproachistoapproximatethefunctionI
(
t,
x+rcos(ϑ),
y+rsin(ϑ))
byaTaylorexpansion:the obtainedmethodisstudiedin[5] and[6].Wenotethatthisprocessisnotefficientinthecaseofnon-constantfunction g2 asshownin[14].Theotherapproachistouseacombinationofinterpolationandnumericalintegration(byusingcubature formulas).Forthepresentstudyweimplementthesecondapproach.WeconsideratwodimensionalcubatureformulaonthediscBδ withpositivecoefficients.ForindexsetJ ⊂N2andfor all
(
i,
j)
∈J,letri∈ [0, δ
]denotethe(
i,
j)
th cubaturepoints’distancefromthecentre point(
x,
y)
∈,and
ϑ
j∈ [0,
2π )
its angle.ThenQ(x,
y)
denotesthesetofcubaturepoints inthedisk Bδ(
x,
y)
parametrizedby polarcoordinates(see[14] or [13]):Q
(
x,
y) =
x+
ricos(ϑ
j),
y+
risin(ϑ
j)
∈
IntBδ(
x,
y), (
i,
j) ∈
J,
where Int denotes the interior of the set. Numerical integration leads then to the following approximation of the term I
(
t,
x,
y)
in(6):T
(
t,
Q(
x,
y)) =
(i,j)∈J
wi,jg1
(
ri)
g2(ϑ
j)
I t,
x+
ricos(ϑ
j),
y+
risin(ϑ
j)
(9)
withsome weights wi,j≥0. Fortheinfectedindividuals beingclosertothe boundary ofthedomain
asthe radius
δ
, theapproximationoftheintegralin(5) needsvaluesofIlyingoutside:forthese,wearegoingtousezerovalues.After theseconsiderationswegetthefollowingsystem,beingstillcontinuousint≥0 and
(
x,
y)
∈:
⎧ ⎪
⎨
⎪ ⎩
∂
tS(
t,
x,
y) = −
S(
t,
x,
y)
T(
t,
Q(
x,
y)) −
c S(
t,
x,
y),
∂
tI(
t,
x,
y) =
S(
t,
x,
y)
T(
t,
Q(
x,
y)) −
bI(
t,
x,
y),
∂
tR(
t,
x,
y) =
c S(
t,
x,
y) +
bI(
t,
x,
y)
(10)
withtheoriginalinitialcondition(7).Wenotethat thereareseveralpossibilities howtochoosethequadratures.Onecan useadirectmethodwhichresultsinauniformcubature,ortransformtheballontoasquare,andusegeneralisedGaussian quadraturesonit.Theresultsin[13] showthatforlessquadraturepoints,theuniformonesresultinasmallererror,while foradenser quadrature,thenon-uniformonesperformbetter. Onecanalsocompute theconvolution in(5) by usingthe FastFourierTransform,whichcanbethedirectionoffurtherresearch.
3.2. Spatialdiscretisation
Inordertodiscretiseproblem(10) inspace,weneedaspatialgridG onthedomain
= [0
,
A]× [0,
B].Tothisendwe choose thearbitrarynumbers K,
L∈Z+ anddefine thegrid resolutions hx:=A/(
K−1) >
0 and hy:=B/(
L−1) >
0 in directionsxand y,respectively.Thenthegriditselfisthefollowingset:G
:=
(
xk,
y) ∈ |
xk= (
k−
1)
hx,
y= ( −
1)
hy,
k=
1, . . . ,
K, = , . . . ,
L.
Forallt≥0,
(
xk,
y)
∈G,and X∈ {S,
I,
R},weconsidertheapproximatenumbersXk,
(
t) ≈
X(
t,
xk,
y)
and Tk,(
t) ≈
T(
t,
Q(
xk,
y)).
(11) InordertodeterminetheformofTk,(
t)
,wefirstproject I(
t,
x,
y)
tothegridG.Notethatthepoints(
xk+ricos(
j),
y+ risin(
j))
mightnotbepartofthegridG,sowecannot assignanyIk, valuestothem.Therefore,weapproximatethem byabilinearinterpolatingmethodusingthenearestknownIk, valuesandpositivecoefficients,resultinginthenotationI.Thenwehave Tk,
(
t) :=
(i,j)∈J
wi,jf1
(
ri)
f2(
j)
I t,
xk+
ricos(
j),
y+
risin(
j)
.
(12)It is worth mentioning that higher order interpolations, like cubic and spline, can be also used. Although they do not preserve non-negativity, fora sufficiently smallspatial grid resolutionthey behave asexpected. It isalso possibleto use other,highorder,non-negativitypreservingmethods,see[13].WenoteherethatifI≥0 holds,thenTk,
(
t)
isnon-negative too,forallt≥0 and(
xk,
y)
∈G.Inordertoeasethenotation,wewillleavethetildethroughoutourcomputations.3.3. Timediscretisation
Themainnoveltyofthepaperisthat(besidesthetraditionaltimediscretisation)weuseanothertimediscretisation-like method:operatorsplitting.Asonecansee,theright-handsideofproblem(10) canbewrittenasasumoftwoterms:one containingtheintegralandonewiththeremainingterms.Theideaofoperatorsplittingisto“split”theproblemintotwo sub-problems withthe correspondingterms alone,andsolve themseparately byusing anappropriate initialcondition to linktheirsolutionstogether.Inthepresentpaperwewillintroduceandstudythesequential,thesequentialweighted,and theStrangsplittingschemes.
Asalreadymentioned,itisnaturaltosplitthespace-discretisedSIRmodel(10) intothesub-problemswithandwithout theintegraltermIspecifyingthespace-dependencyoftheinfectionprocess:
⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
∂
tS[1](
t,
x,
y) = −
c S[1](
t,
x,
y),
∂
tI[1](
t,
x,
y) = −
bI[1](
t,
x,
y),
∂
tR[1](
t,
x,
y) =
bI[1](
t,
x,
y) +
c S[1](
t,
x,
y)
(Sub.1)
and
⎧
⎪ ⎪
⎨
⎪ ⎪
⎩
∂
tS[2](
t,
x,
y) = −
S[2](
t,
x,
y)
I[2](
t,
x,
y),
∂
tI[2](
t,
x,
y) =
S[2](
t,
x,
y)
I[2](
t,
x,
y),
∂
tR[2](
t,
x,
y) =
0(Sub.2)
forallt≥0 and
(
x,
y)
∈.Thelinkbetweenthesub-problemsistheinitialcondition,aswillbeshowninthenextsections.
Forthelateruseweremarkthatsub-problem(Sub.1) canbesolvedexactly:
⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
S[1]
(
t+ τ ,
x,
y) =
e−cτS[1](
t,
x,
y),
I[1](
t+ τ ,
x,
y) =
e−bτI[1](
t,
x,
y),
R[1]
(
t+ τ ,
x,
y) =
R[1](
t,
x,
y) + (
1−
e−cτ)
S[1](
t,
x,
y) + (
1−
e−bτ)
I[1](
t,
x,
y)
(13)
forallt≥0 and
(
x,
y)
∈,where
τ
≥0 isanarbitrarytimedifference.Ontheotherhand,sub-problem(Sub.2) cannotbe solvedexactly.Its approximatesolutioncanbe obtainedbyanother timediscretisationmethod.Forinstance,theuseofthefirst-orderexplicitEulermethodwithtimestep
τ >
0 leadsto⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
S[2]
((
n+
1) τ ,
x,
y) =
S[2](
nτ ,
x,
y) − τ
S[2](
nτ ,
x,
y)
I[2](
nτ ,
x,
y),
I[2]((
n+
1) τ ,
x,
y) =
I[2](
nτ ,
x,
y) + τ
S[2](
nτ ,
x,
y)
I[2](
nτ ,
x,
y),
R[2]((
n+
1) τ ,
x,
y) =
R[2](
nτ ,
x,
y)
(14)
foralln∈Nwith X[2]
(
0,
x,
y)
=X0(
x,
y)
foreach X∈ {S,
I,
R}.Wenotethatwetake0∈N.Theuseofthesecond-orderHeun’smethodinShu–Osherform(whichpreservesthestrongstability,see[7])withtime step
τ >
0 resultsinthefollowingsteps:⎧ ⎪
⎨
⎪ ⎩
S[2]((
n+
1) τ ,
x,
y) =
S[2](
nτ ,
x,
y) − τ
S[2](
nτ ,
x,
y)
I[2](
nτ ,
x,
y),
I[2]((
n+
1) τ ,
x,
y) =
I[2](
nτ ,
x,
y) + τ
S[2](
nτ ,
x,
y)
I[2](
nτ ,
x,
y),
R[2]((
n+
1) τ ,
x,
y) =
R[2](
nτ ,
x,
y),
(15)
⎧ ⎪
⎪ ⎪
⎪ ⎪
⎪ ⎪
⎨
⎪ ⎪
⎪ ⎪
⎪ ⎪
⎪ ⎩
S[2]
((
n+
1) τ ,
x,
y) =
12S[2](
nτ ,
x,
y)
+
12 S[2]((
n+
1) τ ,
x,
y) − τ
S[2]((
n+
1) τ ,
x,
y)
I[2](
n+
1)
nτ ,
x,
y) ,
I[2]((
n+
1) τ ,
x,
y) =
12I[2](
nτ ,
x,
y)
+
12 I[2]((
n+
1) τ ,
x,
y) + τ
S[2]((
n+
1) τ ,
x,
y)
I[2]((
n+
1) τ ,
x,
y) ,
R[2]((
n+
1) τ ,
x,
y) =
12R[2](
nτ ,
x,
y) +
12R[2]((
n+
1) τ ,
x,
y) =
R[2](
nτ ,
x,
y).
(16)
Wedonotplugformulae(15) into(16),becausethemethodwillbemoresuitableforanalysisinitspresentform.
We note that the use of an additional time discretisation inside one time step might also lead to a non-negativity preservingmethod.Thistechniquemaybeappliedonlyforonesub-problem.Thenthetimestepcouldbechosenindepen- dentlyoftheconstraintsbutrelatedtotheaccuracyofthescheme.Formoredetailsonsuchkindofadaptivetimestepping wereferto[4].
4. Qualitativeproperties
When combiningthespaceandtime discretisationmethods presentedintheprevious section, weobtain anumerical method represented by a system of algebraic equations. By denoting the approximation of X
(
nτ ,
xk,
y)
by Xkn, for all X∈ {S,
I,
R},the unknown values Xnk,+1 of these algebraic equations are computed with the help of Xnk, forall n∈N,(
xk,
y)
∈G andX∈ {S,
I,
R}.In what followswe list severalimportant propertieswhich reflect real-lifeexpectations. Inthe presentwork we will studywhetherthenumericalsolutionpossessesthem.
1. Byaddingtheequationsofsystem(10),one obtainsthat the totalsizeofthepopulationremains constantintime at eachspaceposition:
∂
tS(
t,
x,
y) + ∂
tI(
t,
x,
y) + ∂
tR(
t,
x,
y) =
0S
(
t,
x,
y) +
I(
t,
x,
y) +
R(
t,
x,
y) =:
C(
x,
y)
for all t≥
0 and(
x,
y) ∈ .
(17) Thisweexpectfromthenumericalsolutionaswell,i.e.,thatthereexistnumbers Nk, suchthat:Snk,
+
Ink,+
Rnk,=
Nk, for all n∈ N, (
xk,
y) ∈
G.
(P1) 2. SincefunctionsS,
I,
R denotedensities,theirvaluesshouldremainnon-negative:X
(
t,
x,
y) ≥
0 for all t≥
0, (
x,
y) ∈
and X∈ {
S,
I,
R} .
(18) Weexpectthesamefromthenumericalvaluesaswell:Xkn,
≥
0 for all n∈ N, (
xk,
y) ∈
G and X∈ {
S,
I,
R} .
(P2) 3. Sinceinfectedorrecoveredindividualscannotbesusceptibleagain,thefunction Sisnon-increasingintimeS
(
t,
x,
y) ≥
S(
t+ τ ,
x,
y)
for allt, τ ≥
0 and(
x,
y) ∈ .
(19) Thesameshouldholdforthenumericalvaluesaswell:Snk,
≥
Snk+,1 for alln∈ N, (
xk,
y) ∈
G.
(P3)4. Similarly,thedensityofrecoveredindividualscannotdecreaseintime:
R
(
t,
x,
y) ≤
R(
t+ τ ,
x,
y)
for all t, τ ≥
0 and(
x,
y) ∈ .
(20) WhichmeansforthenumericalvaluesthatRnk,
≤
Rnk,+1 for all n∈ N, (
xk,
y) ∈
G.
(P4)In [13] itwas shownthatthe properties(17)–(20) holdforthesystems(6) and (10).So ouraimis toconstructsuch numericalmethodswhichpreservetheirdiscreteversions(P1)–(P4).
5. Technicaltools
Beforethederivationandanalysisofthemethods,wecollectsomenotationsandtechnicaltoolswewilluselateron.
Notation5.1.
(i) Foreach X∈ {S
,
I,
R}weintroducethenotation Xn:= ((
Xnk,)
k=1,...,K,=1,...,L) ∈ R
K L×K L.
(ii) LetM:RK L×K L→RK L×K Ldenotetheboundedlinearoperator(representedbyamatrixinapplications)thatmapsIn toTn bytheruleTn=M
(
In)
.Furthermore,letM
:=
M∞·
S0+
I0+
R0∞inwhich
.
∞meansthemaximummatrixnormtakenelement-wise.Wenotethatcondition(8) impliesM>
0.(iii) LetW−1: [−1
/
e,
0)
→(−∞,
−1]andW0:(−
1/
e,
+∞)→(−
1,
+∞)denotethetwobranchesoftheLambert-W func- tion,thatis,theinverseofthemapx→xex.(iv) Forarbitraryp
,
q>
0,wedefinethesetT
p,q:=
0
, −
1pW0−
pq∪
−
1pW−1−
pq, +∞
⊂ R.
Furthermore,wedefine
T
0,q:= [
0,
1q) ⊂ R.
Thelatternotationmakessensebecauseofthefollowingconsideration.
Lemma5.2.WithNotation5.1,thelimit−1pW0 −pq−−−→p→0 1q holdsforarbitraryq
>
0.Proof. ItsufficestoshowthatW0
(
x)/
x−−→x→0 1 forx= −p/
q<
0.TheL’Hospitalrule,thederivativeoftheinversefunction, andtheidentityW0(
0)
=0 implythatxlim→0 W0
(
x)
x
=
limx→0W0
(
x) =
limx→0
1
eW0(x)
+
W0(
x)
eW0(x)=
1.
Remark5.3.Sincewewilluseitseveraltimesthroughoutthepaper,weanalysethesolutionx
<
0 toequationxex
= μ
(21)forsomeparameter
μ <
0.(i) For
μ <
−1/
e,thereisnosolutiontoequation(21).(ii) For
μ
= −1/
e,thereisonesolution:x1= −1.(iii) For
μ >
−1/
e,therearetwosolutions:x−1=W−1( μ )
andx0=W0( μ )
. Wealsoknowthat x−1≤x1= −1<
x0.Hence,fortheinequalityxex
≥ μ
(22)wehavethefollowingcases.
Fig. 1.Graph of functionx→xex. The horizontal lines indicate theμ-values−0.25 and−1/e.
(i) For
μ <
−1/
e,theinequality(22) holdsforevery x<
0.(ii) For
μ
= −1/
e,theinequality(22) holdsforevery x<
0 (wehavexex=μ
forx= −1).(iii) For
μ >
−1/
e,wehave:x<
x−1=W−1( μ )
orx>
x0=W−1( μ )
. Thegraphoffunctionx→xex isdepictedonFig.1.Inthenextsectionswewillpresenttheconditiononthetimestep
τ
underwhichthequalitativeproperties(P1)–(P4) holdforthevariousoperatorsplittingschemes.Weareespeciallyinterestedinthecaseswhentheapplicationofoperator splittingleadstolesssevereconditionthantheoneobtainedwithoutsplitting.In[14] theauthorsappliedthesamespace discretisationasshowedinSection3.2andtheexplicitEulermethodforthewholesystem(6) withouttakingintoaccount the vaccination(c=0).Theyfound that property (P1) was automaticallysatisfied,andproperties(P2)–(P4) held truefor timestepsτ
satisfyingτ ≤
min 1M
,
1 b.
Thecasec
>
0 wasstudiedin[13],andresultedinasimilarbound,namelyτ ≤
min 1 M+
c,
1b
.
(23)Fromnowon,theupperbound(23) willbeconsideredasareferencevalue,andwewillstudytheconditionsunderwhich theapplicationofoperatorsplittingproceduresleadstoahigherone.
6. Sequentialsplitting
Operatorsplittingisbasedontheideatosimplifytheproblembysplittingitintotwoormoresub-problemswhichare easiertosolveortreatnumerically.Sincethesub-problemsneedtobesolvedseparately,weshouldderiveawaytoconnect their solutions. Dependingon theserules,we distinguish severalsplittingmethods. Themostbasic oneis thesequential splitting (initiatedfirstin [1])when thesub-problems are solvedone afterthe other ona time interval oflength
τ >
0, alwaystakingthesolutionoftheprevioussub-problemasinitialconditionfortheactualone.Aswewillsee,theproperties ofthesequentialsplittingdependontheorderofthesub-problems,therefore,wewilltreatthetwocasesseparately.Another splitting procedure is derived when the solutions of the two types of sequential splittings are weighted by a parameter
∈ [0
,
1].This kind ofmethod is calledweighted sequential splitting,see in [3], and will be discussed in Section7.ThethirdoperatorsplittingtobediscussedinSection 8istheStrangsplitting(derivedin[12] and[11])solving three problems ina single time step: one with thefirst sub-problem over a time interval of lengthτ /
2, then withthe secondsub-problemonanintervaloflengthτ
,andfinallywiththefirstsub-problemagainonaτ /
2 interval.In what follows we analyse the splitting procedures in the light of whetherthey preserve the qualitative properties introducedinSection4.
6.1. Sequentialsplitting1–2
Firstwetreatthesequentialsplittinginthecasewhenthesub-problemsaretakenintheorder(Sub.1)–(Sub.2).Thenthe applicationofthesequentialsplittingmeansthat inasingle time stepwefirst solvesub-problem(Sub.1) whosesolution
(13) servesastheinitialconditiontosub-problem(Sub.2). Moreprecisely,we considerthefollowingiterationstepsforall n∈Nand
(
x,
y)
∈: (Sub.1) for all t
∈ (
nτ , (
n+
1) τ ]
with initial condition X[1]
(
nτ ,
x,
y) =
Xspl(
nτ ,
x,
y)
(Sub.2) for all t∈ (
nτ , (
n+
1) τ ]
with initial condition X[2]
(
nτ ,
x,
y) =
X[1]((
n+
1) τ ,
x,
y)
Xspl((
n+
1) τ ,
x,
y) :=
X[2]((
n+
1) τ ,
x,
y)
(24)
where Xspl
(
0,
x,
y)
=X0(
x,
y)
istheoriginalinitialvaluein(7) foreach X∈ {S,
I,
R}.Afterdiscretisingsub-problem(Sub.2) bytheexplicitEulermethod,anddiscretisinginspacesub-problems(13) and(14),wegetthefollowingtwosub-problems:⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
Sk[1,],n+1
=
e−cτSk[1,],n,
Ik[1,],n+1=
e−bτIk[1,],n,
Rk[1,],n+1
=
R[k1,],n+ (
1−
e−cτ)
S[k1,],n+ (
1−
e−bτ)
Ik[1,],n(25)
and
⎧
⎪ ⎪
⎨
⎪ ⎪
⎩
Sk[2,],n+1
=
S[k2,],n− τ
S[k2,],nTk[2,],n,
Ik[2,],n+1=
I[k2,],n+ τ
S[k2,],nTk[2,],n,
Rk[2,],n+1=
R[k2,],n.
(26)
Bytakingintoaccounttheinitialconditionsasstatedin(24),thesub-problemshavethefollowingformforalln∈N and givenSnk,
,
Ink,,
Rnk,:⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
S[k1,],n+1
=
e−cτSnk,,
I[k1,],n+1=
e−bτIkn,,
R[k1,],n+1
=
Rnk,+ (
1−
e−cτ)
Snk,+ (
1−
e−bτ)
Ink,,
(27)
⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
Snk,+1
=
Sk[1,],n+1− τ
S[k1,],n+1Tk[1,],n+1,
Ink,+1=
I[k1,],n+1+ τ
Sk[1,],n+1Tk[1,],n+1,
Rnk,+1=
Rk[1,],n+1.
(28)
Notation5.1(ii)andthelinearityofoperatorMimplythefollowingrelation:
Tk[1,],n+1
=
M(
Ik[2,],n) =
M(
Ik[1,],n+1) =
M(
e−bτInk,) =
e−bτM(
Ink,) =
e−bτTkn,.
(29) Bycombiningthesub-problems(27)–(28),andtherelation(29),wearriveatthenumericalscheme:⎧ ⎪
⎪ ⎨
⎪ ⎪
⎩
Snk,+1
=
e−cτSnk,(
1− τ
e−bτTkn,),
Ink,+1=
e−bτ(
Ikn,+ τ
e−cτSnk,Tkn,),
Rnk,+1
=
Rnk,+ (
1−
e−cτ)
Snk,+ (
1−
e−bτ)
Ikn,.
(30)
Inwhatfollowsweshowtheconnectionbetweenproperties(P1)–(P4),andinvestigatetheconditionsunderwhichthey arefulfilled.
Proposition6.1.Wehavethefollowingassertions.
(a) Property(P1)holdsforthenumericalmethod(30)withoutanyrestriction.
(b) Property(P3)and(P4)areconsequencesofproperty(P2).
Proof. (a) Property(P1) followsbyaddinguptheequationsofsystem(30).
(b) Since Tnk,≥0 holdsifInk,≥0,ande−bτ
>
0 inthefirstandthethirdequationsofsystem(30),wegetthatproperties (P3) and(P4) alsohold.Thisconcludestheproof.