københavns universitet
københavns universitet
1
System af differensligninger
2
System af lineære differensligninger
3
Autonome differensligninger. Ligevægt og stabilitet
4
Epidemimodel
5
Vært-parasit model
Systemer af differensligninger
Henrik L. Pedersen
Institut for Grundvidenskab og Miljø
Institut for Matematiske Fag
www.matdat.life.ku.dk/mat-model/
16.5.2012
Dias 2/26
16.5.2012
Dias 1/26
københavns universitet
københavns universitet
To (samhørende) differensligninger
xt+1
= f (t, xt , yt )
yt+1
= g (t, xt , yt )
System af n differensligninger
hvor f og g er givne funktioner
= f1 (t, xt,1 , . . . , xt,n )
xt+1,2
= f2 (t, xt,1 , . . . , xt,n )
..
.
xt+1,n
= fn (t, xt,1 , . . . , xt,n )
hvor f1 , f2 , . . . , fn er givne funktioner
Eksempel
xt+1
=
3xt − xt yt + 2 = f (t, xt , yt )
yt+1
=
2xt + 3yt − t = g (t, xt , yt )
t
xt
yt
16.5.2012
Dias 4/26
xt+1,1
0
1
2
1
3
8
2
-13
29
...
...
...
t
??
??
Vektornotation
xt+1 = f(t, xt )
hvor
16.5.2012
Dias 5/26


xt,1


xt =  ... 
xt,n


f1


f =  ... 
fn
københavns universitet
københavns universitet
System af n differensligninger: løsning
• En løsning er en følge af vektorer (xt ) = (x0 , x1 , . . .) som opfylder
xt+1 = f(t, xt ) for all t ≥ 0.
• En løsning er fastlagt ved et sæt af begyndelsesværdier, også kaldet
en startvektor
System af lineære differensligninger
Definition
Et system af lineære differensligninger er på formen
xt+1 = Axt + b = f(t, xt ) = f(xt )


x0,1


x0 =  ...  .
x0,n
Løsningen (xt ) = (x0 , x1 , . . .) kan nemlig findes ved successiv
udregning:
hvor f(x) = Ax + b er en affin afbildning (A er en given n × n matrix og
b ∈ Rn en given vektor)
(Et system af lineære differensligninger er det samme som en affin model)
Eksempel: nationaløkonomi
x1 = f(0, x0 ), x2 = f(1, x1 ), x3 = f(2, x2 ) osv.
16.5.2012
Dias 6/26
Ct+1
It+1
=
a
ac
Ct
b
+
It
bc
16.5.2012
Dias 8/26
københavns universitet
københavns universitet
Autonome differensligninger: ligevægt
Sætning: formel for løsningerne
Hvis A − E har en invers matrix (er regulær), så er den fuldstændige
løsning til ligningen
xt+1 = Axt + b
Definition: system af autonome differensligninger
givet ved
xt = At c − (A − E)−1 b
xt+1 = f(xt )
(c ∈ Rn )
Ingen direkte afhængighed af t på højre side
Sætning: formel vha. egenværdier og -vektorer
Hvis A er diagonaliserbar med egenværdier λ1 , . . . , λn (alle forskellige fra
1) og tilhørende egenvektorer q1 , . . . , qn , så kan den fuldstændige løsning
til
xt+1 = Axt + b
xt = c1 λt1 q1 + . . . + cn λtn qn − (A − E)−1 b
Definition: ligevægt
En vektor x∗ er en ligevægt for et autonomt system af differensligninger
xt+1 = f(xt ) hvis
x∗ = f(x∗ )
dvs. hvis x∗ “ikke ændrer sig fra gang til gang”
også skrives på formen
16.5.2012
Dias 9/26
a
(a − 1)c
(c ∈ Rn )
16.5.2012
Dias 11/26
københavns universitet
Autonome differensligninger: stabilitet
københavns universitet
Stabilitet af ligevægt: system af lineære
differensligninger
Definition: stabil ligevægt
En ligevægt x∗ for et autonomt system xt+1 = f(xt ) er stabil hvis der
gælder følgende:
For alle startvektorer x0 i nærheden af x∗ skal de tilsvarende
vektorer xt konvergere mod x∗ for t → ∞
+ At en ligevægt ikke er stabil eller ustabil betyder altså, at der findes
startvektorer x0 vilkårligt tæt på ligevægten, hvor de tilsvarende
vektorer xt ikke konvergerer mod ligevægten
Selv små ændringer i startvektoren kan altså betyde meget i det
lange løb
16.5.2012
Dias 12/26
københavns universitet
Stabilitet af ligevægt: generelt tilfælde
• Hvordan afgør vi om en ligevægt er stabil eller ustabil?
• For en enkelt differensligning xt+1 = f (xt ) er en ligevægt x ∗ stabil
når |f 0 (x ∗ )| < 1 og ustabil når |f 0 (x ∗ )| > 1 (sætning fra sidst)
• Hvordan generaliseres dette til systemer af differensligninger?
Definition: funktionalmatrix
Lad f være en vektorfunktion af n variable. Matricen
 0

