• Nem Talált Eredményt

* PetraCsomós , BálintTakács Operator splitting for space-dependent epidemic model Applied Numerical Mathematics

N/A
N/A
Protected

Academic year: 2022

Ossza meg "* PetraCsomós , BálintTakács Operator splitting for space-dependent epidemic model Applied Numerical Mathematics"

Copied!
22
0
0

Teljes szövegt

(1)

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/).

(2)

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

,

R00 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 withtheinitialcondition

S

(

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}attimet0.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

)

,wedenotetheinfectiousdomainby

(

x

,

y

)

,beingthecirclewithorigin

(

x

,

y

)

andradius

δ >

0.Tothisend,letr0 denoteanarbitrarypoint’sdistancefrom

(

x

,

y

)

and

ϑ

∈ [0

,

2

π )

itsangle.The mainideaistoreplace I

(

t

)

intheterms±S

(

t

)

I

(

t

)

inthespace-dependentversionofmodel(1) bythefollowingweighted integralontheball

(

x

,

y

)

:

I

(

t

,

x

,

y

) :=

0

2π

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

(

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 maybeusefulformodellingthespreadofdiseaseswhenthereisa

(3)

constant 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

) =

δ 0

2π

0

g1

(

r

)

g2

(ϑ)

I

(

t

,

x

+

rcos

(ϑ),

y

+

rsin

(ϑ))

rd

ϑ

dr (5)

forallt0,

(

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)

forallt0,

(

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,j0. Fortheinfectedindividuals beingclosertothe boundary ofthedomain

asthe radius

δ

, theapproximationoftheintegralin(5) needsvaluesofIlyingoutside

:forthese,wearegoingtousezerovalues.After theseconsiderationswegetthefollowingsystem,beingstillcontinuousint0 and

(

x

,

y

)

:

(4)

⎧ ⎪

⎪ ⎩

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

/(

K1

) >

0 and hy:=B

/(

L1

) >

0 in directionsxand y,respectively.Thenthegriditselfisthefollowingset:

G

:=

(

xk

,

y

) |

xk

= (

k

1

)

hx

,

y

= (

1

)

hy

,

k

=

1

, . . . ,

K

, = , . . . ,

L

.

Forallt0,

(

xk

,

y

)

G,and X∈ {S

,

I

,

R},weconsidertheapproximatenumbers

Xk,

(

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].WenoteherethatifI0 holds,thenTk,

(

t

)

isnon-negative too,forallt0 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)

forallt0 and

(

x

,

y

)

.Thelinkbetweenthesub-problemsistheinitialcondition,aswillbeshowninthenextsections.

Forthelateruseweremarkthatsub-problem(Sub.1) canbesolvedexactly:

⎧ ⎪

⎪ ⎨

⎪ ⎪

S[1]

(

t

+ τ ,

x

,

y

) =

ecτS[1]

(

t

,

x

,

y

),

I[1]

(

t

+ τ ,

x

,

y

) =

ebτI[1]

(

t

,

x

,

y

),

R[1]

(

t

+ τ ,

x

,

y

) =

R[1]

(

t

,

x

,

y

) + (

1

ecτ

)

S[1]

(

t

,

x

,

y

) + (

1

ebτ

)

I[1]

(

t

,

x

,

y

)

(13)

(5)

forallt0 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

) +

12

R[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

) =

0

S

