Numerische Simulation mathematischer Modelle€¦ · Numerische Simulation mathematischer Modelle 5...

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


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.



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



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.



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


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.



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.



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



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



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



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



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



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.



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



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


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



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.



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.



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.



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



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’)



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.



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)‖ <∞.



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



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.



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



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→∞.



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.



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.



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.



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



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.



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.



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



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.



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



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



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.



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.



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)



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)



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.



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?



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)



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.



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



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.



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



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.



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.



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



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



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



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.



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



0 u_1* 1/b_2 10

u_2*

1/b_1

1



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.



(((!!

""

x1 x2 x3

0 < f ′(N∗) < 1



HHaa

aaaH

HHH

QQQ

cc

x1 x2x3

−1 < f ′(N∗) < 0



SS

eee

@@eeeSSTTTAAAA

x1x2 x3

|f ′(N∗)| > 1



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)



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



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



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



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?



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.



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



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



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.



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


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

].



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

].



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.



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



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



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 .



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.



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



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.



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)



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.



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



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



0 1 2 3 4 5 6 7 8 9 10

0

1

2

Warteschlange

t

X(t

)



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



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.



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.



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



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



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



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.



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



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.



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.



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 .



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



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.



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.



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.



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

].



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.



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



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



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.



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.



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.



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.



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



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.



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.



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.

...



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

].



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.

...



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 = π(α).



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



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.



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.



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



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.



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



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



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

].



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



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



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.



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.



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 .



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



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.



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 .



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.



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.



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.



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?



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 π?



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



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‖∞.



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.



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.



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



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



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.



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



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.



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 −−> ∞



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



‖π − π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



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



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

.



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.



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



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.



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

.



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.



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.



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



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)



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



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

.



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

.



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

.



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



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



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



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



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.



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



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



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 .



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(τ).



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

.



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.



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.



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.



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



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



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.



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


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.



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.



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:



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



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



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.



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



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 .



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‖∞.



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.



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



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 .



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.



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



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

).



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



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



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



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



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

.



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.



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



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.



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



Λ(A)

Λ(S1) Λ(S

2) Λ(S

3)



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.



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(ε).



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)



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



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


Numerische Simulation mathematischer Modelle€¦ · Numerische Simulation mathematischer Modelle 5...

Documents

Transcript of Numerische Simulation mathematischer Modelle€¦ · Numerische Simulation mathematischer Modelle 5...