Institut fur Numerische Mathematik und Optimierung
Numerische Simulationmathematischer Modelle
Sommersemester 2014
Michael Eiermann
Numerische Simulation mathematischer Modelle 1
Themen dieser Vorlesung
Die Vorlesung besteht aus zweiTeilen:
In Teil I werden Modelle der Populationsdynamik vorgestellt, die auf
gewohnlichen Differentialgleichungen und Differenzengleichungen basieren.
Teil II ist der Untersuchung von Markoff-Ketten gewidmet, mit denen man
sehr viele Prozesse modellieren kann, bei denen der Zufall eine Rolle spielt. Wir
werden uns unter anderem mit Kapazitatsproblemen bei Netzwerken befassen.
In allen Teilen wird die numerische Behandlung der jeweiligen Modelle betont.
Organisation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 2
Literatur
• Edward Beltrami. Von Krebsen und Kriminellen. Mathematische
Modelle in Biologie und Soziologie. Friedr. Vieweg & Sohn
Verlagsgesellschaft mbH, Braunschweig 1993.
• Martin Braun. Differentialgleichungen und ihre Anwendungen.
Springer-Verlag, Berlin 1979.
• Andrew C. Fowler. Mathematical Models in the Applied Sciences.
Cambridge University Press, Cambridge 1997.
• Franz-Josef Fritz, Bertram Huppert und Wolfgang Willems.
Stochastische Matrizen. Springer-Verlag, Berlin 1979.
• Glenn Fullford, Peter Forrester und Arthur Jones. Modeling with
Differential and Difference Equations. Cambridge University Press,
Cambridge 1997.
Organisation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 3
• Guy Latouche und V. Ramaswami. Introduction to Matrix Analytic
Methods in Stochastic Modeling. SIAM, Philadelphia (PA) 1999.
• C. C. Lin und L. A. Segal. Mathematics Applied to Deterministic
Problems in the Natural Sciences. SIAM, Philadelphia (PA) 1988.
• Douglas D. Mooney und Randall J. Swift. A Course in Mathematical
Modeling. The Mathematical Association of America, 1999.
• James D. Murray. Mathematical Biology, 2nd corrected edition.
Springer-Verlag, Berlin 1993.
• James R. Norris. Markov Chains. Cambridge University Press, Cambridge
1996.
• Thomas Sonar. Angewandte Mathematik, Modellbildung und
Informatik. Friedr. Vieweg & Sohn Verlagsgesellschaft mbH,
Braunschweig 2001.
Organisation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 4
Leistungspunkte und Noten
Die Modulprufung besteht aus einer Klausurarbeit im Umfang von 120
Minuten, die in der Prufungsperiode nach dem Sommersemesters 2014 statt
findet [genauer Termin liegt noch nicht fest]. Eine weitere Klausur findet in der
Prufungsperiode nach dem Wintersemester 2014/2015 statt.
Im Modul werden 6 Leistungspunkte erworben. Die Modulnote ergibt sich aus
der Note der schriftlichen Prufung (Wichtung 1).
Der Zeitaufwand betragt 180 h und setzt sich aus 60 h Prasenzzeit und 120 h
Selbststudium zusammen. Letzteres umfasst die Vor- und Nachbereitung der
Lehrveranstaltungen, Vorbereitung und Bearbeiten der Klausur sowie das
Losen von Ubungsaufgaben.
Organisation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 5
1 Mathematische Modellbildung und
numerische Simulation am Beispiel eines
Wasserkreislaufs
”Simulation ist die Nachbildung eines dynamischen Prozesses in einem
Modell, um zu Erkenntnissen zu gelangen, die auf die Wirklichkeit
ubertragbar sind“ (VDI-Richtlinie 3633).
Zum einen ist die rechnerische Simulation dann unumganglich, wenn reale
Experimente mit den Untersuchungsobjekten undurchfuhrbar sind: Denken Sie
etwa an die Entstehung von Galaxien oder an Untersuchungsobjekte, die erst
geplant sind, also noch gar nicht existieren. Aber auch wenn reale Experimente
moglich sind, ist es oft kostengunstiger und ressourcenschonender, stattdessen
numerische Simulationen einzusetzen.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 6
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 7
A. Physikalische Grundlagen.
Evangelista Torricelli (1608–1647):
Abflussgeschwindigkeit v =√
2gh, g = 9.81 (Gravitationsbeschleunigung),
h = Hohe des Wasserspiegels.
Abflussrate als Funktion des im Behalter befindlichen Wasservolumens V (falls
es sich um einen Zylinder mit Grundflache A handelt)
f = a√
2gV/A = c√V mit c := a
√2g/A .
Der Parameter c kann uber a variiert werden, wenn der Abfluss einen Hahn
besitzt. Wir sprechen von einem Steuerungsparameter.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 8
B. Mathematisches Modell.
U(t), V (t),W (t), R(t) : Wassermengen zur Zeit t in den Behaltern.
f1, . . . , f5 : Abflussfunktionen mit den Steuerungsparametern
c1, . . . , c5.
p = p(t) :”Pumpenfunktion“
Anderungsraten der Wasservolumina: Zuflusse weniger Abflusse, d.h.
U ′(t) = p(t)− f1(U(t))− f2(U(t))
V ′(t) = f1(U(t))− f3(V (t))− f4(V (t))
W ′(t) = f2(U(t)) + f4(V (t))− f5(W (t))
R′(t) = f3(V (t)) + f5(W (t))− p(t).
(1.1)
Diese Gleichungen, sog. (gewohnliche) Differentialgleichungen, heißen die
Kontinuitatsgleichungen unseres Systems.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 9
Anfangszustand: Wassermengen in Behaltern zu einem festen Zeitpunkt,
etwa fur t = 0.
Das Verhalten unseres Systems ist fur alle Zeiten t > 0 durch die obigen
Differentialgleichungen eindeutig bestimmt. Die Aufgabe, eine Losung des
Systems (1.1) zu bestimmen, welche gegebene Anfangsbedingungen erfullt,
nennt man ein Anfangswertproblem.
Addiert man alle Gleichungen, so ergibt sich
U ′(t) + V ′(t) +W ′(t) +R′(t) = 0,
ein globales Erhaltungsprinzip, welches besagt, dass sich die
Gesamtwassermenge in unserer Apperatur nicht verandert. (Es handelt sich
hier um ein geschlossenes System.)
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 10
C. Algorithmus.
Anfangswertprobleme lassen sich nur in Ausnahmefallen geschlossen losen
(reine Mathematik: in unserem Fall gibt es genau eine Losung).
Aufgabe der Numerik: Bereitstellung von Naherungslosungen.
Idee: Wir betrachten die Gleichungen nicht mehr fur jeden beliebigen
Zeitpunkt, sondern nur noch zu bestimmten diskreten Zeitpunkten, etwa fur
t0 = 0, t1 = 1, . . . (Diskretisierung).
Un := U(tn), . . . , Rn := R(tn) (Volumina, die sich zum Zeitpunkt tn in den
Behaltern U, . . . , R befinden).
Anderungsrate U ′(tn) wird durch U(tn+1)− U(tn) = Un+1 − Un approximiert
(wir nahern hier eine Tangentensteigung durch eine Sekantensteigung an).
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 11
Un+1 = Un + pn − f1(Un)− f2(Un)
Vn+1 = Vn + f1(Un)− f3(Vn)− f4(Vn)
Wn+1 = Wn + f2(Un) + f4(Vn)− f5(Wn)
Rn+1 = Rn + f3(Vn) + f5(Wn)− pn.
(1.2)
Diese vier Gleichungen heißen die diskreten Kontinuitatsgleichungen unserer
Kreislaufs (System von vier Differenzengleichungen).
Addition liefert globales Erhaltungsprinzip
Un+1 + Vn+1 +Wn+1 +Rn+1 = Un + Vn +Wn +Rn.
Legt man noch einen Anfangszustand fest (etwa U0 = V0 = W0 = 0 sowie R0
= Gesamtwassermenge = 100) und wahlt geeignete Werte fur die Parameter,
etwa c1 =√
12, c2 = c4 =√
2, c3 = 1, c5 = 2 sowie p = 17, so konnen wir
unser Modell”laufen lassen“.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 12
0 5 10 15 20 25 300
10
20
30
40
50
60
70
80
90
100
Zeit
Was
serm
enge
Un
Vn
Wn
Rn
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 13
D. Gleichgewichtswerte.
Simulation: Jede der Großen Un, Vn, Wn und Rn nahert sich mit
zunehmendem n einem Gleichgewichtswert U∞, V∞, W∞, R∞, wenn wir die
Steuerungsparameter nicht andern.
Bestimme Gleichgewichtswerte ohne (zeitaufwendige) Simulation:
p = f1(U∞) + f2(U∞)
f1(U∞) = f3(V∞) + f4(V∞)
f5(W∞) = f2(U∞) + f4(V∞).
Die vierte Gleichung (p = f3(V∞) + f5(W∞)) ist redundant.
Hier – im Gegensatz zur”Praxis“ – Gleichungen einfach (Dreiecksform).
Vorsicht: Die theoretisch ermittelten Gleichgewichtswerte konnen, aber
mussen nicht im Fassungsbereich der Behalter liegen (Nebenbedingungen).
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 14
E. Steuerung.
Wesentliches Ziel von Simulationen: Optimierung des Systemverhaltens bzw.
Entscheidungshilfen fur die Steuerung des Systems.
In unserem Beispiel etwa: Wie muss man die Steuerungsparameter wahlen,
damit sich ein erwunschter (vorgegebener) Gleichgewichtszustand einstellt
(auch hier: Nebenbedingungen, man kann z.B. die Hahne nicht beliebig weit
offnen).
Fixiert man p, so fuhrt dies in unserem Fall zu drei Bedingungen fur die funf
Parameter c1, . . . , c5:
c1 = −c2 + p/√U∞
c4 = −c3 + c1√U∞/
√V∞
c5 = c2√U∞/
√W∞ + c4
√V∞/
√W∞.
In realen Systemen ist ein solches Steuerungsproblem nicht explizit losbar, man
wird es nur naherungsweise und iterativ losen konnen.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 15
F. Kritik.
Realitat → mathematisches Modell → Algorithmus
→ numerische Simulation der Realitat.
Bei jedem dieser drei Ubergange haben wir Fehler begangen,
— Modellierungsfehler. Unser Modell setzt wirbelfreien Wasserfluss voraus;
in der Realitat werden sich aber Wirbel bilden. Die Torricellische
Ausflussformel ist nur gultig, wenn sich die Spiegelhohe langsam andert
und keine Druckdifferenz zwischen Spiegel und Austrittsoffnung besteht,
Voraussetzungen, die in der Realitat nicht immer erfullt sind.
— Diskretisierungsfehler. Wir haben den stetigen Strom des Wassers durch
”Durchschnittswerte“ (bez. Zeit und Raum) ersetzt.
— Rundungsfehler. Computer”rechnen falsch“.
Modellbildung und numerische Simulation Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 16
2 Populationsdynamik
2.1 Zwei einfache Modelle
Sei N(t) die Population einer Spezies zur Zeit t.
Alle Modelle basieren auf dem Erhaltungssatz
dN
dt= Geburten− Todesfalle + Migration.
a. Das einfachste Modell (T. R. Malthus, 1766–1834)
Vernachlassige Migration. Geburts- und Todesfalle seien proportional zu N ,
d.h.dN
dt= bN − dN = rN mit r = b− d. (2.1)
Die allgemeine Losung dieser gewohnlichen Differentialgleichung (gDG) ist
N(t) = N0 exp(rt) mit N0 = N(0).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 17
Exponentielles Wachstum: N(t) = N0 exp(rt)
t
N
N0
0
r > 0r = 0r < 0
Exponentielles Wachstum ist i. Allg. naturlich ziemlich unrealistisch und nur in
sehr speziellen Situationen zu beobachten (etwa bei Bakterienkulturen mit i.
W. unbeschrankten Ressourcen).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 18
b. Ein modifiziertes Modell (P. F. Verhulst, 1804–1849)
dN
dt= rN
(1− N
K
)mit positiven Konstanten r und K. (2.2)
K heißt Tragerkapazitat und hangt etwa von den verfugbaren
Nahrungsressourcen ab.
Die Reproduktionsrate [hier r(1−N/K)] ist jetzt von N abhangig:
r
(1− N
K
) > 0, N < K,
< 0, N > K.
Die allgemeine Losung von (2.2) ist
N(t) = N0K exp(rt)
K +N0(exp(rt)− 1).
Man spricht von logistischem Wachstum. Offenbar gilt: limt→∞N(t) = K.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 19
Logistisches Wachstum
t
N
K
K/2
N0 > K
K > N0 > K/2
N0 < K/2
Beachte, dass sich die Losungen fur 0 < N0 < K/2 und K/2 < N0 < K
qualitativ unterscheiden.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 20
2.2 Elementares aus der Theorie gDGen
Wir betrachten Systeme von expliziten gDGen erster Ordnung:
y′1 = f1(t, y1, y2, . . . , yn)
y′2 = f2(t, y1, y2, . . . , yn)
... =...
y′n = fn(t, y1, y2, . . . , yn)
(2.3)
mit den n unbekannten Funktionen y1 = y1(t), y2 = y2(t), . . . , yn = yn(t).
Sei I ein Intervall.
Jedes System von n Funktionen y1 = y1(t), . . . , yn = yn(t) ∈ C1(I), das
(2.3) fur alle t ∈ I erfullt, heißt Losung von (2.3) uber I.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 21
Beispiel. Das System
y′1 = 1
y′2 = 2y1
besitzt die Losungen
y1(t) = t+ α, y2(t) = t2 + 2αt+ β (α, β ∈ R)
uber (−∞,∞).
Fur eine eindeutige Losung sind Anfangsbedingungen erforderlich.
Z.B. y1(0) = 1, y2(0) = 2.
Dann ist y1(t) = t+ 1, y2(t) = t2 + 2t+ 2 die einzige Losung.
Allgemein: Das Problem, eine Losung von (2.3) zu finden, die die
Anfangsbedingung
y1(t0) = y0,1, . . . , yn(t0) = y0,n (2.4)
erfullt, heißt Anfangswertproblem (AWP) fur die gDG (2.3).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 22
Mit
y :=
y1
...
yn
, f :=
f1
...
fn
, y0 :=
y0,1
...
y0,n
fuhren wir fur (2.3), (2.4) folgende Kurzschreibweise ein:
y ′ = f (t,y), (2.3’)
y(t0) = y0. (2.4’)
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 23
Satz 2.1 (Existenz- und Eindeutigkeitssatz von Picard-Lindelof)
Gegeben ist das AWP (2.3’), (2.4’).
f sei stetig im ‘Quader’ Q := (t,y) : |t− t0| ≤ a, ‖y − y0‖ ≤ b ⊆ Rn+1
und es sei M := max‖f (t,y)‖ : (t,y) ∈ Q.Außerdem erfulle f in Q die Lipschitz-Bedingung
‖f (t,y)− f (t, y)‖ ≤ L‖y − y‖ ∀ (t,y), (t, y) ∈ Q. (2.5)
Dann besitzt das AWP genau eine Losung uber I := [t0 − α, t0 + α], wobei
α = mina, b/M.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 24
Bemerkungen.
1. Das AWP (2.3), (2.4) besitzt in [t0 − a, t0 + a] eine eindeutige Losung,
wenn f die Lipschitz-Bedingung (2.5) in
Q = (t,y) : |t− t0| ≤ a, ‖y‖ <∞ erfullt.
2. Ist f auf Q bez. y stetig differenzierbar und bezeichnet
f ′y = [∂fi/∂yj ]1≤,i,j≤n die zugehorige Jacobi-Matrix, dann folgt aus dem
Mittelwertsatz, dass die Voraussetzungen von Satz 2.1 erfullt sind mit
L = sup(t,y)∈Q
‖f ′y (t,y)‖ <∞.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 25
Satz 2.2 (Stetige Abhangigkeit von den Daten)
Die Voraussetzungen von Satz 2.1 seien erfullt.
Sind y , z Losungen der AWPe y ′ = f (t,y), y(t0) = y0, bzw. z ′ = g(t, z ),
z (t0) = z0 uber I (g sei stetig in Q) und gilt
‖y0 − z0‖ ≤ γ sowie ‖f (t,y)− g(t,y)‖ ≤ δ (∀ (t,y) ∈ Q),
dann folgt fur t ∈ I
‖y(t)− z (t)‖ ≤ γ eL(t−t0) +δ
L
(eL(t−t0) − 1
).
Die folgenden Seiten behandeln ein wichtiges Beispiel fur gDGen...
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 26
2.2.1 Lineare Systeme von gDGen mit konstanten Koeffizienten
y ′ = Ay + b(t) mit A = [ai,j ] ∈ Rn×n (unabhangig von t). (2.6)
Das System heißt homogen, falls b = 0 , sonst inhomogen.
Es gelten:
• Sind b1(t), b2(t), . . . , bn(t) stetig in I = [t0 − a, t0 + a], dann besitzt (2.6)
mit der Anfangsbedingung y(t0) = y0 genau eine Losung in I.
• Sind y und z zwei Losungen von (2.6), dann lost w = y − z das
homogene System y ′ = Ay .
• Die Losungen von y ′ = Ay bilden einen Vektorraum der Dimension n.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 27
Noch spezieller (n = 2):
y ′ = Ay mit A =
[a b
c d
]∈ R2×2 bzw.
y′1 = ay1 + by2,
y′2 = cy1 + dy2.(2.7)
Die Eigenwerte von A sind
λ1,2 =(a+ d)±
√(a+ d)2 − 4 detA
2(detA = ad− bc).
Losungen (und ihr Verhalten in der Phasenebene = (y1, y2)-Ebene):
Fall 1: A besitzt zwei linear unabhangige Eigenvektoren v1, v2 (z.B. wenn die
Eigenwerte λ1, λ2 von A verschieden sind).
Die allgemeine Losung der gDG (2.7) ist dann
y(t) = α exp(λ1t)v1 + β exp(λ2t)v2 (α, β ∈ R).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 28
Fall 1a: λ1,2 sind verschieden und reell und haben das gleiche Vorzeichen.
Sei etwa λ2 < λ1 < 0. Jede Losung strebt gegen 0 (fur t→∞) und, wenn
α 6= 0,
y(t) =
[y1(t)
y2(t)
]∼ α exp(λ1t)v1 (t→∞).
D.h.: Genugend nahe am Ursprung streben alle Losungen gegen (0, 0) entlang
sign(α)v1 (in der Phasenebene), falls α 6= 0. Wenn α = 0 und β 6= 0, streben
sie auf sign(β)v2 gegen (0, 0). Wir sprechen von einem Knotenpunkt (Typ I).
Ist λ2 < λ1 ≤ 0, so ist der Knotenpunkt stabil, weil alle Kurven
(y1(t), y2(t))0≤t<∞ gegen (0, 0) streben fur t→∞.
Ist λ1 > λ2 > 0, so ist der Knotenpunkt instabil, denn
(y1(t), y2(t))0≤t<∞ → (∞,∞) fur t→∞.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 29
Fall 1b: λ1,2 sind reell und haben verschiedene Vorzeichen.
Sei etwa λ1 < 0 < λ2. Dann exp(λ1t)v1 → 0 fur t→∞ und
exp(λ2t)v2 →∞ fur t→∞. D.h. y(t)→∞ fur t→∞ (entlang v2), falls
β 6= 0. Wenn β = 0 und α 6= 0, dann y(t)→ 0 fur t→∞ auf v1. Man
spricht von einem Sattelpunkt. Er ist immer instabil.
Fall 1c: λ1,2 = λ sind gleich und damit reell (es gibt aber zwei linear
unabhangige Eigenvektoren).
In diesem Fall ist A = λI2 und jeder Vektor v 6= 0 ist ein Eigenvektor von A.
Die Losungen streben fur t→∞ auf Geraden gegen (0, 0), falls λ < 0
(stabiler Fall), bzw. von (0, 0) weg, falls λ > 0 (instabiler Fall). Man spricht
von einem Stern.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 30
Sind λ1,2 = µ1 ± iµ2 (µ1, µ2 ∈ R) konjugiert komplex, so ist die obige
Losungsformel zwar korrekt, liefert aber komplexe (d.h. fur uns unbrauchbare)
Losungen.
Ist v1 = w1 + iw2 (w1,w2 ∈ R2) Eigenvektor von A zum Eigenwert
λ1 = µ1 + iµ2, dann ist v2 = w1 − iw2 Eigenvektor zu λ2 = µ1 − iµ2.
Damit konstruieren wir eine reelle Basis des Losungsraums und erhalten als
allgemeine Losung von (2.7) (mit α, β ∈ R)
y(t) = exp(µ1t) α [cos(µ2t)w1 − sin(µ2t)w2] + β [sin(µ2t)w1 + cos(µ2t)w2] .
Fall 1d: λ1,2 = µ1 ± iµ2 sind konjugiert komplex mit µ1 6= 0.
Man spricht von einem Strudelpunkt, der stabil ist fur µ1 < 0 und instabil fur
µ1 > 0.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 31
Fall 1e: λ1,2 = µ1 ± iµ2 sind konjugiert komplex mit µ1 = 0 (µ2 6= 0).
Die Singularitat ist weder stabil noch instabil. Man nennt sie manchmal
neutral stabil und spricht von einem Zentralpunkt. Die Phasenkurven sind
Ellipsen mit Mittelpunkt (0, 0).
Fall 2: A besitzt nur einen linear unabhangigen Eigenvektor (und folglich nur
einen Eigenwert λ).
Die allgemeine Losung von (2.7) ist in diesem Fall
y(t) = α exp(λ1t)v1 + β exp(λ1t) (tv1 + v2) (α, β ∈ R).
Dabei ist v1 ein Eigenvektor von A (Av1 = λv1) und v2 ein zugehoriger
Hauptvektor (Av2 = λv2 + v1).
Man spricht von einem Knotenpunkt (Typ II), der stabil (λ1 < 0) oder instabil
(λ1 ≥ 0) sein kann.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 32
2.2.2 Phasenportraits
y1
y2
v1
v2
Knotenpunkt (Typ I), hier stabil
y1
y2
v1
v2
Sattelpunkt (immer instabil)
y1
y2
Sternpunkt, hier stabil
y1
y2
Strudelpunkt, hier stabil
y1
y2
Zentralpunkt (immer neutral stabil)
y1
y2
v1
Knotenpunkt (Typ II), hier stabil
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 33
2.2.3 Gleichgewichtslosungen
Beim logistischen Wachstumsmodell dNdt = rN(1− N
K
)sind zwei Losungen
ausgezeichnet: N(t) ≡ 0 und N1(t) ≡ K.
Man nennt sie die Gleichgewichtslosungen (GGL) oder Gleichgewichtspunkte
und erhalt sie durch Losen der Gleichung rN(1− NK ) = 0. N ist eine instabile
GGL (stort man N , so strebt die Losung weg von N). N1 ist eine stabile GGL
(stort man N1, so strebt die Losung wieder gegen N1).
Fur das allgemeinere Populationsmodell
dN
dt= f(N) mit einer (nichtlinearen) Funktion f (2.8)
(man spricht hier von einer autonomen gDG, d.h. f hangt nicht explizit von t
ab) ist N∗ eine GGL, wenn f(N∗) = 0. Sie ist instabil, wenn f ′(N∗) > 0 gilt,
und stabil, wenn f ′(N∗) < 0 gilt.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 34
Allgemein heißt eine beliebige Losung N = N(t) von (2.8) stabil, wenn fur
jede Losung N , die zur Zeit t = 0 genugend nahe bei N(t) startet, fur t→∞gegen N(t) strebt.
Genauer: ∀ε > 0 ∃δ = δ(ε) > 0, so dass fur alle Losungen N von (2.8) mit
|N(0)− N(0)| < δ immer |N(t)− N(t)| < ε ∀t ∈ I folgt.
Satz 2.3 (Stabilitat bei linearen gDGen)
Es gelten:
• Jede Losung von y ′ = Ay ist stabil, wenn alle Eigenwerte vonA negativen
Realteil haben.
• Jede Losung von y ′ = Ay ist instabil, wenn mindestens ein Eigenwert
von A positiven Realteil hat.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 35
2.2.4 Qualitative Analyse in Rezeptform
Wir betrachten die autonome gDG y ′ = f (y). Sie ist i.a. n-dimensional, wir
konzentrieren uns hier auf den Fall n = 2:
du
dt= f(u, v) und
dv
dt= g(u, v). (2.9)
1. Bestimme Gleichgewichtslosungen (Gleichgewichtspunkte), d.h. die
Losungen von (2.9), die nicht mit t variieren.
Lose das Gleichungssystem (in zwei Unbekannten)
f(u, v) = 0 und g(u, v) = 0.
2. Sind die Gleichgewichtspunkte stabil?
Besitzen f und g stetige partielle Ableitungen zweiter Ordnung bez. u und v in
einer Umgebung von (u∗, v∗), dann gibt es eine stetige Funktion H : R2 → R2
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 36
mit limx→0 H(x )/‖x‖ = 0 und[f(u∗ + h, v∗ + k)
g(u∗ + h, v∗ + k)
]=
[f(u∗, v∗)
g(u∗, v∗)
]+A
[h
k
]+H
([h
k
]).
Dabei ist
A = A(u∗, v∗) =
∂f(u∗,v∗)∂u
∂f(u∗,v∗)∂v
∂g(u∗,v∗)∂u
∂g(u∗,v∗)∂v
. (2.10)
Fur uns bedeutet das, dass sich die Losungen von (2.9) in einer Umgebung von
(u∗, v∗) im wesentlichen wie die von y ′ = Ay verhalten.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 37
Satz 2.4 (Stabilitat von Gleichgewichtslosungen)
Unter den obigen Voraussetzungen sei (u∗, v∗) eine Gleichgewichtslosung von
(2.9).
• (u∗, v∗) ist stabil, wenn alle Eigenwerte von A(u∗, v∗) negativen Realteil
haben.
• (u∗, v∗) ist instabil, wenn mindestens ein Eigenwert von A(u∗, v∗) posi-
tiven Realteil hat.
• Im nichtlinearen Fall sind keine Aussagen moglich, wenn alle Eigenwer-
te einen Realteil ≤ 0 besitzen, aber mindestens einer verschwindenden
Realteil besitzt.
Um den Charakter des Gleichgewichtspunktes beim Ubergang von y ′ = Ay
nach y ′ = f (y) zu erhalten (Knotenpunkt, Sattelpunkt etc.), mussen die
Voraussetzungen an f verscharft werden (vgl. die folgenden Beispiele).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 38
3. Phasenportrait.
Bestimme Bahnen oder Trajektorien (u(t), v(t))t≥0, wobei u, v Losungen
von (2.9) sind, uberdv
du=g(u, v)
f(u, v).
Satz 2.5 (Eigenschaften von Bahnen)
Es gilt
• Durch jeden Punkt (u0, v0) lauft eine Bahn (außer durch Gleichgewichts-
punkte, fur die die Bahnen zu einem Punkt degenerieren). Haben zwei
Bahnen einen gemeinsamen Punkt, so sind sie identisch.
• Kehrt eine Bahn nach einer Zeit T > 0 zu ihrem Ausgangspunkt zuruck,
dann ist sie geschlossen (d.h. die zugehorigen Losungen sind periodisch).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 39
Satz 2.6 (Satz von Poincare-Bendixson.)
Bleibt eine Losung (u, v) von (2.9) in einem beschrankten Gebiet der Ebene,
das keine stabilen Gleichgewichtspunkte von (2.9) enthalt, dann dreht sich ihre
Bahn in eine einfach geschlossene Kurve hinein, die Bahn einer periodischen
Losung von (2.9) ist.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 40
2.3 Western spruce budworm (choristoneura occidentalis)
einer der gefahrlichsten Nadelwaldschadlinge in Nordamerika; periodische
Ausbruche (etwa alle 40 Jahre) mit dramatischen Folgen fur die
Forstwirtschaft.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 41
2.4 Modell fur Populationswachstum
dN
dt= rBN
(1− N
KB
)− p(N) mit p(N) =
BN2
A2 +N2. (2.11)
B ist ein Mass fur die Rauber-Effektivitat der Vogel, A heißt Schalterwert.
0 A N 0
B/2
B
p(N)
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 42
Dieses Problem enthalt vier Parameter:
A,KB [Einheit von N ], rB [1/t], B [Einheit von N/t].
Obligatorisch: Transformation auf dimensionslose Großen (dazu spater mehr).
Setze hier: u = NA , r = ArB
B , q = KB
A , τ = BtA . Einsetzen liefert:
du
dτ= ru
(1− u
q
)− u2
1 + u2= f(u; r, q).
Gleichgewichtslosungen:
f(u; r, q) = 0 ⇒ ru(1− uq ) = u2
1+u2 ⇒ u = 0 ∨ r(1− uq ) = u
1+u2 .
0 q u0
r
N1
N2
N3
r/(1−u/q)
u/(1+u2)
0 u
N1
N2
N3
f(u;r,q)
(q,r)=(8,0.4)(q,r)=(8,0.5)(q,r)=(8,0.6)
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 43
Quelle: James D. Murray. Mathematical Biology, 2nd corrected edition.
Springer-Verlag, Berlin 1993
Hysterese Effekt: Fixiere q und variiere r entlang ABCD: Die
Gleichgewichtswerte”springen“ in C. Variiert man r entlang DCBA, so findet
ein Sprung bei B statt.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 44
2.5 Modelle fur zwei PopulationenGenauer: Populationen, deren Großen N und P sich wechselweise beeinflussen.
• Rauber-Beute-Modelle,
• Wettbewerbs-Modelle,
• Symbiose-Modelle.
Umberto d’Ancona stellt 1925 den prozentualen Anteil der Haie am
Gesamtfang (Speisefische und Haie) im Hafen von Triest fest:
1914 1915 1916 1917 1918 1919 1920 1921 1922 19230
5
10
15
20
25
30
35
40
Benachteiligt eingeschrankter Fischfang (1. Weltkrieg) die Speisefische?
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 45
N(t): Beutepopulation zur Zeit t (Speisefische),
P (t): Rauberpopulation zur Zeit t (Haie).
Annahmen (Volterra-Lotka-Modell):
• Ohne Rauber wurde die Beute exponentiell wachsen.
• Rauber reduzieren Beute proportional zu N(t)P (t).
• Ohne Beute wurden die Rauber exponentiell abnehmen.
• Rauber profitieren von der Beute proportional zu N(t)P (t).
dN
dt= aN(t)− bN(t)P (t) (mit Konstanten a, b > 0)
dP
dt= −c P (t) + eN(t)P (t) (mit Konstanten c, e > 0)
(2.12)
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 46
Transformation
τ = at, u(τ) =eN(t)
c, v(τ) =
bP (t)
a, α =
c
a
liefert dimensionslose Form
du
dτ= u(1− v),
dv
dτ= αv(u− 1). (2.13)
Wir bestimmen die Bahnen von (2.13), d.h. die Losungen von
dv
du= α
v(u− 1)
u(1− v)bzw.
1− vv
dv = αu− 1
udu
(gDG mit getrennten Variablen) und erhalten
log(v)− v +H = αu− α log(u) bzw. αu+ v − log(uαv) = H (2.14)
mit einer Konstanten H.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 47
Es gelten:
• Die Gleichung (2.14) hat nur fur H ≥ 1 + α positive Losungen u und v.
• (2.14) beschreibt fur H > 1 + α eine geschlossene Kurve in der
Phasenebene. Folglich sind die Losungen von (2.13) periodisch.
(Beachte auch: u(τ) = 1 ⇒ v′(τ) = 0 und v(τ) = 1 ⇒ u′(τ) = 0.)
• Ist T die gemeinsame Periodenlange von u bzw. v, dann
1
T
∫ T
0
u(τ) dτ = 1 und1
T
∫ T
0
v(τ) dτ = 1.
Gleichgewichtslosungen:
f(u, v) = u(1− v) = 0
g(u, v) = αv(u− 1) = 0⇒ (u, v) = (0, 0) oder (u, v) = (1, 1).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 48
Lineare Stabilitatsanalyse: ∂f∂u
∂f∂v
∂g∂u
∂g∂v
=
[1− v −uαv α(u− 1)
].
Das bedeutet ∂f∂u
∂f∂v
∂g∂u
∂g∂v
(u,v)=(0,0)
=
[1 0
0 −α
]mit Eigenwerten 1 und − α.
∂f∂u
∂f∂v
∂g∂u
∂g∂v
(u,v)=(1,1)
=
[0 −1
α 0
]mit Eigenwerten ±
√αi.
Also ist ein (0, 0) ein (instabiler) Sattelpunkt und (1, 1) ein Zentralpunkt.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 49
Rucktransformation ergibt fur (2.12) die Bahnen
eN + bP − log(N cP a) = K
und die Mittelwerte
1
T
∫ T
0
N(t) dt =c
eund
1
T
∫ T
0
P (t) dt =a
b.
a/b
c/e
N(t)
P(t)
c/e
a/b Mittelwert
Beute
Räu
ber
Phasenebene
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 50
Berucksichtige Fischfang:
dN
dt= aN − bNP − f N = (a− f)N − bNP,
dP
dt= −c P + eNP − f P = −(c+ f)P + eNP (f > 0).
Gleiches System mit neuen Koeffizienten: a→ a− f und c→ c+ f .
Mittelwerte: (c+ f)/e > c/e (Beute), (a− f)/b < a/b (Rauber).
c/e (c+f)/e
(a−f)/b
a/b
Beute
Rae
uber
ohne
mit
Fischfang
Volterras Prinzip:
Angemessener Fischfang
(f < a) steigert die durch-
schnittliche Zahl der Spei-
sefische und reduziert die
durchschnittliche Zahl der
Haie.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 51
Ein realitatsnaheres Rauber-Beute-Modell hat die Form
dN
dt= Nr
(1− N
K
)− kNP
N +D(r,K, k,D > 0),
dP
dt= Ps
(1− hP
N
)(s, h > 0).
Nach Substitution u(τ) = N(t)/K, v(τ) = hP (t)/K, τ = rt, a = k/(hr),
b = s/r, d = D/K ergibt sich
du
dτ= f(u, v) = u(1− u)− auv
u+ dund
dv
dτ= g(u, v) = bv
(1− v
u
).
(2.15)
Der einzige Gleichgewichtspunkt mit positiven Komponenten ist
(u∗, v∗) mit u∗ = v∗ =1− a− d+
√(1− a− d)2 + 4d
2.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 52
Die zugehorige Funktionalmatrix ist ∂f∂u
∂f∂v
∂g∂u
∂g∂v
(u,v)=(u∗,v∗)
=
u∗[
au∗
(u∗+d)2 − 1]−au∗u∗+d
b −b
.(u∗, v∗) ist genau dann instabil, wenn die Parameter (a, b, d) in der skizzierten
Teilmenge des R3 liegen:
00.1
0.20.3
0.40.5
0.60.7
0.8 00.5
11.5
22.5
0
0.2
0.4
0.6
0.8
1
ad
b
d(a)=(a2+4a)1/2−(1+a)
b=1/a b=2a−1
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 53
a = 1, b = 0.05, d = 0.2 liefert z.B. ein instabiles Gleichgewicht
u∗ = v∗ ≈ .36:
0 u* .8 10
v*
.7
1
g<0
g>0
f<0
f>0
0 50 100 150 200 250 3000
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
u
τ=rt
v
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 54
2.6 Ein Konkurrenzmodell
dN1
dt= r1N1
(1− N1
K1
)− a1N1N2 (r1,K1, a1 > 0),
dN2
dt= r2N2
(1− N2
K2
)− a2N1N2 (r2,K2, a2 > 0)
wird dimensionslos durch u1 = N1
K1, u2 = N2
K2, τ = r1t, ρ = r2
r1, b1 = a1
K2
r1,
b2 = a2K1
r2:
du1
dτ= f(u1, u2) = u1 (1− u1)− b1u1u2,
du2
dτ= g(u1, u2) = ρu2 (1− u2)− ρb2u1u2.
Wir betrachten den Fall b1 > 1 und b2 > 1 (es gibt drei weitere Falle, vgl.
Ubungsaufgaben).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 55
1. Gleichgewichtspunkte
u1 (1− u1)− b1u1u2 = 0⇔ u1 = 0 ∨ u2 =1
b1(1− u1) ,
ρu2 (1− u2)− ρb2u1u2 = 0⇔ u2 = 0 ∨ u2 = 1− b2u1.
Das bedeutet: Im ersten Quadranten gibt es vier Gleichgewichtspunkte; (0, 0),
(1, 0), (0, 1) und einen weiteren mit positiven Komponenten (der Schnittpunkt
der Nulllinien), namlich ( b1−1b1b2−1 ,
b2−1b1b2−1 ).
2. Stabilitat
A(u1, u2) =
∂f∂u1
∂f∂u2
∂g∂u1
∂g∂u2
=
1− 2u1 − b1u2 −b1u1
−ρb2u2 ρ(1− 2u2 − b2u1)
.A(0, 0) =
[1 0
0 ρ
], A(1, 0) =
[−1 −b10 ρ(1− b2)
], A(0, 1) =
[1− b1 0
−ρb2 −ρ
].
Also ist (0, 0) instabil, wahrend (1, 0) und (0, 1) stabil sind.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 56
A
(b1 − 1
b1b2 − 1,b2 − 1
b1b2 − 1
)=
1
1− b1b2
[b1 − 1 b1(b1 − 1)
ρb2(b2 − 1) ρ(b2 − 1)
]
mit Eigenwerten λ2 < 0 < λ1.
Es handelt sich also um einen instabilen Sattelpunkt.
Interpretation: Die beiden Spezies konkurrieren so aggressiv, dass eine von
beiden ausstirbt (Ausschluss durch Konkurrenz).
Es gibt eine Bahn (Separatrix), die (theoretisch) gegen das instabile
Gleichgewicht strebt. Alle anderen Punkte liegen im Einzugsbereich eines der
beiden Attraktoren (1, 0) (d.h. Spezies 2 stirbt aus) oder (0, 1) (d.h. Spezies 1
stirbt aus).
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 57
0 u_1* 1/b_2 10
u_2*
1/b_1
1
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 58
2.7 Diskrete Modelle
Fur viele Arten ist es sinnvoll anzunehmen, dass die Population in diskreten
Zeitschritten wachst:
Nt+1 = f(Nt), hier f(Nt) = Nt · F (Nt). (2.16)
(Oft ergibt sich die diskrete Form aber auch durch Diskretisierung einer gDG.)
Beispiele sind:
Nt+1 = rNt (fur r > 0) ⇒ Nt = rtN0 (diskretes exponentielles Wachstum),
Nt+1 = rNt(1− Nt
K
)(fur r > 0, K > 0) (diskretes logistisches Wachstum).
Nachteil des letzten Modells: Nt+1 < 0, wenn Nt > K.
Gleichgewichtspunkte erhalt man als Schnittpunkte N∗ von Nt+1 = f(Nt) mit der
”ersten Winkelhalbierenden“ Nt+1 = Nt.
Der Gleichgewichtspunkt N∗ ist stabil, wenn |f ′(N∗)| < 1 gilt, und instabil, wenn
|f ′(N∗)| > 1 gilt. Im Fall von |f ′(N∗)| = 1 spricht man von einem
Verzweigungspunkt.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 59
(((!!
""
x1 x2 x3
0 < f ′(N∗) < 1
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 60
HHaa
aaaH
HHH
QQQ
cc
x1 x2x3
−1 < f ′(N∗) < 0
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 61
SS
eee
@@eeeSSTTTAAAA
x1x2 x3
|f ′(N∗)| > 1
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 62
2.8 Ein”
chaotisches“ Beispiel
Normiere das diskrete logistische Modell durch ut = Nt
K :
ut+1 = rut(1− ut) (mit r > 0)
und wahle u0 ∈ (0, 1).
Gleichgewichtspunkte (bestimme Losungen u∗ von f(u) = ru(1− u) = u),
Stabilitat (bestimme λ = f ′(u∗)):
u∗1 = 0, λ1 = f ′(0) = r, und u∗2 =r − 1
r, λ2 = f ′
(r − 1
r
)= 2− r. (2.17)
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 63
1. Fur r ∈ (0, 1) ist u∗1 = 0 der einzige realistische (positive)
Gleichgewichtspunkt,
2. Erster Bifurkationspunkt r = 1.
3. Fur r ∈ [1, 3) ist u∗2 ein (realistischer) stabiler Gleichgewichtspunkt.
0 2 4 6 8 10 12 14 16 18 200
0.1
0.2
0.3
0.4
0.5
0.6
0.7
r = 0.5r = 1.0r = 2.0r =2.8
0 10
11 <= r <= 3
ut
ut+1
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 64
4. Zweiter Bifurkationspunkt r = 3.
5. Fur r ∈ [3, r4) existiert ein stabile 2-periodische Gleichtgewichtlosung.
Betrachte ut+2 = f2(ut) = r(rut(1− ut))(1− rut(1− ut)). Diese
Iteration hat vier Gleichgewichtspunkte: u∗1,2 (beide instabil) und
u∗3,4 =(r + 1)± [(r + 1)(r − 3)]1/2
2r> 0
mit |(f2)′(u∗3,4)| < 1.
0 5 10 15 20 25 30 35 400.3
0.4
0.5
0.6
0.7
0.8
0.9
1r = 3.0r = 3.3r = 3.5
0 10
1r = 3.3
ut
ut+2
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 65
6. Der dritte Bifurkationspunkt r4 ist der (kleinste) Wert von r mit
(f2)′(u∗3,4) = −1.
7. Zwischen r4 und r8 gibt es eine stabile 4-periodische Gleichgewichtslosung
u∗5,6,7,8 bis (f4)′(u∗5,6,7,8) = −1.
8. Dieser Prozess (aus einer 2k-periodischen Gleichgewichtslosung wird eine
2k+1-periodische) setzt sich fort bis zum Wert rc.
0 10 20 30 40 50 60 70 800.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
r = rc
0 10
1
r = rc
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 66
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 67
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 68
Ein Modell fur Masernepidemien
Die Kinderkrankheit Masern wird durch Viren ubertragen. Ein typischer
Krankheitsverlauf stellt sich (leicht vereinfacht) wie folgt dar: Ein Kind, das
sich angesteckt hat, durchlauft zunachst eine einwochige Inkubationszeit, in
der es weder sichtbar krank ist noch andere Kinder anstecken kann. Danach
wird es fur eine Woche infektios (d.h. es kann die Krankheit ubertragen).
Danach wird es sichtbar krank (Exanthem), aber nicht mehr (in Wirklichkeit
deutlich weniger) infektios, erholt sich wieder und ist in Zukunft immun gegen
Masern.
Man hat beobachtet, dass es in Industrielandern (U.S.A., U.K.) alle zwei bis
drei Jahre zu epidemischen Ausbruchen dieser Krankheit kommt. In sog.
Dritte-Welt-Landern ist der Zeitraum zwischen zwei Epidemien wesentlich
kurzer. Warum?
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 69
Sei Sn die Anzahl der gesunden, aber nicht immunen Kinder in der Woche n.
Sei In die Anzahl der infektiosen Kinder in der Woche n. Sei B > 0 die Anzahl
der Neugeborenen pro Woche. Sei µ > 0 der Anteil der Kinder (aus Sn), die
pro Woche von einem infektiosen Kind angesteckt werden.
In = f(In−1, Sn−1) = µIn−1Sn−1
Sn = g(In−1, Sn−1) = Sn−1 − µIn−1Sn−1 +B.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 70
I0 = 20, S0 = 30000, µ = 0.3 · 10−4, B = 120 (U.K.):
0 50 100 150 200 250 3002.8
3
3.2
3.4
3.6
3.8x 104
0 50 100 150 200 250 3000
100
200
300
400
500
600
700
Sn
In
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 71
I0 = 20, S0 = 30000, µ = 0.3 · 10−4, B = 360 (Nigeria):
0 50 100 150 200 250 3002.5
3
3.5
4
4.5x 104
0 50 100 150 200 250 3000
500
1000
1500
2000
Sn
In
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 72
Wir berechnen die Gleichgewichtspunkte aus
f(I, S) = I und g(I, S) = S.
Als einziger realistischer Gleichgewichtspunkt ergibt sich (I∗, S∗) = (B, 1/µ),
der wegen ∂f∂I
∂f∂S
∂g∂I
∂g∂S
(I,S)=(I∗,S∗)
=
[1 µB
−1 1− µB
]
ein Zentralpunkt ist (falls µB ≤ 4 sind die Eigenwerte der Funktionalmatrix
konjugiert komplex und vom Betrag 1). Wir beobachten eine periodisches
Verhalten mit Periodenlange ≈ 126 im Fall B = 120 bzw. Periodenlange ≈ 73
im Fall B = 360.
Populationsdynamik Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 73
3 Markoff-Ketten
3.1 Zwei Zufallsprozesse
a. Rekursionen mit Zufallsmatrizen.
Wahle
• k ∈ N, k ≥ 2,
• ein Wahrscheinlichkeitsmaß π = [π1, π2, . . . , πk] auf 1, 2, . . . , k,• k Matrizen A1, A2, . . . , Ak ∈ R2×2,
• k Vektoren b1, b2, . . . , bk ∈ R2,
• und einen Startvektor x0 ∈ R2.
Dann fuhre folgenden Algorithmus aus
for m = 1 : M (M groß, z.B. M = 10000)
wahle j ∈ 1, 2, . . . , k nach π
xm = Ajxm−1 + bjend
Plotte xm, m = m0,m0 + 1, . . . ,M (mit z.B. m0 = 100).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 74
Ein Farnblattk = 2, π = [.2993, .7007],
A1 =
[+.4000 −.3733
+.0600 +.6000
],
b1 =
[+0.3533
+0.0000
],
A2 =
[−.8000 −.1867
+.1371 +.8000
],
b2 =
[+1.1000
+0.1000
].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 75
Sierpinski-Dreieck
k = 3, π = [1, 1, 1]/3,
A1 = 12
[1 0
0 1
],
b1 =
[0
0
],
A2 = A1,
b2 =
[1
0
],
A3 = A1,
b3 = 12
[1
1
].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 76
b. Irrfahrt eines Springers.
Stelle einen Springer in die Ecke eines ansonsten leeren Schachbretts und
wahle zufallig einen der nach den Regeln erlaubten Zuge. Jeder der moglichen
Zuge werde mit der gleichen Wahrscheinlickeit ausgewahlt.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 77
Naturliche Fragen:
1. Wieviele Zuge sind durchschnittlich notwendig, bis der Springer zu seiner
Ausgangsposition zuruckkehrt?
2. Wieviele Zuge sind durchschnittlich notwendig, bis der Springer alle 64
Felder besucht hat?
Um ein Gefuhl fur diese Großen zu entwickeln, simulieren wir 100 solche
Irrfahrten mit dem Zufallszahlengenerator rand von Matlab. In genau 35
dieser Experimente ist der Springer nach weniger als 50 Zugen zu seiner
Ausgangsposition zuruckkehrt, wahrend er in genau einem Experiment dazu
k ∈ [1050, 1100) Zuge benotigt. Die durchschnittliche first return time unserer
Stichprobe ist 191.78, aber wir werden spater beweisen, dass diese Frage ohne
numerische Experimente beantwortet werden kann (der Erwartungswert der
first return time betragt 168).
Eine ahnliche Simulation (mit wieder 100 Experimenten) liefert als Antwort
auf die zweite Frage nach der cover time den Durchschnittswert 556.9 (auch
hier kann die Antwort ohne numerische Experimente gegeben werden).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 78
1 2 3 4 5 6 7 8 9 1002468
10121416182022242628303234
x102 2 3 4 5 6 7 8 9 10 11 120
3
6
9
12
15
x102
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 79
3.2 Elementares aus der Wahrscheinlichkeitstheorie
Wir untersuchen folgendes Spiel: Nach Zahlung eines Einsatzes darf mit drei
Wurfeln gewurfelt werden. Fur jede 4 wird 1 Euro, fur jede 5 bzw. 6 werden 2
bzw. 3 Euro ausgezahlt.
Wie hoch muss der Einsatz sein, damit das Spiel fair ist?
zugehoriger Wahrscheinlichkeitsraum (Ω,A ,W ):
• Ergebnisraum Ω := (w1, w2, w3) : wj ∈ 1, 2, ..., 6(63 = 216 Elementarereignisse),
• Ereignisalgebra A := A : A ⊆ Ω(Potenzmenge von Ω, 2216 ≈ 1065 Ereignisse),
• Wahrscheinlichkeitsmaß W (A) = 1216 card(A), A ∈ A .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 80
Zufallsvariable X : Ω→ R : Wurf 7→ Auszahlung,
X((w1, w2, w3)) =∑wj=4 1 +
∑wj=5 2 +
∑wj=6 3.
Offenbar: X(Ω) = 0, 1, ..., 9.
diskrete Dichtefunktion f (Wahrscheinlichkeitsverteilung)
f(x) := W (X = x) = W ((w1, w2, w3) ∈ Ω : X((w1, w2, w3)) = x),x ∈ 0, 1, ..., 9. Beachte
∑9x=0 f(x) = 1.
Im folgenden steht χ fur eine Augenzahl zwischen 1 und 3.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 81
x W (X = x)
0 (χ, χ, χ) [27|216] 27/216
1 (χ, χ, 4) [27|216] 27/216
2 (χ, 4, 4) [9|216]; (χ, χ, 5) [27|216] 36/216
3 (4, 4, 4) [1|216]; (χ, 4, 5) [18|216]; (χ, χ, 6) [27|216] 46/216
4 (4, 4, 5) [3|216]; (χ, 4, 6) [18|216]; (χ, 5, 5) [9|216] 30/216
5 (4, 4, 6) [3|216]; (4, 5, 5) [3|216]; (χ, 5, 6) [18|216] 24/216
6 (4, 5, 6) [6|216]; (5, 5, 5) [1|216]; (χ, 6, 6) [9|216] 16/216
7 (4, 6, 6) [3|216]; (5, 5, 6) [3|216] 6/216
8 (5, 6, 6) [3|216] 3/216
9 (6, 6, 6) [1|216] 1/216
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 82
Verteilungsfunktion F
F (y) := W (X ≤ y) = W ((w1, w2, w3) ∈ Ω : X((w1, w2, w3)) ≤ y)=∑x≤y f(x), y ∈ R.
Beachte limy→−∞ F (y) = 0, limy→+∞ F (y) = 1, F (y1) ≤ F (y2) fur y1 ≤ y2.
Erwartungswert von X
E(X) =∑9x=0 xf(x) = 3, d.h. 3 Euro Einsatz ⇒ faires Spiel.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 83
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25f(x)= W(X=x)
0 1 2 3 4 5 6 7 8 90
0.2
0.4
0.6
0.8
1
F(x)=W(X <= x)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 84
Warteschlangen
Zu den zufalligen Zeitpunkten t0 = 0, t1, t2, ... treffen Kunden in der
Warteschlange vor einem Bedienungssystem ein. Wir setzen
Sn := tn − tn−1 (n = 1, 2, ...)
(Zwischenankunftszeiten).
Die Zufallsvariablen Sn seien unabhangig und identisch verteilt (iid).
Ihre Verteilung sei F (x) = W (Sn ≤ x). Beliebt ist die Exponentialverteilung
F (x) = Eλ(x) := 1− exp(−λx), x ≥ 0 (mit einem Parameter λ > 0).
Fλ(x) besitzt eine Dichtefunktion, namlich fλ(t) = λ exp(−λx) (fur x ≥ 0)
(Fλ(x) =∫ x
0fλ(t)dt).
Die Exponentialverteilung”besitzt kein Gedachtnis“:
W (X > x+ y |X > x) = W (X > y), x, y ≥ 0,
was spater noch sehr wichtig wird.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 85
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
λ = 2
λ = .75
Exponentialverteilung
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2Dichte der Exponentialverteilung
λ = 2
λ = .75
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 86
Zuruck zur Warteschlange: Die wartenden Kunden werden vom
Bedienungssystem in der Reihenfolge ihres Eintreffens bedient. Durch Tn wird
die Bedienungszeit des Kunden n bezeichnet. Diese Zufallsvariable besitze
die Verteilung G. Auch die Tn seien iie und unabhangig von Sn.
Wir interessieren uns fur die Lange der Warteschlange X(t). Sind auch die Tnexponentialverteilt (mit Parameter µ), spricht man von einem
M/M/1/∞-System (”1“ bedeutet: ein Server,
”∞“ bedeutet: unendlich
grosser Warteraum).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 87
0 1 2 3 4 5 6 7 8 9 10
0
1
2
Warteschlange
t
X(t
)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 88
3.3 Diskrete Markoff-Ketten
Ein System S , das sich in genau einem von n Zustanden Z1, Z2, . . . , Zn (oder
kurzer 1, 2, . . . , n) befinden kann, wird einem stochastischen Prozess A
unterworfen. Bekannt sind die Ubergangswahrscheinlichkeiten
ai,j = W(i
A−→ j)
(1 ≤ i, j ≤ n).
ai,j hange nicht davon ab, wie S in den Zustand i kam. Der zukunftige
Zustand des Prozesses hangt also nur vom gegenwartigen, aber nicht von den
vergangenen ab (Markoff-Eigenschaft).
Das Paar (S ,A ) heißt (zeit-) diskrete homogene (d.h. ai,j ist zeitunabhangig)
Markoff-Kette (mit endlichem Zustandsraum 1, 2, . . . , n). Wir werden nur
solche Markoff-Ketten betrachten.
Die Ubergangsmatrix A = [ai,j ] ∈ Rn×n beschreibt diesen Prozess.
A ist eine stochastische Matrix, d.h. A ist nichtnegativ (ai,j ≥ 0, 1 ≤ i, j ≤ n)
mit Zeilensummen gleich 1 (∑nj=1 ai,j = 1, 1 ≤ i ≤ n, oder Ae = e mit
e = [1, 1, . . . , 1]T ∈ Rn).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 89
Zusammengesetzte Prozesse
Sind A und B stochastische Vorgange im System S mit den
Ubergangsmatrizen A bzw. B, dann besitzt der zusammengesetzte Prozess
C = B A die Ubergangsmatrix C = AB. (Insbesondere ist das Produkt
zweier stochastischer Matrizen wieder stochastisch.)
Hier sollen Prozesse A m mit den Ubergangsmatrizen Am
(m-Schritt-Ubergangswahrscheinlichkeiten) und ihr Verhalten fur m→∞untersucht werden:
SA−→ S
A−→ · · · A−→ SA−→ · · ·
Wird das System S im Zeitpunkt m = 0 durch das Wahrscheinlichkeitsmaß
(die Ausgangsverteilung) π(0) = [π(0)1 , π
(0)2 , . . . , π
(0)n ] beschrieben
(Wahrscheinlichkeitsvektoren sind hier stets Zeilenvektoren), dann beschreibt
π(m) = π(0)Am das System S im Zeitpunkt m ≥ 0.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 90
Beispiel. Eine 1-bit Nachricht (”ja“ oder
”nein“ bzw.
”wahr“ oder
”falsch“
bzw.”0“ oder
”1“) wird z.B. mundlich uber viele Zwischenstationen verbreitet
und dabei moglicherweise verfalscht: Mit Wahrscheinlichkeit p (0 ≤ p ≤ 1)
wird”0“ als
”1“ weitergereicht und folglich ist die Wahrscheinlichkeit 1− p,
dass”0“ korrekt als
”0“ wiedergegeben wird. Die
Verfalschungswahrscheinlichkeit von”1“ sei q (0 ≤ q ≤ 1) und deshalb wird
”1“ mit Wahrscheinlichkeit 1− q als
”1“ weitergegeben.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 91
Unser Ubertragungssystem besteht also aus zwei Zustanden, namlich 0 und 1,
und wird durch die vier Ubergangswahrscheinlichkeiten
W (0→ 0) = 1− p, W (0→ 1) = p,
W (1→ 0) = q, W (1→ 1) = 1− q
beschrieben, die wir in der Ubergangsmatrix
A =
[1− p p
q 1− q
]∈ R2×2
zusammenfassen.
Drei Falle:
1. p = q = 0, dann Am = I2 ∀m und limm→∞Am = I2.
2. p = q = 1, dann A2m+1 = A, A2m = I2 ∀m und limm→∞Am existiert
nicht.
3. 0 < p, q < 1, dann limm→∞Am = 1p+q [ q pq p ].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 92
Im Fall von p = q ∈ (0, 1) (unparteiische Verfalschung) ist
limm→∞
Am =
12
12
12
12
=: A∞.
Nach einer langen Reihe von Zwischentragern kann die Ankunft der
unverfalschten Ausgangsnachricht nur mit der Wahrscheinlichkeit 1/2 erwartet
werden (und zwar unabhangig vom Wert von p = q).
Was”lange“ bedeutet, hangt aber vom Wert von p = q ab:
‖Am −A∞‖∞ = |1− 2p|m,
d.h. die Konvergenz ist umso schneller, je naher p bei 1/2 liegt.
(Interessante Beobachtung: A besitzt die Eigenwerte 1 und 1− 2p.)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 93
Markoff-Ketten konnen alternativ durch eine Folge Xtt=0,1,... von
Zufallsvariablen beschrieben werden. Dabei ist Xt der Zustand, in dem sich
das System zum Zeitpunkt t befindet.
Die Markoff-Eigenschaft ist dann
W (Xt+1 = it+1|Xt = it, Xt−1 = it−1, . . . , X0 = i0) = W (Xt+1 = it+1|Xt = it).
Die Eintrage der Ubergangsmatrix A konnen als bedingte Wahrscheinlichkeiten
ai,j = W (Xt+1 = j|Xt = i)
interpretiert werden.
Unter der Bedingung Xt = i besitzt die Zufallsvariable Xt+1 die
Wahrscheinlichkeitsverteilung W (Xt+1 = j) = ai,jj=1,2,...,n.
Xt+1 ist stochastisch unabhangig von X0, X1, . . . , Xt−1.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 94
Beliebte Beispiele fur Markoff-Ketten sind Irrfahrten wie die des Springers aus
dem 1. Abschnitt. Dort besteht der Zustandsraum aus 64 Elementen, den
Feldern des Schachbretts.
Die folgenden Abbildungen zeigen die Ubergangsmatrix fur zwei
Nummerierungen der Zustande (Felder). Bei der lexikografischen Ordnung
werden die Felder von links unten (a1) nach rechts oben (h8) zeilenweise
durchnummeriert. Bei der”Schwarz/Weiss-Ordnung“ nummeriert man erst
alle schwarzen und danach alle weissen Felder (jeweils lexikografisch).
0 10 20 30 40 50 60
0
10
20
30
40
50
60
nz = 336
Lexikographische Ordnung
0 10 20 30 40 50 60
0
10
20
30
40
50
60
nz = 336
"Schwarz/Weiß"−Ordnung
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 95
Klassifizierung von Zustanden (und Markoff-Ketten)
Erinnerung: Mit Am =[a
(m)i,j
]1≤i,j≤n
ist W(i
A m
−→ j)
= a(m)i,j .
Der Zustand j heißt vom Zustand i aus erreichbar (i→ j), wenn es ein m ∈ Ngibt mit a
(m)i,j > 0. Gilt i→ j und j → i, so kommunizieren die Zustande i
und j (i↔ j).
”↔“ ist eine Aquivalenzrelation auf dem Zustandsraum 1, 2, . . . , n und
zerlegt ihn damit in Aquivalenzklassen (sog. Kommunikationsklassen).
1, 2, . . . , n = K1 ∪K2 ∪ · · · ∪Ks, Kk 6= ∅, Kk ∩K` = ∅ fur k 6= `,
wobei alle Zustande in Kk miteinander kommunizieren und kein Zustand aus
Kk mit einem Zustand kommuniziert, der außerhalb von Kk liegt.
Die Markoff-Kette heißt irreduzibel, wenn alle ihre Zustande miteinander
kommunizieren, d.h. wenn es nur eine Kommunikationsklasse gibt. Sonst heißt
sie reduzibel.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 96
3.4 Reduzible und irreduzible Matrizen
Ist π eine Permutation der Indexmenge 1, 2, . . . , n dann ist die zugehorige
Permutationsmatrix P = Pπ = [pi,j ] ∈ Rn×n durch
pi,j =
1, falls j = π(i),
0, falls j 6= π(i),
definiert. Es gilt P−1π = Pπ−1 = PTπ .
Ist A = [ai,j ] ∈ Cn×n mit den Zeilen zT1 , . . . , zTn und Spalten s1, . . . , sn, so
besitzt PA die Zeilen zTπ(1), . . . , zTπ(n) und APT die Spalten sπ(1), . . . , sπ(n).
Insbesondere ist PAPT = [aπ(i),π(j)].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 97
Eine Matrix A ∈ Cn×n heißt reduzibel (zerlegbar), falls eine
Permutationsmatrix P existiert, so dass
PAPT =
[A1,1 O
A2,1 A2,2
]mit quadratischen Matrizen A1,1, A2,2 der Ordnung p > 0 und q > 0
(p+ q = n). A heißt irreduzibel, wenn A nicht reduzibel ist.
Wir ordnen der Matrix A = [ai,j ] ∈ Cn×n n Knoten N1, N2, . . . , Nn zu. Fur
jedes Element ai,j 6= 0 wird Ni mit Nj durch eine gerichtete Kante verbunden.
Falls ai,i 6= 0 ist, so wird Ni mit einer geschlossenen Schleife versehen.
Die Gesamtheit der Knoten und Kanten bildet den gerichteten Graph G(A)
von A.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 98
A =
0 0 1 0
0 0 1/2 1/2
1/3 1/3 1/3 0
1/3 1/3 0 173
⇒ G(A) :
N2
N1
N3
N4
6
?
6
@@@@@@R@@
@@
@@I
-
Ein gerichteter Graph heißt vollstandig zusammenhangend, wenn es fur jedes
geordnete Paar (Ni, Nj) von Knoten einen gerichteten Weg (d.h. eine Folge
von gerichteten Kanten) gibt, der von Ni nach Nj fuhrt.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 99
Satz 3.1 (Charakterisierung von Irreduzibilitat)
Die folgenden Aussagen sind einander aquivalent:
1. Die Markoff-Kette (S ,A ) ist irreduzibel.
2. Die zugehorige Ubergangsmatrix A ist irreduzibel.
3. Der gerichtete Graph von A ist vollstandig zusammenhangend.
Aquivalent: A = [ai,j ] ist genau dann irreduzibel, wenn es fur jede beliebige
Zerlegung Z = S ∪T der Indexmenge (hier: des Zustandraums)
Z := 1, 2, . . . , n mit S ∩T = ∅, S 6= ∅, T 6= ∅ ein Element ai,j 6= 0 gibt
mit i ∈ S und j ∈ T .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 100
Zu jeder Matrix A ∈ Cn×n gibt es eine Permutationsmatrix P ∈ Rn×n, so dass
PAP> =
A1,1 0 · · · · · · 0
A2,1 A2,2 0 0
A3,1 A3,2 A3,3
......
. . ....
Ak,1 Ak,2 Ak,3 · · · Ak,k
untere Blockdreiecksform besitzt. Die quadratischen Diagonalblocke Aj,j(j = 1, 2, . . . , k) sind dabei entweder irreduzibel oder eindimensionale
Nullmatrizen (Normalform einer reduziblen Matrix).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 101
Sei (S ,A ) eine Markoff-Kette mit Zustandsraum 1, 2, ..., n .
Die Zufallsvariable Tj sei die Zeit, zu der der Prozess zum ersten Mal in den
Zustand j kommt und fi,j(k) = W (Tj = k|X0 = i) sei die Wahrscheinlichkeit,
dass der Prozess, wenn er im Zustand i startet, zur Zeit k zum ersten Mal
nach j kommt. Dann gelten
fi,j(k) =
ai,j , k = 1,∑l 6=j
ai,`f`,j(k − 1), k ≥ 2, und a(k)i,j =
δi,j(Kronecker-δ), k = 1,k∑`=1
fi,j(`)a(k−l)j,j , k ≥ 2,
Aus der ersten Formel folgt fur k ≥ 2
fi,j(k) =∑
ai,`1a`1`2 · · · a`k−1,j ,
wobei uber alle (k − 1)-Tupel (`1, `2, . . . , `k−1) summiert wird mit `s 6= j fur
alle s = 1, 2, . . . k − 1.
Weiter ist f∗i,j =∑∞k=1 fi,j(k) die Wahrscheinlichkeit, dass der Prozess
ausgehend vom Zustand i in endlich vielen Schritten nach j kommt.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 102
Der Zustand i heißt rekurrent, wenn f∗ii = 1, d.h. wenn mit Wahrscheinlichkeit
1 der Prozess ausgehend von i irgendwann wieder nach j zuruckkehrt.
Der Zustand i heisst transient, wenn f∗ii < 1, d.h. es gibt eine positive
Wahrscheinlichkeit (1− f∗i,i > 0) dafur, dass der Prozess nie mehr nach i
zuruckkehrt.
Ein rekurrenter Zustand i heißt periodisch mit Periode p, falls
p = ggTm : a(m)i,i > 0. Ist p = 1, so heißt der Zustand aperiodisch.
Sei Ni die Anzahl, wie oft der Prozess in den Zustand i zuruckkehrt, wenn er
dort gestartet wird (Startpunkt wird mitgezahlt.) Die Zufallsvariable Ni is
geometrisch verteilt:
W (Ni = m) =(f∗i,i)m−1
(1− f∗i,i) (m = 1, 2, . . .).
Das bedeutet:
W (Ni <∞) =∞∑m=1
W (Ni = m) =∞∑m=1
(f∗i,i)m−1
(1− f∗i,i) =
0, f∗i,i = 1,
1, f∗i,i = 0.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 103
Ein transienter Zustand i ist also dadurch charakterisiert, dass der Prozess mit
Wahrscheinlichkeit 1 nur endlich oft nach i zuruckkehrt, wahrend ein
rekurrenter Zustand mit Wahrscheinlichkeit 1 unendlich oft besucht wird.
Satz 3.2 (Klasseneigenschaften bei (endlichen) Markoff-Ketten)
1. Nicht alle Zustande konnen transient sein.
2. Angenommen, die Zustande i und j kommunizieren.
Dann gelten die Solidaritatsprinzipien:
a) Ist i rekurrent, so auch j.
b) Ist i transient, so auch j.
c) Ist i aperiodisch, so auch j.
d) Ist i periodisch mit Periode p ≥ 2, so auch j.
3. Jede (endliche) Markoff-Kette besitzt mindestens eine Kommunikations-
klasse, in der alle Zustande rekurrent sind.
4. Ist die endliche Markoff-Kette irreduzibel, so sind alle Zustande rekurrent.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 104
Wir beschreiben die Struktur der Ubergangsmatrix A:
Der Zustandsraum von (S ,A ) zerfalle in r1 rekurrente und r2 transiente
Kommunikationsklassen.
Dann kann A nach geeigneter Umnummerierung der Zustande in der folgenden
Form dargestellt werden:
A =
A1,1 0 0 0
. . .
0 Ar1,r1 0 0
Ar1+1,1 · · · Ar1+1,r1 Ar1+1,r1+1 0...
. . ....
.... . .
Ar1+r2,1 · · · Ar1+r2,r1 Ar1+r2,r1+1 · · · Ar1+r2,r1+r2
=
[AR,R O
AT,R AT,T
].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 105
Die quadratischen Matrizen Ak,k (1 ≤ k ≤ r1) sind irreduzible stochastische
Matrizen. Insbesondere besitzen sie alle den Spektralradius 1.
Die quadratischen Matrizen Ak,k (r1 + 1 ≤ k ≤ r1 + r2) sind entweder
Nullmatrizen der Dimension 1× 1 oder sie sind irreduzibel und
substochastisch, d.h. ihre Zeilensummen sind ≤ 1, wobei fur mindestens eine
Zeile strikte Ungleichheit gilt. Das bedeutet, dass die Spektralradien ρ(Ak,k)
echt kleiner als 1 sind und daher limm→∞Amk,k = 0 gilt
(r1 + 1 ≤ k ≤ r1 + r2). Daraus folgt: limm→∞AmT,T = 0.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 106
3.5 Stationare Verteilungen
Sei (S ,A ) eine diskrete Markoff-Kette mit Zustandsraum 1, 2, . . . , n und
Ubergangsmatrix A.
Ein Wahrscheinlichkeitsvektor π heißt stationar (invariant) bez. (S ,A ), wenn
π = πA
gilt.
Aquivalent: π ist linker Eigenvektor von A zum Eigenwert 1.
Aquivalent: πT lost das homogene LGS (In −AT )x = 0 .
Aquivalent: πT liegt in N (In −AT ), dem Nullraum von (In −AT ).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 107
• Existiert π = limm→∞ π0Am (fur eine beliebige Anfangsverteilung π0),
dann ist π stationar.
• Existiert limm→∞Am = A∞, dann ist π0A∞ stationar fur jede Verteilung
π0.
• Gilt limm→∞ a(m)i,j = πj (fur alle i und j), m.a.W.: konvergiert die j-te
Spalte von Am gegen πje (m→∞) fur alle j, dann ist
π = [π1, π2, . . . , πn] stationar.
Satz 3.3 (Satz von Perron (1908) und Frobenius (1912))
Ist A ∈ Rn×n nichtnegativ und irreduzibel, so gelten:
1. Der Spektralradius ρ(A) = max|λ| : λ ist Eigenwert von A > 0 ist
ein (algebraisch) einfacher Eigenwert von A.
2. Der (bis auf skalare Vielfache eindeutige) Eigenvektor von A zum Eigen-
wert ρ(A) kann positiv gewahlt werden.
3. Gilt Ax = λx fur ein x ≥ 0 , x 6= 0 , so folgt λ = ρ(A).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 108
Satz 3.4 (Satz von Perron-Frobenius, Teil II)
Ist A = [ai,j ] ∈ Rn×n nichtnegativ und irreduzibel, so gelten:
1. Besitzt A genau p verschiedene Eigenwerte λ1, . . . , λp vom Betrag ρ(A), dann folgt
(nach einer geeigneten Umnummerierung)
λj = ρ(A) exp
(2πj
pi
)(j = 1, . . . , p).
λ1, . . . , λp sind einfache Eigenwerte von A. Außerdem ist das Spektrum von A
( = Menge aller Eigenwerte) invariant unter der Drehung um den Winkel 2π/p.
2. Ist p > 1, so gibt es eine Permutationsmatrix P ∈ Rn×n mit
PAP> =
O A1,2
O A2,3
. . .. . .
O Ap−1,p
Ap,1 O
(die Nulldiagonalblocke sind quadratisch).
3. Ist wenigstens ein Diagonalelement von A positiv, so folgt p = 1.
4. Gibt es i, j ∈ 1, . . . , n mit ai,jaj,i > 0, so folgt p ≤ 2.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 109
Sei A ∈ Rn×n nichtnegativ und irreduzibel. p sei die Anzahl der Eigenwerte
von A vom Betrag ρ(A). Ist p = 1, so heißt A primitiv. Ist p > 1, so heißt A
zyklisch vom Index p.
Eine nichtnegative irreduzible Matrix A ∈ Rn×n ist genau dann primitiv, wenn
es eine Potenz m gibt, so dass Am nur positive Eintrage besitzt.
Sei A zyklisch vom Index p. Liegt A in der Normalform von Satz 3.4 vor, so
folgt:
Ap =
B1
B2
. . .
Bp−1
Bp,
mit Bj = Aj,j+1Aj+1,j+2 · · ·Ap,1A1,2A2,3 · · ·Aj−1,j .
Alle Bj , j = 1, 2, . . . , p, sind quadratisch, nichtnegativ, irreduzibel und
primitiv.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 110
Satz 3.5 (Romanovsky, 1936)
Es sei G(A) der gerichtete Graph einer nichtnegativen irreduziblen Matrix
A ∈ Rn×n. Fur jeden Knoten Nj von G(A) sei Sj die Menge aller gerichteten
Wege, die Nj mit sich selbst verbinden. Die Lange `(w) eines gerichteten Wegs
w sei die Anzahl der beteiligten Kanten. Wir setzen
pj = ggT `(w) : w ∈ Sj.
Dann gilt p1 = p2 = · · · = pn =: p und p ist die Anzahl der Eigenwerte von
A vom Betrag ρ(A), d.h. A ist zyklisch vom Index p.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 111
Synopsis
Sei A = [ai,j ] ∈ Rn×n die Ubergangsmatrix der Markoff-Kette (S ,A ).
1. Immer existiert
P := limm→∞
1
m+ 1
m∑k=0
Ak.
P ist stochastisch. Außerdem ist P eine Projektion, die mit A
kommutiert, genauer: P = P 2 = PA = AP . Es gilt
N (P ) =⊕
λ∈Λ(A)\1
N ((A− λI)n).
π0P ist stationar fur jede Verteilung π0. Insbesondere besitzt jede
Markoff-Kette (unter den hier ublichen Voraussetzungen) eine
stationare Verteilung.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 112
2. Wenn
Q := limm→∞
Am
existiert, so gilt Q = P .
3. Ist A primitiv, dann gibt es genau einen Vektor z = [z1, z2, . . . , zn]> ∈ Rnmit z > 0 , ‖z‖1 = 1 und z>A = z>.
Der Grenzwert Q = limm→∞Am existiert und es gilt
Q =
z1 z2 · · · zn
z1 z2 · · · zn...
......
z1 z2 · · · zn
.
Ist (S ,A ) irreduzibel und aperiodisch, dann existiert genau eine
stationare Verteilung π. Jede Komponente von π ist positiv. Die
Kette strebt fur jede Ausgangsverteilung gegen π.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 113
4. Ist A zyklisch vom Index p ≥ 2, so existiert limm→∞Am nicht.
Ist (S ,A ) irreduzibel und periodisch, dann existiert genau eine
stationare Verteilung π, aber die Kette strebt nicht fur jede
Ausgangsverteilung gegen π.
5. Die (eindeutig bestimmte) stationare Verteilung π einer p-periodischen
Kette lasst sich wie folgt bestimmen:
Ap = diag(B1, B2, . . . , Bp) besitzt Blockdiagonalform (vgl. die Aussage
auf Seite 109). Die Diagonalblocke Bj sind stochastisch und primitiv, d.h.
sie besitzen eine eindeutig bestimmte stationare Verteilung π(j), die nach
Punkt 3 bestimmt werden kann. Es gilt
π = 1p
[π(1),π(2), . . . ,π(p)
].
Jede Komponente von π ist positiv.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 114
6. Ist A reduzibel mit Normalform
SAST =
A1,1
A2,1 A2,2
.... . .
Ak,1 Ak,2 · · · Ak,k
,so existiert Q = limm→∞Am genau dann, wenn jeder Diagonalblock Ai,iprimitiv ist, fur den ρ(Ai,i) = 1 gilt.
limm→∞Am existiert genau dann, wenn jede rekurrente
Kommunikationsklasse aperiodisch ist.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 115
7. Andere Formulierung: Die Ubergangsmatrix A einer reduziblen Kette
besitzt die Normalform
A =
A1,1 0 0 0
. . .
0 Ar1,r1 0 0
Ar1+1,1 · · · Ar1+1,r1 Ar1+1,r1+1 0...
. . ....
.... . .
Ar1+r2,1 · · · Ar1+r2,r1 Ar1+r2,r1+1 · · · Ar1+r2,r1+r2
=
[AR,R O
AT,R AT,T
],
wenn die Kette r1 rekurrente und r2 transiente Kommunikationsklassen
besitzt.
...
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 116
Der Grenzwert limm→∞Am existiert genau dann, wenn alle rekurrenten
Kommunikationsklassen aperiodisch sind, d.h. genau dann, wenn alle Ak,k(1 ≤ k ≤ r1) primitiv sind.
Dann gelten: limm→∞Ak,k = Qk (1 ≤ k ≤ r1), d.h.
limm→∞AmR,R = QR,R := diag(Q1, Q2, . . . , Qr1) und
limm→∞
Am =
[QR,R O
(I −AT,T )−1AT,RQR,R O
].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 117
8. Eine reduzible Markoff-Kette besitzt i.a. (unendlich) viele stationare
Verteilungen. Zerfallt sie in r1 rekurrente und r2 transiente
Kommunikationsklassen und sind π(1), π(2), . . . , π(r1) die (eindeutig
bestimmten) stationaren Verteilungen der rekurrenten Klassen, so sind
π(α) =
α1π(1), α2π
(2), . . . , αr1π(r1), 0 ,0 , . . . ,0︸ ︷︷ ︸
r2 Nullvektoren
(αj ≥ 0,
∑r1j=1 αj = 1) genau die stationaren Verteilungen der gesamten
Kette.
...
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 118
Sind alle rekurrenten Kommunikationsklassen aperiodisch, dann gilt fur
die Anfangsverteilung
π0(α) =
α1π(1)0 , α2π
(2)0 , . . . , αr1π
(r1)0 , 0 ,0 , . . . ,0︸ ︷︷ ︸
r2 Nullvektoren
(π
(k)0 ist eine beliebige Anfangsverteilung in der rekurrenten
Kommunikationsklasse Kk und αj ≥ 0,∑r1j=1 αj = 1) die Grenzbeziehung
limm→∞
π0(α)Am = π(α).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 119
3.6 Absorbierende Ketten und Ubergangszeiten
Erinnerung: Tj ist die Zeit, in der die Kette zum ersten Mal in den Zustand j
kommt und fi,j(k) ist die Wahrscheinlichkeit dafur, dass dies zur Zeit k
geschieht, unter der Bedingung, dass der Prozess im Zustand i gestartet wird.
f∗i,j =∑∞k=1 fi,j(k) ist die Wahrscheinlichkeit mit der ein Ubergang von i
nach j in endlicher Zeit stattfindet.
Sei jetzt A eine Teilmenge des Zustandsraums, d.h. A ⊆ 1, 2, . . . , n. Wir
definieren
TA := minn ≥ 0 : Xn ∈ A sowie f∗i,A := W (TA <∞|X0 = i).
Die Zufallsvariable TA (erstmalige Ubergangszeit, hitting time) gibt an, wann
der Prozess zum ersten Mal einen der Zustande aus A erreicht. f∗i,A ist die
Wahrscheinlichkeit dafur, dass der Prozess — ausgehend von i — in endlich
vielen Schritten A erreicht. Ist A eine rekurrente Kommunikationsklasse, so
nennt man f∗i,A eine Absorptionswahrscheinlichkeit (wenn diese Klasse erreicht
ist, kann sie mit Wahrscheinlichkeit 1 nicht mehr verlassen werden).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 120
Besteht eine rekurrente Klasse aus genau einem Zustand j, so nennt man
diesen absorbierend. Aquivalent: Ein Zustand j ist genau dann absorbierend,
wenn aj,j = 1 gilt, d.h. wenn er mit Wahrscheinlichkeit 1 nicht mehr verlassen
werden kann.
Schließlich ist (falls A von i aus erreichbar ist)
ki,A = E(TA |X0 = i) =∞∑k=0
kW (TA = k |X0 = i)
die durschschnittliche Zeit, die der Prozess benotigt, um von i erstmalig nach
A zu kommen (mittlere (erstmalige) Ubergangszeit). Ist A eine rekurrente
Kommunikationsklasse, nennt man die durchschnittliche Zeit bis zur
Absorption in A auch Kettenlange.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 121
Im Folgenden sei A immer eine rekurrente Kommunikationsklasse. Ist i
rekurrent, so gilt offenbar
f∗i,A =
1, i ∈ A,0, i 6∈ A,
und ki,A = 0, wenn i ∈ A.
Bevor wir diese Großen fur transiente i bestimmen, betrachten wir einen
wichtigen Spezialfall:
Sei zunachst (S ,A ) eine absorbierende Kette, d.h.
• sie besitzt s ≥ 1 absorbierende Zustande,
• die n− s nicht-absorbierenden Zustande sind transient,
• alle transienten Zustande kommunizieren miteinander.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 122
Die Ubergangsmatrix hat (nach eventueller Umnummerierung der Zustande)
die Form
A =
[Is O
C B
](In−s −B ist invertierbar) und es gilt
limm→∞
Am =
[Is O
(In−s −B)−1C O
].
Bezeichnet A = 1, 2, . . . , s die Menge der absorbierenden Zustande, so gilt
f∗i,A = 1 fur alle i. (Die Wahrscheinlichkeit einer endgultigen Absorption ist 1.)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 123
Kanonische Fragestellungen bei absorbierenden Markoff-Ketten
Frage 1. i und j seien transiente Zustande. ti,j bezeichne die
durchschnittliche Anzahl, wie oft der Zustand j erreicht wird, wenn in i
gestartet wird (der Ausgangszustand wird mitgezahlt, falls i = j).
Antwort: Setze T = [ti,j ]s+1≤i,j≤n ∈ R(n−s)×(n−s), dann gilt
T = (In−s −B)−1.
Frage 2. Der Prozess werde im transienten Zustand i gestartet. Wie groß ist
die durchschnittliche Kettenlange ki,A, das ist die durchschnittliche Anzahl der
Schritte bis zur Absorption?
Antwort: Setze e = [1, 1, . . . , 1]T ∈ Rn−s, dann gilt
ki,A =∑nj=s+1 ti,j = [Te ]i (i = s+ 1, s+ 2, . . . , n).
Frage 3. Der Prozess werde im transienten Zustand i gestartet. Wie groß ist
die Wahrscheinlichkeit qi,j fur Absorption im Zustand j (j ≤ s)?
Antwort: Setze Q = [qi,j ]s+1≤i≤n,1≤j≤s ∈ R(n−s)×s, dann gilt
Q = (I −B)−1C.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 124
Beispiel 1 (Vererbung bei reiner Inzucht [Pharaonen]) Eine vererbliche
Eigenschaft werde durch zwei Gene (a und b) bestimmt. Jeder Mensch habe
zwei solche Gene, was 6 Zustande bei der Genverteilung eines Paares ergibt:
1: aa × aa 2: aa × ab 3: aa × bb
4: ab × ab 5: ab × bb 6: bb × bb
Beim Erbvorgang bekommt das Kind je ein Gen von beiden Eltern, und zwar
jedes der Gene mit Wahrscheinlichkeit 1/2. Mit den Mendelschen Gesetzen
ergibt sich
aa × aa 7→ aa aa × ab 7→ 12 aa + 1
2 ab aa × bb 7→ ab
ab × ab 7→ 14 aa + 1
2 ab + 14 bb ab × bb 7→ 1
2 ab + 12 bb bb × bb 7→ bb
Die erstgeborene Tochter und der erstgeborene Sohn bilden das nachste
Pharaonenpaar. (Wir nehmen an, dass es immer mindestens eine Tochter und
einen Sohn gibt.)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 125
Wir erhalten als Ubergangsmatrix
A =
1 0 0 0 0 0
1/4 1/2 0 1/4 0 0
0 0 0 1 0 0
1/16 1/4 1/8 1/4 1/4 1/16
0 0 0 1/4 1/2 1/4
0 0 0 0 0 1
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 126
bzw. nach der Permutation (1, 6, 2, 3, 4, 5)
A =
1 0 0 0 0 0
0 1 0 0 0 0
1/4 0 1/2 0 1/4 0
0 0 0 0 1 0
1/16 1/16 1/4 1/8 1/4 1/4
0 1/4 0 0 1/4 1/2
=
[I2 O
C B
].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 127
1. Durchschnittliche Anzahl der Besuche in j, wenn in i gestartet wird
und
2. Kettenlange `i bei Start in i.
i/j aa × ab aa × bb ab × ab ab × bb `i
aa × ab 83
16
43
23
296 = 4.83 . . .
aa × bb 43
43
83
43
203 = 6.66 . . .
ab × ab 43
13
83
43
173 = 5.66 . . .
ab × bb 23
16
43
83
296 = 4.83 . . .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 128
3. Wahrscheinlichkeit fur Absorption im Zustand j, wenn in i gestartet
wird.
i/j aa × aa bb × bb
aa × ab 34
14
aa × bb 12
12
ab × ab 12
12
ab × bb 14
34
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 129
Beispiel 2 (Bankrott eines Spielers)
Zwei Spieler spielen mit einem Geldvorrat von insgesamt n Euro. Pro Spiel
wird um einen Euro gespielt, der den Besitzer wechselt. Dabei gewinnt Spieler
1 mit Wahrscheinlichkeit p > 0, Spieler 2 mit Wahrscheinlichkeit
q = 1− p > 0. Der Zustand j liegt vor, wenn der erste Spieler genau j Euro
besitzt. Das Spiel endet, wenn einer der Spieler bankrott ist.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 130
Die Ubergangsmatrix ist
A =
1 0 0 0 · · · 0 0 0
q 0 p 0 · · · 0 0 0
0 q 0 p 0 0 0
0 0 q 0 0 0 0...
.... . .
......
0 0 0 0 · · · 0 p 0
0 0 0 0 · · · q 0 p
0 0 0 0 · · · 0 0 1
∈ R(n+1)×(n+1).
Die Zustande 0 und n sind absorbierend, alle anderen sind transient und
kommunizieren miteinander.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 131
Es gilt
limm→∞
Am =
1 0 · · · 0 0
s1 0 · · · 0 1− s1
s2 0 0 1− s2
......
......
sn−1 0 · · · 0 1− sn−1
0 0 · · · 0 1
mit
sj = t
n−1∑k=j
(qp
)kund t
[1 + q
n−1∑k=2
(qp
)k]= q.
1. Fall: p = q = 1/2. Dann folgt sj = n−jn .
2. Fall: p 6= q. Setze r = qp . Dann folgt sj = rn−rj
rn−1 .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 132
Interpretation: Sei Spieler 2 die Bank. Sie sorgt dafur, dass r (leicht) großer als 1 ist
(Roulette). Die Bank beginne mit k Euro, Spieler 1 mit n− k Euro. Die
Wahrscheinlichkeit fur den Ruin von Spieler 1 ist also
sn−k =rn − rn−k
rn − 1>rk − 1
rk.
0 0.5 1 1.5 2 2.5 3 3.5 4
x 104
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1(rk−1)/rk
r=1.0001
r=1.001
r=1.01
Bei gegebenem r > 1 kann die Bank also ihren Geldvorrat k so bestimmen, dass die
Wahrscheinlichkeit fur den Ruin des Spielers deutlich hoher als 1/2 ist (und zwar
unabhangig vom Anfangskapital n− k des Spielers!).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 133
Satz 3.6 (Mittlere Ubergangszeiten)
Der Vektor [ki,A]1≤i≤n der mittleren Ubergangszeiten ist die minimale nicht-
negative Losung de linearen Gleichungssystems
ki,A =
0, i ∈ A1 +
∑j 6=Aai,jkj,A, i 6= A.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 134
Satz 3.7 (Mittlere Ubergangszeiten bei rekurrenten Zustanden)
Sei
A =
[AR,R O
AT,R AT,T
]die Ubergangsmatrix einer reduziblen Kette mit s rekurrenten und t transien-
ten Zustanden (d.h. AR,R ∈ Rs×s und AT,T ∈ Rt,t).
Die mittlere Ubergangszeit von einem transienten Zustand i in die Ver-
einigung A aller rekurrenten Klassen ist ki,A =∑s+tj=s+1 ti,j . Dabei ist
T = [ti,j ]s+1≤i,j≤s+t = (I −AR,R)−1.
Die mittlere Ubergangszeit von einem transienten Zustand i in eine feste re-
kurrente Klasse K ist ki,K = ai,K∑s+tj=s+1 ti,j . Dabei ist ai,K =
∑j∈K ai,j .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 135
Beispiel 3 (Irrfahrt auf Graphen)
Unter einem Graph verstehen wir jetzt einen vollstandig zusammenhangenden
ungerichteten Graphen mit n Knoten 1, 2 . . . , n ohne Schleifen, also ohne
Kanten, die einen Knoten mit sich selbst verbinden. Mit deg(i), dem Grad des
Knoten i, bezeichnen wir die Anzahl der Kanten in i. In der Irrfahrt wechseln
wir von einem Knoten i zu einem der Nachbarknoten uber mit
Wahrscheinlichkeit 1/deg(i), also zu jedem der Nachbarknoten mit gleicher
Wahrscheinlichkeit. Der Prozess befindet sich im Zustand i, wenn die Irrfahrt
im Knoten i angelangt ist.
Ein solcher Markoff-Prozess ist irreduzibel, besitzt also eine eindeutig
bestimmte stationare Verteilung, namlich
π = 1d [deg(1), deg(2), . . . , deg(n)] mit d =
n∑i=1
deg(i) = 2 cardKanten.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 136
Die Irrfahrt des Springers aus Abschnitt 3.1 ist beispielsweise eine Irrfahrt auf
einem Graphen. Die Knoten des Graphen sind die Felder des Schachbretts.
Zwei Felder sind durch eine Kante verbunden, wenn der Springer von einem
dieser Felder zum anderen ziehen kann. (Die zugehorige Markoff-Kette ist
periodisch mit der Periode 2.)
Satz 3.8 (Mittlere Ruckkehrzeit)
Sei π = [π1, π2, . . . , πn] die stationare Verteilung einer irreduzibeln Markoff-
Kette. Die mittlere Ruckkehrzeit in den Zustand i, also∑∞k=1 kW (Ti =
k |X0 = i), hat den Wert 1/πi.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 137
2
2
2
2
3
3
3
3
3
3
3
3
4 4
4
4
4 4
4
4
4 4
4
4
4 4
4
4
4
4
4
4
6 6
6
6
6 6
6
6
6 6
6
6
6 6
6
6
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
8
Die mittlere Ruckkehrzeit fur den Springer in eine der vier Ecken des
Schachbretts ist also∑64i=1 deg(i)/2 = 336/2 = 168.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 138
3.7 Geschwindigkeit der Konvergenz zum Gleichgewicht.
A sei die Ubergangsmatrix einer Markoff-Kette.
reduzibel
limk→∞ P k = ?
irreduzibel
primitiv
∃ limk→∞ P k
irreduzibel
zyklisch
6 ∃ limk→∞ P k
Frage: Angenommen A ist primitiv (wird ab jetzt immer vorausgesetzt), d.h.
A∞ = limm→∞Am = eπ existiert (π bezeichnet die eindeutig bestimmte
stationare Verteilung der Kette), wie schnell konvergiert dann Am gegen A∞?
Oder: Wie schnell strebt die Markoff-Kette gegen ihre stationare Verteilung?
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 139
Irrfahrt auf Cn, dem zyklischen Graph mit n Knoten
Zeit 0: W (X0 = 0) = 1.
Zeit 1: W (X1 = n− 1) = 1/3,
W (X1 = 0) = 1/3,
W (X1 = 1) = 1/3.
Zeit 2: W (X2 = n− 2) = 1/9,
W (X2 = n− 1) = 2/9,
W (X2 = 0) = 3/9,
W (X2 = 1) = 2/9,
W (X2 = 2) = 1/9.
Offensichtlich ist die stationare Verteilung π die Gleichverteilung
π = 1n [1, 1, . . . , 1].
Wie schnell strebt πm = π0Am (π0 = [1, 0, . . . , 0]) gegen π?
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 140
0 20 40 60 80 1000
0.01
0.02
0.03
0.04
0.05n=100, m=100
0 20 40 60 80 1000
0.01
0.02
0.03
0.04
0.05n=100, m=500
0 20 40 60 80 1000
0.01
0.02
0.03
0.04
0.05n=100, m=1000
0 20 40 60 80 1000
0.01
0.02
0.03
0.04
0.05n=100, m=2000
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 141
Sei E := A−A∞. Dann (beachte A∞ = A∞A∞ = AA∞ = A∞A)
Am −A∞ = (A−A∞)m
= Em (m = 1, 2, . . .).
Seien π die stationare Verteilung der Kette, π0 eine beliebige
Anfangsverteilung und πm = π0Am die Verteilung nach Schritt m. Dann
(beachte A∞ = eπ)
πm − π = π0 (Am −A∞) = π0Em (m = 1, 2, . . .).
Messe Abstand von Verteilungen durch
‖πm − π‖ := sup |πm(X)− π(X)| : X ⊆ 1, 2, . . . , n = 12‖πm − π‖1
(total variation distance). Beachte außerdem
‖M‖ = max‖σ‖=1
‖σM‖ = ‖M‖∞.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 142
Das bedeutet: Fur jede Anfangsverteilung π0 ist
‖πm − π‖ = ‖π0Em‖ ≤ ‖π0‖‖Em‖∞ = 1
2‖Em‖∞
und es gibt ein π0 (worst case), so dass hier Gleichheit gilt.
Die Abschatzung
‖πm − π‖ ≤ 12‖E
m‖∞ ≤ 12‖E‖
m∞
ist oft zu grob, um brauchbar zu sein. Bei der Irrfahrt auf Cn gilt
beispielsweise ‖E‖∞ = 2− 6n ≈ 2.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 143
Satz 3.9
Es seien A die Ubergangsmatrix einer aperiodischen Kette mit stationarer
Wahrscheinlichkeitsverteilung π, d.h. limm→∞Am = eπ, und E = eπ −A.
Dann folgt Λ(E) = 0 ∪ (Λ(A) \ 1).Setzt man γ(A) := ρ(E) = max|λ| : λ ∈ Λ(A) \ 1, dann gilt
limm→∞
‖πm − π‖1/m ≤ γ(A),
wobei fur fast alle π0 Gleichheit gilt.
Die Markoff-Kette strebt also asymptotisch linear mit Konvergenzfaktor
γ(A) gegen ihre stationare Verteilung.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 144
Sn =
1
. . .
1
1
= Fn
ω0n
ω1n
. . .
ωn−1n
F−1, F =[ω(i−1)(j−1)n
]1≤i,j≤n
Sn ∈ Rn×n heißt zirkulanter Shift. Fn ∈ Cn×n ist die Fourier-Matrix.
Die Eigenwerte von Sn sind die n-ten Einheitswurzeln:
ωjn = exp(i 2π jn) = cos(2π j
n) + i sin(2π j
n) (j = 0, 1, . . . , n− 1)
ω60=1
2 π/6
ω6 ω
62
ω63=−1
ω64 ω
65
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 145
Eine Matrix A ∈ Cn×n heißt zirkulant, wenn sie die Form
A = p(Sn) = α0In + α1Sn + α2S2n + · · ·+ αn−1S
n−1n
=
α0 α1 · · · αn−2 αn−1
αn−1 α0 αn−3 αn−2
.... . .
...
α2 α3 α0 α1
α1 α2 αn−1 α0
besitzt, d.h. wenn sie ein Polynom in Sn ist.
Die Eigenwerte λj von A sind dann
λj = p(ωjn) = α0 +α1ωjn+α2ω
2jn +· · ·+αn−1ω
(n−1)jn (j = 0, 1, . . . , n−1).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 146
Die Ubergangsmatrix A der Irrfahrt auf Cn ist 13 (I + Sn + Sn−1
n ), besitzt also
die Eigenwerte
λj = 13 (1 + ωjn + ω
(n−1)jn ) = 1
3 (1 + 2 Real(ωjn)) = 13 (1 + 2 cos(2π jn )). Das
bedeutet γ(A) = 1− 83π2
n2 + O( 1n4 ) (n→∞).
Fur die Irrfahrt auf Cn kann man zeigen:
12γ(A)m ≤ ‖πm − π‖ ≤ 1
2
n−1∑j=1
| 13 (1 + 2 cos(2π jn ))|2m1/2
≤ 12
√nγ(A)m.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 147
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
10−3
10−2
10−1
100
log
of to
tal v
aria
tion
dist
ance
m
Irrfahrt auf Cn, n=100
untere SchrankeFehlerasymptotikobere Schrankeschärfere obere SchrankeAbstand
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 148
Irrfahrt auf HN , dem Hyperkubus in N Dimensionen
Zeit 0: W (X0 = (0, 0, 0)) = 1.
Zeit 1: W (X1 = (0, 0, 0)) = 1/4,
W (X1 = (1, 0, 0)) = 1/4,
W (X1 = (0, 1, 0)) = 1/4,
W (X1 = (0, 0, 1)) = 1/4.
Zeit 2: W (X2 = (0, 0, 0)) = 2/8,
W (X2 = (1, 0, 0)) = 1/8,
W (X2 = (0, 1, 0)) = 1/8,
W (X2 = (0, 0, 1)) = 1/8,
W (X2 = (1, 1, 0)) = 1/8,
W (X2 = (1, 0, 1)) = 1/8,
W (X2 = (1, 1, 0)) = 1/8.
Die Kette besitzt n = 2N Zustande. Die stationare Verteilung ist die
Gleichverteilung.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 149
Ein”
cutoff“-Phanomen
0 2 4 6 80
0.2
0.4
0.6
0.8
1N=8, n=256
0 50 1000
0.2
0.4
0.6
0.8
1N=64, n=1.844674e+019
0 500 1000 15000
0.2
0.4
0.6
0.8
1N=512, n=1.340781e+154
0 n log(n)/4 n log(n)/20
0.2
0.4
0.6
0.8
1N −−> ∞
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 150
Semilogarithmisch
0 2 4 6 8
10−1
100
N=8, n=256
0 50 100
10−1
100
N=64, n=1.844674e+019
0 500 1000 1500
10−1
100
N=512, n=1.340781e+154
γ = .969... γ = .777...
γ = .996...
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 151
‖π − πm‖/‖π − πm−1‖
0 2 4 6 80.6
0.7
0.8
0.9
1N=8, n=256
0 50 1000.94
0.95
0.96
0.97
0.98
0.99
1N=64, n=1.844674e+019
0 500 1000 1500
0.993
0.994
0.995
0.996
0.997
0.998
0.999
1N=512, n=1.340781e+154
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 152
Die Anzahl der Zustande = dim(An) = 2N = n wachst zu schnell mit N , um die
Irrfahrt auf HN (selbst fur moderate Werte von N) nach dem klassischen
Strickmuster analysierenzu konnen.
Bilde durch Aggregation neue Kette mit nur (N + 1) neue Zustanden:
Zj := Knoten ∈ HN , deren Koordinaten genau j 1-Komponenten enthalten
(card(Zj) =(Nj
)). Liegt ein Knoten in Zj , so bedeutet das: Er hat j Nachbarn in
Zj−1, keine Nachbarn (außer sich selbst) in Zj und N − j Nachbarn in Zj+1.
Offensichtlich sind von Zj aus nur Ubergange nach Zj−1, Zj und Zj+1 moglich. Die
Ubergangsmatrix der aggregierten Kette ist also
BN+1 =
1N+1
NN+1
1N+1
1N+1
N−1N+1
2N+1
1N+1
N−2N+1
. . .. . .
. . .
N−1N+1
1
N+11
N+1
NN+1
1N+1
∈ R(N+1)×(N+1).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 153
Formale Bestimmung von BN+1, zum Beispiel fur N = 3, d.h. n = 8:
1
4
1 3 0 0
1 1 2 0
0 2 1 1
0 0 3 1
︸ ︷︷ ︸
BN+1
=
1 0 0 0 0 0 0 0
0 13
13
13
0 0 0 0
0 0 0 0 13
13
13
0
0 0 0 0 0 0 0 1
︸ ︷︷ ︸
IN+1n
1
4
1 1 1 1 0 0 0 0
1 1 0 0 1 1 0 0
1 0 1 0 1 0 1 0
1 0 0 1 0 1 1 0
0 1 1 0 1 0 0 1
0 1 0 1 0 1 0 1
0 0 1 1 0 0 1 1
0 0 0 0 1 1 1 1
︸ ︷︷ ︸
An
1 0 0 0
0 1 0 0
0 1 0 0
0 1 0 0
0 0 1 0
0 0 1 0
0 0 1 0
0 0 0 1
︸ ︷︷ ︸
InN+1
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 154
Allgemein: BN+1 = IN+1n AnI
nN+1 mit der Klassen-Inzidenzmatrix
InN+1 = [ki,j ] ∈ 0, 1n×(N+1), ki,j =
1 vi ∈ Zj0 sonst
und der Klassen-Wahrscheinlichkeitsmatrix
IN+1n = [`i,j ] ∈ R(N+1)×n, `i,j =
1/card(Zi) vj ∈ Zi0 sonst
.
Es gelten
IN+1n InN+1 = IN+1 und InN+1I
N+1n = P ∈ Rn×n ist blockdiagonal,
genauer P = diag(P1, P2, . . . , PN+1) mit
Pj =1(Nj
)eeT ∈ R(Nj )×(Nj ).Insbesondere ist P eine Projektion (P 2 = P ), die mit An kommutiert (AnP = PAn)
und PInN+1 = InN+1 sowie IN+1n P = IN+1
n erfullt.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 155
Weiter gelten:
P ist eine Orthogonalprojektion (P = PT ) und R(P ) besteht aus genau den
Vektoren, die in jeder Klasse Zj konstant sind (also im Fall N = 3 die Form
[a, b, b, b, c, c, c, d] haben).
BmN+1 =(IN+1n AnI
nN+1
)m= IN+1
n Amn InN+1,
π(B) = 2−N[(N0
),(N1
), . . . ,
(NN
)].
Ist schließlich π(A)0 ∈ Rn eine Ausgangsverteilung, die auf jeder Klasse Zj
konstant ist (d.h. fur die π(A)0 = π
(A)0 P gilt), und ist
π(B)0 = π
(A)0 InN+1 ∈ RN+1 (beachte, dass π
(B)0 ein Verteilungsvektor ist), so
gilt
‖π(A) − π(A)0 Am‖ = ‖π(B) − π
(B)0 Bm‖
(m = 0, 1, 2, . . .).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 156
Beispiele
Irrfahrt auf einer Geraden mit absorbierenden Randern.
Hier (aber nicht nur hier) treten tridiagonale
Ubergangswahrscheinlichkeitsmatrizen der Form
A =
1
p1 q1 r1
p2 q2 r2
. . .. . .
. . .
pn−1 qn−1 rn−1
1
∈ R(n+1)×(n+1)
auf.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 157
Nach der ublichen Umnummerierung (absorbierende Zustande zuerst) ist
A =
[I2 0
C B
]mit
C =
p1 0
0 0...
...
0 0
0 rn−1
und B =
q1 r1
p2 q2 r2
p3. . .
. . .
. . . qn−2 rn−2
pn−1 qn−1
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 158
A∞ = limm→∞Am existiert beispielsweise, wenn pi > 0, ri > 0 (i = 1, 2, . . . , n− 1)
gilt. Dann ist
A∞ =
[I2 0
S 0
]mit S = (I −B)−1C.
Da A∞ stochastisch ist, besitzt S die Form
S =
s1 1− s1s2 1− s2...
...
sn−1 1− sn−1
.
Ein Koeffizientenvergleich in S = (I −B)−1C bzw. in S −BS = C liefert
s1 = p1 + q1s1 + r2s2,
sj = pjsj−1 + qjsj + rjsj+1 (j = 2, 3, . . . , n− 2),
sn−1 = pn−1sn−2 + qn−1sn−1.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 159
Setzt man jetzt sj =∑n−1k=j tk, d.h. tn−1 = sn−1 und tj = sj − sj+1
(j = 1, 2, . . . , n− 2), so folgt (mit pj + qj + rj = 1)
(1− q1)t1 + p1
∑n−1k=2 tk = p1,
rjtj = pjtj−1 (j = 2, 3, . . . , n− 1),
was fur j = 2, 3, . . . , n− 1 wiederum
tj =pjrjtj−1 =
pjpj−1
rjrj−1tj−2 = · · · = pjpj−1 · · · p2
rjrj−1 · · · r2t1
impliziert. Den Wert von t1 bestimmt man aus
t1
(1− q1 + p1
∑n−1k=2
p2···pkr2···rk
)= p1.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 160
Zwei Spieler spielen mit n ≥ 2 Karten. Im Elementarprozess wird eine dieser
Karten (gleichverteilt) zufallig ausgewahlt, die dann ihren Besitzer wechselt.
Der Zustand j liegt vor, wenn der erste Spieler genau j Karten auf der Hand
hat (der zweite Spieler hat dann also n− j Karten). Das Spiel endet, wenn
einer der Spieler keine Karten mehr besitzt.
Die Ubergangsmatrix ist
A =
1
1/n 0 (n− 1)/n
2/n 0 (n− 2)/n
. . .. . .
. . .
(n− 1)/n 0 1/n
1
∈ R(n+1)×(n+1).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 161
sj = t1∑n−1k=j
2·3···k(n−2)(n−3)···(n−k) mit t1(1 + 1
n
∑n−1k=2
2·3···k(n−2)(n−3)···(n−k) ) = 1
n
ist die Wahrscheinlichkeit dafur, dass der erste Spieler alle seine Karten verliert
unter der Bedingung, dass er das Spiel mit j ≥ 1 Karten beginnt.
Fur große n ist sj , nahezu unabhanig von j, etwa gleich 1/2:
1
2≤ sj ≤
n+ 4
2n+ 2=
1
2+
3
2n+ 2(1 ≤ j < n
2)
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 162
1 2 3 4 5 6 7 8 90.4
0.45
0.5
0.55
0.6
n=10
2 4 6 8 10 12 14 16 180.46
0.48
0.5
0.52
0.54
0.56n=20
5 10 15 20 250.48
0.5
0.52
n=30
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 163
Warteschlangen. Ein Sessellift kann pro Zeiteinheit eine Person
abtransportieren. Die Wahrscheinlichkeit, dass in einer Zeiteinheit j Personen
in der Talstation ankommen, sei pj (∑∞j=0 pj = 1). Der Wartesaal fasse
maximal n Personen. Personen, die einen vollen Wartesaal antreffen, kehren
um. Der Zustand j liegt vor, wenn sich am Ende einer Zeiteinheit, d.h. in dem
Moment, in dem der Lift abfahrt, genau j Personen im Wartesaal aufhalten.
Die Ubergangsmatrix hat Hessenberg-Form
A =
p0 + p1 p2 · · · pn−1 pn∑∞j=n+1 pj
p0 p1 · · · pn−2 pn−1
∑∞j=n pj
p0 · · · pn−3 pn−2
∑∞j=n−1 pj
. . .. . .
...
p0 p1
∑∞j=2 pj
p0
∑∞j=1 pj
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 164
Um die Rechnung einfach zu halten, nehmen wir an, dass pj = qpj
(j = 0, 1, . . .) gilt fur ein p mit 0 < p < 1. (Aus der Normierung∑∞j=0 pj = 1
ergibt sich q = 1− p.) Es folgt∑∞j=k pj = pk, so dass A jetzt eine einfachere
Form besitzt:
A =
q + qp qp2 · · · qpn−1 qpn pn+1
q qp · · · qpn−2 qpn−1 pn
q · · · qpn−3 qpn−2 pn−1
. . .. . .
...
q qp p2
q p
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 165
A ist primitiv, d.h. es gibt eine eindeutig bestimmte stationare Verteilung
π = [π0, π1, . . . , πn]. Ein Komponentenvergleich in π = πA ergibt
π0 = (q + qp)π0 + qπ1,
πj = qpj+1π0 + qpjπ1 + · · ·+ qpπj + qπj+1 (j = 1, 2, . . . , n− 1),
πn = pn+1π0 + pnπ1 + · · ·+ pπn−1 + pπn.
Durch vollstandige Induktion erhalt man
πj = pj+1
qj π0 (j = 1, 2, . . . , n) und π0 =[1 + p2
q + p3
q2 + · · ·+ pn+1
qn
]−1
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 166
0 2 4 6 8 10 12 14 16 18 200.04
0.06
0.08
0.1n=20, p=0.5
0 2 4 6 8 10 12 14 16 18 200.02
0.04
0.06
0.08n=20, p=0.505
0 2 4 6 8 10 12 14 16 18 200.02
0.04
0.06
0.08n=20, p=0.50845395...
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 167
0 2 4 6 8 10 12 14 16 18 200.02
0.04
0.06
0.08n=20, p=0.51
0 2 4 6 8 10 12 14 16 18 200
0.05
0.1
0.15
0.2n=20, p=0.55
pn+1
(1− p)n≤ 1 ⇒ leerer Wartesaal ist am wahrscheinlichsten
pn+1
(1− p)n≥ 1 ⇒ voller Wartesaal ist am wahrscheinlichsten
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 168
Lagerhaltung. Ein Lager fasst n Maschinen. Im Laufe einer Woche werden
mit Wahrscheinlichkeit pj genau j Maschinen abgerufen (j = 0, 1, . . . , n− 1).
Mit Wahrscheinlichkeit pn werden n oder mehr Maschinen nachgefragt
(∑nj=0 pj = 1). Sei m eine feste Zahl zwischen 1 und n. Im Lager wird
folgende Strategie verfolgt. Sind zu Beginn einer Woche m+ 1 oder mehr
Maschinen vorhanden, so werden die im Laufe der Woche eingehenden
Auftrage (soweit es moglich ist) erledigt. Sind am Anfang der Woche aber nur
m oder weniger Maschinen verfugbar, so wird das Lager zunachst auf n
Maschinen aufgefullt und dann die Abrufe (soweit es moglich ist) bedient.
Der Zustand j liege vor, wenn am Ende einer Woche genau j Maschinen
(j = 0, 1, 2, . . . , n) im Lager sind.
Die Ubergangsmatrix hat die Blockstruktur
A =
[A1,1 A1,2
A2,1 A2,2
]∈ R(n+1)×(n+1)
mit
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 169
A1,1 =
pn pn−1 · · · pn−m
pn pn−1 · · · pn−m...
......
pn pn−1 · · · pn−m
=
1
1...
1
[pn · · · pn−m
]∈ R(m+1)×(m+1),
A1,2 =
1
1...
1
[pn−m−1 pn−m−2 · · · p0
]∈ R(m+1)×(n−m),
A2,1 =
Pm+1 pm · · · p1
Pm+2 pm+1 · · · p2...
......
Pn pn−1 · · · pn−m
∈ R(n−m)×(m+1),
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 170
A2,2 =
p0 0 · · · 0
p1 p0 0
. . .. . .
pn−m−1 pn−m−2 · · · p0
∈ R(n−m)×(n−m).
Dabei ist Pj :=∑nk=j pk (j = m+ 1,m+ 2, . . . , n).
Setzt man pj > 0 voraus (j = 1, 2, . . . , n), so ist A irreduzibel und primitiv. Es
existiert eine (eindeutig bestimmte) stationare Verteilung π = [πj ]0≤j≤n, die wir in
zwei Teilvektoren π =[
π(1) π(2)]
(mit π(1) ∈ Rm+1 und π(2) ∈ Rn−m)
zerlegen. Außerdem setzen wir t :=∑m+1j=0 πj = ‖π(1)‖1.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 171
Aus π = πA folgt
π(1) = π(1)A1,1 + π(2)A2,1 = t[pn, . . . , pn−m] + π(2)A2,1,
π(2) = π(1)A1,2 + π(2)A2,2 = t[pn−m−1, . . . , p0] + π(2)A2,2.
Aus der zweiten Gleichung ergibt sich
π(2) = t[pn−m−1, . . . , p0](In−m −A2,2)−1 (beachte ρ(A2,2) = p0 < 1, so dass
In−m −A2,2 invertierbar ist). Einsetzen in die erste Gleichung liefert
π(1) = t[pn, . . . , pn−m] + t[pn−m−1, . . . , p0](In−m −A2,2)−1A2,1. Schließlich
ergibt sich t aus
1− t = π(t)2 e = t[pn−m−1, . . . , p0](In−m −A2,2)−1et =: tγ
(beachte, dass
γ = [pn−m−1, . . . , p0](In−m −A2,2)−1e =∑j=0[pn−m−1, . . . , p0]Aj2,2e eine
positive Zahl ist und dass damit t = 1/(γ + 1) ∈ (0, 1) gilt.)
Fur die Beispielrechnung setzen wir pj = qpj (j = 0, 1, . . . , n).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 172
0 2 4 6 8 10 12 14 16 18 200
0.02
0.04
0.06
0.08n=20, m=5, p=.5
0 2 4 6 8 10 12 14 16 18 200
0.05
0.1n=20, m=10, p=.5
0 2 4 6 8 10 12 14 16 18 200
0.05
0.1
0.15
0.2n=20, m=15, p=.5
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 173
Labyrinthe. Gegeben ist ein Labyrinth, das aus n Kammern besteht. Die Kammer j
besitzt wj > 0 Turen, die in beide Richtungen passierbar sind. Im Labyrinth befindet
sich eine Maus. Der Zustand j liegt vor, wenn sich die Maus in Kammer j aufhalt.
Im Elementarprozess verweilt die Maus mit Wahrscheinlicheit b, 0 ≤ b < 1, in der
Kammer, in der sie gerade ist, oder sie geht (mit gleicher Wahrscheinlichkeit) durch
eine der verfugbaren Turen.
Die Ubergangsmatrix A = [ai,j ] ist
ai,j =
b, i = j,
(1− b)/wi, falls es zwischen i und j eine Tur gibt,
0 sonst.
Ein Labyrinth heißt zusammenhangend, wenn jede Kammer von jeder anderen
Kammer aus erreichbar ist (also genau dann, wenn A irreduzibel ist). A ist dann
entweder primitiv (sicher dann, wenn b > 0) oder zyklisch vom Index 2 (genau dann,
wenn man die Kammern in zwei disjunkte Klassen K1 und K2 zerlegen kann, so dass
von Kammern aus K1 nur direkte Ubergange nach K2 und von Kammern aus K2
nur direkte Ubergange nach K1 moglich sind). Die stationare Verteilung ist
(unabhangig von b) π = [w1 w2 · · · wn] /∑nj=1 wj .
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 174
Mischen von Spielkarten.
Gegeben sind N Karten, der einfach halber seien sie von 1 bis N durchnummeriert.
Als Zustande betrachten wir die n = N ! verschiedenen Reihenfolgen
(Permutationen), in der diese Karten angeordnet werden konnen. Der Zustandsraum
entspricht also der symmetrischen Gruppe SN , d.h. der Menge aller Permutationen
τ1, τ2, . . . , τn von N Objekten (zusammen mit der ublichen Komposition von
Abbildungen).
Ist etwa N = 3, so besteht der Zustandsraum aus den n = 3! = 6 Permutationen
τ1 = [1 2 3], τ2 = [1 3 2], τ3 = [2 1 3], τ4 = [2 3 1], τ5 = [3 1 2], τ6 = [3 2 1].
Sei außerdem p eine Wahrscheinlichkeitsverteilung auf SN , d.h. p : SN → [0, 1] mit∑τ∈SN
p(τ) = 1. (Man kann p = [pi]1≤i≤n, pi = p(τi), als nichtnegativen Vektor
aus Rn mit∑ni=1 pi = 1 auffassen.).
Die Ubergangswahrscheinlichkeit ai,j wird jetzt wie folgt definiert: Es gibt genau eine
Permutation τ ∈ SN , die die Reihenfolge i (reprasentiert durch die Permutation τi)
in die Reihenfolge j (reprasentiert durch die Permutation τj) uberfuhrt, namlich
τ = τjτ−1i . Wir setzen ai,j = p(τ).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 175
Setzt man in unserem Beispiel p = 16[3 1 0 0 1 1], so erhalt man die Ubergangsmatrix
A =1
6
3 1 0 0 1 1
1 3 0 0 1 1
0 1 3 1 1 0
1 0 1 3 0 1
0 1 1 1 3 0
1 0 1 1 0 3
.
A ist offensichtlich (und dies bezieht sich auf den allgemeinen Fall) doppelt
stochastisch, d.h. wenn der Grenzwert limm→∞Am existiert, so ist
limm→∞
Am =1
neeT =
1
n
1 1 · · · 1
1 1 1...
. . ....
1 1 · · · 1
.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 176
Ein Mischverfahren (definiert durch p) heißt fair, wenn limm→∞Am existiert
(und der Mischprozess damit automatisch gegen die Gleichverteilung
konvergiert).
T (p) := τ ∈ SN ; p(τ) > 0 ⊆ SN
heißt Trager von p.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 177
Satz 3.10
Mit den oben eingefuhrten Bezeichungen gilt:
• A ist genau dann irreduzibel, wenn jede Permutation aus SN ein Produkt
von Permutationen aus T (p) ist.
(Man sagt dann, dass T (p) die Gruppe SN erzeugt.)
• Ist A irreduzibel, so ist A entweder primitiv oder zyklisch vom Index 2.
• Ist A zyklisch vom Index 2, so enthalt T (p) nur ungerade Permutationen.
Schließlich sind die folgenden Aussagen einander aquivalent:
• Das Mischverfahren ist fair.
• A ist primitiv.
• T (p) erzeugt SN und enthalt mindestens eine gerade Permutation.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 178
In unserem Beispiel ist
T (p) = [1 2 3], [1 3 2], [3 1 2], [3 2 1].
T (p) erzeugt S3 ([2 1 3] = [3 1 2] [3 2 1], [2 3 1] = [1 3 2] [3 2 1]) und
enthalt gerade Permutationen ([1 2 3] und [3 1 2]). Unsere Mischverfahren ist
daher fair.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 179
Fur ein (etwas realistischeres) Beispiel sei T die Menge aller Permutationen, die wie
folgt entstehen. Wir zerlegen die N Karten 1, 2, . . . , N in zwei oder drei Teilpakete
(stets sei N ≥ 3):
1, . . . , k1︸ ︷︷ ︸1. Teilpaket
, k1 + 1, . . . , N︸ ︷︷ ︸2. Teilpaket
,
oder
1, . . . , k1︸ ︷︷ ︸1. Teilpaket
, k1 + 1, . . . , k2︸ ︷︷ ︸2. Teilpaket
, k2 + 1, . . . , N︸ ︷︷ ︸3. Teilpaket
.
Im Fall von drei Teilen legen wir das oberste (erste) Paket zuunterst, darauf das
mittlere und schließlich darauf das letzte:
k2 + 1, . . . , N︸ ︷︷ ︸3. Teilpaket
, k1 + 1, . . . , k2︸ ︷︷ ︸2. Teilpaket
, 1, . . . , k1︸ ︷︷ ︸1. Teilpaket
.
Im Fall von zwei Teilen wird das untere auf das obere Paket gelegt. Mit anderen
Worten, wir erlauben alle Permutationen der Form
[k2 + 1, . . . , N, k1 + 1, . . . , k2, 1, . . . , k1] oder [k1 + 1, . . . , N, 1, . . . , k1]
mit 1 < k1 < k2 < N (und nur diese).
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 180
Diese Menge T erzeugt SN . Sei 2 ≤ j ≤ N − 2:
[1, . . . , j − 1︸ ︷︷ ︸1. TP
, j︸︷︷︸2. TP
, j + 1, . . . , N︸ ︷︷ ︸3. TP
]
↓ τ1
[j + 1, . . . , N︸ ︷︷ ︸3. TP
, j︸︷︷︸2. TP
, 1, . . . , j − 1︸ ︷︷ ︸1. TP
] = [j + 1, . . . , N, j︸ ︷︷ ︸1. TP
, 1, . . . , j − 1︸ ︷︷ ︸2. TP
]
↓ τ2
[1, . . . , j − 1︸ ︷︷ ︸2. TP
, j + 1, . . . , N, j︸ ︷︷ ︸1. TP
] = [1, . . . , j − 1, j + 1︸ ︷︷ ︸1. TP
, j + 2, . . . , N︸ ︷︷ ︸2. TP
, j︸︷︷︸3. TP
]
↓ τ3
[ j︸︷︷︸3. TP
, j + 2, . . . , N︸ ︷︷ ︸2. TP
, 1, . . . , j − 1, j + 1︸ ︷︷ ︸1. TP
] = [j, j + 2, . . . , N︸ ︷︷ ︸1. TP
, 1, . . . , j − 1, j + 1︸ ︷︷ ︸2. TP
]
↓ τ4
[1, . . . , j − 1, j + 1︸ ︷︷ ︸2. TP
, j, j + 2, . . . , N︸ ︷︷ ︸1. TP
] = [1, . . . , j − 1, j + 1, j, j + 2, . . . , N ].
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 181
Wir haben gezeigt, dass die Vertauschung (j, j+1) = τ4 τ3 τ2 τ1 als Produkt von
Permutationen aus T geschrieben werden kann (die Beweise fur die Vertauschungen
(1, 2) und (N − 1, N) sind leichte Modifikationen des vorgestellten Verfahrens). Die
Vertauschungen der Form (j, j + 1), 1 ≤ j ≤ N − 1, erzeugen aber die Gruppe SN .
Die Permutation τ , die dem Abheben der obersten N − 2 Karten entspricht,
τ = [3, 4, . . . , N, 1, 2] ∈ T , ist gerade.
Bemerkungen.
1. Lasst man zur Menge T alle Permutationen zu, die sich bei der Aufteilung der
Karten in irgendeine Anzahl von Teilpaketen (nach dem oben skizzierten
Verfahren) ergeben, so ist dieser Mischprozess erst recht fair.
2. Ist N ≥ 5, so kommt man — um einen fairen Mischprozess zu definieren — mit
den Permutationen aus, die sich aus der Aufteilung der Karten in drei Pakete
ergeben.
3. Ein Verfahren mit T = [2, 1, 3, . . . , N ], [2, 3, . . . , N, 1] ist fair, wenn N
ungerade ist. Ein Verfahren mit T = [2, 1, 3, . . . , N ], [1, 3, . . . , N, 2] ist fair,
wenn N gerade ist.
Markoff-Ketten Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 182
4 Nichtnegative Matrizen
4.1 Der Ergodensatz
Problem: Sei A ∈ Cn×n. Wann existieren die Grenzwerte
limm→∞
Am und limm→∞
1
m
m−1∑j=0
Aj ?
Bezeichnung. (Jordan-Block)
J = J(λ) =
λ 1
λ 1
. . .. . .
λ 1
λ
∈ Cn×n
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 183
Lemma 1 [Konvergenz von Jordan-Blocken].
(a) Sei |λ| < 1. Dann gilt limm→∞mk[J(λ)]m = O fur jedes k ∈ N0 .
(b) Der Grenzwert limm→∞[J(λ)]m existiert genau dann, wenn die beiden
folgenden Bedingungen erfullt sind:
1. |λ| ≤ 1,
2. Fur |λ| = 1 muss λ = 1 und n = 1 folgen, d.h. J(λ) = 1(∈ R1×1).
Lemma 2 [Permanenz der Cesaro-Mittel].
Sei Amm≥0 eine Folge aus Cn×n. Existiert limm→∞Am, dann existiert
auch limm→∞1m
∑m−1j=0 Aj und es gilt
limm→∞
1
m
m−1∑j=0
Aj = limm→∞
Am.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 184
Lemma 3 [Cesaro-Mittel fur Potenzen von Jordan-Blocken].
(a) Ist |λ| < 1, dann gilt limm→∞1m
∑m−1j=0 [J(λ)]j = O.
(b) Sind |λ| = 1 und n = 1 (d.h. J(λ) ist ein Skalar vom Betrag 1), dann gilt
limm→∞
1
m
m−1∑j=0
[J(λ)]j =
0, falls λ 6= 1,
1, falls λ = 1.
(c) Ist |λ| ≥ 1 und existiert limm→∞1m
∑m−1j=0 [J(λ)]j , dann gelten |λ| = 1
und n = 1, d.h. J(λ) ist ein Skalar vom Betrag 1.
Bezeichnungen. Fur A ∈ Cn×n heißen Λ(A) := λ ∈ C : det(A− λI) = 0das Spektrum von A und ρ(A) := max|λ| : λ ∈ Λ(A) der Spektralradius
von A.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 185
Satz 4 [Ergodensatz, Teil I].
Die Matrix A ∈ Cn×n besitze die Jordansche Normalform
T−1AT =
J1(λ1)
J2(λ2)
. . .
Jr(λr)
mit Jk(λk) =
λk 1
. . .. . .
λk 1
λk
∈ Cnk×nk ,r∑
k=1
nk = n.
Fur λ ∈ Λ(A) sei d(A, λ) := maxnk : λ = λk, k = 1, 2, . . . , r die
Dimension des großten Jordan-Blocks zum Eigenwert λ. Dann gelten:
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 186
(a) Ist ρ(A) < 1, dann folgt fur alle k ∈ N0
limm→∞
mkAm = limm→∞
1
m
m−1∑j=0
Aj = O.
(b) Ist ρ(A) > 1, dann divergieren die beiden Folgen Amm≥0 und
1m
∑m−1j=0 Ajm≥1.
(c) Ist ρ(A) = 1, dann existiert limm→∞Am genau dann, wenn λ = 1 der
einzige Eigenwert von A mit |λ| = 1 ist und daruberhinaus d(A, 1) = 1
gilt.
(d) Ist ρ(A) = 1, dann existiert limm→∞1m
∑m−1j=0 Aj genau dann, wenn
d(A, 1) = 1 fur alle λ ∈ Λ(A) mit |λ| = 1 gilt.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 187
Satz 5 [Ergodensatz, Teil II].
Fur die Matrix A ∈ Cn×n existiere P = limm→∞1m
∑m−1j=0 Aj .
Dann gelten:
(a) P 2 = P = PA = AP , d.h. der Grenwert P ist eine Projektion, die mit A
kommutiert.
(b)
R(P ) = N (A− I) = N ((A− I)m) fur alle m = 1, 2, . . . ,
N (P ) =⊕
λ∈Λ(A)\1
N ((A− λI)n) .
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 188
Lemma 6 [Approximation des Spektralradius durch Matrixnormen].
Sei A ∈ Cn×n. Dann gibt es zu jedem ε > 0 eine Matrixnorm ‖ · ‖ = ‖ · ‖ε auf
Cn×n mit
‖A‖ − ε ≤ ρ(A) ≤ ‖A‖.
Es gibt genau dann eine Matrixnorm ‖ · ‖ auf Cn×n mit ‖A‖ = ρ(A), wenn
d(A, λ) = 1 fur alle λ ∈ Λ(A) mit |λ| = ρ(A) erfullt ist. Das ist insbesondere
dann der Fall, wenn A diagonalisierbar ist.
Satz 7 [Ergodensatz, Teil III].
Bezeichne ‖ · ‖ eine Matrixnorm auf Cn×n und sei A ∈ Cn×n mit ‖A‖ ≤ 1.
Dann gelten:
(a) Der Grenzwert P = limm→∞1m
∑m−1j=0 Aj existiert.
(b) Enthalt die Menge ζ ∈ C : |ζ| = 1 \ 1 keine Eigenwerte von A, dann
existiert auch limm→∞Am.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 189
Lemma 8 [Satz von Gerschgorin].
Fur A = [ai,j ] ∈ Cn×n sind die Gerschgorin-Kreisscheiben durch
Di :=
ζ ∈ C : |ζ − ai,i| ≤n∑
j=1j 6=i
|ai,j |
(i = 1, 2, . . . , n)
definiert. Dann gilt
Λ(A) ⊆n⋃i=1
Di.
Sind k der Gerschgorin-Kreisscheiben, etwa D1,D2, . . . ,Dk, disjunkt von den
restlichen n− k, d.h.(∪ki=1Di
)∩(∪ni=k+1Di
)= ∅, so enthalt ∪ki=1Di genau k
Eigenwerte von A (algebraische Vielfachheiten mitzahlen!).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 190
Satz 9 [Ergodensatz fur stochastische Matrizen].
Es sei A = [ai,j ] ∈ Rn×n eine stochastische Matrix. Dann gelten:
(a) 1 ∈ Λ(A) und ρ(A) = 1.
(b) Die Matrix P = limm→∞1m
∑m−1j=0 Aj existiert und ist stochastisch.
(c) Ist λ = 1 der einzige Eigenwert von A mit |λ| = 1, dann existiert auch
limm→∞Am, und es gilt limm→∞Am = P .
(d) Gilt ai,i > 0 fur i = 1, 2, . . . , n, dann ist λ = 1 der einzige Eigenwert von
A mit |λ| = 1.
(e) Ist N (A− I) eindimensional, dann gilt P = eyT . Dabei ist
e = [1, 1, . . . , 1]T ∈ Rn und y ∈ Rn, y ≥ 0 , ist eindeutig bestimmt durch
yTA = yT und ‖y‖1 = 1.
(f) Ist N (A− I) eindimensional und ist A doppelt-stochastisch (d.h. A
und AT sind stochastisch), dann gilt P = 1nee
T .
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 191
Satz 10 [Erganzung zum Satz von Gerschgorin].
Es sei A ∈ Cn×n irreduzibel. Mit D1,D2, . . . ,Dn werden die zugehorigen
Gerschgorin-Kreisscheiben bezeichnet.
Liegt ein Eigenwert λ von A auf dem Rand der Vereinigung aller
Gerschgorin-Kreisscheiben, so liegt λ auf dem Rand jeder
Gerschgorin-Kreisscheibe:
λ ∈ Λ(A), λ ∈ δ
n⋃j=1
Dj
⇒ λ ∈ δDj ∀ j = 1, 2, . . . , n.
Korollar. Die Matrix A = [ai,j ] ∈ Cn×n heißt irreduzibel diagonaldominant,
wenn
(a) A irreduzibel ist und
(b)∑nj=1,j 6=i |ai,j | ≤ |ai,i| fur alle i = 1, 2, . . . , n gilt, wobei fur
mindestens ein i Ungleichheit vorliegt.
Fur eine solche Matrix gilt dann ρ(A) ‖A‖∞.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 192
4.2 Die Perron-Frobenius-Theorie nichtnegativer Matrizen
Bezeichnungen. Fur A = [ai,j ], B ∈ Cm×n seien
A ≥ O :⇔ ai,j ≥ 0 (1 ≤ i ≤ m, 1 ≤ j ≤ n),
A > O :⇔ ai,j > 0 (1 ≤ i ≤ m, 1 ≤ j ≤ n),
A ≥ B :⇔ A−B ≥ O, A > B :⇔ A−B > O.
Außerdem setzen wir |A| = [|ai,j |] ∈ Rm×n (Betrag von A).
Lemma 11. Es seien A,B,C,D ∈ Cm×n. Dann gelten:
(a) |A| = O ⇔ A = O.
(b) |αA| = |α||A| ≥ O ∀ α ∈ C.
(c) |A+B| ≤ |A|+ |B|.(d) A,B ≥ O und α, β ≥ 0 ⇒ αA+ βB ≥ O.
(e) A ≥ B und C ≥ D ⇒ A+ C ≥ B +D.
(f) A ≥ B und B ≥ C ⇒ A ≥ C.
Lemma 12. Es seien A,B,C,D ∈ Cn×n und x ,y ∈ Cn. Dann gelten:
(a) |Ax | ≤ |A||x |.(b) |AB| ≤ |A||B|, insbesondere |Am| ≤ |A|m.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 193
(c) O ≤ A ≤ B, O ≤ C ≤ D ⇒ O ≤ AC ≤ BD,
insbesondere O ≤ Am ≤ Bm.
(d) A > O, x ≥ 0 , x 6= 0 ⇒ Ax > 0 .
(e) A ≥ O, x > 0 , Ax = 0 ⇒ A = O.
(f) |A| ≤ |B| ⇒ ‖A‖F ≤ ‖B‖F .
(g) ‖A‖F = ‖|A|‖F .
Satz 13. Es seien A,B ∈ Cn×n mit |A| ≤ B. Dann gilt
ρ(A) ≤ ρ(|A|) ≤ ρ(B).
Korollar. Ist B eine beliebige Hauptuntermatrix der nichtnegativen Matrix
A = [ai,j ] ∈ Rn×n, so folgt ρ(B) ≤ ρ(A). Insbesondere gilt ai,i ≤ ρ(A) fur
i = 1, 2, . . . , n.
Sind A,B ∈ Rn×n mit O ≤ A < B, so folgt ρ(A) < ρ(B).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 194
Lemma 14. Sei A ∈ Rn×n, A ≥ O. Sind die Zeilensummen von A alle
identisch, dann gilt ‖A‖∞ = ρ(A). Sind die Spaltensummen von A alle
identisch, dann gilt ‖A‖1 = ρ(A).
Satz 15. Sei A = [ai,j ] ∈ Rn×n nichtnegativ. Dann gelten
min1≤i≤n
n∑j=1
ai,j ≤ ρ(A) ≤ max1≤i≤n
n∑j=1
ai,j und
min1≤j≤n
n∑i=1
ai,j ≤ ρ(A) ≤ max1≤j≤n
n∑i=1
ai,j .
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 195
Satz 16 [Quotientensatz von Collatz].
Seien A = [ai,j ] ∈ Rn×n nichtnegativ und x = [x1, . . . , xn]> ∈ Rn positiv.
Dann gelten
min1≤i≤n
∑nj=1 ai,jxj
xi= min
1≤i≤n
[Ax ]ixi
≤ ρ(A) ≤ max1≤i≤n
[Ax ]ixi
= max1≤i≤n
∑nj=1 ai,jxj
xi
und
min1≤j≤n
xj
∑nj=1 ai,j
xi≤ ρ(A) ≤ max
1≤j≤nxj
∑nj=1 ai,j
xi.
Lemma 17. Ist A ∈ Rn×n nichtnegativ und irreduzibel, dann ist (I +A)n−1
positiv.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 196
Lemma 18. Es seien A ∈ Rn×n nichtnegativ und irreduzibel und
M := x ∈ Rn : x 6= 0 , x ≥ 0. Fur x ∈M definieren wir
g(x ) := minxi 6=0
∑nj=1 ai,jxj
xi= minxi 6=0
[Ax ]ixi
.
Dann existiert r = maxx∈M g(x ).
Satz 19. [Satz von Perron (1908) und Frobenius (1912)]
Ist A ∈ Rn×n nichtnegativ und irreduzibel, so gelten:
(a) Die Zahl r = maxx∈M g(x ) (vgl. Lemma 18) ist positiv. Außerdem ist r
ein (algebraisch) einfacher Eigenwert von A und es gilt r = ρ(A).
(b) Gilt r = g(y) fur ein y ≥ 0 , so folgt y > 0 und Ay = ry . Der (bis auf
skalare Vielfache eindeutige) Eigenvektor von A zum Eigenwert r kann
also positiv gewahlt werden.
(c) Gilt Ax = λx fur ein x ≥ 0 , x 6= 0 , so folgt λ = r = ρ(A).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 197
Korollar. Seien A ∈ Rn×n nichtnegativ und irreduzibel und B ∈ Cn×n mit
|B| ≤ A. Gilt ρ(A) = ρ(B), so gibt es einen Winkel ϕ ∈ [0, 2π) und eine
Diagonalmatrix D = diag(δ1, . . . δn) ∈ Cn×n mit |δj | = 1 (j = 1, . . . , n), so
dass
B = eiϕDAD−1
gilt (insbesondere folgt |B| = A).
Satz 20. Ist A = [ai,j ] ∈ Rn×n nichtnegativ und irreduzibel, so gelten:
(a) Besitzt A genau k verschiedenen Eigenwerte λ1, . . . , λk vom Betrag
ρ(A), dann folgt (nach einer geeigneten Umnummerierung)
λj = ρ(A) exp
(2πj
ki
)(j = 1, . . . , k).
λ1, . . . , λk sind einfache Eigenwerte von A. Außerdem ist das Spektrum
von A invariant unter der Drehung um den Winkel 2π/k,
Λ(A) = Λ(A) exp
(2π
ki
).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 198
(b) Ist k > 1, so gibt es eine Permutationsmatrix P ∈ Rn×n mit
PAP> =
O A1,2
O A2.3
. . .. . .
O Ak−1,k
Ak,1 O
(die Diagonalblocke sind quadratisch).
(c) Ist wenigstens ein Diagonalelement von A positiv, so folgt k = 1.
(d) Gibt es i, j ∈ 1, . . . , n mit ai,jaj,i > 0, so folgt k ≤ 2.
Sei A ∈ Rn×n nichtnegativ und irreduzibel. k sei die Anzahl der Eigenwerte
von A vom Betrag ρ(A). Ist k = 1, so heißt A primitiv. Ist k > 1, so heißt A
zyklisch vom Index k.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 199
Satz 21 [Romanovsky, 1936].
Es sei G(A) der gerichtete Graph einer nichtnegativen irreduziblen Matrix
A ∈ Rn×n. Fur jeden Knoten Pj von G(A) sei Sj die Menge aller gerichteten
Wege, die Pj mit sich selbst verbinden. Die Lange `(w) eines gerichteten Wegs
w sei die Anzahl der beteiligten Kanten. Wir setzen
kj = großter gemeinsamer Teiler von `(w) : w ∈ Sj.
Dann gilt k1 = k2 = · · · = kn =: k und k ist die Anzahl der Eigenwerte von A
vom Betrag ρ(A).
Lemma 22. Seien A ∈ Rn×n nichtnegativ und m ∈ N. Der gerichtete Graph
G(Am) von Am enthalt die Kante Pi Pj genau dann, wenn es in G(A) einen
gerichteten Weg der Lange m gibt, der von Pi auf Pj fuhrt.
Lemma 23. Sei A ∈ Rn×n nichtnegativ. A ist genau dann primitiv, wenn es
ein m ∈ N gibt mit Am > O. (Die kleinste Zahl m = m(A) mit Am > O heißt
Primitivitatsindex von A.)
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 200
Bemerkung: m(A) ≤ n2 − 2n+ 2 und diese Abschatzung ist scharf.
Satz 24. Sei A ∈ Rn×n primitiv. Seien z = [z1, . . . , zn]> und
y = [y1, . . . , yn]> positive linke bzw. rechte Eigenvektoren von A zum
Eigenwert r = ρ(A) mit y>z = 1 (Az = rz und y>A = ry>). Dann folgen:
(a) P = limm→∞ r−mAm existiert und es gilt P = zy> =
[z1y1 ··· z1yn
......
zny1 ··· znyn
].
(b) Fur beliebige nichtnegative Vektoren v ,w ∈ Rn \ 0 gelten
limm→∞
[Amv ]j[Amw ]j
=y>v
y>w(j = 1, 2, . . . , n)
und
limm→∞
[Amv ]i[Amv ]j
=yiyj
(i, j = 1, 2, . . . , n).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 201
4.3 Stochastische Komplemente
Gegeben ist eine Markoff-Kette mit n = k1 + k2 + · · ·+ k` Zustanden, die wirzu ` Aggregaten zusammenfassen,
z1, . . . , zk1︸ ︷︷ ︸Z1
, zk1+1, . . . , zk2︸ ︷︷ ︸Z2
, . . . zk`−1+1, . . . , zk`︸ ︷︷ ︸Z`
,
und der irreduziblen Ubergangswahrscheinlichkeitsmatrix A ∈ Rn×n.Als Beispiel betrachten wir im folgenden immer
A =
.85 .0 .149 .0009 .0 .00005 .0 .00005
.1 .65 .249 .0 .0009 .00005 .0 .00005
.1 .8 .0996 .0003 .0 .0 .0001 .0
.0 .0004 .0 .7 .2995 .0 .0001 .0
.0005 .0 .0004 .399 .6 .0001 .0 .0
.0 .00005 .0 .0 .00005 .6 .2499 .15
.00003 .0 .00003 .00004 .0 .1 .8 .0999
.0 .00005 .0 .0 .00005 .1999 .25 .55
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 202
mit den drei Aggregaten Z1 = 1, 2, 3, Z2 = 4, 5, Z3 = 6, 7, 8.
Aggregation fuhrt innerhalb der Aggregate auf ` Ketten mit UWM’n
Sj := A(Kj ,Kj) +A(Kj , Lj)(I −A(Lj , Lj))−1A(Lj ,Kj),
wobei Kj := [kj−1 + 1 : kj ] (k0 := 0) und Lj := [1 : n] \Kj (j = 1, 2 . . . , `)
(stochastisches Komplement von A(Kj ,Kj) in A).
In unserem Beispiel
S1 =
0.8503 0.0004 0.1493
0.1003 0.6504 0.2493
0.1001 0.8002 0.0997
, S2 =
[0.7003 0.2997
0.3995 0.6005
],
S3 =
0.6000 0.2499 0.1500
0.1000 0.8000 0.0999
0.1999 0.2500 0.5500
.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 203
Satz 25.
Sj , j = 1, 2, . . . , `, ist stochastisch und irreduzibel. Sj besitzt daher eine
(eindeutig bestimmte) stationare Wahrscheinlichkeitsverteilung
π(j) =[π
(j)kj−1+1, π
(j)kj−1+2, . . . , π
(j)kj
].
Bemerkung. Ist A primitiv, so ist Sj nicht notwendig primitiv. Ist aber
A(Kj ,Kj) primitiv (z.B. wenn mindestens ein Hauptdiagonalelement von 0
verschieden ist), so ist auch Sj primitiv.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 204
Zwischen den Aggregaten wird eine Markoff-Kette mit UWM
C =
π(1)A(K1,K1)e π(1)A(K1,K2)e . . . π(1)A(K1,K`)e
π(2)A(K2,K1)e π(2)A(K2,K2)e . . . π(2)A(K2,K`)e...
. . ....
π(`)A(K`,K1)e π(`)A(K`,K2)e . . . π(`)A(K`,K`)e
(Koppelungsmatrix) induziert. In unserem Beispiel
C =
9.9911 · 10−1 7.9083 · 10−4 1.0000 · 10−4
6.1433 · 10−4 9.9929 · 10−1 1.0000 · 10−4
5.5556 · 10−5 4.4444 · 10−5 9.9990 · 10−1
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 205
Satz 26 [Koppelungssatz].
C ist stochastisch und irreduzibel. C besitzt daher eine (eindeutig bestimmte)
stationare Wahrscheinlichkeitsverteilung
γ = [γ1, γ2, . . . , γ`] .
Die stationare Wahrscheinlichkeitsverteilung von A ist durch[γ1π
(1), γ2π(2), . . . , γ`π
(`)]
gegeben.
Abweichung von der vollstandigen Reduzibilitat:
δ(A) := 2 max1≤j≤`
‖A(Kj , Lj)‖∞
In unserem Beispiel: δ(A) = 2 max10−3, 10−3, 10−4 = 2 · 10−3.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 206
1 2 3 4 5 6 7 80
0.05
0.1
0.15
0.2
0.25
0.3
0.35
1 2 3 4 5 6 7 80
0.2
0.4
0.6
0.8
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 207
Λ(A)
Λ(S1) Λ(S
2) Λ(S
3)
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 208
Satz 27. Sei
S := diag(S1, S2, . . . , S`) ∈ Rn×n
mit Eigenwerten
λ1 = λ2 = · · · = λ` = 1 > |λ`+1| ≥ |λ`+2| ≥ · · · ≥ |λn|.
Dabei seien alle Si primitiv und S sei diagonalisierbar
(T−1ST = diag(λ1, λ2, . . . , λn)). Seien τ (0) eine beliebige Anfangsverteilung
und τ (m) = τ (0)Am (m = 1, 2, . . . ). Außerdem sei
σ :=
∑j∈K1
τ(0)j
π(1),
∑j∈K2
τ(0)j
π(2), . . . ,
∑j∈K`
τ(0)j
π(`)
= lim
m→∞τ (0)Sm.
Dann gilt fur alle m = 1, 2, . . .
‖τ (m) − σ‖1 ≤ mδ(A) + cond∞(T )|λ`+1|m.
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 209
Dynamik einer fast vollstandig reduziblen Kette (δ(A) ‘klein’)
Kurzfristig: Falls |λ`+1|m 1 und m 1/δ(A):
τ (m) ≈ σ =[β1π
(1), β2π(2), . . . , β`π
(`)]
mit βj :=∑ν∈Kj
τ (0)ν .
Kurzfristiges Stabilitatsintervall (Niveau ε):
I(ε) :=m : ‖τ (m) − σ‖1 < ε
wird geschatzt durch E(ε) := m : f(m) < ε mit
f(x) := δ(A)x+ cond∞(T )|λ`+1|x:m :
log(ε/(2cond∞(T )))
log |λ`+1|< m <
ε
2δ(A)
⊆ E(ε) ⊆ I(ε).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 210
1 2 3 4 5 6 7 80.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
0.65
0.7τ(m)./σ, (m=0,2,...,10)
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 211
Mittelfristig: Sei ε > 0 mit I(ε) 6= ∅. Dann gibt es zu jedem
m ≥ m0 := k : k ∈ I(ε) Konstanten α1(m), α2(m), . . . , α`(m) mit∥∥∥τ (m) −[α1(m)π(1), α2(m)π(2), . . . , α`(m)π(`)
]∥∥∥1< ε.
Dies gilt etwa fur
αj(m) =∑ν∈Kj
τ (m−m0)ν (j = 1, 2, . . . , `).
Langfristig: Falls A primitiv,
limm→∞
αj(m) = γj (j = 1, 2, . . . , `)
(linear mit asymptotischem Konvergenzfaktor max|λ| : λ ∈ Λ(A) \ 1= 0.9998 in unserem Beispiel).
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg
Numerische Simulation mathematischer Modelle 212
1 2 3 4 5 6 7 80.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
n=∞n=10 n=100 n=1000 n=10000
Nichtnegative Matrizen Technische Universitat Bergakademie Freiberg