(

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-increasingintime

S

(

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)

(6)

4. Similarly,thedensityofrecoveredindividualscannotdecreaseintime:

R

(

t

,

x

,

y

)

R

(

t

+ τ ,

x

,

y

)

for all t

, τ

0 and

(

x

,

y

).

(20) Whichmeansforthenumericalvaluesthat

Rnk,

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,let

M

:=

M

·

S0

+

I0

+

R0

inwhich

.

∞meansthemaximummatrixnormtakenelement-wise.Wenotethatcondition(8) impliesM

>

0.

(iii) LetW1: [−1

/

e

,

0

)

(−∞,

1]andW0:

(−

1

/

e

,

+∞)→

(−

1

,

+∞)denotethetwobranchesoftheLambert-W func- tion,thatis,theinverseofthemapxxex.

(iv) Forarbitraryp

,

q

>

0,wedefinetheset

T

p,q

:=

0

,

1pW0

pq

1pW1

pq

, +∞

⊂ R.

Furthermore,wedefine

T

0,q

:= [

0

,

1q

) ⊂ R.

Thelatternotationmakessensebecauseofthefollowingconsideration.

Lemma5.2.WithNotation5.1,thelimit1pW0pq−−−→p0 1q holdsforarbitraryq

>

0.

Proof. ItsufficestoshowthatW0

(

x

)/

x−−→x0 1 forx= −p

/

q

<

0.TheL’Hospitalrule,thederivativeoftheinversefunction, andtheidentityW0

(

0

)

=0 implythat

xlim→0 W0

(

x

)

x

=

lim

x0W0

(

x

) =

lim

x0

1

eW0(x)

+

W0

(

x

)

eW0(x)

=

1

.

Remark5.3.Sincewewilluseitseveraltimesthroughoutthepaper,weanalysethesolutionx

<

0 toequation

xex

= μ

(21)

forsomeparameter

μ <

0.

(i) For

μ <

1

/

e,thereisnosolutiontoequation(21).

(ii) For

μ

= −1

/

e,thereisonesolution:x1= −1.

(iii) For

μ >

1

/

e,therearetwosolutions:x1=W1

( μ )

andx0=W0

( μ )

. Wealsoknowthat x1x1= −1

<

x0.Hence,fortheinequality

xex

μ

(22)

wehavethefollowingcases.

(7)

Fig. 1.Graph of functionxxex. The horizontal lines indicate theμ-values0.25 and1/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

<

x1=W1

( μ )

orx

>

x0=W1

( μ )

. Thegraphoffunctionxxex 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

1

M

,

1 b

.

Thecasec

>

0 wasstudiedin[13],andresultedinasimilarbound,namely

τ

min

1 M

+

c

,

1

b

.

(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

(8)

(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

=

ecτSk[1,],n

,

Ik[1,],n+1

=

ebτIk[1,],n

,

Rk[1,],n+1

=

R[k1,],n

+ (

1

ecτ

)

S[k1,],n

+ (

1

ebτ

)

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

=

ecτSnk,

,

I[k1,],n+1

=

ebτIkn,

,

R[k1,],n+1

=

Rnk,

+ (

1

ecτ

)

Snk,

+ (

1

ebτ

)

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

(

ebτInk,

) =

ebτM

(

Ink,

) =

ebτTkn,

.

(29) Bycombiningthesub-problems(27)–(28),andtherelation(29),wearriveatthenumericalscheme:

⎧ ⎪

⎪ ⎨

⎪ ⎪

Snk,+1

=

ecτSnk,

(

1

τ

ebτTkn,

),

Ink,+1

=

ebτ

(

Ikn,

+ τ

ecτSnk,Tkn,

),

Rnk,+1

=

Rnk,

+ (

1

ecτ

)

Snk,

+ (

1

ebτ

)

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,andebτ

>

0 inthefirstandthethirdequationsofsystem(30),wegetthatproperties (P3) and(P4) alsohold.

Thisconcludestheproof.

Ábra

Fig. 1. Graph of function x → xe x . The horizontal lines indicate the μ -values − 0
Fig. 2. Graph of function V , b ( τ ) = τ ( 1 − ( 1 − e − b τ )) for  = 0 . 95 and b = 0
Fig. 4. The numerical solutions S n k , , I n k , , R n k , shown in columns, respectively, at time levels t = 0, t = 5, t = 10, t = 15, t = 30, for the sequential splitting (30).
Fig. 5. Order plot for the relative global error of the four splitting schemes identified with colours

Hivatkozások

KAPCSOLÓDÓ DOKUMENTUMOK

For this LPV model identification step, two LPV model structures (a black-box and a gray box structure, respectively) are handled, as explained in Sub-Section 5.1. As far as

parapsilosis is able to induce inflammasome activation with sub- sequent IL-1β and IL-18 release in both human THP-1 macrophages and MDMs, but the induction of cytokine

The protocol works in two steps, where at step 1 sub-channels are approximately evenly distributed to the terminals and at step 2 terminals within in a sub- channel will contend

Pancreatic elastase-1 in stools, a marker of exocrine pancreas function, correlates with both residual b -cell secretion and metabolic control in type 1 diabetic sub-

The extreme loads special presentation of two extreme vehicle described in the following two sub-heading, which will be presented as an example of the supposedly

Positive feedback loops and the merged controller sub-network in lesional psoriatic skin Individual 730. positive feedback loops with 2, 3 or 4 nodes

A kollaboratív művészeti tevékenység eredményeképp létrejött sub tentorium / sátor alatt című helyspecifikus műegyüttest most, az Alkalmazott Művészeti Intézet

The Independent Strongly Nondecisive Clause Rule states that if a clause set contains an independent strongly nondecisive clause, then it is satisfiable and a sub- model generated