close

Anmelden

Neues Passwort anfordern?

Anmeldung mit OpenID

9. Qualitative Analyse mathematischer Modelle In - comedia.co.at

EinbettenHerunterladen
9. Qualitative Analyse mathematischer Modelle
In diesem Abschnitt zeigen wir, wie sich Information unabh¨angig von den genauen
Zahlenwerten aus einem mathematischen Modell ableiten l¨aßt, wenn es einfach genug f¨
ur
eine analytische Behandlung (mit Bleistift und Papier) ist. Im Allgemeinen sind solche
Modelle, genau wie unsere Beispiele, sehr grobe Verallgemeinerungen, eher Karikaturen
als getreue Abbilder der Wirklichkeit. Gerade darum kann aber die Arbeit an solchen
Modellen zum Verst¨andnis der Prinzipien beitragen, die in einem System wirksam sind.
9.1. Die Gleichgewichte des SIR-Modells.
Wir betrachten das Infektionsmodell (SIR-Modell) aus Beispiel 8.5. Die Gleichungen und
Bedeutung der Parameter sind in Tabelle 8.1 zu finden. Wir erkl¨aren kurz die Bedeutung
der Gleichungen:
Die Mengenbilanzen werden durch vier Vorg¨ange beeinflußt: Zuwanderung, Tod,
Infektion und Heilung. F¨
ur die drei Bev¨olkerungsgruppen ergeben sich daraus die
folgenden Bilanzen:
Zuwachs=Zuwanderung
d
S(t)
dt
d
I(t)
dt
d
R(t)
dt
=
=
=
β
–Tod ± Infektion ± Heilung
−µS(t) −λS(t)I(t)
−µI(t) +λS(t)I(t) −γI(t)
−µR(t)
+γI(t)
Beispiel 9.1. Bestimmen Sie die Gleichgewichtslagen des SIR-Modells.
L¨
osung. Die statische Mengenbilanz erh¨alt man durch Nullsetzen der Ableitung in
der dynamischen Mengenbilanz: Im Gleichgewicht ver¨andern sich die Best¨ande nicht,
daher sind ihre Zuwachsraten gleich Null. Die Zustandsgr¨oßen S, I, R h¨angen dann nicht
mehr von der Zeit ab.
(9.1)
0 = β − µS − λSI,
(9.2)
0 = −µI + λSI − γI,
(9.3)
0 = −µR + γI.
Aus der letzten Gleichung (9.3) ergibt sich R, sobald I bekannt ist:
γ
R = I.
µ
Wir untersuchen jetzt die Gleichung (9.2), indem wir I herausheben:
0 = I(λS − (µ + γ)).
Diese Gleichung wird in zwei F¨allen gel¨ost: Fall 1: I = 0 oder Fall 2: λS − (µ + γ) = 0.
Wir untersuchen die beiden F¨alle separat. Zur Unterscheidung der beiden
Gleichgewichtslagen indizieren wir die Werte von S, I, R durch S1 , I1 , R1 bzw. S2 , I2 , R2 .
Fall 1: I1 = 0 bedeutet, daß in der Population gar keine infizierten Individuen pr¨asent
sind, die Infektion ist ausgestorben, die Bev¨olkerung ist gesund. Wir wissen dann auch:
γ
R1 = I1 = 0.
µ
99
100
Zur Berechnung von S berufen wir uns auf Gleichung (9.1):
0 = β − µS1 − λS1 0
β
S1 = .
µ
Im “gesunden Gleichgewicht” stellen sich also die folgenden Werte ein:
β
S1 = , I1 = 0, R1 = 0.
µ
Fall 2: Aus λS2 − (µ + γ) = 0 folgt sofort
S2 =
µ+γ
.
λ
Einsetzen in (9.1) liefert
µ+γ
µ+γ
−λ
I2
λ
λ
β
µ
I2 =
− .
µ+γ λ
0=β−µ
Daraus folgt sofort
γ
βγ
γ
I2 =
− .
µ
µ(µ + γ) λ
Insgesamt erhalten wir die Zahlenwerte
µ+γ
β
µ
βγ
γ
S2 =
, I2 =
− , R2 =
− .
λ
µ+γ λ
µ(µ + γ) λ
R2 =
Allerdings sind diese Werte nur sinnvoll, wenn sie nichtnegativ sind. Negative Infizierte
oder Geheilte kann es ja nicht geben. Die Gleichgewichtslage I = 0 wurde bereits
untersucht. Daher kann eine zweite Gleichgewichtslage nur auftreten, wenn gilt I2 > 0,
also
β
µ
(9.4)
> .
µ+γ
λ
Es gibt bei diesem Gleichgewicht eine positive Anzahl Infizierter in der Population, die
Krankheit hat sich eingenistet, sie ist endemisch.
Zusammenfassend sagen wir: Es gibt jedenfalls das Gleichgewicht, in dem die Krankheit
gar nicht vorkommt, mit den Werten S1 , I1 , R1 . Wenn die Ungleichung (9.4) gilt, so
existiert zus¨atzlich ein Gleichgewicht, in dem die Krankheit endemisch ist, mit den
Werten S2 , I2 , R2 .
Bevor wir in der Untersuchung fortfahren, machen wir eine Beobachtung, die uns viel
Arbeit ersparen kann: Es gibt keine R¨
uckwirkung von R auf die Zustandsgr¨oßen S und
I. Das heißt, S und I k¨onnen zun¨achst ohne Ber¨
ucksichtigung von R mit Hilfe der ersten
beiden Differentialgleichungen untersucht werden. Wenn man u
¨ber S und I Bescheid
weiß, kann man leicht R mit Hilfe der dritten Gleichung behandeln.
9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE
101
9.2. Stabilit¨
at.
In nat¨
urlichen Systemen k¨onnen wir nicht erwarten, daß eine Gleichgewichtslage
eingenommen und f¨
ur alle Zeit festgehalten wird: Jedes nat¨
urliche System ist gewissen
geringf¨
ugigen Schwankungen und St¨orungen ausgesetzt. Wenn sich ein Gleichgewicht
¨
eingespielt hat, und es tritt anschließend eine kleine, vor¨
ubergehende Anderung
der
Systemumgebung auf, so kann das System auf drei Arten auf die St¨orung reagieren:
ubergehend aus dem Gleichgewicht ausgelenkt, es erholt
• Das System wird vor¨
sich aber und strebt wieder auf das alte Gleichgewicht zu. Das Gleichgewicht
heißt asymptotisch stabil. Gilt dieses Verhalten nur f¨
ur kleine St¨orungen, so
heißt das System lokal asymptotisch stabil. Wird von jeder beliebigen Lage aus
letztlich immer dasselbe Gleichgewicht angestrebt, so heißt das Gleichgewicht
global asymptotisch stabil.
• Das System wird aus dem Gleichgewicht ausgelenkt und strebt nicht wieder
zur¨
uck, aber wenn die St¨orung klein genug war, entfernt sich das System nicht
weit vom alten Gleichgewicht. Das Gleichgewicht heißt stabil.
• Auch bei kleinen St¨orungen bricht das alte Gleichgewicht zusammen, und das
System steuert weit entfernte Zust¨ande an. Das Gleichgewicht ist instabil.
In der Natur sind instabile Gleichgewichte nur sehr selten zu beobachten: Sie brechen ja
bei der kleinsten St¨orung zusammen. Ein instabiles Gleichgewicht erzielt man zum
Beispiel, wenn man einen Bleistift auf seiner Spitze ausbalanciert. Will man selbst ein
System in einer bestimmten Gleichgewichtslage halten, so muß man Bedingungen
schaffen, die bewirken, daß diese Gleichgewichtslage stabil, am besten asymptotisch
stabil ist.
Wir geben jetzt eine mathematisch exakte Definition der Stabilit¨atsbegriffe:
Definition 9.2. Wir betrachten die Differentialgleichung dtd x(t) = f (x(t)). (Dabei kann
x(t) ein Vektor sein, das System kann also mehrere Zustandsgr¨oßen umfassen. f ist eine
gegebene Funktion.) Sei x0 eine Gleichgewichtslage, also f (x0 ) = 0.
• Das Gleichgewicht x0 heißt global asymptotisch stabil, wenn gilt: F¨
ur jede
L¨osung x(t) der Differentialgleichung ist der Grenzwert
lim x(t) = x0 .
t→∞
• Das Gleichgewicht x0 heißt lokal asymptotisch stabil, wenn es einen Radius
δ > 0 gibt, sodaß gilt: Jede L¨osung x(t), die h¨ochstens im Abstand δ von x0
entfernt beginnt, strebt gegen x0 :
Wenn |x(0) − x0 | ≤ δ, dann gilt lim x(t) = x0 .
t→∞
• Das Gleichgewicht x0 heißt stabil, wenn gilt: Zu jeder beliebig vorgegebenen
Toleranz > 0 l¨aßt sich ein Radius δ finden, sodaß sich L¨osungen, die innerhalb
des Abstandes δ von x0 starten, nie weiter als von x0 entfernen:
Zu jedem
> 0 gibt es ein δ > 0 sodaß gilt:
Wenn |x(0) − x0 | ≤ δ, dann gilt f¨
ur alle t : |x(t) − x0 | ≤ .
• Ist das Gleichgewicht nicht stabil, so heißt es instabil.
102
Eine genauere Diskussion des Begriffes Stabilit¨at h¨oren Sie in einer Vorlesung u
¨ber
Differentialgleichungen.
9.3. Stabilit¨
at der Gleichgewichte des SIR-Modells.
Beispiel 9.3. Untersuchen Sie, ob die Gleichgewichtslage S = βµ , I = 0 des SIR-Modells
stabil ist.
Wir verwenden den folgenden Lehrsatz:
Satz 9.4 (von der linearisierten Stabilit¨at). Wir betrachten die Differentialgleichung
d
x(t) = f (x(t)), in Komponenten ausgeschrieben:
dt

 

f1 (x1 (t), · · · , xn (t))
x1 (t)
x (t)   f (x (t), · · · , xn (t)) 
d 
,
 2.  =  2 1
..

dt  ..  
.
fn (x1 (t), · · · , xn (t))
xn (t)
und eine Gleichgewichtslage x0 . Sei J die Jacobimatrix von f an der Stelle x0 , also
 ∂

∂
f
(x
)
·
·
·
f
(x
)
1
0
1
0
∂x1
∂xn