0
f1,x1 (x) · · · f1,x
(x)
n


..
..
..
Df(x) = 

.
.
.
0
0
fn,x
(x)
·
·
·
f
(x)
n,xn
1
kaldes funktionalmatricen for vektorfunktionen
Sætning: stabilitet (repetition fra Modul 1)
Hvis alle egenværdierne for A har modulus mindre end 1, så har modellen
xt+1 = Axt + b netop en ligevægt x∗ og den er stabil: Der gælder
xt → −(A − E)−1 b = x∗
for t → ∞ uanset startvektoren x0
16.5.2012
Dias 13/26
københavns universitet
Stabilitet af ligevægt: generelt tilfælde
Sætning
Lad x∗ være en ligevægt for et system af autonome differensligninger. Da
gælder
• Hvis alle egenværdier for Df(x∗ ) har modulus < 1, så er ligevægten
x∗ stabil.
• Hvis blot én af egenværdierne for Df(x∗ ) har modulus > 1, så er
ligevægten x∗ ustabil.
+ Man skal indsætte ligevægten x∗ på x’s plads i funktionalmatricen
 0

0
f1,x1 (x∗ ) · · · f1,x
(x∗ )
n


..
..
..
Df(x∗ ) = 

.
.
.
0
∗
0
∗
fn,x
(x
)
·
·
·
f
(x
)
n,xn
1
og dernæst bestemme egenværdierne for denne matrix.
16.5.2012
Dias 14/26
16.5.2012
Dias 15/26
københavns universitet
københavns universitet
Eksempel (stabilitet)
Eksempel (ligevægte)
• Partielle afledede
• System af autonome differensligninger
xt+1
=
yt+1
=
1
2 xt
1
4 xt
fx0 (x, y ) =
− 13 xt yt
gx0 (x, y ) =
+ 12 yt
fy0 (x, y ) = − 13 x
− 13 y
gy0 (x, y ) =
1
2
• Stabilitet af ligevægten (x ∗ , y ∗ ) = (0, 0):
• Udtrykkene på højresiderne
f (x, y )
=
g (x, y )
=
1
2x
1
4x
1 ∗
1 ∗ ∗
2x − 3x y
1 ∗
1 ∗
4x + 2y
1
2
1
4
− 13 xy
Df(0, 0) =
+ 12 y
• Ligevægte (x ∗ , y ∗ ) findes ved at løse ligningerne
= x∗
= y
0
1
2
har to egenværdier med modulus < 1
• Ustabilitet af ligevægten (x ∗ , y ∗ ) = (−3, − 32 ):
∗
Df(−3, − 32 )
• Resultat (x ∗ , y ∗ ) = (0, 0) eller (x ∗ , y ∗ ) = (−3, − 32 )
=
1
1
1
4
1
2
har en egenværdi med modulus > 1
16.5.2012
Dias 17/26
16.5.2012
Dias 16/26
københavns universitet
københavns universitet
Eksempel (succesiv udregning med R)
xt+1
=
yt+1
=
1
2 xt
1
4 xt
− 13 xt yt
Epidemimodel
Vi vil udbygge modellen
+ 12 yt
>
>
>
>
>
>
f <- function(x,y){x/2-x*y/3};
g <- function(x,y){x/4+y/2};
v0<-c(1,3) # begyndelsesvektoren v0=(1,3)
V<-v0; v<-v0;
for (t in (1:9)){v<-c(f(v[1],v[2]),g(v[1],v[2]));V<-cbind(V,v)}
V
V
v
v
v
v
v
v
[1,] 1 -0.50 0.04166667 0.01041667 0.003870081 0.001683082 0.0007862104
[2,] 3 1.75 0.75000000 0.38541667 0.195312500 0.098623770 0.0497326557
v
v
v
[1,] 0.0003800718 0.0001868607 9.264386e-05
[2,] 0.0250628804 0.0126264582 6.359944e-03
St+1
St
= (1 − a)St + bSt 1 −
N
hvor
• St angiver antallet af smittede og smittende individer i en
population af N individer
• a angiver raskhedsraten (bedre på engelsk: recovery rate)
• b angiver infektionsraten
til en såkaldt SIR model.
SIR model – warning about english and danish
English abbriviation: Susceptible, Infected,
Dansk forkortelse:
Modtagelig, Syg,
Udregningerne bekræfter, at (0, 0) er en stabil ligevægt
16.5.2012
Dias 18/26
1
2
1
4
16.5.2012
Dias 20/26
Recovered
Immun
københavns universitet
københavns universitet
SIR model
SIR model
En population på N individer inddeles i tre grupper, modtagelige,
smittede og immune:
Antagelserne leder til
• Mt antallet af individer, som er modtagelige over for smitten
Mt+1
• St antallet af individer, som er smittede og smittende
St+1
• It antallet af individer, som er (midlertidigt) immune over for smitten
Følgende parametre indgår
• a recovery rate
• b infektionsraten
• c angiver raten for tab af immunitet
It+1
b
St Mt + cIt
N
b
= (1 − a)St + St Mt
N
= (1 − c)It + aSt
=
Mt −
Elimination af en af gruppebetegnelserne
Bemærk, at
Mt+1 + St+1 + It+1 = Mt + St + It = N
Antagelser
• Individer, der netop er blevet raske, er midlertidigt immune overfor
smitten
• Antallet af nye tilfælde er proportional med antal modtagelige og
er (konstant) og lig med populationens størrelse.
Der gælder dermed Mt = N − St − It og derfor kan antallet af
modtagelige elimineres fra modellen.
antal smittede
16.5.2012
Dias 21/26
16.5.2012
Dias 22/26
københavns universitet
københavns universitet
SIR model
Vært-parasit model
(Kursorisk stof)
Efter elimination opstilles følgende system af differensligninger
St + It
St+1 = (1 − a)St + bSt 1 −
N
It+1
=
(1 − c)It + aSt
Fortsættelse følger i miniprojektet!
16.5.2012
Dias 23/26
• Parasitter lægger æg i værtsinsekter
• En angreben vært bliver ædt indefra af parasitlarven, som derefter
udvikler sig til en parasit året efter
• I hver vært kan der kun blive udviklet én parasitlarve
• Værter og parasitter dør efter deres første leveår
• Vt :
værtspopulationens størrelse til tiden t målt i år
• Pt :
parasitpopulationens størrelse til tiden t
16.5.2012
Dias 25/26
københavns universitet
• α: det gennemsnitlige antal overlevende æg der lægges af et
•
•
•
•
•
ikke-angrebent værtsinsekt (α > 1)
λt : det gennemsnitlige antal parasitæg, der i år t lægges pr. vært
(afhænger både af Vt og Pt )
Antager at sandsynligheden for at en vært ikke er angrebet af et
parasitæg er givet ved e −λt (antallet af parasitæg i en vært er
Poissonfordelt)
Vt+1 = antal overlevende æg af ikke-angrebne værtsinsekter i år t
Pt+1 = antal angrebne værtsinsekter i år t
Dette giver
Vt+1
= αVt e −λt
Pt+1
= Vt (1 − e −λt )
Fx: α = 2 og λt = 0.1 Pt0.5 :
16.5.2012
Dias 26/26
0.5
Vt+1
= 2Vt e −0.1 Pt
Pt+1
= Vt (1 − e −0.1 Pt )
0.5