Toepassingen op discrete dynamische systemen Een discreet dynamisch systeem is een proces van de vorm xk+1 = Axk ,
k = 0, 1, 2, . . .
met A een vierkante matrix. Een Markov-proces is een speciaal geval van een discreet dynamisch systeem. Als de (n × n)-matrix A diagonaliseerbaar is, dan bestaat er een basis {v 1 , . . . , v n } van Rn bestaande uit eigenvectoren van A, dat wil zeggen : Av i = λi v i ,
i = 1, 2, . . . , n.
Dat betekent dat elke startvector x0 ∈ Rn geschreven kan worden als lineaire combinatie van de basisvectoren v 1 , . . . , v n , zeg x0 = c1 v 1 + . . . + cn v n . Hieruit volgt x1 = Ax0 = c1 Av 1 + . . . + cn Av n = c1 λ1 v 1 + . . . + cn λn v n , x2 = Ax1 = c1 λ1 Av 1 + . . . + cn λn Av n = c1 λ21 v 1 + . . . + cn λ2n v n , enzovoorts. Dus : xk = c1 λk1 v 1 + . . . + cn λkn v n ,
k = 0, 1, 2, . . . .
Voorbeeld 1. Beschouw het discreet dynamisch systeem xk+1 = Axk , k = 0, 1, 2, . . . met 1/2 1/4 A = . Merk op dat A een stochastische matrix (of Markov-matrix) is. Nu 1/2 3/4 volgt : 1/2 − λ 1/4 = λ2 − 5 λ + 1 = (λ − 1)(λ − 1 ). |A − λI| = 1/2 3/4 − λ 4 4 4 De eigenwaarden van A zijn dus λ1 = 1 en λ2 = 1/4. Verder volgt : −1/2 1/4 −2 1 1 λ1 = 1 : ∼ =⇒ v 1 = 2 1/2 −1/4 0 0 en
λ2 = 1/4 :
1/4 1/4 1/2 1/2
∼
1 1 0 0
=⇒
v2 =
De oplossing kan dus geschreven worden als k 1 1 1 k k xk = c1 λ1 v 1 + c2 λ2 v 2 = c1 + c2 , 2 −1 4
1 −1
.
k = 0, 1, 2, . . . ,
waarbij c1 en c2 bepaald kunnen worden uit c1 v 1 + c2 v 2 = x0 . We zien dat lim xk bestaat k→∞
voor elke startvector x0 :
lim xk = c1
k→∞
1
1 2
,
1/2 waarbij c1 afhankelijk is van de keuze van de startvector x0 . Als we bijvoorbeeld x0 = 1/2 kiezen, dan is c1 = 1/3. Dit kunnen we snel vinden omdat het een Markov-proces betreft. In het algemeen vinden we c1 en c2 door : 1 1 1/2 c1 v 1 + c2 v 1 = x0 ⇐⇒ c1 + c2 = 2 −1 1/2 en dus
1 1 1 1 1/2 ∼ 2 −1 1/2 3 0
1/2 3 0 ∼ 1 0 3
1 1/2
=⇒
c1 = 1/3
c2 = 1/6.
Voorbeeld 2. Beschouw het voorbeeld uit de tweede week : van het autoverhuurbedrijf 0, 85 0, 05 0, 10 xk+1 = Axk met k = 0, 1, 2, . . . en A = 0, 05 0, 85 0, 05 . Dan volgt (tel de tweede en 0, 10 0, 10 0, 85 de derde rij op bij de eerste) : 0, 85 − λ 0, 05 0, 10 1 − λ 1−λ 1 − λ 0, 85 − λ 0, 05 = 0, 05 0, 85 − λ 0, 05 |A − λI| = 0, 05 0, 10 0, 10 0, 85 − λ 0, 10 0, 10 0, 85 − λ 1 0 0 = (1 − λ)(0, 80 − λ)(0, 75 − λ). 0 = (1 − λ) 0, 05 0, 80 − λ 0, 10 0 0, 75 − λ De eigenwaarden van A zijn dus λ1 = 1, λ = 0, 80 en λ3 = 0, 75. Verder volgt : −0, 15 0, 05 0, 10 −3 1 2 7 0, 05 ∼ 1 −3 1 =⇒ E1 = Span{ 5 }, λ1 = 1 : 0, 05 −0, 15 0, 10 0, 10 −0, 15 2 2 −3 8 zoals we eerder hebben gezien, 0, 05 0, 05 0, 10 1 1 2 1 λ2 = 0, 80 : 0, 05 0, 05 0, 05 ∼ 1 1 1 =⇒ E0,80 = Span{ −1 }, 0, 10 0, 10 0, 05 2 2 1 0 en λ3 = 0, 75 :
0, 10 0, 05 0, 10 2 1 2 1 0, 05 0, 10 0, 05 ∼ 1 2 1 =⇒ E0,75 = Span{ 0 }. 0, 10 0, 10 0, 10 2 2 2 −1
De algemene oplossing van het dynamische systeem xk+1 = Axk met k = 0, 1, 2, . . . is dus 7 1 1 xk = c1 5 + c2 (0, 80)k −1 + c3 (0, 75)k 0 , k = 0, 1, 2, . . . , 8 0 −1 2
waarbij c1 , c2 en c3 gevonden kunnen worden uit : 7 1 1 c1 5 + c2 −1 + c3 0 = x0 . 8 0 −1
0, 25 Als bijvoorbeeld x0 = 0, 25 , dan volgt (ga na !) : 0, 50
7 1 1 0, 25 1 0 0 5 −1 0 0, 25 0 1 0 ∼ 8 0 −1 0, 50 0 0 1
0, 05 0, 00 −0, 10
In dat geval is de oplossing dus gelijk aan 7 1 xk = 0, 05 5 − 0, 10 (0, 75)k 0 , 8 −1
=⇒
c1 = 0, 05 c2 = 0 c3 = −0, 10.
k = 0, 1, 2, . . . .
Duidelijk is dan dat lim xk bestaat en dat k→∞
7 0, 35 lim xk = 0, 05 5 = 0, 25 . k→∞ 8 0, 40 Merk op dat we nu in staat zijn de oplossing expliciet te bepalen. Dat betekent dat we de situatie voor elke waarde van k ∈ {0, 1, 2, . . .} kunnen berekenen. We zien nu ook dat de limiet lim xk daadwerkelijk bestaat. Voorheen konden we de limiet alleen bepalen als we k→∞
aannemen dat deze bestaat. Met de methode die we in voorbeeld 2 gebruikten om het karakteristieke polynoom van de matrix A te bepalen kunnen we nu ook aantonen dat elke stochastische matrix (of Markovmatrix) een eigenwaarde 1 heeft. Beschouw namelijk |A − λI| en tel de tweede tot en met de laatste rij op bij de eerste rij. Alle elementen in de eerste rij zijn dan gelijk aan 1 − λ, omdat in de Markov-matrix de som van alle elementen in ´e´en kolom gelijk is aan 1. Dit betekent dat |A − λI| deelbaar is door 1 − λ en dat λ = 1 dus een eigenwaarde van A is. Dat betekent dat elke Markov-matrix dus een evenwichtstoestandsvector q (met Aq = q) heeft (zie : Lay, pag. 286). De oplossing van een discreet dynamisch systeem xk+1 = Axk met k = 0, 1, 2, . . . kan geschreven worden in de vorm xk = Ak x0 , k = 0, 1, 2, . . . , waarbij x0 de startvector is. Als A diagonaliseerbaar is, dan is A = P DP −1 voor zekere inverteerbare matrix P en diagonaalmatrix D. Dan is Ak = P Dk P −1 voor k = 0, 1, 2, . . .. 3
Stel dat A diagonaliseerbaar is en dat A = P DP −1 voor zekere inverteerbare matrix P en diagonaalmatrix D. Stel vervolgens xk = P y k voor k = 0, 1, 2, . . ., dan volgt : xk+1 = Axk
⇐⇒
P y k+1 = AP y k = P DP −1 P y k = P Dy k ,
k = 0, 1, 2, . . . .
Aangezien P inverteerbaar is kunnen we links en rechts vermenigvuldigen met P −1 . Dan volgt dus : xk+1 = Axk ⇐⇒ y k+1 = Dy k , k = 0, 1, 2, . . . . c1 λk1 Als D = diag(λ1 , . . . , λn ), dan volgt hieruit dat y k = ... voor k = 0, 1, 2, . . .. Met cn λkn P = v 1 . . . v n volgt dan :
xk = P y k = v 1 . . . v n
c1 λk1 . k k .. = c1 λ1 v 1 + . . . + cn λn v n , cn λkn
k = 0, 1, 2, . . . .
In het geval dat A een (2 × 2)-matrix is kunnen we de vectoren {xk }∞ k=0 die voldoen aan 2 xk+1 = Axk voor k = 0, 1, 2, . . . ’plotten’ in het platte vlak R . De grafiek van {xk }∞ k=0 wordt een baan van het discrete dynamische systeem xk+1 = Axk genoemd. Merk op dat zo’n baan bestaat uit allemaal ’losse’ punten in R2 en dat deze alleen afhankelijk is van de keuze van de startvector x0 . Als voor de eigenwaarden λ1 en λ2 van A geldt dat |λ1 | < 1 en |λ2 | < 1, dan gaan alle oplossingen xk = c1 λk1 v 1 + c2 λk2 v 2 voor k → ∞ naar de oorsprong O. Men zegt dan dat de oorsprong O een aantrekker (attractor) is van het dynamische systeem xk+1 = Axk . Als voor de eigenwaarden λ1 en λ2 van A geldt dat |λ1 | > 1 en |λ2 | > 1, dan gaan alle oplossingen xk = c1 λk1 v 1 + c2 λk2 v 2 voor k → ∞ naar oneindig (weg van de oorsprong O). Men zegt dan dat de oorsprong O een afstoter is van het dynamische systeem xk+1 = Axk . Als |λ1 | < 1 en |λ2 | > 1, dan noemt men de oorsprong wel een zadelpunt van het dynamische systeem xk+1 = Axk . Sommige banen bewegen dan in de richting van de oorsprong, terwijl andere naar oneindig gaan (afhankelijk van de startvector x0 ). In het geval van niet-re¨ele (en dus complex geconjugeerde) eigenwaarden hebben de banen de vorm van een spiraal (zie : Lay, pag. 344). Als λ1,2 p= α ± iβ met α, β ∈ R en β 6= 0, dan zijn de banen naar depoorsprong toe gericht als |λ| = α2 + β 2 < 1 en juist van de oorsprong af gericht als |λ| = α2 + β 2 > 1.
4
Numerieke methoden voor het bepalen van eigenwaarden We bestuderen twee iteratieve methoden voor het benaderen van de eigenwaarden van een matrix. De eerste methode is de machtmethode of ook wel powermethode en de tweede is de inverse machtmethode. De tweede methode is afgeleid van de eerste. De machtmethode of powermethode kan worden gebruikt om van een (n×n)-matrix A de strikt dominante eigenwaarde λ1 te vinden. Een eigenwaarde λ1 heet (strikt) dominant als deze in absolute waarde (strikt) groter is dan de absolute waarden van de andere eigenwaarden. Neem aan dat A diagonaliseerbaar en dat er dus een basis {v 1 , . . . , v n } van Rn bestaat bestaande uit eigenvectoren van A behorende bij de eigenwaarden λ1 , . . . , λn respectievelijk, waarbij |λ1 | > |λ2 | ≥ |λ3 | ≥ . . . ≥ |λn |. Dan is dus λ1 de strikt dominante eigenwaarde van A. Omdat {v 1 , . . . , v n } een basis van Rn is kan elke vector x ∈ Rn geschreven worden in de vorm x = c1 v 1 + . . . + cn v n . Omdat Av i = λi v i voor i = 1, 2, . . . , n volgt dan : Ak x = c1 λk1 v 1 + . . . + cn λkn v n ,
k = 0, 1, 2, . . . .
Delen door λk1 geeft dan 1 k A x = c1 v 1 + c2 λk1
λ2 λ1
k
v 2 + . . . + cn
λn λ1
k vn,
k = 0, 1, 2, . . . .
Nu is |λi /λ1 | < 1 voor alle i = 2, 3, . . . , n en dus 1 k A x → c1 v 1 λk1
voor k → ∞.
Dit geeft aanleiding tot het volgende algoritme (details laten we achterwege) : • Kies een startvector x0 waarvan de (in absolute waarde) grootste co¨ordinaat gelijk aan 1 is. • Voor k = 0, 1, 2, . . . – Bereken Axk . – Kies voor µk de in absolute waarde grootste co¨ordinaat van de vector Axk . – Bereken xk+1 = (1/µk )Axk . • Voor bijna elke keuze van de startvector x0 zal de rij {µk }∞ k=0 naderen naar de strikt ∞ dominante eigenwaarde van A en de rij {xk }k=0 naar een bijbehorende eigenvector.
5 2 Voorbeeld 3. De matrix A = heeft de eigenwaarden λ1 = 6 en λ2 = 3 (ga na !). 1 4 Er geldt dus : |λ1 | > |λ2 | (dus : λ1 is een strikt dominante eigenwaarde van A). 5
0 1
2 4
Kies nu (bijvoorbeeld) x0 = Ax0 =
5 2 1 4
0 1
. Dan volgt :
=
=4
1/2 9/2 = = Ax1 = 1 9/2 5 2 1 7 Ax2 = = =7 1 4 1 5 5 2 1 4
5 2 1 4
Ax3 =
1 5/7
=
2/4 1
1 1
9 2
1 5/7
45/7 27/7
Ax4 =
5 2 1 4
1 3/5
=
µ0 = 4
45 = 7
31/5 17/5
=⇒
en
=⇒
µ2 = 7
1 3/5
31 = 5
1 17/31
en x3 =
en
1/2 1
1 1
1 5/7
1 3/5
x2 =
x4 =
,
, ,
,
31 µ4 = = 6, 2 5
=⇒
en x1 =
9 µ1 = 2
45 µ3 = ≈ 6, 42857 7
=⇒
=⇒
en
x5 =
1 17/31
.
De convergentie is niet erg snel, maar er geldt : lim µk = 6
k→∞
en
lim xk =
k→∞
1 1/2
.
De convergentie is sneller als de startvector beter in de richting van de eigenvector 1 gekozen wordt. Kiezen we bijvoorbeeld x0 = , dan volgt : 0 Ax0 =
5 2 1 4
Ax1 =
1 0
5 2 1 4
5 1
1 1/5
=
=5
=
1 1/5
27/5 9/5
Ax2 =
5 2 1 4
1 1/3
=
17/3 7/3 =⇒
=⇒
27 = 5
µ0 = 5
1 1/3
17 = 3
1 7/17
,
en
x2 =
1 1/3
,
17 µ2 = ≈ 5, 66667 3 6
en x1 =
1 1/5
27 µ1 = = 5, 4 5
=⇒
2 1
en x3 =
1 7/17
,
Ax3 =
5 2 1 4
1 7/17
=
99/17 45/17 =⇒
Ax4 =
5 2 1 4
1 5/11
=
65/11 31/11 =⇒
99 1 = 17 5/11 99 µ3 = ≈ 5, 82353 17
en x4 =
1 5/11
,
65 1 = 11 31/65 65 1 µ4 = ≈ 5, 90909 en x5 = . 31/65 11
Door dezelfde methode los te laten op een andere matrix kunnen we ook andere eigenwaarden van een matrix A benaderen. Stel dat α ∈ R een redelijke benadering is van ´e´en van de eigenwaarden λ1 , . . . , λn van A. Stel nu B = (A − αI)−1 , dan geldt dat 1 , λ1 − α
1 , λ2 − α
...,
1 λn − α
de eigenwaarden zijn van B. De in absolute waarde grootste eigenwaarde is die met de in absolute waarde kleinste noemer λp − α. Passen we de machtmethode of powermethode toe op de matrix B dan wordt die strikt dominante eigenwaarde 1/(λp − α) benaderd. Aangezien we de waarde van α kennen kunnen we hiermee de waarde van de eigenwaarde λp benaderen. Dit wordt de inverse machtmethode of inverse powermethode genoemd. Hierbij treedt echter een klein probleem op als we Bxk = (A − αI)−1 xk willen berekenen. Daarvoor hebben we immers de inverse van A − αI nodig. In plaats van het berekenen van y k = (A − αI)−1 xk lossen we echter de vergelijking (A − αI)y k = xk op (bijvoorbeeld via een LU -decompositie). Dit geeft aanleiding tot het volgende algoritme : • Kies een geschikte waarde voor α (voldoende dicht bij de gezochte eigenwaarde λ). • Kies een startvector x0 waarvan de (in absolute waarde) grootste co¨ordinaat gelijk aan 1 is. • Voor k = 0, 1, 2, . . . – Los (A − αI)y k = xk op. – Kies voor µk de in absolute waarde grootste co¨ordinaat van de vector y k . – Bereken νk = α + (1/µk ). – Bereken xk+1 = (1/µk )y k . • Voor bijna elke keuze van de startvector x0 zal de rij {νk }∞ k=0 naderen naar de eigenwaarde λ van A en de rij {xk }∞ naar een bijbehorende eigenvector. k=0 Voorbeeld 4. Als we α =2 kiezen als een (grove) benadering van λ2 = 3 de eigenwaarde 5 2 3 2 van de matrix A = , dan geldt dus : A − αI = A − 2I = . 1 4 1 2 7
Kiezen we nu als startvector x0 = (A − 2I)y 0 = x0 : =⇒
µ0 =
(A − 2I)y 1 = x1 : =⇒
(A − 2I)y 2 = x2 : =⇒
µ2 =
3 2 1 2
3 2 1 2
3 2 1 2 =⇒
3 2 1 2
y2 =
y4 = µ4 =
ν4 = 2 +
683 ≈ 0, 99854, 684
A
=
=
172 ≈ 3, 00585 171
=⇒ y 5 = ν5 = 2 +
684 ≈ 3, 00146 683
5 2 1 4
1 −1
11 = 12
1 −21/22 1 x3 = , −21/22
43 44
en
171 = 172
683/684 −455/456
Het lijkt er op dat νk → 3 en dat xk → 1 −1
43/44 −85/88
171/172 −341/344
11/12 −7/8
en
44 ν3 = 2 + ≈ 3, 02326 43
=⇒ y 4 =
µ5 =
12 ≈ 3, 09091 11
y5 =
=⇒ y 2 =
=⇒ y 3 =
1 −341/342
1 −5/6
ν2 = 2 +
171 ≈ 0, 99419, 172
1 −1/2
1 −21/22
1 −85/86
1 0
3 1 3/4 y1 = =⇒ y 1 = = −5/8 4 −5/6 4 1 ν1 = 2 + = 3, 33333 en x2 = , −5/6 3
43 µ3 = ≈ 0, 97727, 44
=⇒
3 2 1 2
11 ≈ 0, 91667, 12
y3 =
=⇒
1 1 1/2 y0 = =⇒ y 0 = = −1/4 2 −1/2 1 1 ν0 = 2 + = 4 en x1 = , −1/2 1/2
3 = 0, 75, 4
µ1 =
, dan volgt :
3 2 1 2
1 = 0, 5, 2
1 0
1 −85/86
1 −1365/1366 1 . x6 = −1365/1366
683 = 684
voor k → ∞. Er geldt inderdaad
8
=
3 −3
=3
,
1 −341/342 1 , x5 = −341/342
1 −1
en
en
1 −85/86 x4 =
1 −1
.