..
..
J =
.
.
.
∂
f (x )
∂x1 n 0
···
∂
f (x )
∂xn n 0
Dann gilt:
(1) Wenn die Realteile aller Eigenwerte von J kleiner als Null sind, dann ist die
Gleichgewichtslage x0 lokal asymptotisch stabil.
(2) Wenn J mindestens einen Eigenwert besitzt, dessen Realteil gr¨oßer als Null ist,
dann ist die Gleichgewichtslage x0 instabil.
(3) (Wenn der gr¨oßte Realteil der Eigenwerte genau Null ist, dann k¨onnen wir mit
diesem Kriterium keine Schl¨
usse ziehen.)
L¨
osung von Beispiel 9.3. In diesem Beispiel lautet die Differentialgleichung (nach
Weglassen von R, wie oben erkl¨art):
d
dt
S(t)
I(t)
=
β − µS(t) − λS(t)I(t)
.
−(γ + µ)I(t) + λS(t)I(t)
Wir bilden die Jacobimatrix, indem wir die beiden Komponenten der rechten Seite erst
nach S und dann nach I partiell differenzieren:
J=
=
∂
(β − µS − λSI)
∂S
∂
(−(γ + µ) + λSI)
∂S
∂
(β
∂I
− µS − λSI)
+ µ)I + λSI)
∂
(−(γ
∂I
−µ − λI
−λS
λI
−(γ + µ) + λS
9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE
Wir setzen dann die Gleichgewichtswerte S1 =
J=
β
µ
103
und I1 = 0 ein:
−µ
− λβ
µ
0 −(γ + µ) +
λβ
µ
.
Zur Berechnung der Eigenwerte bilden wir die Matrix zE − J, wobei E die
Einheitsmatrix ist, und bilden die Determinante. Dieser Ausdruck ist das
charakteristische Polynom:
1 0
det(zE − J) = det z
−
0 1
= det
−µ
− λβ
µ
0 −(γ + µ) +
λβ
z+µ
µ
0
z + (γ + µ) −
= (z + µ) (z + γ + µ −
λβ
µ
λβ
µ
λβ
).
µ
Die Eigenwerte sind die Nullstellen des charakteristischen Polynoms, also die L¨osungen
von:
λβ
(z + µ)(z + γ + µ −
) = 0.
µ
Damit das Produkt Null wird, muß mindestens einer der Faktoren Null sein. Die beiden
L¨osungen sind daher
λβ
z1 = −µ und z2 = −γ − µ +
.
µ
Der Realteil von z1 (das ist z1 selbst) ist auf jeden Fall negativ. Alle Stabilit¨at h¨angt
daher an z2 . Der Eigenwert z2 ist genau dann negativ, wenn gilt
−γ−µ+
λβ
< 0, d.h.
µ
λβ
<γ+µ
µ
λ
γ+µ
<
µ
β
µ
β
>
.
λ
γ+µ
Wir erhalten daher die Antwort: Die “gesunde” Gleichgewichtslage ist lokal asymptotisch
stabil, wenn gilt
β
µ
>
.
λ
γ+µ
Sie ist instabil, wenn z2 > 0 gilt, also
µ
β
<
.
λ
γ+µ
Diese Ungleichung ist genau Ungleichung (9.4). Das heißt aber: Wenn die endemische
Gleichgewichtslage existiert, dann ist die gesunde Gleichgewichtslage instabil.
104
Wir k¨onnten jetzt ebenfalls die endemische Gleichgewichtslage nach dem Prinzip der
linearisierten Stabilit¨at untersuchen (wir m¨
ussen nur S2 und I2 in die Jacobimatrix
einsetzen). Das Ergebnis einer l¨angeren Rechnung w¨are: Wenn die endemische
Gleichgewichtslage existiert, dann ist sie lokal asymptotisch stabil.
Das Prinzip der linearisierten Stabilit¨at ist eine Routinemethode zur Untersuchung
lokaler Stabilit¨at. Das globale Verhalten eines dynamischen Systems ist viel schwieriger
zu untersuchen, und obwohl es verschiedene Methoden gibt, gibt es kein Patentrezept.
Im Fall des SIR-Modells f¨
uhrt eine globale Untersuchung zu folgenden Schl¨
ussen:
Wenn die Ungleichung (9.4) erf¨
ullt ist, dann gibt es zwei Gleichgewichtslagen, die
gesunde und die endemische. Jede L¨osung, die mit einer positiven Anzahl von Infizierten
beginnt, konvergiert gegen das endemische Gleichgewicht. Nur wenn die Infektion in der
Population gar nicht auftritt, steuert die Population auf das gesunde Gleichgewicht zu.
Wenn die Ungleichung (9.4) nicht erf¨
ullt ist, dann gibt es nur das gesunde Gleichgewicht,
und es ist global asymptotisch stabil. Gleichg¨
ultig, von welchem Zustand man ausgeht,
die L¨osung steuert immer auf das gesunde Gleichgewicht zu, die Epidemie stirbt aus.
Betrachten Sie jetzt noch einmal die Simulationen aus Abbildung 8.11. Die linke
Abbildung zeigt eine L¨osung, die auf das endemische Gleichgewicht zustrebt. In der
rechten Abbildung sind die Parameter so gew¨ahlt, daß nur das gesunde Gleichgewicht
¨
existiert, und wir sehen eine L¨osung, bei der die Infektion ausstirbt. Uberzeugen
Sie sich
durch Einsetzen der jeweiligen Modellparameter in Ungleichung (9.4), daß das Verhalten
der L¨osungen mit den Voraussagen der Theorie u
¨bereinstimmt.
9.4. Andere qualitative Fragestellungen.
Definition 9.5. Alle Zustandsgr¨oßen eines Systems fassen wir zum Zustandsvektor
zusammen. Der Raum aller Zustandsvektoren ist der Zustandsraum.
Wenn wir aus dem SIR-Modell wie oben die Zustandsgr¨oße R weglassen, erhalten wir als
S(t)
Zustandsvektor zur Zeit t den Vektor
. Der Zustandsraum f¨
ur dieses SI-Modell ist
I(t)
also zweidimensional, eine Ebene.
Gleichgewichte und Stabilit¨at sind Begriffe, die sich mit geometrischen Vorstellungen gut
verkn¨
upfen lassen: Sie beschreiben, ob L¨osungen im Zustandsraum stets an derselben
Stelle stehen bleiben, ob sie sich auf Gleichgewichte zu bewegen, oder ob sie sich weit
entfernen.
Wir z¨ahlen noch einige andere Begriffe auf, die mit geometrischen Vorstellungen
verbunden sind, und sich im Fall einfacher dynamischer Systeme manchmal mit Bleistift
und Papier untersuchen lassen.
Positivit¨at von L¨osungen: Die L¨osungen des SIR-Modells sind nur sinnvoll, wenn
S(t), I(t) und R(t) nicht negativ sind. Wenn das Modell die Wirklichkeit ad¨aquat
9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE
105
wiedergeben soll, muß es also positive Zustandsgr¨oßen liefern, wenn es mit positiven
Anfangswerten gestartet wird. F¨
ur das SIR-Modell kann man diese Eigenschaft beweisen.
Wenn man dagegen zeigen k¨onnte, daß die L¨osungen nicht positiv bleiben, m¨
ußte man
das Modell als unrealistisch ablehnen.
Invariante Teilmengen: Eine Teilmenge des Zustandsraumes, aus der die L¨osung nicht
entweichen kann, wenn sie darin begonnen hat, heißt eine invariante Teilmenge. Zum
Beispiel haben wir eben festgestellt, daß eine L¨osung des SIR-Modells, die im positiven
Quadranten (S ≥ 0, I ≥ 0) startet, den positiven Quadranten nie mehr verl¨aßt. Der
positive Quadrant ist also eine invariante Teilmenge. Wenn man die invarianten
Teilmengen des Zustandsraumes f¨
ur ein System kennt, kann man sich bereits ein
wesentlich besseres Bild u
¨ber das Verhalten machen.
Zum Beispiel k¨onnten die L¨osungen in einer beschr¨ankten Teilmenge gefangen bleiben,
dann k¨onnte man folgern, daß sie nicht gegen unendlich streben k¨onnen. Wenn die
L¨osungen eines Populationsmodells in einer invarianten Teilmenge gefangen sind, die die
Null nicht enth¨alt, sagt das Modell voraus, daß die Population nicht ausstirbt.
Periodische L¨osungen: Eine L¨osung, die nach endlicher Zeit T zu ihrem Anfangszustand
zur¨
uckkehrt, wiederholt denselben Verlauf alle T Zeiteinheiten (dies gilt f¨
ur autonome
Systeme x(t)
˙
= f (x(t)), deren rechte Seite nicht explizit von t abh¨angt). Man sagt, die
L¨osung ist periodisch mit Periode T . Periodische L¨osungen durchlaufen den
Zustandsraum in Form von geschlossenen Kurven. F¨
ur manche Systeme, z.B.
R¨auber-Beute Modelle kann man mit analytischen Methoden entscheiden, ob sie
periodische L¨osungen besitzen.
Beispiel 9.6. Zeigen Sie, daß das SIR-Modell außer den beiden Gleichgewichtslagen
keine periodischen L¨osungen besitzen kann.
L¨
osung. Nat¨
urlich sind die Gleichgewichtslagen selbst periodische L¨osungen mit
jeder beliebigen Periode T . Wir haben festgestellt, daß jede L¨osung gegen eine der
Gleichgewichtslagen konvergiert. Wenn eine L¨osung also von einem anderen Punkt als
einem Gleichgewichtspunkt ausgeht, kann sie nach einiger Zeit nicht mehr zum
Ausgangspunkt zur¨
uckkehren, da dieser zu weit von den Gleichgewichten entfernt w¨are.
Also kann sie nicht periodisch sein.
Viele Naturph¨anomene zeigen periodisches Verhalten, zum Beispiel mechanische
Schwingungen oder der circadiane Rhythmus in der Physiologie. Wegen der Rotation der
Erde um ihre Achse und um die Sonne werden in Modellen ¨okologischer Systeme die
tages- und jahreszeitlichen Bedingungen periodisch angesetzt. Das Modell ist dann
nicht-autonom, x(t)
˙
= f (t, x(t)), wobei aber f periodisch in t ist. Auch in solchen
Modellen kann es periodische L¨osungen geben.
Attraktoren: Eine Teilmenge des Zustandsraumes, auf die jede L¨osung zul¨auft, heißt ein
Attraktor. Ein global asymptotisch stabiles Gleichgewicht ist ein Attraktor, der nur aus
106
Abbildung 9.1. Periodische L¨osung im Zustandsraum als Attraktor
1.5
1
y
0.5
0
−0.5
−1
−1.5
−1.5
−1
−0.5
0
x
0.5
1
1.5
Durchgezogen: periodische L¨osung. Strichliert: zwei weitere L¨osungen.
Die L¨osungen werden im Gegenuhrzeigersinn durchlaufen.
Kreis: Gleichgewicht.
einem Punkt besteht. Aber auch die Kurve einer periodischen L¨osung k¨onnte ein
Attraktor sein, in diesem Fall spielt sich jede L¨osung im Lauf der Zeit mehr und mehr
auf einen Zyklus ein, der letztlich durch die periodische L¨osung wiedergegeben wird. Der
Attraktor ist dann eine geschlossene Kurve im Zustandsraum. Chaotische Systeme
besitzen oft Attraktoren von viel komplizierterer Gestalt.
Beispiel 9.7. Abbildung 9.1 zeigt eine periodische L¨osung der Gleichung
d
x(t) = (1 − x2 + y 2 )x(t) − 5y(t),
dt
d
y(t) = (1 − x2 + y 2 )y(t) + 5x(t).
dt
Der Zustandsraum ist hier zweidimensional. Die Kurve der periodischen L¨osung bildet
einen Attraktor f¨
ur alle L¨osungen außer der Gleichgewichtslage bei Null. Zwei andere
L¨osungen, die sich dem Attraktor n¨ahern, sind eingezeichnet. Eine einzige L¨osung wird
nicht angezogen: Die Gleichgewichtsl¨osung x = y = 0. (Nach dem Satz von
Poincar´e-Bendixson uml¨auft jede periodische L¨osung in einem zweidimensionalen
Zustandsraum mindestens eine Gleichgewichtsl¨osung.) Abbildung 9.2 zeigt dieselben
L¨osungen, diesmal werden die beiden Zustandsgr¨oßen gegen die Zeit aufgetragen.
Chaos: Auch streng deterministische Systeme k¨onnen ein Verhalten zeigen, das
Voraussagen praktisch unm¨oglich macht, n¨amlich dann, wenn die L¨osung so empfindlich
9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE
107
Abbildung 9.2. L¨osungen aus Abbildung 9.1 in anderer Darstellung
1.5
1
0.5
x
0
−0.5
−1
−1.5
0
0.5
1
1.5
t
2
2.5
3
0.5
1
1.5
t
2
2.5
3
2
y
1
0
−1
−2
0
Durchgezogen: periodische L¨osung. Strichliert: zwei weitere L¨osungen.
Strichpunkt: Gleichgewicht.
¨
von den Anfangsdaten abh¨angt, daß kleinste Anderungen
der Anfangswerte zu total
verschiedenen L¨osungen f¨
uhren k¨onnen. Die L¨osungskurven solcher Systeme zeigen
bizarre, scheinbar v¨ollig unregelm¨aßige Formen.
In der Mathematik spricht man von Chaos, wenn das System drei Eigenschaften zeigt:
(1) Es gibt einen Mindestabstand M > 0, sodaß f¨
ur jedes Paar von zwei
verschiedenen L¨osungen x, y ein Zeitpunkt t existiert, zu dem der Abstand
zwischen x(t) und y(t) gr¨oßer als M ist (auch wenn die Anfangswerte der
L¨osungen sehr nahe aneinander liegen).
(2) Zu jedem Punkt x des Zustandsraumes und zu jedem beliebig klein
vorgegebenen Abstand > 0 gibt es periodische L¨osungen, die n¨aher als an x
vorbeilaufen. (Es gibt ungeheuer viele periodische L¨osungen mit den
unterschiedlichsten Perioden.)
(3) Zu je zwei Punkten x, y im Zustandsraum und jedem beliebig klein vorgegebenen
Abstand > 0 gibt es eine L¨osung, die zuerst n¨aher als an x, und dann n¨aher
als an y vorbeil¨auft. (Beinahe, n¨amlich wenn man von beliebig kleinen
Toleranzen absieht, k¨onnte man von jedem Punkt des Systems aus jeden anderen
erreichen: Daher ist der Verlauf der L¨osung v¨ollig unvorhersagbar, wenn die
Anfangsdaten nicht absolut genau vorgegeben sind: Sie kann u
¨berall hingehen.)
Diskrete chaotische Systeme gibt es schon im eindimensionalen Fall. Kontinuierliche
Systeme k¨onnen erst ab Dimension 3 chaotisch sein.
Beispiel 9.8. Abbildung 9.3 zeigt eine Vorrichtung, bei der ein Pendel mit einer
Eisenkugel u
¨ber vier gleichstarken Magneten h¨angt, die in einem Quadrat angeordnet
108
Abbildung 9.3. Pendel aus Beispiel 9.8
sind. Die Pendelbewegung wird von der Erdanziehungskraft und den Magnetfeldern
beeinflußt. Wir halten das Pendel etwas abseits von der Ruhelage fest und lassen es dann
aus. Abbildung 9.4 zeigt die Pendelbewegung, die sich ergibt, von oben gesehen. Verfolgen
Sie die Bewegung.
Beschreibung der Pendelbewegung: Wir beginnen die Bewegung bei den
Koordinaten x = 0.6, y = 2. Das Pendel beginnt zu schwingen. Abgelenkt von den
Magneten der rechten Seite, schwingt es aber nicht aufs Zentrum zu, sondern zun¨achst
zwischen diesen beiden Magneten. Im Zuge der dritten Schwingung ger¨at es aber nahe
genug in den Einflußbereich der linken Magneten, daß es anschließend den Magneten
links hinten umschwingt. Auch die vierte und f¨
unfte Schwingung sind weit nach links
ausgelenkt.
Abbildung 9.4. Pendelbewegung aus Beispiel 9.8
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−2
−1.5
−1
−0.5
0
0.5
1
1.5
Start der Bewegung an den Koordinaten x = 0.6, y = 2.
Kreise: Magneten. Stern: Zentrum.
2
9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE
109
Wir verstehen leicht, warum dieses System sich chaotisch verh¨alt. Die Schwerkraft und
vier Magneten ziehen in verschiedene Richtungen am Pendel. Wenn sich das Pendel in
Bereichen bewegt, in denen mehrere Magneten ungef¨ahr gleich stark wirken, entscheiden
kleinste Unterschiede der Lage und Bewegung dar¨
uber, welcher Magnet die Oberhand
gewinnt und f¨
ur die n¨achste Schwingung das Pendel am st¨arksten anzieht.
¨
Ubrigens
widerspricht dieses System nicht unserer Behauptung, dass Chaos erst ab
Dimension 3 auftritt. Zu den zwei Freiheitsgraden der Lage kommen ja noch dieselben
zwei Freiheitsgrade der Geschwindigkeit. Erst wenn wir Lage und Geschwindigkeit des
Pendels im Augenblick kennen, k¨onnen wir die weitere Bewegung vorausberechnen. Der
Zustandsvektor, der Lage- und Geschwindigkeitskoordinaten umfaßt, ist
vierdimensional.
110
10. Verteilte Parameter
In unseren bisherigen Modellen haben wir stets die Einfl¨
usse der r¨aumlichen
Gegebenheiten vernachl¨assigt. Zum Beispiel wurde im Punktreaktormodell angenommen,
daß die Neutronen, instabile Isotopen und das spaltbare Material u
¨ber den ganzen
Reaktor gleichm¨aßig verteilt sind. Beim Entwurf eines Reaktors ist aber gerade die
r¨aumliche Anordnung der verschiedenen Elemente ein wesentliches Kriterium, daher gibt
das Punktreaktormodell nur eine grobe Bilanz, aber keine Auskunft u
¨ber die spezifischen
Vor- und Nachteile von Konstruktionsdetails.
Als ein anderes Beispiel k¨onnten wir die Stickstoffbelastung von Grundwasser durch
D¨
ungung heranziehen. Nat¨
urlich k¨onnen wir wieder eine grobe Bilanz aufstellen, indem
wir absch¨atzen, wie groß die ged¨
ungte Fl¨ache ist, und welcher Anteil des eingetragenen
Stickstoffs ins Grundwasser absickert. Wenn aber vorhergesagt werden soll, welche
Brunnen besonders gef¨ahrdet sein werden, m¨
ussen wir auch einkalkulieren, an welchen
Stellen der Boden besonders durchl¨assig f¨
ur Sickerwasser ist, und wohin der Schadstoff
durch die Grundwasserstr¨omungen transportiert wird, die in verschiedenen Schichten
verschieden schnell und in verschiedene Richtungen laufen k¨onnen. Je nachdem, wohin
die Str¨omung den Schadstoff tr¨agt, werden manche Bereiche besonders stark durch
Stickstoff belastet werden.
¨
Wenn wir einen Stahltr¨ager belasten, treten Spannungen auf, die bei Uberlastung
zum
Bruch f¨
uhren. Das Gesamtgewicht, das der Tr¨ager aushalten muß, ist nur ein Faktor in
der Beurteilung, ob eine Belastung zul¨assig ist. Ebenso wichtig ist, wie sich die Last auf
den Tr¨ager verteilt. Auch die Spannungen sind nicht im ganzen Tr¨ager gleich. An den
¨
Stellen mit den st¨arksten Spannungen ist bei Uberlastung
der Bruch zu erwarten.
Wir werden uns also jetzt Modellen zuwenden, in denen die Zustandsgr¨oßen nicht nur
von der Zeit, sondern auch vom Ort abh¨angen. Wir sollten von (r¨aumlich) verteilten
Zust¨anden sprechen, aber es hat sich die Phrase “verteilte Parameter” eingeb¨
urgert.
10.1. Ein Modell fu
¨ r die Schadstoffbelastung eines Flusses.
Beispiel 10.1. An einem Flußlauf befindet sich eine punktf¨ormige Schadstoffquelle (etwa
der Abfluß einer Fabrik oder Kl¨aranlage), die zu bestimmten Zeitpunkten organische
Stoffe in das Wasser einleitet. Durch Diffusion (zuf¨allige Molekularbewegung) und
Dispersion (zuf¨allige Bewegung der Wassertr¨opfchen) verteilt sich der Schadstoff im
Wasser, vor allem aber treibt ihn die Str¨omung flußabw¨arts. Im Lauf der Zeit wird der
Stoff durch Oxidation abgebaut, und dadurch belastet er die Wasserqualit¨at, denn die
Sauerstoffkonzentration im Wasser sinkt. An der Wasseroberfl¨ache kann das Wasser
Sauerstoff aus der Luft aufnehmen und sich dadurch regenerieren. (Siehe
Abbildung 10.1.)
Gesucht ist ein Modell, das in der Zeit nach einer Schadstoffaussch¨
uttung den Verlauf
der Schadstoff- und Sauerstoffkonzentration an verschiedenen Stellen des Flusses
vorhersagt. Insbesondere sollen die Stellen erfaßt werden, an denen die
Sauerstoffverknappung besonders drastisch versp¨
urt wird.
10. VERTEILTE PARAMETER
111
Abbildung 10.1. Schadstoffbelastung eines Flusses
Wir nehmen gleich vorweg, daß der gr¨oßte Engpaß f¨
ur die Sauerstoffversorgung erst ein
St¨
uck stromabw¨arts von der Schadstoffquelle auftritt. Auf dem Weg dahin wird n¨amlich
Schadstoff abgebaut, und erst dadurch Sauerstoff aufgebraucht.
Der Zustand unseres Systems besteht aus den Konzentrationen von Schadstoff und
Sauerstoff, in ihrer r¨aumlichen Verteilung u
uck des Flußlaufs. Diese
¨ber ein St¨
Konzentrationen k¨onnen und werden sich mit der Zeit a¨ndern, wenn nicht die
Schadquelle zu allen Zeiten dieselbe Menge an Schadstoff ausst¨oßt.
Folgende Faktoren bestimmen die Struktur des Systems:
• Quellen: Das Modell dreht sich um die Auswirkungen einer punktf¨ormigen
Quelle, den Fabriksabfluß. Aber auch das von weiter oben ankommende
Flußwasser kann mit Schadstoff vorbelastet sein. Entlang der Ufer kann, zum
Beispiel als Folge landwirtschaftlicher Aktivit¨aten, weiteres organisches Material
eindringen, wodurch sich zus¨atzlich eine verteilte Schadstoffquelle ergibt.
• Advektion: Das ist Transport durch die Str¨omung. Der treibende Faktor ist also
die Wasserstr¨omung, und die Richtung ist vorgegeben, n¨amlich flußabw¨arts.
• Diffusion und Dispersion: Das sind zuf¨allige Bewegungen, die f¨
ur sich allein
dazu f¨
uhren w¨
urden, das sich das Material gleichf¨ormig u
ber
das
Wasser verteilt.
¨
Die Bewegung des einzelnen Molek¨
uls oder Tr¨opfchens ist rein zuf¨allig, aber in
der Gesamtschau bewirken diese Vorg¨ange, daß die Mehrheit der Teilchen von
Bereichen großer Konzentration zu Bereichen kleiner Konzentration abwandern.
• Chemische und physikalisch-chemische Reaktionen: An jedem Ort im Wasser
findet die Oxidation des Schadstoffes statt. Treibende Kraft sind die
Konzentration an Schadstoff und Sauerstoff. Die Reaktion bewirkt, daß letztlich
aller Schadstoff abgebaut und dabei Sauerstoff aufgebraucht wird. An der
Wasseroberfl¨ache findet ein Austausch zwischen dem gel¨osten Sauerstoff im
Wasser und dem Luftsauerstoff statt. Treibende Kraft ist der Partialdruck des
Luftsauerstoffs im Vergleich zur Konzentration des gel¨osten Sauerstoffs. Je
weniger Sauerstoff gel¨ost ist, desto mehr wird durch die Wasseroberfl¨ache
aufgenommen.
112
Abbildung 10.2. Einteilung in Zellen
• Nicht ber¨
ucksichtigt in unserer Systembeschreibung sind zum Beispiel die
Wasserzufuhr von Seitenarmen (was zur Verd¨
unnung der
Schadstoffkonzentration f¨
uhren w¨
urde), Versickerung und Verdunstung von
Wasser, zeitliche Schwankungen des Wasserpegels, Ablagerung des Schadstoffes
am Boden und Ufer des Flusses, und vieles mehr.
10.2. Modellierung durch Zellen.
Die Idee dieses einfachen Modellansatzes besteht darin, den Flußlauf (willk¨
urlich) in
einzelne Zellen aufzuteilen, und f¨
ur jede Zelle eine Schadstoff- und Sauerstoffbilanz
aufzustellen. Wir werden dieses Modell unter vereinfachten Annahmen erstellen, wir
nehmen n¨amlich an, daß der Schadstoff u
¨ber die Tiefe und Breite des Flusses gleichm¨aßig
verteilt ist, sodaß Konzentrationsunterschiede nur in der Richtung stromauf–stromab
auftreten. Außerdem vernachl¨assigen wir die Effekte der Dispersion und Diffusion.
Wir teilen also den Flußlauf der L¨ange nach in Zellen, die wir stromabw¨arts mit den
Zahlen 1, 2, 3, · · · , N durchnumerieren (siehe Abbildung 10.2). Wir machen die
vereinfachende (aber falsche) Annahme, daß in jeder Zelle Schadstoff und Sauerstoff
gleichm¨aßig verteilt sind. Je kleiner die Zellen, desto eher entspricht diese Annahme der
Wirklichkeit, aber desto umfangreicher wird das Modell.
Die Zellen mit Nummern 1, 2, · · · , N haben also zur Zeit t ihre Schadstoffkonzentrationen
L1 (t), L2 (t), · · · , LN (t) und ihre Sauerstoffkonzentrationen C1 (t), C2 (t), · · · , CN (t). Das
sind die gesuchten Zustandsgr¨oßen. Es bew¨ahrt sich, die Schadstoffmasse nicht in Mol
Schadstoff, sondern indirekt in Sauerstoff¨aquivalenten anzugeben, also durch die Menge
von Sauerstoff, die ben¨otigt wird, um den Schadstoff v¨ollig zu oxidieren.
Als Parameter hat jede Zelle ihr Wasservolumen V1 , · · · , VN und ihre Wasseroberfl¨ache
(Phasengrenzfl¨ache zur Luft) a1 , · · · , aN . In jede Zelle sickert pro Zeiteinheit aus
verteilten Quellen die Schadstoffmenge S1 (t), · · · , SN (t) ein. S1 (t) beinhaltet auch den
10. VERTEILTE PARAMETER
113
Schadstoff, der aus dem Fabriksabfluß in den Fluß in Zelle 1 eingeleitet wird. Die
gesamte Wassermenge, die den Flußquerschnitt in einer Zeiteinheit durchsetzt, sei F ,
und weil wir keine Seitenarme, Versickerung und Verdunstung und auch keine zeitlichen
Pegelschwankungen ber¨
ucksichtigen, muß diese Zahl das gesamte betrachtete Flußst¨
uck
entlang konstant sein: Dieselbe Menge, die ganz oben zufließt, muß ganz unten wieder
abfließen. Die Konzentrationen an Schadstoff und Sauerstoff im Wasser, das auf Zelle 1
zufließt, bezeichen wir mit L0 (t) und C0 (t). Weitere Parameter werden eingef¨
uhrt
werden, wenn wir die einzelnen Prozesse modellieren.
Wir betrachten jetzt eine einzelne Zelle, die Zelle mit Nummer i, und stellen eine
Mengenbilanz auf. Zum Zeitpunkt t befindet sich in der Zelle die Menge Vi Li (t) an
Schadstoff und Vi Ci (t) an Sauerstoff. (Das gilt deshalb, weil in Vi m3 Wasser Sauerstoff
mit einer Konzentration von Ci (t) kmol
gel¨ost ist.)
m3
Modellierung der Advektion: Oberhalb der Zelle i befindet sich die Zelle i − 1 mit einer
Schadstoffkonzentration von Li−1 (t). Pro Sekunde fließen F Liter des Wassers aus der
oberen Zelle in die Zelle i. Das liefert einen Schadstoffeintrag von F Li−1 (t) und einen
Sauerstoffeintrag von F Ci−1 (t) kmol pro Sekunde in die Zelle i. Vom Wasser der Zelle i
(mit Schadstoffkonzentration Li (t)) fließen pro Zeiteinheit F m3 ab, und nehmen
Schadstoff mit sich. Das ergibt einen Schadstoffabgang von F Li (t) kmol pro Sekunde.
Ebenso errechnet sich der Sauerstoffabgang F Ci (t). Insgesamt ergibt sich als Bilanz:
Schadstoffzuwachs durch Advektion: F Li−1 (t) − F Li (t),
Sauerstoffzuwachs durch Advektion: F Ci−1 (t) − F Ci (t).
Beachten Sie, daß die Namen L0 (t) und C0 (t) f¨
ur die Konzentrationen des auf Zelle (1)
zufließenden Wassers so gew¨ahlt wurden, daß diese Bilanzgleichungen auch in der ersten
Zelle gelten.
Modellierung der Oxidation: Die treibenden Faktoren der chemischen Reaktion sind die
Konzentrationen von Schadstoff und Sauerstoff. (Wir vernachl¨assigen die Abh¨angigkeit
von der Konzentration verschiedener Katalysatoren, oder vom Bestand an
Mikroorganismen bei Oxidation infolge von biologischen Prozessen.) Die Reaktionsrate
einer Reaktion erster Ordnung ist proportional dem Produkt der Schadstoff- und
Sauerstoffkonzentration. Wir machen aber die vereinfachende Annahme, daß trotz der
Sauerstoffverknappung immer noch soviel Sauerstoff vorhanden ist, daß die Oxidation
nicht wesentlich beeintr¨achtigt wird: Der limitierende Faktor ist also die
Schadstoffkonzentration. Auf diese Weise bekommen wir ein lineares Modell. Mit kS
bezeichnen wir die Reaktionsgeschwindigkeit. Pro kmol gel¨osten Schadstoff (in
Sauerstoff¨aquivalenten) werden in der Sekunde kS kmol Schadstoff oxidiert, und dabei
dieselbe Menge Sauerstoff verbraucht. Selbstverst¨andlich h¨angt kS stark von der
Temperatur ab. Weil in der Zelle i insgesamt die Schadstoffmenge Vi Li (t) gel¨ost ist,
ergibt sich als Bilanz:
Schadstoffbilanz der Oxidation: − kS Vi Li (t),
Sauerstoffbilanz der Oxidation: − kS Vi Li (t).
114
Modellierung der Sauerstoffaufnahme aus der Luft: Wir nehmen an, daß der Partialdruck
von Sauerstoff in der Luft konstant bleibt, was durchaus realistisch scheint. Nach dem
Henryschen Gesetz gibt es eine S¨attigungskonzentration cS von Sauerstoff im Wasser, die
diesem Partialdruck das Gleichgewicht h¨alt. Wird das Gleichgewicht gest¨ort, weicht also
Ci (t) von cS ab, so geht pro m2 Phasengrenzfl¨ache und pro Sekunde kO (cS − Ci (t)) kmol
Sauerstoff aus der Luft in L¨osung. Ist Ci (t) gr¨oßer als die S¨attigungskonzentration cS , so
ergibt sich ein negatives Vorzeichen: Sauerstoff entweicht aus der L¨osung in die Luft. Wir
bemerken, daß auch kO nat¨
urlich stark temperaturabh¨angig ist. Die Zelle i hat eine
Phasengrenzfl¨ache von ai Quadratmetern. Diese Zahl ergibt sich nicht einfach aus einer
geod¨atischen Vermessung des Flusses, es m¨
ussen auch die Str¨omungseigenschaften des
Wassers in Rechnung gestellt werden: Die feinen Tr¨opfchen u
¨ber einem Wasserfall haben
zum Beispiel eine gewaltige Oberfl¨ache. Die Bilanz der Sauerstoffaufnahme aus der Luft
ist daher
Sauerstoffzuwachs aus der Luft: kO ai (cS − Ci (t)).
Modellierung der Quellen: Jede Zelle hat eine Schadstoffzufuhr von Si (t) kmol pro
Sekunde (in Sauerstoff¨aquivalenten):
Schadstoffzuwachs aus Quellen: Si (t).
Summieren wir die Bilanzen aller Prozesse, so erhalten wir die Zuwachsrate f¨
ur den
gesamten Schadstoff und des gesamten Sauerstoff in Zelle i. Wir erinnern uns, das der
gesamte Schadstoff durch Vi Li (t) und der gesamte Sauerstoff durch Vi Ci (t) gegeben sind:
d
(Vi Li (t)) = F Li−1 (t) − F Li (t) − kS Vi Li (t) + Si (t),
dt
d
(Vi Ci (t)) = F Ci−1 (t) − F Ci (t) − kS Vi Li (t) + kO ai (cS − Ci (t)).
dt
Division durch die Konstante Vi liefert die Konzentrationsbilanz der Zelle i:
d
F
F
Si (t)
Li (t) = −( + kS )Li (t) + Li−1 (t) +
,
dt
Vi
Vi
Vi
d
F + kO ai
F
kO ai cS
Ci (t) = −
Ci (t) − kS Li (t) + Ci−1 (t) +
.
dt
Vi
Vi
Vi
Diese Gleichungen gelten f¨
ur jede der Zellen i = 1, · · · , N . Mit 10 Zellen w¨
urden wir ein
System von 20 Differentialgleichungen f¨
ur die 20 Zustandsgr¨oßen
L1 (t), · · · , L10 (t), C1 (t), · · · , C10 (t) erhalten. Mit Simulationssoftware f¨
ur
Differentialgleichungen l¨aßt sich dann das Zellenmodell auf dem Computer nachbilden.
Je kleiner die einzelnen Zellen, desto genauer trifft das Modell die Realit¨at, aber desto
gr¨oßer ist auch der Rechenaufwand und damit die Anf¨alligkeit auf Rundungsfehler.
Abbildung 10.3 zeigt eine M¨oglichkeit, den Fluß in Zellen einzuteilen, wenn die
vereinfachende Annahme fallen gelassen werden soll, daß die Schadstoffkonzentration
u
¨ber den ganzen Querschnitt des Flußbettes konstant sein soll. W¨ahrend das oben
entwickelte Zellenmodell sich mit einer Raumdimension begn¨
ugt, n¨amlich der
10. VERTEILTE PARAMETER
115
Tabelle 10.1. Eindimensionales Zellenmodell zu Beispiel 10.1
Gr¨oße
t
N
Li (t)
Ci (t)
L0 (t)
C0 (t)
Si (t)
Vi
ai
F
kS
kO
cS
Einheit
sec
1
kmol/m3
kmol/m3
kmol/m3
kmol/m3
kmol/sec
m3
m2
3
m /sec
1/sec
m/sec
kmol/m3
Modellgr¨
oßen
Benennung
Zeit
Anzahl der Zellen
Schadstoffkonzentration in Zelle i
Sauerstoffkonzentration in Zelle i
Schadstoffkonzentration Zufluß
Sauerstoffkonzentration Zufluß
Schadstoffzufluß zu Zelle i (Sauerstoff¨aquivalente)
Volumen in Zelle i
Wasser-Luft-Oberfl¨ache in Zelle i
Str¨omung
Reaktionsgeschwindigkeit Oxidation
Austauschgeschwindigkeit Wasser—Luft
S¨attigungskonzentration O2
Modellgleichungen
d
F
F
Si (t)
Li (t) = −( + kS )Li (t) + Li−1 (t) +
,
dt
Vi
Vi
Vi
d
F + kO ai
F
kO ai cS
Ci (t) = −
Ci (t) − kS Li (t) + Ci−1 (t) +
.
dt
Vi
Vi
Vi
(Die Gleichungen gelten f¨
ur jede Zellen-Nummer i = 1 · · · N .)
Abbildung 10.3. Dreidimensionales Zellenmodell
L¨angsrichtung des Flusses, bildet Abbildung 10.3 die Grundlage f¨
ur ein
dreidimensionales Modell. Die Modellierung w¨
urde wieder darin bestehen, f¨
ur jede Zelle
eine Mengenbilanz aufzustellen. Die Details werden komplizierter, weil jede Zelle mehr
als nur 2 Nachbarn hat. Vor allem m¨
ußte man die Str¨omungsverh¨altnisse nicht nur in
L¨angsrichtung, sondern auch quer und vertikal, und nach Tiefe und Abstand vom Ufer
116
aufgeschl¨
usselt kennen. Der Rechenaufwand f¨
ur dreidimensionale Modelle ist sehr viel
gr¨oßer als f¨
ur eindimensionale. Wenn man zur Erh¨ohung der Genauigkeit die
Seitenl¨angen der Zellen zehntelt, erh¨alt man im eindimensionalen Modell zur
¨
Uberdeckung
desselben Flußabschnittes zehnmal so viele Zellen, aber im
dreidimensionalen Modell die 1000-fache Anzahl der Zellen.
10.3. Das eindimensionale Modell in Matlab.
Dieser Abschnitt soll zeigen, wie man das Zellenmodell recht einfach und kompakt
programmieren kann. In Matlab kommt uns dabei die M¨oglichkeit, mit Vektoren zu
arbeiten, sehr zugute. Zun¨achst schreiben wir ein Matlab-Script, genannt zellen.m, das
die Modell-Parameter festlegt, einen Differentialgleichungsl¨oser (hier ode45) aufruft und
die Simulationsergebnisse zu vier Zeitpunkten plottet:
% zellen.m: script zur Fluss / Schadstoff Simulation mit dem Zellenmodell
% Bezeichnungen und Einheiten wie im Skriptum
global N V a F kS kO cS tein L_0 taus
N=12;
Laenge=10000;
v=2; % Stroemungsgeschwindigkeit in Zelle 1
Breite=10*ones(N,1); % Spalten-Vektor der Zellenbreiten und der
Tiefe=3*ones(N,1);
% Zellentiefen (hier alle gleich 10 bzw 3 gewaehlt)
V=Laenge/N*Breite.*Tiefe; % (hier alle Zellenlaengen gleich
a=Breite.*Laenge/N;
%
Laenge/N gewaehlt)
F=v*Breite(1).*Tiefe(1);
kO=3.98e-8; % Richtwerte fuer die Belueftungs-Parameter
cS=6.44e-4; % stammen aus der Literatur (Snape et al.)
kS=5e-4; % kS zur Verschoenerung angepasst
Lanf=0*ones(N,1); % Anfangskonzentrationen L_i(0) = 0
Canf=cS*ones(N,1); %
und C_i(0) = cS
L_0=5e-3; % voruebergehende Schadstoffkonzentration oberhalb
tein=0;
% der ersten Zelle, Beginn und
taus=360; % Ende der Schadstoffzufuhr
tend=3600; % Ende der Simulation
[t,y]=ode45(’reseze’,[0:tend/4:tend],[Lanf;Canf]);
for i=2:5, subplot(4,1,i-1)
plot([1:N],y(i,1:N),’*-’,[1:N],y(i,N+1:end),’o-’)
axis tight, title([’t = ’ num2str(t(i)) ’:’])
end
¨
Die Sicherung dieser Zeilen in einem file hat den Vorteil, dass man zur Anderung
der
Parameter nur den file editieren und wieder aufrufen muss. Um z.B. N auf 120 zu setzen,
braucht man nur eine Null anh¨angen, wie man im Programmierer-Jargon sagt. Dann l¨ost
10. VERTEILTE PARAMETER
117
das Programm eben ein System von 240 Differentialgleichungen, was f¨
ur dieses Beispiel
auf einem herk¨ommlichen PC auch nicht l¨anger als eine halbe Sekunde dauert. Die
Variablen Breite, Tiefe, V, a, Lanf und Canf sind Vektoren der L¨ange N. Die
Punkt-Operatoren .* und ./ bewirken komponentenweise Multiplikation und Division
von Matrizen (hier Spaltenvektoren) derselben Gr¨oße. Mit [a;b] werden die
Spaltenvektoren a und b zu einem Spaltenvektor aneinandergf¨
ugt. Wenn wir die
Zellenabmessungen nicht f¨
ur alle Zellen gleich w¨ahlen, k¨onnen wir dennoch den Rest des
Programms unver¨andert verwenden: das Programm ist wie das Modell f¨
ur N
Zellenvolumina Vi formuliert.
Die rechte Seite der Modellgleichungen aus Tabelle 10.1 wird von ode45 als function
reseze aufgerufen, die in folgendem file reseze.m enthalten ist:
function y=reseze(t,x)
% rechte Seite des Fluss / Schadstoff Zellenmodells
global N V a F kS kO cS tein L_0 taus
L=x(1:N);
C=x(N+1:2*N);
yL = -(F./V+kS).*L + F./V.*[L0(t); L(1:N-1)];
yC = -(F+kO*a)./V.*C - kS*L + F./V.*[cS; C(1:N-1)] + kO*a*cS./V;
y=[yL;yC];
%---------------------------------------------------------------function schadst=L0(t)
% Schadstoffeintrag in Zelle 1 zur Zeit t
global N V a F kS kO cS tein L_0 taus
if t >= tein && t <= taus
schadst = L_0;
else
schadst = 0;
end
Der Einfachheit halber haben wir die Terme Si (t) weggelassen - der Schadstoffeintrag
erfolgt nur in Zelle 1 und zwar u
¨ber die function L0(t). Die Schadstoffkonzentration,
die in Zelle 1 einstr¨omt, wird zwischen den Zeitpunkten tein und taus konstant auf L 0
gesetzt, sonst auf 0. Wir w¨ahlen L 0=5e-3 um h¨
ubsche Graphiken zu erhalten, dieser
Wert ist aber im Vergleich zu cS = 6.44 × 10−4 unrealistisch hoch (und liefert bei
genauerer Rechnung (z.B. N = 120) negative Ci ). Den Sauerstoffzustrom in Zelle 1
setzen wir konstant auf cS .
In zellen.m w¨ahlen wir das Simulationsintervall [0, 3600] und den Belastungszeitraum
[0, 360]. Dann geben wir im command window den Befehl zellen und erhalten die vier
118
Abbildung 10.4. Simulationsergebnisse mit 12 Zellen zu den Zeitpunkten
t=900, 1800, 2700, 3600. * = Schadstoffkonzentration, o = Sauerstoffkonzentration
−4
t = 900:
x 10
8
6
4
2
−4
1
x 10
2
3
4
5
6t = 1800:7
8
9
10
11
12
−4
1
x 10
2
3
4
5
6t = 2700:7
8
9
10
11
12
1
−4
x 10
2
3
4
5
6t = 3600:7
8
9
10
11
12
1
2
3
4
5
6
8
9
10
11
12
6
4
2
6
4
2
6
4
2
7
plots in Abbildung 10.3. Man sieht, wie das Maximum der Schadstoff-Belastung
flussabw¨arts str¨omt, und der Sauerstoffmangel dem Schadstoffabbau entsprechend
mitwandert. Das im Beispiel 10.1 gesuchte Minimum an Sauerstoff liegt in Abbildung
10.3 bei t = 2700 und zwar in Zelle 7 und betr¨agt etwa 1.55 × 10−4 . Nach Beendigung
der Belastung und durch die Aufnahme von Sauerstoff aus der Luft erholt sich die
Sauerstoffkonzentration; nach gen¨
ugend langer Zeit (nicht gezeichnet) wird der
Ausgangszustand wieder angen¨ahert, n¨amlich Schadstoff Li = 0 und Sauerstoff am
S¨attigungswert Ci = cS , i = 1, . . . , 12.
Es f¨allt auf, dass die zur Verdeutlichung mit Geraden verbundenen Zellen-Werte ziemlich
glatte r¨aumliche Verteilungen suggerieren. Zum Beispiel ist die Schadstoffverteilung bei
t = 900 schon bis zur Zelle 7 eindeutig positiv, wogegen die Verteilung durch die
Str¨omung physikalisch h¨ochstens bis x = 900v = 1800, also bis Zelle 3 reichen d¨
urfte
(wegen 1800/Zellenl¨ange = 1800/(10000/12) = 2.16). Dies ist ein charakteristischer
Fehler in Modellen mit gew¨ohnlichen Differentialgleichungen f¨
ur r¨aumlich gemittelte
¨
Zellen, weil durch die Koppelung der Gleichungen sich Anderungen
in einer Zelle
mathematisch sofort auf alle Zellen auswirken. Das Modell verschmiert also die
r¨aumliche Verteilung, was auch eine Verringerung des Maximums von Li und eine
10. VERTEILTE PARAMETER
119
Vergr¨oßerung des Minimums von Ci bewirkt. Immerhin ist wenigstens die Wanderung
des Maximums von Li selbst mit gut 2 Zellen pro plot (Abstand 900s) gleich der
Str¨omungsgeschwindigkeit v =2m/s, wie es sein soll.
10.4. Partielle Ableitungen.
Bevor wir weitere Modelle entwickeln, wiederholen wir kurz, was man unter einer
partiellen Ableitung versteht.
Definition 10.2. Ist w eine Funktion zweier Ver¨anderlicher t und x, so ist die partielle
Ableitung nach x an der Stelle (t, x) gegeben durch
∂
w(t, x + h) − w(t, x)
w(t, x) = lim
.
h→0
∂x
h
Das bedeutet, dass man sich eine der Ver¨anderlichen (hier t) festgehalten denkt, und von
der so entstandenen Funktion w(t, ·) einer Ver¨anderlichen (hier x) die Ableitung wie
u
¨blich definiert. Der Wert der partiellen Ableitung an der Stelle (t, x) ist der Anstieg der
Tangentialebene im Punkt t, x, w(t, x) an die von w(·, ·) gebildete Fl¨ache u
¨ber der
(t, x)-Ebene, und zwar in Richtung der entsprechenden Koordinate.
Definition 10.3. Sei f (x1 , x2 , . . . , xn ) eine Funktion, die von n Variablen abh¨angt. Die
partielle Ableitung von f nach der Variablen xj errechnet man, in dem man in f alle
Variablen außer xj wie Konstante ansieht, und f nach der Variablen xj differenziert.
Man schreibt f¨
ur die partielle Ableitung
∂
f (x1 , · · · , xn ).
∂xj
Die partiellen Ableitungen h¨oherer Ordnung
∂k
∂xjk ∂xjk−1 · · · ∂xj1
f (x1 , · · · , xn )
errechnet man, indem man f hintereinander nach den Variablen xj1 , xj2 , · · · , xjk
differenziert.
Beispiel 10.4. Berechnen Sie die partiellen Ableitungen erster und zweiter Ordnung von
f (x, y) = 3x4 y 2 + x2 + 5y.
120
L¨
osung.
∂
(3x4 y 2 + x2 + 5y) = 12x3 y 2 + 2x
∂x
∂
(3x4 y 2 + x2 + 5y) = 6x4 y + 5
∂y
∂2
∂
(3x4 y 2 + x2 + 5y) =
(12x3 y 2 + 2x) = 36x2 y 2 + 2,
2
∂x
∂x
∂2
∂
(3x4 y 2 + x2 + 5y) =
(12x3 y 2 + 2x) = 24x3 y
∂y∂x
∂y
2
∂
∂
(3x4 y 2 + x2 + 5y) =
(6x4 y + 5) = 24x3 y
∂x∂y
∂x
2
∂
∂
(3x4 y 2 + x2 + 5y) =
(6x4 y + 5) = 6x4
2
∂y
∂y
2
2
∂
∂
Beachten Sie, daß die Reihenfolge bei den partiellen Ableitungen ∂x∂y
und ∂y∂x
keine
Rolle spielt. Das gilt f¨
ur alle Funktionen, deren partielle Ableitungen 2. Ordnung stetig
sind.
10.5. Advektion und Reaktion in einer Dimension.
Die Schw¨ache des Zellenmodells besteht in der Annahme, daß die Konzentrationen
innerhalb einer Zelle, also u
¨ber einen gewissen Bereich des Flußlaufes, konstant sind.
Erst im Grenz¨
ubergang zu “unendlich vielen, unendlich kleinen” Zellen wird das Modell
realistisch. Die mathematische Formulierung eines solchen Grenz¨
uberganges f¨
uhrt zu
einer Differentialgleichung, die partielle Ableitungen der Zustandsgr¨oßen enth¨alt, einer
sogenannten partiellen Differentialgleichung. W¨ahrend die Theorie und auch die
numerische L¨osung von partiellen Differentialgleichungen auf dem Computer ziemlich
kompliziert ist, gen¨
ugt zum Erstellen und Verst¨andnis der Modelle oft ein bißchen
naturwissenschaftlicher Hausverstand.
Wir nehmen wieder an, daß in jedem Querschnitt quer zur Str¨omungsrichtung alle
Konzentrationen konstant sind, sodaß wir uns mit einem eindimensionalen Modell
begn¨
ugen k¨onnen. Wir f¨
uhren daher nur eine Raumkoordinate ein: Den Abstand x eines
Flußquerschnittes vom Fabriksauslaß. Wir betrachten den Flußabschnitt vom Auslaß bis
l Meter stromabw¨arts, also liegt x im Intervall von 0 bis l.
Zu jedem Zeitpunkt t herrscht am Flußquerschnitt mit Raumkoordinate x eine
Schadstoffkonzentration von L(t, x) kmol/m3 Sauerstoff¨aquivalenten, und eine
Sauerstoffkonzentration von C(t, x) kmol/m3 . Die Zustandsgr¨oßen C und L sind also
Funktionen, die von zwei Variablen, n¨amlich Zeit und Raum, abh¨angen.
Die Idee der Modellbildung besteht darin, eine kleine Zelle im Fluß, n¨amlich den Bereich
zwischen den Raumkoordinaten x und x + ∆x zu betrachten, wobei ∆x eine L¨ange ist,
10. VERTEILTE PARAMETER
121
die wir uns sehr klein vorstellen. F¨
ur diese Zelle stellen wir eine Mengenbilanz auf, und
betrachten dann, was geschieht, wenn wir die L¨ange ∆x gegen Null gehen lassen.
(Achtung: Entsprechend einer verbreiteten Konvention haben wir hier eine Gr¨oße,
n¨amlich die L¨ange der Zelle, mit zwei Buchstaben ∆x bezeichnet. Das ist eine einzige
Gr¨oße, und nicht etwa ein Produkt einer Zahl ∆ mit x.)
Geometrie und Hydrodynamik des Flusses: Mit A(x) bezeichnen wir die
Querschnittsfl¨ache des Flusses an der Stelle mit Abstand x stromabw¨arts vom Auslaß.
Die Str¨omungsgeschwindigkeit des Wassers wird nat¨
urlich nahe am Ufer und am Grund
des Flusses langsamer sein als in der Mitte. Wir bilden den Mittelwert der
Str¨omungsgeschwindigkeit u
¨ber den ganzen Querschnitt an der Stelle x und benennen
ihn mit u(x). Die Wassermenge, die dann pro Sekunde diesen Querschnitt durchsetzt, ist
A(x)u(x). Wir nehmen an, daß A und u nicht von der Zeit abh¨angen, also daß der Fluß
im beobachteten Zeitraum immer die gleiche Menge Wasser f¨
uhrt. Zur Modellierung der
Sauerstoffaufnahme ben¨otigen wir wieder die Wasser-Luft-Grenzfl¨ache. Wir nehmen an,
daß an der Stelle x auf ein kurzes L¨angenst¨
uck ∆x die Oberfl¨ache b(x)∆x kommt. Die
Gr¨oße b(x) hat die Dimension einer L¨ange. Sie h¨angt von der Breite des Flusses und dem
Str¨omungsverhalten ab: Wellen und schwebende Tr¨opfchen vergr¨oßern die
Phasengrenzfl¨ache.
Schadstoff- und Sauerstoffbestand in der Zelle: Wir gehen davon aus, daß die L¨ange ∆x
der Zelle ziemlich klein ist. Wenn sich der Querschnitt A entlang dieses kurzen St¨
uckes
nicht stark a¨ndert, enth¨alt die Zelle ungef¨ahr ein Wasservolumen von A(x)∆x. In diesem
Wasservolumen sind Sauerstoff und Schadstoff mit Konzentrationen C(t, x) und L(t, x)
gel¨ost. Daher ergibt sich
Schadstoffgehalt in der Zelle: L(t, x)A(x)∆x,
Sauerstoffgehalt in der Zelle: C(t, x)A(x)∆x.
Modellierung der Advektion: Von stromaufw¨arts dringen pro Sekunde A(x)u(x) m3
Wasser durch den Querschnitt am Zellanfang x in die Zelle ein. Dieses Wasser hat eine
Sauerstoffkonzentration von C(t, x) und eine Schadstoffkonzentration von L(t, x), daher
treten mit dem Wasser in jeder Sekunde A(x)u(x)C(t, x) kmol Sauerstoff und
A(x)u(x)L(t, x) kmol Schadstoff (in Sauerstoff¨aquivalenten) ein. Am unteren Ende der
Zelle, an der Stelle x + ∆x, werden nach derselben Schlußrechnung jede Sekunde
A(x + ∆x)u(x + ∆x)C(t, x + ∆x) kmol Sauerstoff und A(x + ∆x)u(x + ∆x)L(t, x + ∆x)
kmol Schadstoff mit dem abfließenden Wasser davongetragen. Das ergibt die folgende
Bilanz
Schadstoffbilanz der Advektion:
A(x)u(x)L(t, x) − A(x + ∆x)u(x + ∆x)L(t, x + ∆x),
Sauerstoffbilanz der Advektion:
A(x)u(x)C(t, x) − A(x + ∆x)u(x + ∆x)C(t, x + ∆x).
122
Modellierung der verteilten Quellen: Wir nehmen an, daß an der Stelle x auf ein kurzes
L¨angenst¨
uck der L¨ange ∆x pro Sekunde die Schadstoffmenge s(t, x)∆x vom Ufer her
zufließt. Wir erhalten daher f¨
ur unsere Zelle:
Schadstoffzufluß aus verteilten Quellen: s(t, x)∆x.
Auf die Modellierung der punktf¨ormigen Quelle bei x = 0 kommen wir sp¨ater zu
sprechen. Wenn 0 nicht gerade innerhalb unserer sehr kleinen Zelle zu liegen kommt,
spielt diese Quelle f¨
ur die Zelle ohnehin nur indirekt u
¨ber die Advektion eine Rolle.
Modellierung der Oxidation: Wie im Zellenmodell beschreiben wir die Oxidation als eine
Reaktion erster Ordnung in Abh¨angigkeit von der Schadstoffkonzentration, das heißt, wir
gehen davon aus, daß das Sauerstoffangebot ausreichend ist. Mit kS bezeichnen wir
wieder die Reaktionsgeschwindigkeit der Oxidation. In der Zelle ist die Schadstoffmenge
A(x)∆xL(t, x) zur Oxidation vorhanden, davon oxidiert in der Sekunde
kS A(x)∆xL(t, x). Daraus ergibt sich die folgende Bilanz:
Schadstoffbilanz der Oxidation: − kS A(x)∆xL(t, x),
Sauerstoffbilanz der Oxidation: − kS A(x)∆xL(t, x).
Modellierung der Sauerstoffaufnahme aus der Luft: Als Phasengrenzfl¨ache stehen b(x)∆x
Quadratmeter zur Verf¨
ugung, das Konzentrationsgef¨alle von S¨attigungskonzentration cS
zur tats¨achlichen Sauerstoffkonzentration ist cS − C(t, x). Mit einer
Austauschgeschwindigkeit kO ergibt sich die Bilanz
Sauerstoffzuwachs aus der Luft: kO b(x)∆x(cS − C(t, x)).
¨
Mengenbilanz: Die Mengenbilanz f¨
ur die Zelle gibt nun die Anderungsraten
der
Schadstoff- und Sauerstoffmengen in der Zelle an. Dazu werden die Beitr¨age von
Advektion, Quellen, Oxidation und Reaeration aus der Luft summiert. Die
¨
Anderungsraten
sind, wie gewohnt, die Ableitungen der Mengen nach der Zeit. Da jetzt
außer der Zeit noch die Raumvariable zur Verf¨
ugung steht, schreiben wir die Ableitungen
als partielle Ableitungen:
∂
(A(x)L(t, x)∆x) = A(x)u(x)L(t, x) − A(x + ∆x)u(x + ∆x)L(t, x + ∆x)
∂t
+ s(t, x)∆x − kS A(x)∆xL(t, x),
∂
(A(x)C(t, x)∆x) = A(x)u(x)C(t, x) − A(x + ∆x)u(x + ∆x)C(t, x + ∆x)
∂t
− kS A(x)∆xL(t, x) + b(x)∆x(cS − C(t, x)).
10. VERTEILTE PARAMETER
123
Wir k¨
urzen beide Gleichungen durch ∆x und erhalten
∂
A(x)u(x)L(t, x) − A(x + ∆x)u(x + ∆x)L(t, x + ∆x)
(A(x)L(t, x)) =
∂t
∆x
+ s(t, x) − kS A(x)L(t, x),
∂
A(x)u(x)C(t, x) − A(x + ∆x)u(x + ∆x)C(t, x + ∆x)
(A(x)C(t, x)) =
∂t
∆x
− kS A(x)L(t, x) + kO b(x)(cS − C(t, x)).
Wenn wir jetzt in der Mengenbilanz der Zelle die L¨ange ∆x gegen Null gehen lassen,
wenden wir die Definition 10.2 der partiellen Ableitung f¨
ur h = ∆x und
w(t, x) = A(x)u(x)L(t, x) bzw. A(x)u(x)C(t, x) an. Wir erhalten im Grenzwert:
∂
∂
L(t, x) = − (A(x)u(x)L(t, x)) + s(t, x) − kS A(x)L(t, x),
∂t
∂x
∂
∂
A(x) C(t, x) = − (A(x)u(x)C(t, x)) − kS A(x)L(t, x) + kO b(x)(cS − C(t, x)).
∂t
∂x
A(x)
Das Modell f¨
uhrt also auf ein System von zwei partiellen Differentialgleichungen, in
denen sowohl die Raum- als auch die Zeitableitungen der Zustandsgr¨oßen L(t, x) und
C(t, x) auftreten.
Modellierung der Zufl¨
usse und der punktf¨ormigen Quelle: Wie schon im Zellenmodell,
bezeichnen wir mit C0 (t) die Sauerstoffkonzentration des zufließenden Wassers
unmittelbar oberhalb des Fabriksauslasses, und mit L0 (t) die Schadstoffkonzentration
des zufließenden Wassers. Mit S(t) bezeichnen wir den Schadstoffausstoß der Fabrik, in
kmol Sauerstoff¨aquivalenten pro Sekunde. Wir nehmen an, daß wir diese Gr¨oßen kennen.
An der Stelle x = 0, unmittelbar nach dem Fabriksauslaß, betr¨agt also die
Sauerstoffkonzentration
C(t, 0) = C0 (t).
Unmittelbar unter dem Auslaß fließen pro Sekunde u(0)A(0) m3 Wasser. Wenn wir
annehmen, daß der Beitrag des Auslasses zur Wassermenge verschwindend klein ist,
fließt fast dieselbe Menge von oberhalb des Auslasses zu. Sie tr¨agt pro Sekunde von
flußaufw¨arts L0 (t)u(0)A(0) kmol Schadstoff (Sauerstoff¨aquivalente) heran. Dazu kommen
jede Sekunde S(t) kmol Schadstoff von der Fabrik. Daher fließen unmittelbar unter dem
Auslaß jede Sekunde S(t) + u(0)A(0)L0 (t) kmol Schadstoff, gel¨ost in u(0)A(0) m3
Wasser. Das ergibt die Konzentration (die wir der K¨
urze halber mit L1 (t) bezeichnen):
L(t, 0) = L0 (t) +
S(t)
= L1 (t).
u(0)A(0)
Die Information u
¨ber die Zuflußkonzentrationen beschreibt die beiden Zustandsgr¨oßen L
und C an einem Randpunkt des betrachteten Gebietes, n¨amlich an x = 0. Sie sind
Randbedingungen.
124
Anfangsbedingungen: Um die Entwicklung der Schadstoffkonzentration im Fluß
vorhersagen zu k¨onnen, brauchen wir letztlich noch Informationen u
¨ber den
Anfangszustand. Wir nehmen an, daß wir die Schadstoffkonzentration La und die
Sauerstoffkonzentration Ca zum Zeitpunkt t = 0 im ganzen untersuchten Flußgebiet
kennen.
L(0, x) = La (x),
C(0, x) = Ca (x).
Diese Bedingungen legen die Zustandsgr¨oßen zu Anfang des betrachteten Zeitraumes
fest, sie sind Anfangsbedingungen.
Insgesamt erhalten wir ein Modell aus zwei partiellen Differentialgleichungen, mit zwei
Anfangsbedingungen, und zwei Randbedingungen. Das vollst¨andige Modell ist in
Tabelle 10.2 zusammengefaßt.
Merksatz 10.5. Eine partielle Differentialgleichung beschreibt eine unbekannte
Funktion, die von mehreren Variablen abh¨angt, durch eine Beziehung zwischen der
Funktion und ihren partiellen Ableitungen. In der Modellierung verteilter Zust¨ande treten
partielle Differentialgleichungen zusammen mit Anfangsbedingungen und
Randbedingungen auf. Anfangsbedingungen geben die Werte der unbekannten
Funktion(en) zu Anfang des Beobachtungszeitraumes an. Randbedingungen
charakterisieren die unbekannten Funktionen am Rand des beobachteten Gebietes.
Bemerkung. In der soeben durchgef¨
uhrten Modellbildung haben wir alle betrachteten
Teil-Bilanzen (f¨
ur die Str¨omung, die Oxidation, den Sauerstoffeintrag und den
Schadstoffzutritt aus den Ufern) zusammengefasst. Falls man nur die Advektion
(=Str¨omung) ber¨
ucksichtigt und alle anderen Effekte einmal wegl¨asst, dann steht in der
durch ∆x dividierten Mengenbilanz f¨
ur die ∆x-Zelle rechts nur der Differenzenquotient
(A(x)u(x)L(t, x) − A(x + ∆x)u(x + ∆x)L(t, x + ∆x))/∆x. Dann folgt aus der
Inkompressibilit¨at des Wassers A(x)u(x) = F mit einer von x und t unabh¨angigen
Konstanten F . Diese kann also, nach Grenz¨
ubergang ∆x → 0, aus der partiellen
Ableitung nach x herausgezogen werden. K¨
urzt man die gesamte Gleichung dann durch
A(x), so erh¨alt man eine partielle Differentialgleichung der Gestalt
∂
∂
L(t, x) + u(x) L(t, x) = 0, t > 0, x ∈ (0, xL ).
∂t
∂x
Dies ist die klassische Advektionsgleichung f¨
ur die Bewegung einer inkrompressiblen
Fl¨
ussigkeit mit der r¨aumlich und zeitlich ver¨anderlichen Eigenschaft L(t, x). Die
Str¨omungsgeschwindigkeit u(x) ist umgekehrt proportional zur Breite des
Str¨omungskanals A(x), n¨amlich u(x) = F/A(x). F¨
ur u > 0 ben¨otigt man
Randbedingungen am Anfang des Kanals bei x = 0, f¨
ur u < 0 kann hingegen nur eine
Randbedingung am Ende bei x = xL vorgegeben sein. Eine Geschwindigkeitsverteilung
u(x) mit wechselndem Vorzeichen ist wegen der Inkompressibilit¨at unm¨oglich (bei
Kompressibilit¨at verwendet man verallgemeinerte Konzepte f¨
ur zul¨assige ’L¨osungen’).
10. VERTEILTE PARAMETER
125
Tabelle 10.2. Reaktions-Advektionsgleichung zu Beispiel 10.1
Gr¨oße
Einheit
t
sec
x
m
L(t, x)
kmol/m3
C(t, x)
kmol/m3
L0 (t)
kmol/m3
C0 (t)
kmol/m3
S(t)
kmol/sec
s(t, x) kmol/(m sec)
La (x)
kmol/m3
Ca (x)
kmol/m3
kS
1/sec
kO
m/sec
cS
kmol/m3
A(x)
m2
u(x)
m/sec
b(x)
m
Modellgr¨
oßen
Benennung
Zeit
Distanz vom Auslaß
¨
Schadstoffkonzentration (O2 -Aquivalente)
Sauerstoffkonzentration
Schadstoffkonzentration Zufluß
Sauerstoffkonzentration Zufluß
Schadstoffzufluß von Auslaß
Schadstoffzufuhr verteilte Quellen pro L¨ange
Schadstoffkonzentration zu Anfang
Sauertoffkonzentration zu Anfang
Reaktionsgeschwindigkeit Oxidation
Austauschgeschwindigkeit Wasser—Luft
S¨attigungskonzentration O2
Querschnittsfl¨ache
Str¨omungsgeschwindigkeit
Phasengrenzfl¨ache pro L¨ange
Partielle Differentialgleichungen
∂
∂
L(t, x) = − (A(x)u(x)L(t, x))
∂t
∂x
+ s(t, x) − kS A(x)L(t, x),
∂
∂
A(x) C(t, x) = − (A(x)u(x)C(t, x))
∂t
∂x
− kS A(x)L(t, x) + kO b(x)(cS − C(t, x)).
A(x)
(10.1)
Randbedingungen
L(t, 0) = L1 (t) = L0 (t) +
S(t)
,
A(0)u(0)
C(t, 0) = C0 (t).
Anfangsbedingungen
L(0, x) = La (x), C(0, x) = Ca (x).
10.6. Charakteristiken.
Die Hauptstoßrichtung dieses Kapitels ist die Modellbildung, und die ist nun
abgeschlossen, sofern wir Diffusion und Dispersion nicht ber¨
ucksichtigen wollen. Auf
Diffusion und Dispersion kommen wir sp¨ater zur¨
uck. Wir werden jetzt die vorliegenden
Modellgleichungen l¨osen. Das Modell in Tabelle 10.2 sieht zwar abschreckend aus, aber es
ist ein sehr einfaches Beispiel einer partiellen Differentialgleichung. Partielle
126
Differentialgleichungen k¨onnen außerordentlich schwierig sein, sowohl bei der
theoretischen Behandlung als auch bei der numerischen L¨osung auf dem Computer.
Um die Darstellung der Berechnungen einfacher zu gestalten, machen wir die folgenden
Annahmen: Querschnitt A(x) = A, Str¨omungsgeschwindigkeit u(x) = u, und die
Phasenoberfl¨ache pro L¨angeneinheit b(x) = b sind konstant, also unabh¨angig von x.
(Statt von einem Fluß sprechen wir wohl eher von einem regulierten Kanal.) Die einzige
Schadstoffquelle ist der Fabrikauslaß: s(t, x) = 0.
Der Beobachterstandpunkt bei der Modellbildung war fest: Wir haben aus dem Flußlauf
eine kurze Zelle herausgegriffen, und berechnet, wieviel Schadstoff und Sauerstoff das
Wasser im Lauf der Zeit durch diese Zelle hindurchtr¨agt. Wir k¨onnten uns vorstellen,
daß ein Beobachter an einer festen Stelle x am Ufer steht, und an einer Angel seine
Meßsonden ausgeh¨angt hat. Er mißt damit den Verlauf der Konzentrationen L(t, x) und
C(t, x) im Lauf der Zeit. Wenn sich diese Konzentrationen im Lauf der Zeit ¨andern, kann
das zwei Gr¨
unde haben: Einerseits verschieben die chemischen Reaktionen die
Konzentrationen. Andererseits kann der Schadstoffausstoß der Quellen zu verschiedenen
Zeiten unterschiedlich sein, und dadurch kann Wasser, das erst auf die Meßtstelle
zufließt, eine andere Vorgeschichte haben und andere Schadstoffwerte liefern als Wasser,
das eben vorbeigeflossen ist. Weil der Beobachter stets am selben Standort x steht, ist x
¨
f¨
ur ihn eine Konstante. Die Anderungsraten
von L und C, die er wahrnimmt, sind die
Ableitungen der Konzentrationen, die man erh¨alt, wenn man nach der Zeit differenziert
∂
∂
und x konstant l¨aßt: die partiellen Ableitungen ∂t
L(t, x) und ∂t
C(t, x).
Die Idee des folgenden L¨osungswegs besteht in einer geschickten Ver¨anderung des
Beobachterstandpunktes: Wir stellen uns vor, daß der Beobachter zu einem gegebenen
Zeitpunkt t0 neben dem Fabrikauslaß in ein Boot steigt und sich mit der Str¨omung
treiben l¨aßt. Vom Boot aus h¨angt er seine Meßsonden in das Wasser, sodaß sie mit der
Str¨omung mittreiben. Weil die Str¨omung v¨ollig von Zuf¨alligkeiten frei ist (laut
Modellannahme), sind die Sonden immer vom selben Wasser umgeben, und alle
¨
Anderungen
der Konzentrationen, die dieser Beobachter wahrnimmt, sind auf die
chemischen Reaktionen zur¨
uckzuf¨
uhren.
Weil der Beobachter zum Zeitpunkt t = t0 an der Stelle x = 0 ablegt, und mit der
Str¨omungsgeschwindigkeit u fortgetragen wird, befindet er sich zur Zeit t an der Stelle
x(t) = (t − t0 )u. An dieser Stelle mißt er die Konzentrationswerte L(t, (t − t0 )u) und
¨
C(t, (t − t0 )u). Mit der Zeit k¨onnen sich diese Werte ¨andern, und die Anderungsraten
d
d
bezeichnen wir mit dt L(t, (t − t0 )u) und dt C(t, (t − t0 )u). Das sind nicht mehr die
partiellen Ableitungen nach der Zeit, denn die Raumkoordinate x ver¨andert sich mit der
Zeit: das Boot treibt weiter. Wir sprechen von totalen Ableitungen: Sie resultieren aus
der Ver¨anderung von Raum und Zeit zugleich.
Um die totale Ableitung mathematisch zu beschreiben, halten wir uns die Bedeutung der
∂
partiellen Ableitungen vor Augen: Die partielle Ableitung ∂x
L(t, x) ist die
10. VERTEILTE PARAMETER
127
¨
Anderungsrate
bezogen auf L¨ange, die sich ergibt, wenn man die Konzentration L
gleichzeitig an der Stelle x und ein “unendlich kleines” St¨
uck ∆x weiter stromabw¨arts
∂
¨
mißt. Die partielle Ableitung ∂t L(t, x) ist die Anderungsrate bezogen auf Zeit, die sich
ergibt, wenn man die Konzentration an derselben Stelle zur Zeit t und eine “unendlich
kleine” Zeitspanne ∆t sp¨ater mißt:
∂
L(t, x),
∂x
∂
L(t + ∆t, x) − L(t, x) ≈ ∆t L(t, x).
∂t
Wenn der Beobachter vom Boot aus kurz hintereinander mißt, mißt er an der Stelle
(t − t0 )u zur Zeit t, und knapp danach an der Stelle (t − t0 )u + u∆t zur Zeit t + ∆t.
¨
Allein durch den zeitlichen Unterschied der Messungen w¨
urde sich eine Anderung
der
Meßwerte von
∂
∆t L(t, x)
∂t
¨
ergeben. Außerdem erfolgt aber eine Anderung des Standortes um den Betrag ∆x = u∆t
und schl¨agt sich in einer zus¨atzlichen Verschiebung der Meßwerte von
L(t, x + ∆x) − L(t, x) ≈ ∆x
u∆t
∂
L(t, x)
∂x
¨
nieder. Daher ist die gesamte Anderung
∆t
∂
∂
L(t, x) + u∆t L(t, x).
∂t
∂x
¨
Die Anderungsrate,
bezogen auf das Zeitintervall ∆t, ist dann die totale Ableitung:
(10.2)
d
L(t, u(t − t0 )) =
dt
d
C(t, u(t − t0 )) =
dt
∂
∂
L(t, u(t − t0 )) + u L(t, u(t − t0 )),
∂t
∂x
∂
∂
C(t, u(t − t0 )) + u C(t, u(t − t0 )).
∂t
∂x
Merksatz 10.6. Ein Beobachter durchl¨auft den Raum entlang einer Kurve
x1 (t), · · · , xn (t). Die totale Ableitung einer Funktion f (x1 (t), · · · , xn (t)) entlang dieser
Kurve errechnet sich durch
∂
d
∂
d
d
f=
f · x1 + · · · +
f · xn .
dt
∂x1
dt
∂xn
dt
¨
Die totale Ableitung beschreibt die Anderungsrate
von f bezogen auf die Zeit, die ein
Beobachter feststellt, der jeweils zum Zeitpunkt t von seinem augenblicklichen Standort
mit den Raumkoordinaten x1 (t), · · · , xn (t) mißt.
Wir ordnen die partiellen Differentialgleichungen (10.1) aus dem Modell in Tabelle 10.2
etwas um und heben den konstanten Faktor A(x) = A heraus. Unter Beachtung, daß
128
auch u konstant und s = 0 ist, erhalten wir
∂
∂
L(t, x) + u L(t, x) = −kS L(t, x),
∂t
∂x
∂
∂
kO b
C(t, x) + u C(t, x) = −kS L(t, x) +
(cS − C(t, x)).
∂t
∂x
A
Wir setzen jetzt in die totalen Ableitungen (10.2) der Schadstoff- und
Sauerstoffkonzentration die partiellen Differentialgleichungen (10.1) aus dem Modell in
Tabelle 10.2 ein:
d
∂
∂
L(t, (t − t0 )u) = L(t, (t − t0 )u) + u L(t, (t − t0 )u)
dt
∂t
∂x
= −kS L(t, (t − t0 )u),
d
∂
∂
C(t, (t − t0 )u) = C(t, (t − t0 )u) + u C(t, u(t − t0 ))
dt
∂t
∂x
kO b
= −kS L(t, (t − t0 )u) +
(cS − C(t, (t − t0 )u)
A
Alle partiellen Ableitungen sind verschwunden, u
¨brig bleiben nur mehr die totalen
¨
Ableitungen, die Anderungsraten,
die der Beobachter im Boot feststellt. Und diese
totalen Ableitungen resultieren aus den Konzentrations¨anderungen durch Oxidation und
Sauerstoffaufnahme aus der Luft. Die Effekte der Str¨omung werden vom Boot aus nicht
wahrgenommen, weil es mitschwimmt. Statt zwei partiellen Differentialgleichungen hat
der Beobachter im Boot nur zwei gew¨ohnliche Differentialgleichungen zu l¨osen, um seine
Meßwerte zu erkl¨aren:
(10.3)
d
L(t, (t − t0 )u) = −kS L(t, (t − t0 )u),
dt
d
kO b
C(t, (t − t0 )u) = −kS L(t, (t − t0 )u) +
(cS − C(t, (t − t0 )u))
dt
A
Das Boot bewegt sich auf dem Weg x(t) = (t − t0 )u durch den Raum. Diese Kurve heißt
in der Mathematik eine Charakteristik der partiellen Differentialgleichungen (10.1): Eine
Charakteristik verl¨auft so im Raum, daß die totale Ableitung der L¨osung entlang einer
Charakteristik gerade alle partiellen Ableitungen aus der Differentialgleichung
zusammenfaßt. Die L¨osung wird entlang der Charakteristik also durch eine
Differentialgleichung beschrieben, die nur eine Ableitung enth¨alt, n¨amlich die totale
Ableitung nach der Zeit.
Wir wiederholen: Die Konstruktion von Charakteristiken, und damit die Umwandlung
der partiellen Differentialgleichungen in gew¨ohnliche Differentialgleichungen, war nur
m¨oglich, weil die Str¨omungsverh¨altnisse a priori klar und frei von Zuf¨alligkeiten waren.
Die Methode scheitert, wenn Diffusions- und Dispersionseffekte zum Tragen kommen.
10. VERTEILTE PARAMETER
129
10.7. L¨
osung der Reaktions-Advektions-Gleichungen.
Die L¨osung der gew¨ohnlichen Differentialgleichungen 10.3 ist nun eine Standardaufgabe
f¨
ur eine Vorlesung aus Differentialgleichungen. Zur Abk¨
urzung schreiben wir
kO b
.
A
k1 =
Wir erhalten die L¨osungen
L(t, (t − t0 )u) = e−kS (t−t0 ) L1 (t0 ),
C(t, (t − t0 )u) = e−k1 (t−t0 ) C0 (t0 ) + cS [1 − e−k1 (t−t0 ) ]
(10.4)
−
kS L1 (t0 ) −kS (t−t0 )
[e
− e−k1 (t−t0 ) ].
k1 − kS
Der Vollst¨andigkeit halber geben wir den L¨osungsweg an.
Wir l¨osen zuerst die Differentialgleichung f¨
ur L, denn sie ist von C entkoppelt: Die
Sauerstoffkonzentration kommt in der Gleichung gar nicht vor:
d
L(t, (t − t0 )u) = −kS L(t, (t − t0 )u)
dt
Das ist eine homogene lineare Differentialgleichung
d
y(t) = −kS y(t)
dt
mit y(t) = L(t, (t − t0 )u). Die L¨osung ergibt die Exponentialfunktion
y(t) = e−kS (t−t0 ) y(t0 ),
also
L(t, (t − t0 )u) = e−kS (t−t0 ) L(t0 , 0) = e−kS (t−t0 ) L1 (t0 ).
Wir setzen in die Gleichung f¨
ur die Sauerstoffkonzentration ein:
d
kO b
C(t, (t − t0 )u) = −kS L(t, (t − t0 )u) +
(cS − C(t, (t − t0 )u)
dt
A
= −k1 C(t, (t − t0 )u) + k1 cS − kS L1 (t0 )e−kS (t−t0 ) .
Diese Gleichung ist eine inhomogene lineare Differentialgleichung der Form
d
y(t) = −k1 y(t) + g(t),
dt
mit y(t) = C(t, (t − t0 )u) und g(t) = k1 cS − kS L1 (t0 )e−kS (t−t0 ) . Ihre L¨osung ist nach der “Variation der
Konstanten”-Formel
y(t) = e−k1 (t−t0 ) y(0) +
t
e−k1 (t−s) g(s) ds,
t0
das heißt also
C(t, (t − t0 )u) = e−k1 (t−t0 ) C(t0 , 0) +
t
e−k1 (t−s) [k1 cS − kS L1 (t0 )e−kS (s−t0 ) ] ds
t0
= e−k1 (t−t0 ) C0 (t0 ) + cS [1 − e−k1 (t−t0 ) ]
−
kS L1 (t0 ) −kS (t−t0 )
[e
− e−k1 (t−t0 ) ].
k1 − kS
130
Um die L¨osung L(t, x) zu bestimmen, m¨
ussen wir jetzt nur den Zeitpunkt t0 so w¨ahlen,
daß (t − t0 )u = x gilt, also
x
t0 = t − .
u
Dann ist
L(t, x) = L(t, (t − t0 )u) = e−kS (t−t0 ) L1 (t0 )
x
= e−kS x/u L1 (t − ).
u
Ebenso errechnet man
−k1 x/u
C(t, x) = e
kS L1 (t − ux ) −kS x/u
x
−k1 x/u
C0 (t − ) + cS 1 − e
e
− e−k1 x/u .
−
u
k1 − kS
Das sind die Streeter-Phelps-Gleichungen. Sie gelten f¨
ur solche (t, x), f¨
ur die t − ux ≥ t0 ,
andernfalls k¨onnten die Funktionen der Randbedingungen L1 und C0 nicht ausgewertet
werden. F¨
ur x ≥ u(t − t0 ) muß man die Anfangsbedingungen statt der Randbedingungen
verwenden. Die Notwendigkeit dieser Fallunterscheidung ist anschaulich verst¨andlich:
F¨
ur grosse x und kleine t sind die Anfangsbedingungen maßgeblich, die Information aus
dem linken Rand kommt erst nach gen¨
ugend langer Zeit nach x. F¨
ur grosses t und
kleines x sind die Anfangsbedingungen schon zur G¨anze vorbeigestr¨omt, nur mehr die
Zustr¨omung aus dem linken Rand ist maßgeblich.
Als Anwendung k¨onnte man das Minimum von C (Nullsetzen der Ableitung nach x),
also f¨
ur festes t jenen Punkt im Fluß bestimmen, an dem die Sauerstoffkonzentration
minimal ist, und damit das maximale Sauerstoffdefizit vorhersagen. Daraus kann man
dann umgekehrt berechnen, wie groß L1 h¨ochstens sein darf, damit die
Sauerstoffkonzentration einen gewissen vorgegebenen Standard nicht unterschreitet.
Beobachten Sie auch einige qualitative Aussagen der Formeln: Die Voraussagen f¨
ur
C(t, x) und L(t, x) h¨angen nur vom gesamten Schadstoffzufluß L1 (dadurch indirekt vom
Schadstoffausstoß L0 ) zur Zeit t − x/u ab. Die Zeitspanne x/u ist genau die Zeit, die der
Schadstoff braucht, um mit der Str¨omung vom Auslaß zur Stelle x abzutreiben. Folgt
man dem Schadstoff auf seinem Weg stromabw¨arts, so ergibt sich als Grenzwert
lim L(t, (t − t0 )u) = 0,
t→∞
lim C(t, (t − t0 )u) = cS .
t→∞
Weit stromabw¨arts regeneriert sich das Wasser auf die volle Sauerstoffs¨attigung, und
aller Schadstoff wird abgebaut.
10. VERTEILTE PARAMETER
131
Abbildung 10.5. Zufallsbewegung
10
10
8
6
5
4
2
0
0
−2
−5
−5
0
5
10
−4
−5
0
5
10
Linkes Bild: Verlauf der Bewegung. Rechtes Bild: Anfangs- und Endpositionen.
Kreise: Anfangsposition, Kreuze: Endposition
10.8. Modellierung von Diffusion und Dispersion.
Diffusion und Dispersion sind Zufallsbewegungen. Diffusion ist Bewegung der einzelnen
Molek¨
ule in der L¨osung, Dispersion ist Bewegung der Tr¨opfchen. Obwohl das
physikalisch zwei v¨ollig verschiedene Vorg¨ange sind, bewirken beide eine Durchmischung
der L¨osung. Im Lauf der Zeit verteilt sich jeder gel¨oste Stoff gleichm¨aßig in der L¨osung.
Wir werden im Folgenden nur von Diffusion reden, aber Dispersion wird genauso
modelliert. Nur die Zahlenwerte der Parameter sind verschieden.
Abbildung 10.5 ist das Ergebnis einer Simulation einer Zufallsbewegung. 20 Teilchen
stehen anfangs in einem Rechteck in der Mitte der gezeigten Fl¨ache angeordnet und
unterliegen zuf¨alligen Bewegungen. Die Bewegung des einzelnen Teilchens ist rein
zuf¨allig, aber insgesamt ist das Ergebnis, daß sich die Teilchen mehr oder minder
gleichm¨aßig u
¨ber die Fl¨ache ausbreiten. Von den Regionen, wo viele Teilchen sitzen,
k¨onnen sich naturgem¨aß mehr Teilchen wegbewegen als von dort, wo nur wenige sitzen:
Daher ergibt sich mit der Zeit ein Ausgleich der Teilchenkonzentrationen.
Die Modellierung der Diffusion betrachtet nat¨
urlich nicht die einzelnen Molek¨
ule oder
Tr¨opfchen, sondern sie geht davon aus, daß die Gesamtwirkung der Diffusion vieler
Teilchen einen Teilchenstrom (Fluß) von den Gebieten hoher Konzentration zu den
Gebieten niederer Konzentration ergibt. Mit Φ(t, x) bezeichnen wir die Teilchenmenge,
die pro Fl¨acheneinheit und Sekunde einen Querschnitt durchsetzt. Der einfachste
Modellansatz ist das Ficksche Gesetz. Wir formulieren es nur f¨
ur den eindimensionalen
Fall.
Merksatz 10.7 (Ficksches Gesetz). Sei L(t, x) die Konzentration eines diffundierenden
Stoffes und Φ der Teilchenfluß auf Grund einer reinen Zufallsbewegung. Positives Φ
bedeutet Fluß in Richtung auf gr¨oßere x, negatives Φ bedeutet Fluß in Gegenrichtung.
Dann ist
∂
(10.5)
Φ(t, x) = −D L(t, x).
∂x
132
Abbildung 10.6. Die Wirkung der Diffusion
Durchgezogen: Konzentration zu Beginn. Strichliert: Konzentration nach einiger Zeit.
Die Pfeile symbolisieren St¨arke und Richtung des Teilchenflusses.
Dabei ist D eine Konstante, der Diffusionskoeffizient. Unter Umst¨anden kann D(x) auch
vom Ort abh¨angen.
Das Ficksche Gesetz gibt zumindest qualitativ das typische Verhalten der Diffusion
wieder: Wird die Konzentration in positive Richtung, also f¨
ur steigende x, gr¨oßer, so ist
∂
die partielle Ableitung ∂x
L positiv. Damit wird Φ negativ, was einem Fluß in negative
Richtung, also in Richtung auf kleinere x, entspricht. Je gr¨oßer das
Konzentrationsgef¨alle, also je gr¨oßer die Ableitung der Konzentration, desto st¨arker ist
auch der Fluß. Abbildung 10.6 gibt eine schematische Darstellung dieses Sachverhaltes.
Mit Methoden der Wahrscheinlichkeitsrechnung l¨aßt sich das Ficksche Gesetz fundierter
¨
begr¨
unden. Ubrigens
gilt dasselbe Gesetz auch f¨
ur die W¨armeausbreitung durch
W¨armeleitung in einem K¨orper.
Wir komplettieren jetzt unser Schadstoffmodell durch Einbau der Diffusion oder
Dispersion. Seien DS (x) und DO (x) die Diffusions- beziehungsweise
Dispersionskonstanten f¨
ur Schadstoff und Sauerstoff. Die Diffusionskonstanten k¨onnen
auf Grund der verschiedenen Beweglichkeit der Molek¨
ule verschieden sein. Die
Dispersionskonstante w¨are in jedem Falle dieselbe, denn beide Stoffe werden mit
denselben zuf¨allig schwirrenden Tr¨opfchen transportiert. Zus¨atzlich zum Transport durch
die Str¨omung erhalten wir auf Grund der Zufallsbewegung f¨
ur beide gel¨oste Stoffe die
folgenden Fl¨
usse:
∂
L(t, x),
∂x
∂
ΦO (t, x) = −DO (x) C(t, x).
∂x
ΦS (t, x) = −DS (x)
10. VERTEILTE PARAMETER
133
Positives Φ zeigt nach unserer Konvention in Richtung auf gr¨oßere x. Wir betrachten
wieder die “unendlich kleine” Zelle zwischen x und x + ∆x. Am oberen Zellrand, also an
der Stelle x, durchsetzt der Fluß ΦS (t, x) den Querschnitt A(x) in Richtung auf h¨ohere
x, also ins Innere der Zelle. Das ergibt einen Schadstoffzufluß von A(x)ΦS (x) kmol
Schadstoff pro Sekunde. Am unteren Ende, an der Stelle x + ∆x, durchsetzt der Fluß
ΦS (x + ∆x) den Querschnitt A(x + ∆x) in Richtung aus der Zelle heraus. Das ergibt
einen Schadstoffabfluß von A(x + ∆x)ΦS (x + ∆x). Insgesamt erhalten wir:
Schadstoffbilanz durch Diffusion:
A(x)ΦS (x) − A(x + ∆x)ΦS (x + ∆x)
∂
∂
= −A(x)DS (x) L(t, x) + A(x + ∆x)DS (x + ∆x) L(t, x + ∆x),
∂x
∂x
Sauerstoffbilanz durch Diffusion:
A(x)ΦO (x) − A(x + ∆x)ΦO (x + ∆x)
∂
∂
= −A(x)DO (x) C(t, x) + A(x + ∆x)DO (x + ∆x) C(t, x + ∆x).
∂x
∂x
Dieser Term muß in der Modellbildung zu den anderen Mengenbilanzen der Zelle addiert
werden. Im Grenz¨
ubergang zu unendlich kleinen Zellen haben wir die Mengenbilanz
durch ∆x dividiert und dann den Grenzwert f¨
ur ∆x → 0 ausgewertet. Der Beitrag des
Diffusionsterms ist dabei
∂
∂
A(x + ∆x)DS (x + ∆x) ∂x
L(t, x + ∆x) − A(x)DS (x) ∂x
L(t, x)
lim
∆x→0
∆x
∂
∂
=
A(x)DS (x) L(t, x) .
∂x
∂x
Dieser Term wird den anderen Termen der Advektions-Reaktions-Gleichung hinzugef¨
ugt.
Die Diffusion des Sauerstoffs wird ebenso modelliert.
Das Schadstoffmodell unter Einbeziehung der Diffusion besteht also aus den beiden
partiellen Differentialgleichungen
∂
∂
∂
∂
A(x) L(t, x) =
A(x)DS (x) L(t, x) −
A(x)u(x)L(t, x)
∂t
∂x
∂x
∂x
+s(t, x) − kS A(x)L(t, x),
A(x)
∂
∂
C(t, x) =
∂t
∂x
∂
∂
C(t, x) −
A(x)u(x)C(t, x)
∂x
∂x
−kS A(x)L(t, x) + kO b(x)(cS − C(t, x)).
A(x)DO (x)
In den neuen Gleichungen werden L und C zweimal nach x partiell differenziert:
Diffusion f¨
uhrt auf partielle Ableitungen zweiter Ordnung in den Raumvariablen. Um die
Gleichung zu behandeln, w¨
urden wir auch je eine weitere Randbedingung f¨
ur Schadstoff
und Sauerstoff brauchen, zum Beispiel die Konzentrationen am unteren Ende des
betrachteten Flußst¨
uckes. Durch die Diffusion kann ja ein kleiner Teil des Schadstoffes
wieder stromaufw¨arts wandern, sodaß auch vom unteren Ende her theoretisch
geringf¨
ugige Einfl¨
usse auf die Schadstoffkonzentration im betrachteten Flußst¨
uck zu
134
Tabelle 10.3. Reaktions-Diffusions-Advektionsgleichung
∂
A(x)L(t, x)
∂t
∂
∂x
=
∂
DS (x)A(x) ∂x
L(t, x)
∂
− ∂x
A(x)u(x)L(t, x)
+s(t, x)
−kS A(x)L(t, x)
¨
Zeitliche Anderung
partielle Ableitung nach t
Diffusion
partielle Ableitungen zweiter Ordnung nach x
Advektion
partielle Ableitung nach x von
Konzentration × Querschnitt
× Str¨omungsgeschwindigkeit
Beitrag der Quellen
Beitrag der chemischen Reaktion
Gestalt h¨angt von der Art der Reaktion ab
erwarten sind. Weil man diese Konzentrationen ja nicht weiß, kann man sich in der
Praxis behelfen, indem man ein ausreichend langes Flußst¨
uck berechnet und am Ende
die Schadstoffkonzentration mit 0 und die Sauerstoffkonzentration mit cS ansetzt. Dabei
gehen wir davon aus, daß weit stromabw¨arts von der Schadstoffquelle die Wasserqualit¨at
weitgehend regeneriert ist. Das Festlegen der Werte der gesuchten Funktionen am Rand
nennt man Dirichlet Randbedingungen. Ein anderer h¨aufiger Typ sind sogenannte
Neumann Randbedingungen, die die Ableitung der gesuchten Funktionen am Rand
vorschreiben. Wir k¨onnten zum Beispiel verlangen, dass die Ableitung der
Konzentrationen am unteren Ende null sind,
∂
∂
(10.6)
L(t, xL ) =
C(t, xL ) = 0,
∂x
∂x
was bedeutet, dass sich die Konzentrationen an dieser Stelle r¨aumlich nicht ¨andern sollen.
Tabelle 10.3 zeigt zur Wiederholung noch einmal die Gestalt einer der
Reaktions-Diffusions-Advektionsgleichungen.
Durch die Einf¨
uhrung der Diffusion ¨andert sich das qualitative Verhalten deutlich.
Stellen wir uns vor, daß die Fabrik zum Zeitpunkt 0 einen Augenblick lang eine gewisse
Schadstoffmenge in den Fluß l¨aßt, und dann ihren Abfluß wieder schließt. An der Stelle x
lassen wir einen Beobachter mit seinen Meßsonden am Ufer warten. Das
Advektionsmodell ohne Diffusion sagt voraus, daß der Schadstoff zur Zeit x/u
(Wegstrecke durch Str¨omungsgeschwindigkeit) die Stelle x erreicht. Bis zu diesem
Zeitpunkt nimmt der Beobachter nichts wahr. Zur Zeit x/u schlagen seine Instrumente
kurz aus, einen Augenblick sp¨ater ist der Schadstoff vorbeigetrieben und der Beobachter
nimmt nichts mehr wahr.
Diffusion und Dispersion bewirken, daß manche Schadstoffteilchen der großen Masse
vorauseilen, manche nachhinken. Das mathematische Modell w¨
urde voraussagen, daß der
Beobachter unmittelbar nach der Zeit 0 bereits einen minimalen Anstieg des Schadstoffes
10. VERTEILTE PARAMETER
135
wahrnehmen w¨
urde. Die Konzentration w¨
urde weiter ansteigen, bis die große Masse
vorbeigetrieben ist, und wieder langsam gegen Null absinken. Je weiter stromabw¨arts der
Beobachter steht, desto langsamer gehen Anstieg und Abfall der Konzentration vor sich.
Die Methode der Charakteristiken scheitert am Diffusionsmodell. Die Grundidee dieser
Methode war, daß ein Beobachter, der mit der Str¨omung mittreibt, die Str¨omungseffekte
nicht wahrnimmt: Er ist immer vom selben Wasser und denselben Schadstoffteilchen
umgeben. Die Durchmischung durch Diffusion und Dispersion macht aber gerade diesen
Effekt zunichte. Wenn man versucht, die totalen Ableitungen der beobachteten
Konzentrationen zu berechnen, heben sich nicht mehr alle partiellen Ableitungen weg.
Bemerkung. Wenn man bei der Modellbildung nur die Diffusion (oder Dispersion)
ber¨
ucksichtigt und alle anderen Effekte wegl¨asst, kommt man (f¨
ur konstanten
Querschnitt) zur klassischen Diffusionsgleichung der Gestalt
∂
∂
∂
L(t, x) =
D(x) L(t, x) , t > 0, x ∈ (0, xL ),
∂t
∂x
∂x
f¨
ur eine Eigenschaft L, die sich durch kleine Zufallsbewegungen ¨andert. Dies ist zum
Beispiel auch f¨
ur die W¨armeverteilung der Fall, man nennt diese Gleichung daher auch
W¨armeleitungsgleichung, wobei dann L die Temperatur ist. Um unter den L¨osungen der
partiellen Differentialgleichung eine eindeutig auszuw¨ahlen, m¨
ussen Anfangsbedingungen
f¨
ur t = t0 und Randbedingungen bei x = 0, x = xL gegeben sein.
10.9. Numerische L¨
osung von partiellen Differentialgleichungen.
Partielle Differentialgleichungen sind schwierig zu l¨osen, und ben¨otigen vor allem bei
zwei- und dreidimensionalen Problemen viel Speicherplatz und Rechenzeit. Wir haben
gesehen, daß sie ganz verschiedene Ph¨anomene beschreiben k¨onnen und ihre L¨osungen
qualitativ ganz unterschiedliches Verhalten zeigen. Jede partielle Differentialgleichung
hat sozusagen ihren ganz eigenen Charakter. Daher gibt es auch ganz verschiedene
Ans¨atze zur numerischen L¨osung, und Programmpakete werden oft f¨
ur spezielle
Anwendungen erstellt. So gibt es zum Beispiel Programme, die auf die Berechnung von
Str¨omungs- und Diffusionsvorg¨angen im Grundwasser spezialisiert sind, Programme zur
Berechnung von Spannungen und Dehnungen in elastischen Strukturen, Programme zur
Simulation der Luftstr¨omung um eine Tragfl¨ache, und viele andere.
Eines der bekanntesten und erfolgreichsten Verfahren ist die Methode der finiten
Elemente. Sie baut darauf auf, die partielle Differentialgleichung durch eine Art
Zellenmodell n¨aherungsweise darzustellen. Ein Element besteht aus einem Raumst¨
uck,
sozusagen einer Zelle, und einer Funktion auf dieser Zelle, die noch durch geeignete
Koeffizienten modifiziert werden kann. Das gesamte untersuchte Gebiet wird in Elemente
zerlegt. Als N¨aherungsl¨osungen werden nur Funktionen zugelassen, die auf jeder
einzelnen Zelle durch ein Element dargestellt werden k¨onnen. Auf diese Weise wird jede
zul¨assige N¨aherungsl¨osung durch endlich viele Koeffizienten beschrieben. Gegebenenfalls
k¨onnen noch zus¨atzliche Bedingungen die Auswahl der N¨aherungsl¨osung einschr¨anken.
136
Abbildung 10.7. Finite Elemente in einer Dimension
Wir erkl¨aren das an einem Beispiel f¨
ur ein eindimensionales Problem: Eine partielle
Differentialgleichung auf dem Intervall [0, l] soll gel¨ost werden. Als Element betrachten
wir ein Intervall mit einer Geradengleichung y = kx + d als Funktion. Zwei Koeffizienten,
n¨amlich k und d, stehen zur Anpassung zur Verf¨
ugung. Wir teilen das gesamte Gebiet,
also das Intervall [0, l], zum Beispiel in 20 kleine Teilintervalle. Als N¨aherungsl¨osungen
werden nur Funktionen zugelassen, die auf jedem der Teilintervalle Geradenst¨
ucke sind.
Daher wird jede N¨aherungsl¨osung durch 40 Koeffizienten, n¨amlich ki und di f¨
ur jede
einzelne Zelle (mit Nummer i = 1, · · · , 20), beschrieben. Als Zusatzbedingung fordern
wir, daß die gesamte N¨aherungsl¨osung stetig ist, also die einzelnen Geradenst¨
ucke sich
stetig treffen, wo die Zellen zusammenstoßen. Abbildung 10.7 zeigt ein einzelnes Element
und eine zul¨assige N¨aherungsl¨osung. Statt Geradenst¨
ucken k¨onnte man auch zum
Beispiel quadratische oder vor allem kubische Parabeln heranziehen (sogenannte
kubische Splines). Zweidimensionale Gebiete k¨onnen zum Beispiel in Rechtecke oder
Dreiecke zerlegt werden, dreidimensionale Gebiete in Quader, Tetraeder oder Prismen
(wie in Abbildung 10.3).
Die exakte L¨osung der partiellen Differentialgleichung ist nat¨
urlich normalerweise nicht
genau aus Elementen aufgebaut, aber wenn die Zellen ausreichend klein gew¨ahlt sind,
gibt es zul¨assige N¨aherungsl¨osungen, die sich der exakten L¨osung recht genau
anschmiegen. Unter den zul¨assigen N¨aherungsl¨osungen sucht man eine, die die partielle
Differentialgleichung wenigstens bis auf einen m¨oglichst kleinen Fehler genau erf¨
ullt. Das
f¨
uhrt rechnerisch auf die Bestimmung von endlich vielen Koeffizienten als L¨osung eines
Minimumproblems. Je kleiner die Zellen, desto genauer wird die N¨aherung, aber desto
mehr Koeffizienten sind zu bestimmen, und desto gr¨oßer wird auch der Rechenaufwand.
Die mathematische Durchf¨
uhrung dieses Konzeptes im Detail ist nat¨
urlich ziemlich
kompliziert.
Die Methode der finiten Elemente ist ein sogenanntes semidiskretes Verfahren, weil nur
die r¨aumliche Verteilung mit endlich vielen Elementen diskretisiert wird; dadurch
¨
entsteht f¨
ur eine partielle Differentialgleichung, in der auch zeitliche Anderungen
modelliert sind, ein großes System von gew¨ohnlichen Differentialgleichungen, in dem nur
mehr Ableitungen nach der Zeit vorkommen. Wie dieses System dann numerisch gel¨ost
wird, ist im Prinzip ein unabh¨angiges Problem, man sucht sich je nach Anforderungen
einen passenden ODE-solver.
10. VERTEILTE PARAMETER
137
Ausgefeilte Programmpakete suchen nicht nur die Koeffizienten der N¨aherungsl¨osung, sie
bestimmen auch automatisch eine Einteilung des Gebietes in Zellen, bieten geeignete
ODE-solver zur Auswahl an und stellen die L¨osung graphisch dar.
Anders als bei semidiskreten Verfahren erfolgt die Diskretisierung bei Verfahren der
finiten Differenzen vollst¨andig, indem die Zeit in Abstimmung mit den Raumkoordinaten
diskretisiert wird. Solche Verfahren kommen vor allem dann zur Anwendung, wenn es auf
¨
genaue Berechnung beinahe unstetiger Anderungen
der Zust¨ande ankommt (z.B.
Schockwellen). Durch die Anpassung der Diskretisierung an spezifische Eigenheiten des
Modells (etwa entlang von Charakteristiken) kann nicht nur die Genauigkeit der
N¨aherungsl¨osungen verbessert werden, sondern auch die Effizienz ihrer numerischen
Berechnung.
10.10. Lo
¨sungen des Reaktions-Advektions-Diffusions-Systems in Femlab.
Femlab ist ein Softwarepaket, das zur bequemen Behandlung von partiellen
Differentialgleichungen in Matlab verwendet wird (FEM=finite element method). Es
kann Systeme von PDGen mit bis zu 3 Raumkoordinaten bearbeiten. Es enth¨alt Module
f¨
ur spezielle Probleme, etwa W¨armeleitung, St¨omungen, Elastizit¨at, Elektromagnetismus,
etc. Wir verwenden f¨
ur unser Schadstoff-Modell ein allgemeines Modul f¨
ur eine lineare
Gleichung mit Zeitableitungen erster Ordnung und Raumableitungen erster und zweiter
Ordnung.
Wenn wir in Matlab den Befehl femlab eingeben, erhalten wir die aus etlichen Fenstern
und klickbaren Menus bestehende graphische Arbeitsumgebung von Femlab. Hier k¨onnen
wir den Modell-Typ w¨ahlen und die r¨aumliche Geometrie auf recht anschauliche Weise
eingeben (f¨
ur zwei oder drei Raumdimensionen mit Hilfe vorbereiteter Bausteine). Ein
Gitter-Generator diskretisiert das Raumgebiet mit w¨ahlbar Feinheit. Die Koeffizienten
der PDG und die gew¨
unschten Anfangs- und Randbedingungen werden durch Eintragen
von Matlab-Formeln in Formulare programmiert, die je nach PDG-Typ gestaltet sind. Es
steht eine Auswahl von finiten Elementen und ODE-solvern zur Verf¨
ugung. Nachdem
Alles vollst¨andig ausgew¨ahlt und konsistent eingegeben wurde, startet die Simulation auf
Knopfdruck. Zur Visualisierung der Ergebnisse werden verschiedene variierbare
Darstellungs-M¨oglichkeiten angeboten.
F¨
ur unserer Schadstoffmodell haben wir nur eine r¨aumliche Dimension, wir klicken 1D.
Dann legen wir fest, dass es 2 Zustandsgr¨oßen gibt, die hier als Spaltenvektor u mit den
Komponenten u1 und u2 bezeichnet werden. Wir w¨ahlen als Geometrie ein
zusammenh¨angendes Intervall, 0 < x < xL = 10000. Als Einheiten verwenden wir Meter,
Kilogramm, Sekunde, also ist der betrachtete Flußabschnitt zehn Kilometer lang. Beim
Ausf¨
ullen der Formulare f¨
ur die (immer passend angezeigten) Koeffizienten, legen wir
uns darauf fest, dass mit u1 durchgehend die Schadstoffkonzentration und mit u2
durchgehend die Sauerstoffkonzentration bezeichnet wird. Dann lassen sich die in den
Formularen vorgegebenen Bezeichnungen eindeutig auf unsere Bezeichnungen abbilden
(f¨
ur die Str¨omungsgeschwindigkeit weichen wir auf die Bezeichnung v aus). Es werden
138
auch Terme angeboten, die in unserem Modell garnicht vorkommen, deren Koeffizienten
setzen wir einfach 0.
Zum Vergleich mit den Zellen-Simulationen (Abschnitt 10.3) w¨ahlen wir die Parameter
wie dort, insbesondere zeit- und ortsunabh¨angig. Es ist zweckm¨aßig, die Werte der
ben¨otigten Konstanten in die Add/Edit Constants Tabelle einzutragen, sodass ihre
Namen dann zur Bildung der Koeffizienten in den Formularen verwendet werden k¨onnen.
Die Schadstoff-Zufuhr modellieren wir der Einfachheit halber hier nicht als
Randbedingung, sondern als kurzzeitigen Eintrag vom Ufer am oberen Ende. Dies
geschieht mit der im Modell mit s(t, x) bezeichneten Inhomogenit¨at: Wir w¨ahlen eine
Konstante L0 und setzen s(t, x) = L0 /(Av) f¨
ur 0 < x < 100, 0 < t < 360, und s(t, x) = 0
sonst. Um orts- und zeitabh¨angige Koeffizienten zu programmieren, m¨
ussen wir in den
Formeln die Variablen x und t verwenden. In die erste Zeile f¨
ur die Inhomogenit¨at
schreiben wir also die Matlab Formel f¨
ur dieses s(t, x), n¨amlich
L 0*(0<x&x<100)*max(0,sign(360-t))/A/v. In der zweiten Zeile steht kO*b*cS (dies
ist die Inhomogenit¨at auf der rechten Seite der C-Gleichung). Durch Probieren w¨ahlen
wir die Paramter L0 und DS , DO so, dass h¨
ubsche, mit Abschnitt 10.3 vergleichbare
Simulationsergebnisse entstehen. Dies f¨
uhrt uns auf die Wahl L0 = 0.2 und DS = 60,
DO = 20.
Als Anfangsbedingungen setzen wir wieder L(0, x) = 0, C(0, x) = cS . Als
Randbedingungen nehmen wir L(t, 0) = 0, C(t, 0) = cS bei x = 0 und (10.6) bei x = xL .
Wir legen den Simulationszeitraum 0 ≤ t ≤ 3600 fest, alles Andere lassen wir zun¨achst
wie voreingestellt (z.B. ode15s als ODE-solver). Dann dr¨
ucken wir gespannt auf die
Taste Solve und hoffen, dass nur wenige Fehlermeldungen kommen. Nach deren
Beseitigung l¨auft die Simulation dann endlich und zeigt uns die Ergebnisse zu den von
uns gew¨ahlten Zeitpunkten an, sagen wir t = 900, 1800, 2700, 3600, wie beim
Zellenmodell. Wir k¨onnen die Zeitpunkte aber auch dichter w¨ahlen, sodass wir dann mit
der Taste Animate einen Film sehen, der die zeitliche Entwicklung je einer der
Komponenten der L¨osung zeigt. Weitere graphische M¨oglichkeiten sind perspektivische
plots der Fl¨achen L(t, x), C(t, x) u
¨ber der (t, x)-Ebene, oder die F¨arbung einer
Komponente je nach der Gr¨oße der anderen an derselben Stelle (t, x).
Abbildung 10.10 zeigt die Ergebnisse f¨
ur die gerade geschilderten Parameter. Um glatte
Kurven zu erhalten, haben wir 120 Elemente (f¨
ur jede der beiden Komponenten)
gew¨ahlt. Ausserdem haben wir ausprobiert, mit welchem der w¨ahlbaren Verfahren f¨
ur die
L¨osung der linearen Gleichungssysteme (f¨
ur die Koeffizienten der N¨aherungsl¨osung) die
besten Rechenzeiten erzielt werden: Diese betragen auf einem Standard PC insgesamt
etwa 6 Sekunden.
Die Kurven mit den Maxima sind die Schadstoffkonzentrationen L(ti , x), jene mit den
Minima die zugeh¨origen Sauerstoffkonzentrationen C(ti , x). Die negativen Werte von
C(900, x) sind nat¨
urlich chemisch unsinnig und stammen daher, dass f¨
ur so grosse
¨
Schadstoffmengen das lineare Reaktionsmodell nicht zutreffend ist. Im Ubrigen
schauen
die Simulationsergebnisse aber vern¨
unftig aus und entsprechen unseren Erwartungen: Da
der Schadstoffzutritt nur f¨
ur kurze Zeit am linken Ende stattfindet, bildet sich ein
zun¨achst gut lokalisiertes L-Paket am linken Ende. Ohne Diffusion und Dispersion (kurz
10. VERTEILTE PARAMETER
Abbildung 10.8. Simulationsergebnisse
Diffusions-Modells in Femlab
des
139
Reaktions-Advektions-
−4
x 10
t=900
15
10
5
0
−4
x 10
t=1800
6
4
2
0
−4
x 10
t=2700
6
4
2
0
−4
x 10
t=3600
6
4
2
0
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
D&D) w¨are die L¨ange des Pakets 360 ∗ v Meter, bei Str¨omungsgeschwindigkeit v = 2 also
720m. Durch D&D verbreitert sich das Paket, seine H¨ohe nimmt dadurch ab.
Gleichzeitig verschiebt sich das Paket wegen der Str¨omung nach rechts. Zus¨atzlich nimmt
die L-Konzentration wegen der Oxidation ab, die auch C verbraucht. Das
durchh¨angende C-Paket verbreitert sich durch seine Kopplung an L und mit seiner
eigenen D&D Konstanten DO . Wegen DS = 60 > 20 = DO bleibt das Maximum der
L-Kurve st¨arker hinter der Str¨omung zur¨
uck als das Minimum der C-Kurve.
Die Verbreiterung der Pakete kann im Modell mit D&D im Gegensatz zu dem
Zellen-Modell (das nur Advektion modellierte) also plausibel gemacht werden. Trotzdem
r¨
uhrt die Verbreiterung der Pakete in den Simulationsergebnissen nicht nur aus den
D-Termen, sondern sie beruht auch auf systematischen Problemen mit finiten
Elementen: Wie im Zellenmodell kommt es zu einer Verschmierung von lokal
¨
konzentrierten Anderungen.
Dies k¨onnen wir zeigen, indem wir im selben Femlab-Modell
die Diffusionskoeffizienten einfach 0 setzen, DS = DO = 0. Dann m¨
ussten die Fronten des
Schadstoff-Pakets genau mit der Str¨omungsgeschwindigkeit wandern, also wegen der
speziellen Wahl von s(t, x) die vordere zum Zeitpunkt t genau die Stelle x = 100 + v ∗ t
passieren, und die hintere Front x = v ∗ (t − 360) (vgl. die exakte L¨osung der
Advektionsgleichung im Abschnitt u
¨ber Charakteristiken). Bei genauerer Untersuchung
140
der zugeh¨origen Simulationsergebnisse stellen wir jedoch fest, dass die Paket-L¨angen
gr¨osser als 720 sind und mit der Zeit zunehmen, selbst bei DS = DO = 0.
Es w¨are naiv, sich ein Simulations-Programm zu w¨
unschen, das alle PDGen zur vollen
Zufriedenheit im Griff hat, zu unterschiedlich sind die Eigenschaften der verschiedenen
Gleichungstypen. Legt man Wert auf beinahe unstetige Wellenfronten, sind angepasste
finite Differenzen Verfahren sicher genauer. Andererseits ist software wie Femlab f¨
ur
Anwender, die sich nicht mit der Numerik partieller Differentialgleichung auseinander
setzen wollen, und die m¨oglichst anschaulich ihr PDG-Modell programmieren wollen,
eine große Hilfe. Die Gefahr, falsche Simulationsergebnisse nicht als solche zu erkennen,
¨
besteht aber immer, und kann nur durch analytische Uberlegungen
und systematische
numerische Tests verringert werden.
11. DISKRETE DETERMINISTISCHE MODELLE
141
11. Diskrete deterministische Modelle
Angenommen, der Zustand eines Systems wird durch n Zustandsgr¨oßen beschrieben, die
in einem Zustandsvektor x = (x1 , . . . , xn ) zusammengefasst werden. Falls sich das System
mit der Zeit t kontinuierlich ver¨andert, und wir dies mathematisch beschreiben wollen,
brauchen wir Funktionen der Zeit x1 (t), . . . , xn (t) f¨
ur t ≥ t0 , mit einem Anfangszeitpunkt
t0 . Bisher haben wir nur Modelle betrachtet, deren Zustandsvektor x(t) koninuierlich von
der Zeit abh¨angt (im Falle r¨aumlich verteilter Zust¨ande sind f¨
ur jedes t die
Komponenten von x(t) selbst Funktionen, die den Raumkoordinaten Werte zuordnen).
Es kann aber vorteilhaft und in vielen Situationen durchaus angebracht sein, den
Zustand des Systems nur zu diskreten Zeitpunkten t0 < t1 < t2 < . . . zu modellieren. Ein
mathematisches Modell, das dies leistet, k¨onnte zum Beispiel von der folgenden Gestalt
sein:
x0 gegeben
(11.1)
xk+1 = F (xk ), k = 0, 1, 2, . . .
Dabei ist y = F (x) eine Funktion, die einem n-Vektor x einen n-Vektor y zuordnet und
xk ist der Zustandsvektor zum Zeitpunkt tk . Es ist denkbar (und in manchen F¨allen
unabdingbar), dass auf der rechten Seite von (11.1) die Zeit explizit vorkommt und dass
zur Berechnung von xk+1 auch fr¨
uhere Zust¨ande xk−1 , xk−2 , . . . ben¨otigt werden. In (11.1)
haben wir uns aber auf Modelle beschr¨ankt, in denen xk+1 aus dem Vorg¨anger-Zustand
allein berechenbar ist. Insbesondere in der Form xk+1 = xk + G(xk ) nennt man so ein
Modell eine Differenzengleichung erster Ordnung. Die Funktion G beschreibt die
¨
Anderung
des Zustands zwischen diskreten Zeitpunkten. Ausgehend vom Anfangszustand
x0 bestimmt das Modell sukzessive die Zust¨ande x1 , x2 , . . . zu den diskreten Zeitpunkten
t1 , t2 , . . . . Die Werte der Zustandsgr¨oßen sind im Allgemeinen nicht diskret, dh. die
Komponenten von xk sind im Allgemeinen keine ganzen Zahlen, sondern beliebig reell.
In diesem Kapitel besch¨aftigen wir uns mit Modellen, in denen der Zustand xk+1
eindeutig durch den Zustand xk bestimmt ist, das sind also deterministische Modelle.
Diese sind einer numerischen Simulation besonders zug¨anglich, man braucht nur F zu
programmieren und wiederholt anzuwenden. Es ist aber auch m¨oglich, eine qualitative
Analyse der durch x0 und F erzeugten Folge x1 , x2 , . . . durchzuf¨
uhren; es geht dabei um
Gleichgewichtslagen, deren Stabilit¨at, periodische L¨osungen, asymptotisches Verhalten
bei k → ∞, etc. Im Besonderen ist es chaotische Dynamik, die an Hand von einfachen
diskreten deterministischen Modellen studiert werden kann: das ber¨
uhmteste Beispiel
daf¨
ur ist die Differenzengleichung xk+1 = αxk (1 − xk ) mit der logistischen Funktion
F (x) = αx(1 − x) f¨
ur α > 1 und x ∈ R.
11.1. Ein diskretes deterministisches Epidemie-Modell.
Ein fr¨
uhes und wegweisendes Epidemie-Modell, das Modell von Kermack-McKendrick,
unterteilt eine Population von N Individuen in die Klasse x der gesunden Ansteckbaren,
die Klasse y der angesteckten Ansteckenden und die Klasse z der geheilten, nicht
ansteckenden Immunen (wie im SIR-Modell des ersten Semesters). Wir betrachten die
Population zu diskreten Zeitpunkten tk in gleichbleibenden Abst¨anden von zwei Wochen,
sagen wir. Mit xk bezeichnen wir die Anzahl der Ansteckenden zum Zeitpunkt tk , f¨
ur die
142
Anzahl der Individuen in den Klassen y, z schreiben wir yk , zk . Wir nehmen an, dass die
Gesamtgr¨oße der Population konstant ist, also
(11.2)
xk + yk + zk = N.
Als vereinfachende Idealisierung lassen wir zu, dass die Zahlen xk , yk , zk keine ganzen
Zahlen sind. Dies ist insbesondere bei großen Klassen zul¨assig, weil es dann nicht auf
einzelne Individuen ankommt. Vielmehr k¨onnen wir uns in diesem Fall die Werte von
xk , yk , zk als Anteile an der Gesamtbev¨olkerung vorstellen, die in geeigneten Einheiten
angegeben sind.
Sei 0 < p < 1 die gleichbleibende Wahrscheinlichkeit, dass ein ansteckbares Individuum
effektiven Kontakt mit einem ansteckenden Individuum hat, dh. dass es w¨ahrend eines
Beobachtungsintervalls angesteckt wird. Dann ist die Wahrscheinlichkeit, dass es keinen
effektiven Kontakt mit einem der ansteckenden Individuen hat, gleich 1 − p, und die
Wahrscheinlichkeit, dass es den effektiven Kontakt mit allen yk Ansteckenden vermeidet
ist gleich (1 − p)yk . Wir setzen daher
xk+1 = (1 − p)yk xk .
Dies ist f¨
ur sich allein als Modell f¨
ur xk+1 nicht brauchbar, weil ja yk ben¨otigt wird. Zur
Abk¨
urzung f¨
uhren wir den Parameter a ein, und zwar a = −ln(1 − p) > 0. Dann ist
e−a = 1 − p. Der Anteil der w¨ahrend des Intervalls neu Infizierten, n¨amlich (1 − e−ayk )xk
ist zum Zeitpunkt tk+1 in der Klasse der Angesteckten y. Aufgrund der Heilung verlassen
ein Teil der Angesteckten y und kommen zur Klasse z. W¨ahrend [tk , tk+1 ] betreffe dies
(1 − b)yk Individuen, wobei wir einen konstanten Parameter 0 < b < 1 ansetzen. Damit
erhalten wir
xk+1 = e−ayk xk
yk+1 = (1 − e−ayk )xk + byk
(11.3)
zk+1 = zk + (1 − b)yk
Die Klasse z ist strikt an die Klasse y gekoppelt und spielt f¨
ur die weitere Entwicklung
ufen wir die Eigenschaft (11.2)
der Epidemie keine Rolle mehr. Zur Kontrolle u
¨berpr¨
indem wir die drei Zeilen summieren: xk+1 + yk+1 + zk+1 = xk + yk + zk , die Summe
bleibt im Modell unver¨andert.
Das Modell enth¨alt zwei positive reelle Zahlen, a und b, als Parameter. Obwohl bei der
Erstellung des Modells mit Wahrscheinlichkeiten argumentiert wurde, ist es kein
probabilistisches Modell. Der Parameter a ist fest vorgegeben, das Modell bestimmt
xk+1 , yk+1 , zk+1 eindeutig aus xk , yk , zk und ist daher deterministisch.
In einer interaktiven Programmiersprache wie Matlab kann, ausgehend von Startwerten
f¨
ur x, y, z im Vektor u=[u(1),u(2),u(3)], durch wiederholte Exekution der Zeile
u = [exp(−a ∗ u(2)) ∗ u(1), (1 − exp(−a ∗ u(2))) ∗ u(1) + b ∗ u(2), u(3) + (1 − b) ∗ u(2)]
sehr einfach eine Simulation ausgef¨
uhrt werden. Wir wollen dies hier nicht weiter
verfolgen, stellen aber eine qualitative Frage: Wird die Zahl der Angesteckten
zunehmen ? Mit anderen Worten: Sagt das Modell den Ausbruch einer Epidemie voraus ?
11. DISKRETE DETERMINISTISCHE MODELLE
143
Man kann die Frage auch so formulieren: wenn y0 = 1 (dies ist durch Wahl der Einheiten
immer erreichbar), wird dann y1 gr¨oßer als 1 sein ? Laut Modell gilt f¨
ur die Differenz
y1 − y0 = (1 − e−a )x0 + b − 1.
Dies ist gr¨oßer als 0 genau dann, wenn
x0 > (1 − b)/(1 − e−a ) = x∗ .
Der kritische Wert x∗ ist ein Schwellenwert f¨
ur die Gr¨oße der Klasse der Ansteckbaren.
Die Signifikanz eines solchen Schwellenwerts (threshold) ist unabh¨angig von
mathematischen Modellen aus Felddaten bekannt, und wird als Test f¨
ur die G¨
ute des
Modells verwendet. Wir sehen, dass der Zusammenhang zwischen den Parametern und
dem Schwellenwert wenigstens plausibel ist: je kleiner (1 − b), desto mehr Ansteckende
u
¨berdauern den Beobachtungszeitraum und desto kleiner ist x∗ . Je gr¨oßer (1 − e−a ),
desto gr¨oßer ist die Wahrscheinlichkeit der Ansteckung p und desto kleiner ist x∗ . Ob
allerdings die Gr¨oßenordnungen passen und mit realistischen, m¨oglichst unabh¨angig zu
eruierenden Parametern erreichbar sind, ist eine andere Frage. Ein weiteres Kriterium,
das vom Modell erf¨
ullt wird, ist die Forderung nach einer wohlbestimmten Bedeutung
der Parameter: a und b sind anschaulich interpretierbar und k¨onnen als Systemgr¨oßen
von der Sanit¨atsverwaltung gemessen und beeinflusst werden. Sind die Parameter einmal
festgelegt, kann das Modell im Prinzip dazu verwendet werden, die St¨arke der Epidemie
vorherzusagen (bei allem Vorbehalt, der solch einfachen Modellen entgegenzubringen ist).
Zusammenhang mit dem kontinuierlichen SIR Modell. Bei der Modellbildung zu
Kermack-McKendrick und zum SIR Modell werden dieselben Argumente verwendet, also
¨
sollte es Uberg¨
ange von einem zum anderen geben. Wenn wir beim diskreten Modell
(11.3) beginnen, suchen wir nach kontinuierlichen Funktionen S(t), I(t), R(t) so, dass
xk = S(kh), yk = I(kh) und zk = R(kh), zumindest in guter N¨aherung f¨
ur kleine
Zeitintervalle h = tk+1 − tk .
Im Modell (11.3) machen wir den folgenden linearen Ansatz f¨
ur die Parameter a und b:
a = λh,
b = 1 − γh.
Dies ist plausibel, denn f¨
ur h → 0 sollten die Populationen gleich bleiben. Mit t = kh
sieht (11.3) dann so aus:
S(t + h) = e−λhI(t) S(t),
I(t + h) = (1 − e−λhI(t) )S(t) + (1 − γh)I(t),
R(t + h) = R(t) + (1 − (1 − γh))I(t),
also
S(t + h) − S(t)
I(t + h) − I(t)
=
=
=
R(t + h) − R(t) =
(e−λhI(t) − 1)S(t) = −λhI(t)S(t) + O(h2 ),
−γhI(t) + (1 − e−λhI(t) )S(t)
−γhI(t) + λhI(t)S(t) + O(h2 ),
γhI(t).
Dabei haben wir die Reihenentwicklung f¨
ur die Exponentialfunktion verwendet und f¨
ur
deren Terme h¨oherer Ordnung das O-symbol geschrieben (insbesondere gilt O(h2 )/h → 0
144
f¨
ur h → 0). Wenn wir durch h dividieren erhalten wir links Differenzenquotienten. Lassen
wir dann h gegen 0 gehen, verschwinden die O(h2 )-Terme und es bleibt:
d
S(t) = −λI(t)S(t),
dt
d
I(t) = λI(t)S(t) − γI(t),
dt
d
R(t) = γI(t).
dt
Dies ist das von f¨
uher bekannte SIR-Modell ohne Zuwanderung (β = 0) und ohne
Sterblichkeit (µ = 0). Zur Kontrolle u
ufen wir die Eigenschaft (11.2) indem wir die
¨berpr¨
d
drei Zeilen summieren: dt (S(t) + I(t) + R(t)) = 0, die Summe bleibt im Modell
unver¨andert.
Wegen des engen Zusammenhangs von gew¨ohnlichen Differentialgleichungen und
Differenzengleichungen, hat der/die Modellierer/in unter gegebenen Umst¨anden die
M¨oglichkeit, den Modell-Typ zu w¨ahlen, je nach Fragestellung, Vorlieben und
Verf¨
ugbarkeit von Methoden.
11.2. Zellular-Automaten.
Wenn wir einen Vorgang in Zeit und Raum diskret beschreiben m¨ochten, dann m¨
ussen
wir auch das betrachtete Raumgebiet in endlich viele Teile zerlegen, oder, was auf
dasselbe hinausl¨auft, nur endlich viele seiner Punkte ber¨
ucksichtigen. Betrachten wir
zun¨achst den eindimensionalen Fall, also ein Raumintervall [0, L] der L¨ange L. In diesem
Intervall w¨ahlen wir N + 1 ¨aquidistante St¨
utzstellen, n¨amlich
xn = nL/N, n = 1, . . . , N − 1 im Inneren und die beiden R¨ander x0 = 0, xN = L.
Viele dynamische Prozesse sind durch lokale Regeln bestimmt, dh. dass nur r¨aumlich
benachbarte Teile des Systems einander beeinflussen (dies ist eine enorm generelle
Aussage und sollte u
¨berdacht werden). Denken wir zum Beispiel an die Str¨omung einer
Fl¨
ussigkeit in einem Kanal [0, L] mit Str¨omungsgeschwindigkeit v von 0 nach L. Die
Fl¨
ussigkeit besitze eine zeitlich und r¨aumlich ver¨anderliche Eigenschaft C(t, x) (z.B. die
Konzentration einer gel¨osten Substanz, wir lassen hier aber chemische Reaktionen oder
Diffusion weg und betrachten nur die reine Advektion). Wir geben eine diskrete
raum-zeitliche Beschreibung der Entwicklung von C, indem wir den Wert Cnk von C zu
jedem Zeitpunkt tk = k∆t an jeder Stelle xn = n∆x angeben. Stimmen wir ∆t und ∆x
aufeinander ab, und zwar so, dass ∆x = v∆t. Dann ist der Wert von C zur Zeit tk+1 an
der Stelle xn genau der Wert von C zur Zeit tk an der Stelle xn−1 . Die Entwicklung der
diskreten Werte von C folgt also folgender lokalen Regel:
k
Cnk+1 = Cn−1
,
k = 0, 1, . . . ,
n = 1, . . . , N.
ussen vorgegeben sein.
Die Anfangswerte Cn0 zur Zeit t0 sowie die Randwerte C0k bei x0 m¨
Dann ist durch die lokale Regel der Wert von C f¨
ur alle Zeiten tk an allen Stellen xn
eindeutig bestimmt. Das Zeit-Raum-Gitter (tk , xn ) zusammen mit den Randbedingungen
und den Regeln zur Berechnung der Zust¨ande Cnk+1 (Update) aus den Zust¨anden Cnk
bezeichnet man als Zellular-Automaten (engl. cellular automaton, dt. auch zellul¨arer
Automat). Die r¨aumliche Diskretisierung (die wir gerade mit Hilfe von St¨
utzstellen
11. DISKRETE DETERMINISTISCHE MODELLE
145
beschrieben haben) denkt man sich dabei als Unterteilung des Raumgebiets in Zellen
und Cnk als Mittelwerte von C w¨ahrend [tk , tk+1 ) in [xn−1 , xn ).
Es gibt ZA, in denen auch die Zust¨ande nur eine endliche Anzahl diskreter Werte
annehmen; wenn es nur zwei m¨ogliche Zust¨ande pro Zelle gibt, spricht man von
boolschen Automaten. Ein ZA ist gepr¨agt durch die Regeln f¨
ur das Update, wobei der
Zustand einer Zelle sich typischerweise aus den Zust¨anden ihrer Nachbarzellen ergibt
(und zwar eindeutig im deterministischen Fall). Dazu wird insbesondere auch festgelegt,
welche Zellen zur Nachbarschaft geh¨oren. In unserem Str¨omungsbeispiel ben¨otigten wir
nur die links liegende, westliche Nachbarzelle. Andere eindimensionale Prozesse
ben¨otigen die Ber¨
ucksichtigung beider Nachbarzellen (westlich und o¨stlich). Denkbar und
durchf¨
uhrbar sind auch gr¨ossere Nachbarschaften, etwa je zwei westliche und ¨ostliche
Zellen. Im Fall von zwei-oder dreidimensionalen r¨aumlichen Gittern gibt es noch mehr
geometrische M¨oglichkeiten f¨
ur Nachbarschaften.
Die lokalen Regeln zum Update k¨onnen oft eindeutig und leicht verst¨andlich formuliert
¨
werden, auch ohne mathematische Formalisierung. Die Ubersetzung
in eine globale
Funktion F wie in (11.1) ist vergleichsweise viel komplizierter, bringt wenig Erkenntnis
und ist f¨
ur eine numerisch-graphische Implementierung am Computer gar nicht
notwendig. Ein deterministischer ZA gilt als vollst¨andig gegeben, wenn zu beliebigen
Anfangsbedingungen geh¨orende Simulationsergebnisse eindeutig festgelegt sind.
Der bekannteste ZA ist zweifellos Game of Life, der 1970 von Conway erfunden und
seither intensiv numerisch-experimentell aber auch analytisch untersucht wurde. Es
handelt sich um eim stilisiertes r¨aumlich zweidimensional verteiltes Populations-Modell,
¨
¨
in dem der Tod auf lokale Unter- und Uberbev¨
olkerung folgt. Uberleben
und Geburt
setzen sehr eingeschr¨ankte lokale Bev¨olkerungsdichten voraus. Betrachten wir ein
Rechteck, das in viele kleine gleich große rechteckige Zellen unterteilt ist. Jede der Zellen
kann zu jedem Zeitpunkt tk , k = 0, 1, . . . entweder tot oder lebendig sein. In einer
graphischen Darstellung des Zustands zu einem festen Zeitpunkt werden tote Zellen leer
gelassen, lebendige eingef¨arbt oder mit einem Symbol besetzt. Ausgehend von einer
gegebenen Ausgangsverteilung wird der Zustand einer Zelle zur Zeit tk+1 aus ihrem und
ihrer Nachbarn Zustand zur Zeit tk bestimmt. Die Regel lautet:
eine Zelle ist zur Zeit tk+1 genau dann lebendig, wenn zur Zeit tk
sie selbst und zwei ihrer Nachbarn oder drei ihrer Nachbarzellen lebendig waren.
Zu den Nachbarzellen einer Zelle z¨ahlen die 8 n¨achstliegenden Zellen. Also: Hat eine
lebendige Zelle zu wenige oder zu viele lebendige Nachbarn, stirbt sie. Eine lebendige
Zelle mit 2 oder 3 lebendigen Nachbarn u
¨berlebt. Eine tote Zelle mit 3 lebendigen
Nachbarn wird lebendig.
Zeichnen Sie auf kariertem Papier einige Beispiele oder eine Sequenz von Schritten. Ein
Beispiel ist in Abbildung 11.2 zu sehen, der sogenannte glider, eine Konfiguration, die
sich in 4 Schritten um je eine Zelle nach S¨
uden und Westen verschiebt (¨
uberpr¨
ufen Sie,
ob die Regel richtig angewendet wurde). Durch Drehungen um 90◦ erhalten wir glider,
die sich nach NW, NO oder SO bewegen. Es gibt auch einfache periodische und statische
(sich nicht ver¨andernde) Konfigurationen.
146
Abbildung 11.1. Ein glider bewegt sich in vier Zyklen um einen Schritt
nach S¨
ud-West
Da wir mit einem endlich großen Rechteck arbeiten, gibt es Randzellen, die zun¨achst
weniger als 8 Nachbarzellen haben. Eine m¨ogliche Vervollst¨andigung besteht in
periodischen Randbedingungen: f¨
ur das Update der Zellen in der obersten Zeile
verwendet man die unterste Zeile als n¨ordlichen Nachbarn, f¨
ur die rechteste Spalte wird
die linkeste Spalte als ¨ostlicher Nachbar verwendet, und analog f¨
ur die unterste Zeile und
die linkeste Spalte (durch dieses Zusammenkleben der R¨ander macht man aus dem
Rechteck einen Schlauchring). Ein glider, der den s¨
udlichen Rand u
¨berschreitet kommt
dann im Norden wieder herein, usw.
Wenn zum Beispiel ein glider in eine statische/periodische Konfiguration eindringt
k¨onnen die erstaunlichsten Entwicklungen ausgel¨ost werden.... Viel besser als durch eine
Beschreibung lernen Sie Game of Life durch Experimentieren kennen: siehe Matlab
>>demo. Nat¨
urlich finden Sie mit einer Suchmaschine im Internet jede Menge sites mit
Game of Life, etwa www.bitstorm.org/gameoflife. Nur als Andeutung zeigen wir in
Abbildung 11.2 eine Momentaufnahme eines Game of Life mit 50 mal 50 inneren Zellen.
Die Kreise sind die lebendigen Zellen des Anfangszustands, die Sterne sind die
lebendigen Zellen nach 100 Zyklen.
Die Faszination, die Game of Life auf viele seiner Betrachter und Erforscher aus¨
ubt,
beruht wohl darauf, dass eine einfache, in einem Satz formulierbare Regel, derart
komplexe und reichhaltige Entwicklungsszenarien in sich birgt. So gibt es z.B. shooter,
das sind periodische, ortsfeste Konfigurationen, die alle 30 Zyklen einen glider ausstossen.
Daraf aufbauend wurde gezeigt, dass es sich reproduzierende Konfigurationen gibt: ist
das Universum gen¨
ugend groß, dann entstehen aus so einer Konfiguration zwei identische
solche. Nat¨
urlich wurden als Varianten zu Game of Life andere Regeln untersucht und
Regeln systematisch klassifiziert oder mit Parametern stufenlos ge¨andert. Bei der ersten,
legend¨aren Tagung zu diesen Themen wurde daf¨
ur der Begriff “artificial life” eingef¨
uhrt.
Bemerkung. In diesem Abschnitt des Skriptums hat sich - unangek¨
undigt - eine
Erweiterung des Begriffs “Mathematisches Modell” ereignet. W¨ahrend wir bisher unter
einem mathematischen Modell eine vereinfachende, idealisierende Abbildung eines
Teilbereichs der Realit¨at verstanden, erheben die Zellular Automaten des artificial life
nicht den Anspruch, in der Natur existierende Systeme abzubilden oder m¨oglichst gut
nachzuahmen. Ausserdem haben wir das Modell gar nicht vollst¨andig mathematisch
ausformuliert: f¨
ur Simulationen gen¨
ugt die programmierbare Beschreibung eines
11. DISKRETE DETERMINISTISCHE MODELLE
147
Abbildung 11.2. o lebendige Ausgangszellen, * lebendige Zellen nach 100 Updates
0
5
10
15
20
25
30
35
40
45
50
0
5
10
15
20
25
30
35
40
45
50
Algorithmus und die graphische Darstellung der Simulation. Vielleicht sollte man solche
Modelle einfach Simulations-Modelle nennen.
Die M¨oglichkeit, mit solchen Modellen zu arbeiten, er¨offnete sich der Wissenschaft seit es
leistungsstarke Computer gibt. Zus¨atzlich zu Theorie und Experiment mit nat¨
urlich
existierenden Objekten gibt es nun Computer-Simulationen in Welten unserer
sprachlichen Phantasie. Dies ist auch f¨
ur sehr grundlegende Fragen von großer
wissenschaftlicher Bedeutung. So wird an Hand von Simulationen untersucht, wie in
einer stilisierten Ursuppe chemischer Bestandteile Makromolek¨
ule als Vorl¨aufer des
Lebens entstehen k¨onnen. Prozesse, die der biologischen Evolution ¨ahnlich sind, finden
auch in recht einfachen virtuellen, digitalen Universen statt. Soziologische und kulturelle
Entwicklungen werden mit vereinfachten Gesellschaften simuliert (multi agent models).
Scheinbar (?) komplexe und komplizierte Ph¨anomene wie Lernen, Anpassung oder
Selbstorganisation k¨onnen mit Simulations-Modellen immitiert werden, und dies kann
wesentlich zum Verst¨andnis der Grundlagen solcher Prozesse beitragen.
148
12. Diskrete probabilistische Modelle
Probabilistische Modelle - immer h¨aufiger auch stochastische Modelle genannt beinhalten Elemente des Zufalls, in der einen oder anderen Form. Dies kann aus vielerlei
Gr¨
unden angebracht sein. Zum Einen kann es sich um die Modellierung von Systemen
handeln, in denen der Zufall explizit eine grosse Rolle spielt. Dies ist etwa bei
¨
Gl¨
ucksspielen der Fall, der Ubertragung
von Krankheitserregern oder der Mutation
genetischer Information. Zum Anderen kann es sich um Modelle handeln, die die
Beschreibung deterministischer Details zugunsten der Ber¨
ucksichtigung mittlerer Effekte
aufgibt und diese wie vom Zufall gesteuert behandelt. Dies ist etwa in der Gasdynamik
der Fall oder bei der Ber¨
ucksichtigung unregelm¨aßiger Schwankungen von Umweltgr¨oßen
wie Niederschlag, Temperatur etc.
In probabilistischen Modellen sind die Zustandsgr¨oßen nur gem¨aß einer gewissen
Wahrscheinlichkeitsverteilung bestimmt. Bei Simulationen werden diese Gr¨oßen mittels
Zufallsgeneratoren festgelegt, sodass auch die Simulationsergebnisse vom Zufall
abh¨angen und von Simulation zu Simulation verschieden sind. Im Rahmen von
sogenannten Monte-Carlo Simulationen f¨
uhrt man eine Serie von Einzelsimulationen
durch und macht statistische Aussagen u
¨ber die Ergebnisse.
Es gibt auch stochastische Modelle, in denen zeitlich kontinuierliche Modellgr¨oßen vom
Zufall abh¨angen und ein Kontinuum von Werten annehmen k¨onnen. In dieser Vorlesung
werden aber keine kontinuierlichen Zufallsprozesse behandelt: Unsere Modellgr¨oßen
werden nur endlich viele diskrete Werte annehmen, und wir werden nur zeitlich diskrete
¨
Modelle behandeln. Der Ubergang
zwischen den Zust¨anden wird mit
Wahrscheinlichkeiten modelliert. Das Paradebeispiel f¨
ur ein System, das so modelliert
werden kann, ist ein Spielw¨
urfel, der vor bzw. nach jedem Wurf in einem von 6 m¨oglichen
¨
Zust¨anden ist. F¨
ur einen perfekten W¨
urfel ist die Ubergangswahrscheinlichkeit
zwischen
1
allen Zust¨anden gleich, n¨amlich 6 .
12.1. Markov-Ketten.
Eine Markov-Kette oder Markov-Prozess n-ter Ordnung ist ein zeitlich diskretes
dynamisches Modell bestehend aus einer Menge {S1 , . . . , Sn } von Zust¨anden und n2
¨
Ubergangswahrscheinlichkeiten
pij , i = 1, . . . , n, j = 1, . . . , n. Der Prozess kann zu jedem
Zeitpunkt k nur in einem der n Zust¨ande sein. Ist der Prozess zum Zeitpunkt k im
Zustand Si , dann ist pij die Wahrscheinlichkeit, dass der Prozess zum Zeitpunkt k + 1 im
¨
Zustand Sj ist. Im einfachsten Fall h¨angen die Ubergangswahrscheinlichkeiten
nicht von
der Zeit k oder vom Zustand zur Zeit k ab (lineare MK mit konstanten
¨
Ubergangswahrscheinlichkeiten).
Der Prozess startet in einem gegebenen
Anfangszustand.
Beispiel. Wir wollen das Wetter in einer bestimmtem Stadt mit Hilfe dreier m¨oglicher
Zust¨ande S, W und R, n¨amlich Sonne, Wolken, Regen beschreiben. Je einer dieser
Zust¨ande soll f¨
ur jeden Tag k zutreffen. In Abbildung 12.1 sind die (f¨
ur diese Stadt
¨
typischen) Ubergangswahrscheinlichkeiten mit Hilfe von Pfeilen und dazugeschriebenen
Wahrscheinlichkeiten angegeben. Dh. zB. falls es heute wolkig ist, dann gibt es eine
12. DISKRETE PROBABILISTISCHE MODELLE
149
¨
Abbildung 12.1. Wettermodell mit Ubergangswahrscheinlichkeiten
Chance von 25%, dass es morgen wieder wolkig ist, eine Chance von 25%, dass es morgen
regnet und eine Chance von 50%, dass es morgen sonnig ist. Wir k¨onnen die
¨
Ubergangswahrscheinlichkeiten
zu einer Tabelle zusammenfassen:
S
W R
S 1/2 1/2 0
W 1/2 1/4 1/4
R 0 1/2 1/2
Zu einem heutigen Wetter finden wir in der entsprechenden Zeile die
Wahrscheinlichkeiten f¨
ur das morgige. Die Tabelle muss nicht symmetrisch sein: z.B.
W → R besitzt die Wkt 1/4, aber R → W besitzt die Wkt 1/2. Ausserdem setzt das
Modell die Wahrscheinlichkeiten, dass Sonne und Regen aufeinander folgen, gleich Null.
Beispiel. Zwei Spieler A und B haben insgesamt eine Anzahl von n − 1 Jetons. In jeder
Runde k des Spiels gewinnt A mit der Wkt p einen Jeton von B, oder aber A verliert
einen Jeton an B, n¨amlich mit der Wkt q = 1 − p (dh. in jeder Runde wechselt genau ein
Jeton den Besitzer). Sobald einer der beiden Spieler alle Jetons hat, kann sich nichts
mehr ¨andern. Betrachten wir ein Spiel mit 4 Jetons: es gibt 5 Zust¨ande, dir wir einfach
mit der Anzahl der Jetons bezeichnen, die in Besitz von A sind. Es gibt also die
Zust¨ande 0,1,2,3 und 4. Das Pfeile-Diagramm sieht wie in Abbildung 12.1 aus, die
zugeh¨orige Tabelle ist
0
1
2
3
4
0
1
q
0
0
0
1
0
0
q
0
0
2
0
p
0
q
0
3
0
0
p
0
0
4
0
0
0
p
1
Es ist klar, wie das gleiche Modell f¨
ur eine andere Anzahl von Jetons ausssieht.
150
¨
Abbildung 12.2. Spielermodell mit Ubergangswahrscheinlichkeiten
Allgemein geh¨ort zu einer Markov Kette n-ter
¨
Ubergangswahrscheinlichkeiten

p11 p12
 p21 p22
P =
..
 ...
.
pn1 pn2
Ordnung eine n × n-Matrix P von
···
···
···

p1n
p2n 
.. 
. 
pnn
Alle Elemente der Matrix sind reelle Zahlen zwischen 0 und 1, 0 ≤ pij ≤ 1. Ferner sind
alle Zeilensummen gleich 1, denn jeder Zustand muss in irgendeinen der m¨oglichen
Zust¨ande u
¨bergehen. Eine quadratische Matrix mit diesen beiden Eigenschaften heisst
eine stochastische Matrix. Und ein Vektor x ∈ Rn heisst Wahrscheinlichkeitsvektor, wenn
alle seine Komponenten nichtnegativ kleiner gleich 1 sind und ihre Summe 1 ist. Wenn
der Zeilenvektor x ein Wahrscheinlichkeitsvektor ist und P eine stochastische Matrix,
dann ist das Vektor mal Matrix Produkt xP wieder ein Wahrscheinlichkeitsvektor:
n
n
n
(xP )j =
j=1
n
n
xi (P )ij =
j=1 i=1
n
xi (P )ij =
i=1 j=1
n
xi
i=1
n
pij =
j=1
xi = 1.
i=1
Daher sind stochastische Matrizen die nat¨
urlichen linearen Transformationen zwischen
Wahrscheinlichkeitsvektoren. Wir haben, wie bei Markov Ketten u
¨blich, die Vektoren als
Zeilen verwendet und multiplizieren diese von rechts mit quadratischen Matrizen
(genausogut k¨onnten wir die f¨
ur MK un¨
ubliche Schreibweise mit Spaltenvektoren und
den transponierten Matrizen verwenden). F¨
ur das Element einer Matrix A in der i-ten
Zeile und in der j-ten Spalte schreiben wir wie u
¨blich (A)ij .
¨
Wenn ein Markov Prozess mit Ubergangsmatrix
(P )ij = pij im Zustand Si startet, dann
ist die Wkt, dass nach einem Schritt der Zustand Sj angenommen wird gleich pij . Dh.
die Wkten sind die Komponenten des Zeilenvektores (pi1 · · · pin ). Dies ist ein
Wahrscheinlichkeitsvektor, die i-te Zeile von P , und entsteht aus P durch Multiplikation
mit dem Vektor
ei = (0 · · · 010 · · · 0)
(2)
von links; in ei steht der Einser an der i-ten Stelle. Wie gross ist die Wkt pij , dass Si in
zwei Schritten nach Sj u
¨bergeht ? Falls dazwischen sicher Sk l¨age, w¨are die Wkt
(2)
pij = pkj . Aber die Wkt, dass Sk dazwischen liegt ist nur pik , und wir m¨
ussen die
12. DISKRETE PROBABILISTISCHE MODELLE
151
Summe u
¨ber alle m¨oglichen Wege bilden,
n
(2)
pij
pik pkj = (P 2 )ij .
=
k=1
2
Hier bezeichnet P das Quadrat der Matrix P , das dadurch entsteht, dass man die n × n
Matrix P mit sich selbst multipliziert, P 2 := P · P (dies ist nicht die Matrix, die durch
¨
Quadrieren der einzelnen Elemente entsteht). Also ist P 2 die 2-Schritt Ubergangsmatrix,
2
deren Element (P )ij die Wkt angibt, dass Si in zwei Schritten nach Sj u
¨bergeht. Wenn
¨
wir diese Uberlegung
induktiv fortsetzen, sehen wir dass die Matrix P m , die dadurch
¨
entsteht, dass man P m mal mit sich selbst multipliziert, die m-Schritt Ubergangsmatrix
m
ist: (P )ij ist die Wkt, in m Schritten von Si nach Sj zu kommen.
¨
Anstatt die m-Schritt Ubergangsmatrix
auszurechnen, k¨onnen wir auch einen iterativen
Prozess als diskretes dynamisches System formulieren: sei x(k) ein Zeilenvektor, dessen
j-te Komponente die Wkt angibt, dass im kten Schritt der Zustand Sj vorliegt. Dann
ergibt sich der Wktvektor x(k + 1) durch Multiplikation mit P von rechts:
(12.1)
x(k + 1) = x(k)P,
k = 0, 1, 2, . . .
Zu gegebenen Anfangsvektor (zB. x(0) = ei , wenn der Anfangszustand sicher Si ist) legt
dieses diskrete dynamische System die Wahrscheinlichkeitsverteilung x(k) f¨
ur alle Zeiten
k fest. Achtung, x(k) ist nicht der k-te Zustand des Markov Prozesses, m¨ogliche
Zust¨ande sind nur S1 , . . . , Sn . x(k) enth¨alt aber die Wahrscheinlichkeiten, dass die MK
zur Zeit k einen der Zust¨ande annimmt. Zur Simulation einer Zustandsentwicklung
m¨
usste man in jedem Schritt einen Zustand gem¨aß den Wkten in x(k) ausw¨ahlen (siehe
Monte Carlo Simulation im n¨achsten Abschnitt).
F¨
ur eine andere Interpretation von x(k), die vielleicht etwas anschaulicher ist, stellen wir
uns eine große Anzahl N gleichartiger Systeme vor, die alle mit dem gleichen Modell
beschrieben werden. Sei zB. N eine Anzahl von St¨adten, f¨
ur die wir das gleiche
Wetter-Modell anwenden. Jede Stadt wird mit der gleichen MK modelliert, wobei der
Markov Prozess in jeder Stadt unabh¨angig von den anderen St¨adten abl¨auft. Auch wenn
alle St¨adte im gleichen Anfangszustand starten, werden sie sich im Lauf der Zeit
unterschiedlich entwickeln. Das Modell liefert Erwartungswerte f¨
ur eine Verteilung der
Zust¨ande zu sp¨ateren Zeitpunkten: zum Zeitpunkt k sind x(k)j N St¨adte im Zustand Sj
zu erwarten. Obwohl das Modell keine verl¨assliche Aussage u
¨ber den Zustand in einer
einzelnen Stadt macht, sagt es f¨
ur die ganze Menge von St¨adten eine Verteilung voraus.
Von diesem Standpunkt aus gesehen sind MK mit unseren fr¨
uher diskutierten Modellen
verwandt: Etwa in den Epidemie Modellen oder den Modellen zur chemischen Kinetik
bewegten sich auch Gruppen von Individuen oder Mole k¨
ulen von einer Klasse in eine
andere.
¨
Abgesehen von der Berechnung der m-Schritt Ubergangswahrscheinlichkeiten
in P m oder
der Iteration (12.1) ist vor Allem die qualitative Analyse von Markov Ketten von
Interesse. Folgende Fragen liegen auf der Hand und k¨onnen f¨
ur lineare Markov Ketten
¨
mit konstanten Ubergangswahrscheinlichkeiten beantwortet werden: Langzeitverhalten,
Grenzzustand, mittlere Anzahl von Schritten, um spezielle Zust¨ande zu erreichen, Wkt,
152
dass ein spezieller Zustand je erreicht wird, etc. Ein wichtiges Resultat in diese Richtung
formulieren wir als Satz:
¨
Satz 12.1. Sei P die Ubergangsmatrix
einer regul¨aren Markov Kette, das heisst, dass f¨
ur
ein m > 0 alle Elemente der Matrix P m positiv sind. Dann gilt
a) Es gibt einen eindeutig bestimmten Wahrscheinlichkeitsvektor pˆ mit der
Eigenschaft
pˆP = pˆ.
b) F¨
ur beliebigen Anfangszustand x(0) = ei gilt
lim ei P m = pˆ,
m→∞
insbesondere ist die Grenzverteilung pˆ unabh¨angig vom Anfangszustand.
c)


lim P m = Pˆ := 

m→∞
pˆ
pˆ
..
.


.

pˆ
Der Beweis dieses Satzes ist eine Anwendung der Linearen Algebra, der Vektor pˆ ist ein
linker Eigenvektor von P zum Eigenwert 1. Die Aussagen b) und c) sind Aussagen u
¨ber
das Langzeitverhalten der Kette. Wenn wir uns wieder vorstellen, dass N Ketten
gestartet werden, dann konvergiert nach vielen Schritten die Anzahl der Ketten, die sich
in Zustand Sj befinden, nach pˆj N , unabh¨angig von der Anfangsverteilung der Zust¨ande.
Beispiel. F¨
ur die obige Wetter Matrix P sind alle Elemente von P 2 positiv, wie Sie
leicht nachrechnen k¨onnen. Also handelt es sich um eine regul¨are MK. Durch einfache
interaktive Iterationen zB. in Matlab kommt man schnell zu der Vermutung, dass


0.4 0.4 0.2
Pˆ =  0.4 0.4 0.2  .
0.4 0.4 0.2
¨
Uberzeugen
Sie sich, dass pˆ = (0.4 0.4 0.2) tats¨achlich ein linker Eigenvektor von P
zum Eigenwert 1 ist. Die Aussage des Satzes bedeutet folgende Langzeit-Statistik laut
Modell: 40% der Tage sonnig, 40% wolkig, 20% regnerisch.
Wir gehen hier nicht auf andere qualitative Fragestellungen in Bezug auf Markov Ketten
ein, ebensowenig auf Markov Ketten mit unendlich vielen diskreten Zust¨anden.
Die beiden einfachen Beispiele dieses Abschnitts machen vielleicht den Eindruck, dass
Markov Ketten triviale mathematische Modelle f¨
ur Entwicklungen sind, die auch ohne
mathematische Modellierung vorhergesagt werden k¨onnen. Doch gibt es sehr wohl MK
¨
Modelle mit konstanten Ubergangswahrscheinlichkeiten,
die wesentliche Beitr¨age zum
Verst¨andnis von Zufallsprozessen geleistet haben: Jahrelang spekulierten Experten der
Gentechnologie, warum die einem Bakterium zus¨atzlich zum Typ A implantierten
Plasmide vom Typ B im Lauf der Generationen wieder separiert werden. Plasmide
tragen zus¨atzliche Erbinformation ausserhalb des Chromosoms des Bakteriums. Bei der
nat¨
urlichen Zellteilung werden die N vorhandenen Plasmide zun¨achst kopiert und dann
12. DISKRETE PROBABILISTISCHE MODELLE
153
zuf¨allig auf die beiden Tochterzellen verteilt, und zwar erh¨alt jede Tochter N Plasmide.
Im Lauf von vielen Zellteilungen entstehen aus einer Population, in der jedes der
Bakterien gleichviele Plasmide vom Typ A wie B besitzt, zwei Unterpopulationen von
Bakterien, die jeweils nur einen Typ Plasmide besitzen. Als man schließlich die
Zellteilungen als Markov Kette modellierte, wobei die Wahrscheinlichkeiten streng
kombinatorisch hergeleitet wurden, stellte sich heraus, dass die Plasmid Separation ein
Prozess ist, der sich von selbst einstellt, allein auf Grund binomialer
Wahrscheinlichkeitsverteilungen (random genetic drift). Irgendwelche Spekulationen
bzgl. Plasmid-Inkompatibilit¨aten u. dgl. haben sich mit diesem Modell er¨
ubrigt. Das
Modell w¨
urde gut in diese Vorlesung passen und ist nicht besonders aufw¨andig, wir
betrachten nun aber eine andere, nichtlineare Markov Kette mit zustandsabh¨angigen
¨
Ubergangswahrscheinlichkeiten.
12.2. Infektion einer Familie.
Wir betrachten nun wieder eine Population ansteckbarer (S), ansteckender (I) und
geheilter (R) Individuen, wie in Abschnitt 11.1. Die Reihenfolge von Infektion und
Heilung ist durch S → I → R gegeben. Hier wollen wir nun eine kleine Population mit
ganzzahligen Klassen modellieren und die Zufallsprozesse bis ins Detail verfolgen. Seien
k = 0, 1, . . . die Zeitpunkte, zu denen das System beobachtet wird, der Abstand zwischen
zwei Beobachtungen soll gleichbleibend lang sein, sagen wir 2 Wochen. Mit Sk , Ik , Rk
bezeichnen wir die Anzahl der Individuen in den drei Klassen zum Zeitpunkt k. Wieder
nehmen wir an, dass die Gr¨oße der Gesamtpopulation konstant N bleibt:
Sk + Ik + Rk = N.
Die m¨oglichen Zust¨ande des Modells sind also alle jene Vektoren mit nichtnegativen
ganzzahligen Komponenten, deren Summe N ist. Zum Beispiel f¨
ur N = 4 gibt es 10
solche Zust¨ande. Die Modellbildung besteht darin, Formeln f¨
ur die Berechnung der
¨
Ubergangswahrscheinlichkeiten
zwischen diesen Zust¨anden zu entwickeln.
Sei p die Wkt des effektiven Kontakts eines der S mit einem der I w¨ahrend eines
Beobachtungszeitraums. Dann ist 1 − p die Wkt, dass ein S den effektiven Kontakt mit
einem der I vermeidet. Wenn es Ik Infizierte gibt, dann ist also die Wkt, dass ein S den
Kontakt mit allen Infizierten vermeidet gleich
qk = (1 − p)Ik ,
und qkn die Wkt, dass n Ansteckbare den Kontakt vermeiden. Ferner ist (1 − qk ) die Wkt,
dass ein S effektiven Kontakt mit einem der I hat und somit ist
(1 − qk )Sk −n
die Wkt, dass Sk − n effektiven Kontakt haben. Wenn Sk und Ik gegeben ist, kann die
¨
Wahrscheinlichkeit, dass beim Ubergang
k → k + 1 gerade n der Ansteckbaren jeden
Kontakt mit Infizierten vermeiden und gerade Sk − n aber effektiven Kontakt haben als
das Produkt
qkn (1 − qk )Sk −n
154
berechnet werden. Aus der Kombinatorik ist bekannt, wieviele M¨oglichkeiten es gibt, aus
Sk Indiviuen jeweils n auszuw¨ahlen. Diese Anzahl ist gegeben durch den
Binomialkoeffizienten
Sk !
Sk
=
,
n
n!(Sk − n)!
wobei f¨
ur jede nichtnegative ganze Zahl K, K Fakult¨at gegeben ist durch
K! = K · (K − 1) · (K − 2) . . . 3 · 2 · 1. Somit ist die Wkt f¨
ur Sk+1 = n gleich
Sk
n
qkn (1 − qk )Sk −n .
Der Einfachheit halber nehmen wir nun an, dass die Infektionsdauer in jedem Fall mit
Sicherheit gleich dem Beobachtungszeitraum ist. Das heisst, jedes infizierte Individuum
ist genau 2 Wochen ansteckend und wechselt dann in die Klasse der nicht mehr
ansteckenden und nicht mehr ansteckbaren Geheilten. Zusammenfassend erhalten wir
Sk
(12.2)
Pr(Sk+1 = n|Sk und Ik ) =
qkn (1 − qk )Sk −n
n
Falls Sk+1 = n, dann Ik+1 = Sk − n und Rk+1 = Rk + Ik .
Hier haben wir eine Schreibweise aus der Wahrscheinlichkeitsrechnung verwendet: mit
der Probabilit¨at Pr(A|B) bezeichnet man die Wahrscheinlichkeit von A unter der
Bedingung B, dh. die Wahrscheinlichkeit f¨
ur das Eintreten von A unter der
Voraussetzung, dass B erf¨
ullt ist. Die Ausdr¨
ucke auf der rechten Seite von (12.2) sind
nur sinnvoll, wenn n ≤ Sk . Andere F¨alle brauchen auch nicht modelliert zu werden, weil
die Anzahl an Ansteckbaren nur abnehmen oder h¨ochstens gleichbleiben kann.
Diese Gleichungen heissen das Reed-Frost Modell und geben die
¨
Ubergangswahrscheinlichkeiten
f¨
ur eine Markov Kette. Die diskreten Zust¨ande, die das
Modell annehmen kann sind Vektoren mit drei nichtnegativen ganzzahligen
Komponenten, deren Summe N ist. Zu jedem Zeitpunkt k ist der Zustand so ein Vektor;
¨
zB. f¨
ur N = 4 gibt es 10 m¨ogliche Zust¨ande. Die Ubergangswahrscheinlichkeiten
sind
hier aber keine Konstanten, sondern h¨angen stark nichtlinear vom Zustand ab, und zwar
von Sk und, u
¨ber qk , von Ik . Eine qualitative Analyse des Reed-Frost Modells ist daher
nicht m¨oglich. Wohl aber k¨onnen Beipiele der Entwicklung simuliert werden und die
Auswertung vieler Beispiele f¨
uhrt zu statistischen Aussagen u
¨ber das bzw. mit Hilfe des
Modells.
Monte Carlo Simulationen. Wegen der Annahme, dass Infizierte nach zwei Wochen
geheilt sind, ist es sicher, dass die Infektion der Familie nach sp¨atestens 2N Wochen
vor¨
uber ist. Ist n¨amlich einmal Ik = 0, dann kann im Modell kein Krankheitsfall mehr
auftauchen. Wie lange wird die Infektion aber wirklich dauern, wieviele
Familienmitglieder werden krank werden und wieviele werden die Sache gesund
u
¨berstehen? Schauen wir uns ein Beispiel an, das mit Reed-Frost Wahrscheinlichkeiten
arbeitet.
Die Familie hat N = 5 Mitglieder und zum Zeitpunkt k = 0 ist genau eines davon
infiziert, also S0 = 4, I0 = 1, R0 = 0. Die Krankheit sei wenig ansteckend, sagen wir
p = 0.05. Mit Hilfe der Formel (12.2) berechnen wir die Wahrscheinlichkeiten f¨
ur alle
12. DISKRETE PROBABILISTISCHE MODELLE
155
m¨oglichen Nachfolgezust¨ande (unm¨oglich ist zB. S1 = 5 oder R1 = 1) und tragen diese in
eine Tabelle ein:
Beobachtung Wahrscheinlichkeit
S1 I1 R1 dieser Beobachtung
a 4 0
1
0.814
b 3 1
1
0.172
c 2 2
1
0.0135
d 1 3
1
0.000495
e 0 4
1
0.000005
Da ja irgendeiner der m¨oglichen Nachfolgezust¨ande angenommen wird, muss die Summe
der Wahrscheinlichkeiten 1 sein, was tats¨achlich der Fall ist. Zur Veranschaulichung
zeichnen wir eine (nicht ganz maßstabgetreue) den Wahrscheinlichkeiten entsprechende
Unterteilung des Intervalls [0,1], und benennen die Teilintervalle wie die Zeilen der
Tabelle: Nun ziehen wir eine Zufallszahl z ∈ [0, 1], das heisst wir verwenden einen
a
b
c
d
e
1
0.999995
0.9995
0.986
0.814
0
Mechanismus oder eine Methode, zum Beispiel einen Zufallszahlengenerator, der unter
den Zahlen in [0,1] eine ausw¨ahlt, wobei alle der Zahlen die gleiche Chance haben,
ausgew¨ahlt zu werden. L¨age z zB. im Intervall c, dann w¨are S1 = 2 ausgew¨ahlt. Durch
diese Vorgangsweise wird erreicht, dass genau eine der Zeilen ausgew¨ahlt wird und zwar
mit der Wahrscheinlichkeit, die der letzten Spalte der Tabelle entspricht. (Wir k¨onnten
die Teilintervalle auch in anderer Reihenfolge in der Zeichnung anordnen, wichtig ist nur,
dass ihre L¨angen gleich den Wahrscheinlichkeiten sind, dass alle vorkommen, und dass
[0,1] l¨
uckenlos u
¨berdeckt ist.) Wegen der geringen Ansteckungswahrscheinlichkeit p ist es
bei Weitem am wahrscheinlichsten, dass I1 = 0, womit die Infektion schon nach dem dem
ersten Beobachtungszeitraum beendet ist.
Angenommen wir ziehen z = 0.96, unser z liegt also im Teilintervall b, und das bedeutet
die Auswahl
S1 = 3 und I1 = 1.
Dies war der erste Simulationsschritt und die Simulation kann nun analog fortgesetzt
werden: wir berechnen die Wahrscheinlichkeiten Pr(S2 = n|S1 = 3 und I1 = 1) f¨
ur n ≤ 3
nach (12.2):
Beobachtung Wahrscheinlichkeit
S2 I2 R2 dieser Beobachtung
a 3 0
2
0.857
b 2 1
2
0.135
2
0.0075
c 1 2
d 0 3
2
0.0005
Die dazupassende Unterteilung von [0,1] in Teilintervalle der L¨ange a, b, c, d hat die
Trennungspunkte bei 0.875, 0.992, 0.9995. Wir ziehen eine Zufallszahl, sagen wir z = 0.9.
156
Damit ist b ausgew¨ahlt, also
S2 = 2 und I2 = 1.
Von hier ausgehend finden wir folgende Wahrscheinlichkeiten
Pr(S3 = n|S2 = 2 und I2 = 1) f¨
ur n ≤ 2:
Beobachtung Wahrscheinlichkeit
S3 I3 R3 dieser Beobachtung
a 2 0
3
0.902
3
0.095
b 1 1
c 0 2
3
0.003
Wir ziehen wieder eine Zufallszahl z ∈ [0, 1], sagen wir z = 0.25, und erhalten damit
S3 = 2
und
I3 = 0.
Wegen I3 = 0 ist die Simulation nun beendet. Die Infektion dauerte 6 Wochen, 3
Mitglieder der Familie sind erkrankt, 2 haben die Ansteckungsgefahr gesund
u
¨berstanden.
Dies war eine Monte Carlo Simulation. Klarerweise arbeitet man am Computer und
schreibt sich ein passendes Programm zur Berechnung der Wahrscheinlichkeiten. Es ist
effizienter, in jedem Schritt die Zufallszahl z zuerst zu bestimmen (wenig
Rechenaufwand) und die Wahrscheinlichkeiten nur soweit zu berechnen, bis erstmals z
kleiner ist als die Summe der Wahrscheinlichkeiten.
Man kann nun viele Beispiele rechnen – gleichwertig damit ist die Denkweise, dass man
viele Familien mit dem gleichen Modell modelliert – und dann eine Statistik der
Ergebnisse ziehen. Auf diese Weise k¨onnen wir etwa die durchschnittliche Dauer der
Krankheit, oder die zu erwartende Anzahl an Infizierten bestimmen. Zur Illustration
betrachten wir die Krankheit in Familien mit 5 Mitgliedern unter drei verschiedenen
Annahmen u
¨ber die Ansteckungswahrscheinlichkeit: p = 0.1, p = 0.3 und p = 0.5. In
jedem der drei F¨alle f¨
uhren wir 40 Simulationen durch, die alle mit S0 = 4, I0 = 1
starten. Wir bewerten die St¨arke der Infektion einer Familie mit Hilfe von S∞ , das ist die
Anzahl von (noch) Ansteckbaren zu dem Zeitpunkt k, zu dem die Krankheit aus der
Familie verschwindet, also Ik = 0 wird (dies ist sp¨atestens nach 5
Beobachtungszeitr¨aumen der Fall). In obigem ausf¨
uhrlichen Beispiel mit p = 0.05 hatten
wir S∞ = 2. Ergebnisse solcher Simulationen sind in folgender Tabelle aufgelistet, wobei
die Anzahl von Beipielen angegeben ist, die mit S∞ endeten.
S∞ p=0.1 p=0.3 p=0.5
0
1
12
23
1
1
8
12
5
3
1
2
3
11
3
1
22
14
3
4
Die Spaltensumme ist jeweils 40, das ist die Anzahl der gerechneten Beispiele (Familien)
f¨
ur jedes der p. F¨
ur p = 0.1 endete nur eine der 40 Simulationen ohne verbliebene
Ansteckbare, aber in 22 Beispielen wurde kein weiters Familienmitglied angesteckt. F¨
ur
p = 0.5 wurden in 23 der Beispiele alle Familienmitglieder angesteckt und in nur 3
Beispielen war die Gefahr nach dem ersten Beobachtungszeitraum vorbei. Interessant ist
12. DISKRETE PROBABILISTISCHE MODELLE
157
¨
die Anderung
der Verteilung auf die f¨
unf Zeilen, wenn die Ansteckungswahrscheinlichkeit
zunimmt: das Maximum ist f¨
ur p = 0.1 eindeutig bei S∞ = 4, das heisst, dass keine
weitere Erkrankung in einer Familie zu erwarten ist, die mit einem Infizierten startet.
F¨
ur p = 0.5 ist das Maximum bei S∞ = 0, das heisst, dass am wahrscheinlichsten die
ganze Familie erkrankt. Aber f¨
ur p = 0.3 macht das Modell keine eindeutige Aussage: es
gibt zwei Maxima und die Anzahl der 12 F¨alle, in denen alle erkranken ist etwa gleich
groß wie die Anzahl der 14 F¨alle, in denen nur ein Mitglied erkrankt. Gibt es in der N¨ahe
von p = 0.3 vielleicht so etwas wie einen Schwellwert ? Dies konnten wir im
deterministischen Modell in 11.1 analytisch beantworten, im probabilistischen Modell
m¨
ussen wir uns mit Simulationen und Statistik begn¨
ugen.
Implementierung in Matlab. Es ist nicht schwierig, das Reed-Frost Modell in
Matlab umzusetzen und damit Simulations-Reihen auszuf¨
uhren. Der Befehl
Y=rand(m,n) erzeugt eine m×n Matrix Y deren Elemente im Intervall (0,1) gleichverteilt
sind (rand steht f¨
ur random, es gibt auch die Funktion randn f¨
ur Normalverteilung). F¨
ur
unsere Zwecke gen¨
ugt es rand immer wieder ohne Argument aufzurufen, jedesmal wird
ein Skalar erzeugt und die Folge der so erzeugten Skalare ist de facto gleichverteilt im
Intervall (0,1).
Die folgende function im file ReedFrost.m f¨
uhrt f¨
ur die Argumente p, S0, I0,
num of fam in einer Schleife num of fam Monte Carlo Simulationen durch und
retourniert S∞ .
function S_00=ReedFrost(p,S0,I0,num_of_fam)
% Monte Carlo Simulations for the Reed-Frost model
S_00=zeros(1,S0+1);
for i=1:num_of_fam
S=S0; I=I0;
for k=1:S0+I0
q=(1-p)^I;
for n=0:S
Pr=nchoosek(S,n)*q^n*(1-q)^(S-n);
a(S+1-n)=Pr;
end
b(1)=a(1);
for j=2:S+1
b(j)=b(j-1)+a(j);
end
z=rand;
for j=S+1:-1:1
if z < b(j)
I=j-1;
end
end
S=S-I;
end
158
S_00(S+1)=S_00(S+1)+1;
end
Die Binomialkoeffizienten werden mit Hilfe der Matlab function nchoosek berechnet, die
Wahrscheinlichkeiten im Vektor a gespeichert und die Unterteilungspunkte von [0,1]
stehen in b. Die Zufallszahlen werden in der Zeile z=rand generiert, erst nach der
Berechnung aller Wahrscheinlichkeiten: da die 120 Simulationen f¨
ur die S∞ Tabelle nur
0.3 Sekunden Rechenzeit verbrauchen, k¨
ummern wir uns hier nicht um Effizienz. Wenn
man dieses Programm im command window je drei mal f¨
ur p=.1, p=.3, p=.5 aufruft,
erh¨alt man zB. folgende Ergebnisse:
>> S_00
S_00 =
1
>> S_00
S_00 =
0
>> S_00
S_00 =
0
>> S_00
S_00 =
10
>> S_00
S_00 =
10
>> S_00
S_00 =
12
>> S_00
S_00 =
30
>> S_00
S_00 =
30
>> S_00
S_00 =
27
>>
= ReedFrost(.1,4,1,40)
0
3
9
27
= ReedFrost(.1,4,1,40)
1
2
11
26
= ReedFrost(.1,4,1,40)
2
3
8
27
= ReedFrost(.3,4,1,40)
12
5
3
10
= ReedFrost(.3,4,1,40)
9
0
6
15
= ReedFrost(.3,4,1,40)
7
5
7
9
= ReedFrost(.5,4,1,40)
3
1
1
5
= ReedFrost(.5,4,1,40)
5
2
2
1
= ReedFrost(.5,4,1,40)
9
1
1
2
Bemerkung. Da der output eines numerischen Programms eindeutig vorbestimmt ist,
ist es prinzipiell unm¨oglich, einen echten Zufallsgenerator mit Hilfe deterministischer
mathematischer Operationen zu programmieren. Es gibt zwar software, die auf
physikalisches weisses Rauschen im Computer zugreift, f¨
ur viele Simulationszwecke
gen¨
ugen aber (im Prinzip vorhersagbare) Generatoren wie rand in Matlab. Um immer
wieder andere Zahlenfolgen zu bekommen, kann man den Zufallszahlengenerator am
Beginn einer jeden Folge immer wieder anders initialisieren, etwa mit Hilfe von Datum
und Uhrzeit.
12. DISKRETE PROBABILISTISCHE MODELLE
159
12.3. Ein probabilistischer Zellular-Automat.
Als Beispiel f¨
ur ein r¨aumlich verteiltes System, dessen Elemente zuf¨alligen
Entwicklungen ausgesetzt sind, denken wir an einen Waldbrand. Ein einfaches aber
qualitativ reichhaltiges Modell dazu ist der ZA Forest Fire. Das Rechteck aller Zellen
stellt einen Wald dar und jede der Zellen kann zu den diskreten Zeitpunkten in einem
von drei m¨oglichen Zust¨anden sein:
T : Baum (tree)
F : brennender Baum (fire)
E: leerer Platz (empty site)
Diese Zust¨ande k¨onnen nur in der Reihenfolge T → F → E → T → . . . durchlaufen
werden, wobei der Zustand F immer nur einen Zeitschritt lang anh¨alt. Die Verweildauer
einer Zelle in E und T h¨angt u
¨ber vorgegebene Wahrscheinlichkeiten p und f vom Zufall
ab. Dies sind die vier lokalen Regeln, die den ZA steuern:
• Jeder brennende Baum wird zu einem leeren Platz.
• Jeder leere Platz wird mit der Wahrscheinlichkeit p zu einem Baum.
• Befindet sich kein brennender Baum in der Nachbarschaft, dann wird ein Baum
mit der Wahrscheinlichkeit f zu einem brennenden Baum (Blitzschlag).
• Befindet sich ein brennender Baum in seiner Nachbarschaft, dann wird ein
Baum zu einem brennenden Baum (Feuerfront).
Dazu m¨
ussen wir noch festlegen, was wir unter einer Nachbarschaft verstehen. Anders als
bei Game of Life z¨ahlen wir hier nur jene 4 Zellen zur Nachbarschaft einer Zelle, die an
die Seiten angrenzen; jene mit nur einem gemeinsamen Eckpunkt z¨ahlen wir nicht dazu.
Als Randbedingung betrachten wir die Situation, dass das Rechteck von leeren Pl¨atzen
umgeben ist, auf denen keine B¨aume wachsen k¨onnen.
Die folgende Matlab function realisiert den ZA Forest Fire:
function waldbrand(T,F,p,f,maxcount)
% Zellular-Automat Waldbrand: Beschreibung siehe forestfire.m oder Skriptum
[m,n]=size(T);
no=[m+2 1:m+1]; e=[2:n+2 1]; s=[2:m+2 1]; w=[n+2 1:n+1];
count=0;
while count < maxcount
count=count+1;
FE
FE
FE
FN
NFF
NFN
NFL
E
NT
=
=
=
=
=
=
=
=
=
[zeros(m,1) F zeros(m,1)];
[zeros(1,n+2); FE; zeros(1,n+2)];
FE(no,:)+FE(:,e)+FE(s,:)+FE(:,w); % +FE(no,e)+FE(s,e)+FE(s,w)+FE(no,w);
FE(2:m+1,2:n+1);
(T & FN);
~FN;
(T & NFN & (rand(m,n) <= f));
(~T & ~F);
(E & rand(m,n) <= p);
160
OT = T-NFF-NFL;
T = OT+NT;
F = NFF+NFL;
X = T-F; % -1=F, 0=E, 1=T
image(X+2);
axis image off;
drawnow
end
Die function erh¨alt als Argumente die Ausgangsmatrizen T und F, die ausser Nullen nur
Einser an den entsprechenden Pl¨atzen enth¨alt, ferner die Wahrscheinlichkeiten p und f ,
sowie die Anzahl der zu simulierenden Zyklen. Alle Operationen werden als
Matrix-Gleichungen geschrieben, dies ist wesentlich effizienter als Schleifen f¨
ur die
Indizes. Die Zeilen- und Spaltenzahl der Matrizen sind m,n. In den Vektoren no, e, s,
w werden die Indizes der n¨ordlichen und s¨
udlichen Zeilen jeder Zeile und der ¨ostlichen
und westlichen Spalten jeder Spalte einer m + 2, n + 2 Matrix gespeichert, und zwar
periodischen Randbedingungen entsprechend. Dies wird dann zur Berechnung der
Anzahl brennender Nachbarn eines jeden Platzes in FE verwendet, wobei FE die
Vergr¨oßerung von F mit leeren Pl¨atzen am Rand ist. Die mit einem Kommentarzeichen
% unsch¨adlich gemachte Summe +FE(no,e)+FE(s,e)+FE(s,w)+FE(no,w) w¨
urde auch
die Ecknachbarn dazuz¨ahlen. FN enth¨alt also die Anzahl brennender B¨aume in der
Nachbarschaft eines jeden der m × n Pl¨atze. Mit Hilfe des logischen und-Operators &
erh¨alt NFF Einser an den Pl¨atzen mit brennenden Nachbarn wo ausserdem ein Baum
steht, dieser beginnt nun zu brennen. Mit Hilfe der Negation˜enth¨alt NFN jene Pl¨atze, in
deren Nachbarschaft kein brennender Baum ist. Falls dort ein Baum steht wird dieser
mit Wahrscheinlichkeit f zu brennen beginnen, was in NFL mit Hilfe einer in jedem
Zeitschritt neu gebildelten Zufallsmatrix rand(m,n) und der Kleiner-Gleich-Relation <=
realisiert ist. Leere Pl¨atze E befinden sich genau dort wo kein Baum und kein brennender
Baum steht. Mit Hilfe einer zweiten Zufallsmatrix und der Wahrscheinlichkeit p wird
festgelegt, ob dort ein neuer Baum NT w¨achst. Alte u
¨berlebende B¨aume OT sind jene, in
deren Nachbarschaft kein brennender Baum ist und die nicht gerade vom Blitz getroffen
werden. Insgesamt gibt es zum n¨achsten Zeitpunkt die B¨aume OT+NT und die brennenden
B¨aume NFF+NFL. Der Rest ist Graphik.
Am besten lernen Sie die Dynamik dieses ZA durch eigene Experimente kennen, von der
Seite teaching propst.html k¨onnen Sie eine Matlab function forestfire.m
herunterladen. Als matten Abglanz der dynamischen Entwicklungen zeigt Abbildung
12.3 drei Momentaufnahmen, die mit diesem Programm erzeugt wurden. Zu sehen ist
jeweils ein Wald mit 100 × 100 Zellen, wobei die weissen Zellen leere Pl¨atze sind, dunkle
(rote) Zellen brennende B¨aume und hellgraue (gr¨
une) Zellen nicht brennende B¨aume
darstellen. Fig. a zeigt einen Ausgangszustand gleichverteilter T und F , wobei 10% der
Pl¨atze mit T besetzt sind, T=rand(m,n) <= 0.1, und 10% der Pl¨atze brennen, sofern sie
nicht schon von einem T besetzt sind, F=(rand(m,n) <= 0.1) & ˜ T. Wir w¨ahlen
p = 0.05 und f = 0, dh. es gibt zun¨achst keinen Blitzschlag. In diesem Fall sind drei
verschiedene Szenarien m¨oglich: Erstens, die Anfangsverteilung ist so, dass die F schon
nach wenigen Schritten verschwinden, ab dann f¨
ullt sich der Wald nur mehr mit gr¨
unen
Zellen. Zweitens, die Anfangsverteilung ist so, dass sich eine oder mehrere Feuerfronten
12. DISKRETE PROBABILISTISCHE MODELLE
161
Abbildung 12.3. Zust¨ande des Zellular-Automaten Waldbrand
a
b
c
d
ausbilden, und zwar dann und dort, wo der Baumbestand schon dicht genug daf¨
ur ist. In
diesem Szenario kommt es typischerweise zu großen Ringwellen, das sind Karo-f¨ormige
Ketten dunkler Zellen, die sich in Richtung Rand ausdehnen und eine Zone leerer Pl¨atze
hinter sich lassen, die sich nur langsam wieder mit gr¨
unen B¨aumen f¨
ullt. Fig. b zeigt die
Situation nach etwa 100 Zyklen, ausgehend von a. Zwei gr¨oßer werdende Ringwellen sind
in der Mitte zusammengewachsen und im Begriff das Rechteck u
¨ber den Rand zu
verlassen. Wenn die zur¨
uckbleibenden kleinen Br¨ande verl¨oschen, w¨achst nach der
Feuersbrunst wieder einen dichter werdender Baumbestand heran. F¨
ur p = 0.05 und der
10%igen Anfangsverteilung ist es aber, drittens, h¨aufig der Fall, dass sich nach der
großen Anfangswelle viele kleinere Br¨ande bilden und nie mehr ein feuerfreier Zustand
eintritt. Es entstehen wandernde Muster aus B¨aumen, Feuern und leeren Pl¨atzen,
¨ahnlich wie in Fig. c. Da sich diese Muster (insbesondere auch in b) aus einem
zufallsverteilten Zustand entwickeln, liegen Beispiele f¨
ur die Selbstorganisation eines
Systems vor. In ¨ahnlich gebauten deterministischen Automaten kommt es zur
Ausbildung periodischer Raum-Zeit-Strukturen (Spiralwellen), die hier aber durch
Zufallsprozesse gest¨ort sind. Zu Spiralwellen, also sich nicht geradlinig oder konzentrisch
fortpflanzenden Strukturen, kommt es dann, wenn jede Zelle angeregt und damit die
162
Nachbarn anregend sein kann (F ), oder ruhend (T ), wobei dazwischen ein refrakt¨arer
Erholungs-Zustand (E) angenommen werden muss, der nicht angeregt werden kann.
Fig. c zeigt eine Momentaufnahme f¨
ur p = 0.05 und f = 0.0001, die am Ende einer
l¨angeren Simulation gemacht wurde. In diesem Stadium entstehen und vergehen
Strukturen verschiedener Gr¨oßenordnungen. Obwohl sich die Verteilungen dauernd
¨andern spricht man von einem quasistation¨aren Zustand, weil das Gesamt-Muster
qualitativ unver¨andert bleibt. Im Unterschied zum Fall f = 0 drittens, scheint dieses
Stadium f¨
ur f = 0.0001 aber etwas unruhiger und unregelm¨aßiger, und ausserdem kann
das Feuer wegen anhaltender Blitzgefahr nie endg¨
ultig verschwinden. Erh¨ohen wir
schließlich die Blitzschlag-H¨aufigkeit auf f = 0.01, dann stellen sich im Lauf der Zeit
bewegte Verteilungen wie in Fig. d ein (wieder p = 0.05 ). Die Strukturen verschwinden,
das Nachwachsen der B¨aume ist zu sehr gest¨ort, die Gleichverteilung von leeren Pl¨atzen
verhindert die Ausbildung von Feuerfronten.
Wie in Abschnitt 11.2, beruht die Faszination an ZA des Typs Forest Fire nicht auf der
Anwendbarkeit auf reale Waldbr¨ande, sondern auf der M¨oglichkeit, Mechanismen der
Musterbildung zu erkunden, insbesondere ihre Abh¨angigkeit von stufenlos einstellbaren
Parametern. So spricht man bei der Ausbildung quasistation¨arer Strukturen (c und,
unter Umst¨anden, b) von selbstorganisierter Kritizit¨at. Und bei einer grunds¨atzlichen
¨
Anderung
des qualitativen Verhaltens im Lauf der Zeit spricht man von
Phasen¨
uberg¨angen: Erstarrung im Fall b, wenn alle F verschwinden, oder Verdampfung
im Fall d, wenn alle Strukturen aufgel¨ost werden. Lebendige strukturelle Vielfalt gibt es
nur in einem kleinen Bereich aufeinander abgestimmter Parameter, dort entstehen sie
allerdings von selbst, auch aus strukturlosen Ausgangsverteilungen.
13. NACHWORT
163
13. Nachwort
Dieses Skriptum sollte Ihnen einen Einblick in verschiedene Methoden und vor allem in
die Denkweise der mathematischen Modellierung geben. Wichtiger als die
Rechentechniken ist der Prozeß, reale Systeme durch mathematische Gleichungen
sinnvoll wiederzugeben. Wir sind im Eiltempo vom einfachen Anpassen weniger
Parameter an einen Datensatz bis zum Aufstellen von partiellen Differentialgleichungen
und komplexen Computersimulationsmodellen fortgeschritten. Sie haben viel gesehen.
Klar, daß das Ergebnis nicht sein kann, daß Sie jetzt das ganze Repertoire von Methoden
systematisch beherrschen, an dem Sie vorbeigekommen sind. Auch klar, daß es viele
andere Aspekte gibt, die ebenso wichtig sind wie die vorgef¨
uhrten, und die wir in der
kurzen Zeit nicht behandeln konnten. (Wir haben uns zum Beispiel nicht mit
Feedback-Stabilisierung, dynamischer Optimierung, oder mit kontinuierlichen
stochastischen Systemen auseinandergesetzt.)
Ganz bewußt habe ich oben das Wort “sinnvoll” statt “richtig” verwendet: Jedes Modell
ist unvollkommen, und ein und dasselbe System kann durch ganz verschiedene Modelle
treffend wiedergegeben werden. Mathematik ist eine Denkweise, die man der
Wirklichkeit u
ur diesen Vorgang
¨berwirft. Die Erfahrung von vielen Wissenschaftern hat f¨
viele Tricks und Standardans¨atze bereit, von denen Sie einige gesehen haben. Letztlich
ist aber bis auf den heutigen Tag die mathematische Beschreibung von “Wirklichkeit”
eine sch¨opferische Leistung, getragen vom Menschen und seiner Intuition, nicht vom
Formalismus.
Wie soll man Intuition lernen? Ich selbst weiß keinen Weg als durch Beispiele. Daher ist
dieses Skriptum keine Sammlung von Rezepten, sondern darauf aufgebaut, einige
Beispiele durchzuziehen. Lernen Sie bitte nicht auswendig. Versuchen Sie, dem
Gedankengang zu folgen. Wie werden die einzelnen Schritte begr¨
undet? Welche Schritte
u
ur windig? Wenn Sie ein Kapitel durchgearbeitet
¨berzeugen Sie, welche halten Sie f¨
haben: K¨onnen Sie das Modell nacherz¨ahlen, seinen Zweck, die Ableitung der
Gleichungen und Diagramme, die Schl¨
usse, die aus dem Modell gezogen werden? Welche
Voraussetzungen m¨
ussen gemacht werden, damit das Modell gilt? In welchen Situationen
k¨onnte man ¨ahnliche Modelle verwenden? Welche mathematischen Methoden werden
eingef¨
uhrt, wozu dienen Sie, und wie wendet man sie an? Was waren die wesentlichen
Aussagen des Kapitels? — Ich weiß, Lernen mit Kopf ist sehr viel Arbeit.
Gerade wenn man die komplexen o¨kologischen Vorg¨ange beschreiben will, trifft man
schnell auf sehr schwierige Probleme. Und nicht viel weiter st¨oßt man auf die Grenzen
dessen, was heute, und vielleicht u
¨berhaupt, einer wissenschaftlichen Behandlung
zug¨anglich ist. Aus dieser Einsicht resultiert eine skeptische Haltung gegen¨
uber der
angewandten Mathematik. Ich habe versucht, diese Skepsis nicht zu vertuschen. Der
modernen Ablehnung des rationalen Denkens schließe ich mich trotzdem nicht an. Wer
alle Mittel wegwirft, die er f¨
ur unvollkommen h¨alt, dem bleibt letztlich keines u
¨brig.
164
Zu den unvollkommenen Mitteln geh¨ort auch dieses Skriptum. Vielleicht ist es zynisch,
Ihnen viel Freude bei der Lekt¨
ure zu w¨
unschen, aber es soll Ihnen immerhin eine Hilfe
beim Lernen sein. Wenn es bei manchen von Ihnen auch etwas Neugier und Interesse
weckt, hat sich der Aufwand gelohnt.
PS. Wenn Sie Fehler oder Ungereimtheiten entdecken, teilen Sie mir diese bitte mit. An
Meinungen, Ideen und Vorschl¨agen in Bezug auf die Vorlesung und das Skriptum bin ich
sehr interessiert.
Document
Kategorie
Seele and Geist
Seitenansichten
32
Dateigröße
550 KB
Tags
1/--Seiten
melden