9. qualitative analyse mathematischer modelle 103...• das gleichgewicht x0 heißt stabil, wenn...

9. QUALITATIVE ANALYSE MATHEMATISCHER MODELLE 103

9. Qualitative Analyse mathematischer Modelle

In diesem Abschnitt zeigen wir, wie sich Information unabhangig von den genauenZahlenwerten aus einem mathematischen Modell ableiten laßt, wenn es einfach genug fureine analytische Behandlung (mit Bleistift und Papier) ist. Im Allgemeinen sind solcheModelle, genau wie unsere Beispiele, sehr grobe Verallgemeinerungen, eher Karikaturenals getreue Abbilder der Wirklichkeit. Gerade darum kann aber die Arbeit an solchenModellen zum Verstandnis der Prinzipien beitragen, die in einem System wirksam sind.

Fur die Definition der Begriffe Zustandsraum und Trajektorie verweisen wir auf das ersteSemester.

9.1. Die Gleichgewichte des SIR-Modells.

Wir betrachten das Infektionsmodell (SIR-Modell) aus Beispiel 8.5. Die Gleichungen undBedeutung der Parameter sind in Tabelle 8.1 zu finden. Wir erklaren kurz die Bedeutungder Gleichungen:Die Mengenbilanzen werden durch vier Vorgange beeinflußt: Zuwanderung, Tod,Infektion und Heilung. Fur die drei Bevolkerungsgruppen ergeben sich daraus diefolgenden Bilanzen:

Zuwachs=Zuwanderung –Tod ± Infektion ± Heilungddt

S(t) = β −µS(t) −λS(t)I(t)ddt

I(t) = −µI(t) +λS(t)I(t) −γI(t)ddt

R(t) = −µR(t) +γI(t)

Beispiel 9.1. Bestimmen Sie die Gleichgewichtslagen des SIR-Modells.

Losung. Die statische Mengenbilanz erhalt man durch Nullsetzen der Ableitung inder dynamischen Mengenbilanz: Im Gleichgewicht verandern sich die Bestande nicht,daher sind ihre Zuwachsraten gleich Null. Die Zustandsgroßen S, I, R hangen dann nichtmehr von der Zeit ab.

0 = β − µS − λSI,(9.1)

0 = −µI + λSI − γI,(9.2)

0 = −µR + γI.(9.3)

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 Fallen gelost: Fall 1: I = 0 oder Fall 2: λS − (µ + γ) = 0.Wir untersuchen die beiden Falle separat. Zur Unterscheidung der beidenGleichgewichtslagen 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 prasentsind, die Infektion ist ausgestorben, die Bevolkerung ist gesund. Wir wissen dann auch:

R1 =γ

µI1 = 0.

104

Zur Berechnung von S berufen wir uns auf Gleichung (9.1):

0 = β − µS1 − λS10

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

0 = β − µµ + γ

λ− λ

µ + γ

λI2

I2 =β

µ + γ− µ

λ.

Daraus folgt sofort

R2 =γ

µI2 =

βγ

µ(µ + γ)− γ

λ.

Insgesamt erhalten wir die Zahlenwerte

S2 =µ + γ

λ, I2 =

β

µ + γ− µ

λ, R2 =

βγ

µ(µ + γ)− γ

λ.

Allerdings sind diese Werte nur sinnvoll, wenn sie nichtnegativ sind. Negative Infizierteoder Geheilte kann es ja nicht geben. Die Gleichgewichtslage I = 0 wurde bereitsuntersucht. 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, dieKrankheit hat sich eingenistet, sie ist endemisch.Zusammenfassend sagen wir: Es gibt jedenfalls das Gleichgewicht, in dem die Krankheitgar nicht vorkommt, mit den Werten S1, I1, R1. Wenn die Ungleichung (9.4) gilt, soexistiert zusatzlich ein Gleichgewicht, in dem die Krankheit endemisch ist, mit denWerten S2, I2, R2. ¤

Bevor wir in der Untersuchung fortfahren, machen wir eine Beobachtung, die uns vielArbeit ersparen kann: Es gibt keine Ruckwirkung von R auf die Zustandsgroßen S undI. Das heißt, S und I konnen zunachst ohne Berucksichtigung von R mit Hilfe der erstenbeiden Differentialgleichungen untersucht werden. Wenn man uber S und I Bescheidweiß, kann man leicht R mit Hilfe der dritten Gleichung behandeln.


9.2. Stabilitat.

In naturlichen Systemen konnen wir nicht erwarten, daß eine Gleichgewichtslageeingenommen und fur alle Zeit festgehalten wird: Jedes naturliche System ist gewissengeringfugigen Schwankungen und Storungen ausgesetzt. Wenn sich ein Gleichgewichteingespielt hat, und es tritt anschließend eine kleine, vorubergehende Anderung derSystemumgebung auf, so kann das System auf drei Arten auf die Storung reagieren:

• Das System wird vorubergehend aus dem Gleichgewicht ausgelenkt, es erholtsich aber und strebt wieder auf das alte Gleichgewicht zu. Das Gleichgewichtheißt asymptotisch stabil. Gilt dieses Verhalten nur fur kleine Storungen, soheißt das System lokal asymptotisch stabil. Wird von jeder beliebigen Lage ausletztlich immer dasselbe Gleichgewicht angestrebt, so heißt das Gleichgewichtglobal asymptotisch stabil.

• Das System wird aus dem Gleichgewicht ausgelenkt und strebt nicht wiederzuruck, aber wenn die Storung klein genug war, entfernt sich das System nichtweit vom alten Gleichgewicht. Das Gleichgewicht heißt stabil.

• Auch bei kleinen Storungen bricht das alte Gleichgewicht zusammen, und dasSystem steuert weit entfernte Zustande an. Das Gleichgewicht ist instabil.

In der Natur sind instabile Gleichgewichte nur sehr selten zu beobachten: Sie brechen jabei der kleinsten Storung zusammen. Ein instabiles Gleichgewicht erzielt man zumBeispiel, wenn man einen Bleistift auf seiner Spitze ausbalanciert. Will man selbst einSystem in einer bestimmten Gleichgewichtslage halten, so muß man Bedingungenschaffen, die bewirken, daß diese Gleichgewichtslage stabil, am besten asymptotischstabil ist.

Wir geben jetzt eine mathematisch exakte Definition der Stabilitatsbegriffe:

Definition 9.2. Wir betrachten die Differentialgleichung ddt

x(t) = f(x(t)) fur t ≥ 0.Dabei ist x(t) ∈ R

n ein Vektor mit n ≥ 1 Komponenten und f : Rn → R

n ist einegegebene Funktion. Sei x0 ∈ R

n eine Gleichgewichtslage, also f(x0) = 0.

• Das Gleichgewicht x0 heißt global asymptotisch stabil, wenn gilt: Fur jedeLosung x(t) der Differentialgleichung ist der Grenzwert

limt→∞

x(t) = x0.

• Das Gleichgewicht x0 heißt lokal asymptotisch stabil, wenn es einen Radiusδ > 0 gibt, sodaß gilt: Jede Losung x(t), die hochstens im Abstand δ von x0

entfernt beginnt, strebt gegen x0:

Wenn |x(0) − x0| ≤ δ, dann gilt limt→∞

x(t) = x0.

• Das Gleichgewicht x0 heißt stabil, wenn gilt: Zu jeder beliebig vorgegebenenToleranz ǫ > 0 laßt sich ein Radius δ finden, sodaß sich Losungen, die innerhalbdes 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 fur alle t ≥ 0 : |x(t) − x0| ≤ ǫ.

• Ist das Gleichgewicht nicht stabil, so heißt es instabil.

106

Eine genauere Diskussion des Begriffes Stabilitat horen Sie in einer Vorlesung uberDifferentialgleichungen. Beachten Sie, dass in der Definition sinnvollerweise nurautonome Dgl. betrachtet werden, das sind solche, deren rechte Seite f nicht explizit vont abhangt, sondern nur indirekt uber den Zustand x(t).

9.3. Stabilitat 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 Stabilitat). Wir betrachten die Differentialgleichungddt

x(t) = f(x(t)), in Komponenten ausgeschrieben:

d

dt

x1(t)x2(t)

...xn(t)

=

f1(x1(t), · · · , xn(t))f2(x1(t), · · · , xn(t))

...fn(x1(t), · · · , xn(t))

,

und eine Gleichgewichtslage x0. Sei J die Jacobimatrix von f an der Stelle x0, also

J =

∂∂x1

f1(x0) · · · ∂∂xn

f1(x0)...

...∂

∂x1

fn(x0) · · · ∂∂xn

fn(x0)

.

Dann gilt:

(1) Wenn die Realteile aller Eigenwerte von J kleiner als Null sind, dann ist dieGleichgewichtslage x0 lokal asymptotisch stabil.

(2) Wenn J mindestens einen Eigenwert besitzt, dessen Realteil großer als Null ist,dann ist die Gleichgewichtslage x0 instabil.

(3) (Wenn der großte Realteil der Eigenwerte genau Null ist, dann konnen wir mitdiesem Kriterium keine Schlusse ziehen.)

Losung von Beispiel 9.3. In diesem Beispiel lautet die Differentialgleichung (nachWeglassen von R, wie oben erklart):

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 erstnach S und dann nach I partiell differenzieren:

J =

(

∂∂S

(β − µS − λSI) ∂∂I

(β − µS − λSI)∂

∂S(−(γ + µ) + λSI) ∂

∂I(−(γ + µ)I + λSI)

)

=

(

−µ − λI −λSλI −(γ + µ) + λS

)


Wir setzen dann die Gleichgewichtswerte S1 = βµ

und I1 = 0 ein:

J =

(

−µ −λβµ

0 −(γ + µ) + λβµ

)

.

Zur Berechnung der Eigenwerte bilden wir die Matrix zE − J , wobei E dieEinheitsmatrix ist, und bilden die Determinante. Dieser Ausdruck ist dascharakteristische Polynom:

det(zE − J) = det[

z

(

1 00 1

)

−(

−µ −λβµ

0 −(γ + µ) + λβµ

)

]

= det

(

z + µ λβµ

0 z + (γ + µ) − λβµ

)

= (z + µ) (z + γ + µ − λβ

µ).

Die Eigenwerte sind die Nullstellen des charakteristischen Polynoms, also die Losungenvon:

(z + µ)(z + γ + µ − λβ

µ) = 0.

Damit das Produkt Null wird, muß mindestens einer der Faktoren Null sein. Die beidenLosungen sind daher

z1 = −µ und z2 = −γ − µ +λβ

µ.

Der Realteil von z1 (das ist z1 selbst) ist auf jeden Fall negativ. Alle Stabilitat hangtdaher an z2. Der Eigenwert z2 ist genau dann negativ, wenn gilt

− γ − µ +λβ

µ< 0, d.h.

λβ

µ< γ + µ

λ

µ<

γ + µ

β

µ

λ>

β

γ + µ.

Wir erhalten daher die Antwort: Die “gesunde” Gleichgewichtslage ist lokal asymptotischstabil, wenn gilt

µ

λ>

β

γ + µ.

Sie ist instabil, wenn z2 > 0 gilt, also

µ

λ<

β

γ + µ.

Diese Ungleichung ist genau Ungleichung (9.4). Das heißt aber: Wenn die endemischeGleichgewichtslage existiert, dann ist die gesunde Gleichgewichtslage instabil. ¤

108

Wir konnten jetzt ebenfalls die endemische Gleichgewichtslage nach dem Prinzip derlinearisierten Stabilitat untersuchen (wir mussen nur S2 und I2 in die Jacobimatrixeinsetzen). Das Ergebnis einer langeren Rechnung ware: Wenn die endemischeGleichgewichtslage existiert, dann ist sie lokal asymptotisch stabil.

Das Prinzip der linearisierten Stabilitat ist eine Routinemethode zur Untersuchunglokaler Stabilitat. Das globale Verhalten eines dynamischen Systems ist viel schwierigerzu untersuchen, und obwohl es verschiedene Methoden gibt (etwa die Suche nach einemsogenannten Lyapunov-Funktional), gibt es kein Patentrezept. Im Fall des SIR-Modellsfuhrt eine globale Untersuchung zu folgenden Schlussen:Wenn die Ungleichung (9.4) erfullt ist, dann gibt es zwei Gleichgewichtslagen, diegesunde und die endemische. Jede Losung, die mit einer positiven Anzahl von Infiziertenbeginnt, konvergiert gegen das endemische Gleichgewicht. Nur wenn die Infektion in derPopulation gar nicht auftritt, steuert die Population auf das gesunde Gleichgewicht zu.Wenn die Ungleichung (9.4) nicht erfullt ist, dann gibt es nur das gesunde Gleichgewicht,und es ist global asymptotisch stabil. Gleichgultig, von welchem Zustand man ausgeht,die Losung steuert immer auf das gesunde Gleichgewicht zu, die Epidemie stirbt aus.

Betrachten Sie jetzt noch einmal die Simulationen aus Abbildung 8.11. Die linkeAbbildung zeigt eine Losung, die auf das endemische Gleichgewicht zustrebt. In derrechten Abbildung sind die Parameter so gewahlt, daß nur das gesunde Gleichgewichtexistiert, und wir sehen eine Losung, bei der die Infektion ausstirbt. Uberzeugen Sie sichdurch Einsetzen der jeweiligen Modellparameter in Ungleichung (9.4), daß das Verhaltender Losungen mit den Voraussagen der Theorie ubereinstimmt.

Im Vergleich zu den Simulationen sind die Aussagen der qualitativen Analyse desmathematischen Modells viel starker: sie sind mathematisch abgeleitet und betreffenganze Familien von Modellen mit gleicher Struktur und mit Parametern, die bestimmteRelationen erfullen. Je nach Modell kann seine qualitative Analyse aufwandiger (wennnicht gar unmoglich) oder weniger aufwandig sein als quantitative Simulationen. Diesesind jedoch nur Beipiel-Rechnungen: mit Hilfe systematischer Simulationen konnenbestenfalls Trends der Auswirkung von Parameter-Variationen zum Vorschein gebrachtwerden. Erstrebenswert ist es, die beiden Methoden (qualitative Analyse des Modellsund quantitative Simulationen) zu kombinieren. Jede der beiden kann auf die anderebefruchtenden Einfluss haben, ein Vergleich der Ergebnisse kann Analyse- und/oderSimulations-Fehler aufdecken.

9.4. Andere qualitative Fragestellungen.

Gleichgewichte und Stabilitat sind Begriffe, die sich mit geometrischen Vorstellungen gutverknupfen lassen: Sie beschreiben, ob Losungen, das sind Trajektorien imZustandsraum, stets an derselben Stelle stehen bleiben, ob sie sich auf Gleichgewichte zubewegen, oder ob sie sich weit entfernen.


Wir zahlen noch einige andere Begriffe auf, die mit geometrischen Vorstellungenverbunden sind, und sich im Fall einfacher dynamischer Systeme manchmal mit Bleistiftund Papier untersuchen lassen.

Positivitat von Losungen Die Losungen des SIR-Modells sind nur sinnvoll, wennS(t), I(t) und R(t) nicht negativ sind. Wenn das Modell die Wirklichkeit adaquatwiedergeben soll, muß es also positive Zustandsgroßen liefern, wenn es mit positivenAnfangswerten gestartet wird. Fur das SIR-Modell kann man diese Eigenschaft beweisen.Wenn man dagegen zeigen konnte, daß die Losungen nicht positiv bleiben, mußte mandas Modell als unrealistisch ablehnen.

Invariante Teilmengen. Eine Teilmenge des Zustandsraumes, aus der die Losung nichtentweichen kann, wenn sie darin begonnen hat, heißt eine invariante Teilmenge. ZumBeispiel haben wir eben festgestellt, daß eine Losung des SIR-Modells, die im positivenQuadranten (S ≥ 0, I ≥ 0) startet, den positiven Quadranten nie mehr verlaßt. Derpositive Quadrant ist also eine invariante Teilmenge. Wenn man die invariantenTeilmengen des Zustandsraumes fur ein System kennt, kann man sich bereits einwesentlich besseres Bild uber das Verhalten machen.Zum Beispiel konnten die Losungen in einer beschrankten Teilmenge gefangen bleiben,dann konnte man folgern, daß sie nicht gegen unendlich streben konnen. Wenn dieLosungen eines Populationsmodells in einer invarianten Teilmenge gefangen sind, die dieNull nicht enthalt, sagt das Modell voraus, daß die Population nicht ausstirbt.

Periodische Losungen. Eine Losung, die nach endlicher Zeit T zu ihrem Anfangszustandzuruckkehrt, wiederholt denselben Verlauf alle T Zeiteinheiten (dies gilt fur autonomeSysteme x(t) = f(x(t)), deren rechte Seite nicht explizit von t abhangt). Man sagt, dieLosung ist periodisch mit Periode T . Periodische Losungen durchlaufen denZustandsraum in Form von geschlossenen Kurven. Fur manche Systeme, z.B.Rauber-Beute Modelle kann man mit analytischen Methoden entscheiden, ob sieperiodische Losungen besitzen.

Beispiel 9.5. Zeigen Sie, daß das SIR-Modell außer den beiden Gleichgewichtslagenkeine periodischen Losungen besitzen kann.

Losung. Naturlich sind die Gleichgewichtslagen selbst periodische Losungen mitjeder beliebigen Periode T . Wir haben festgestellt, daß jede Losung gegen eine derGleichgewichtslagen konvergiert. Wenn eine Losung also von einem anderen Punkt alseinem Gleichgewichtspunkt ausgeht, kann sie nach einiger Zeit nicht mehr zumAusgangspunkt zuruckkehren, da dieser zu weit von den Gleichgewichten entfernt ware.Also kann sie nicht periodisch sein. ¤

110

Abbildung 9.1. Periodische Losung als Attraktor im Zustandsraum

−1.5 −1 −0.5 0 0.5 1 1.5−1.5

−1

−0.5

0

0.5

1

1.5

x

y

Durchgezogen: periodische Losung. Strichliert: zwei weitere Losungen.

Die Losungen werden im Gegenuhrzeigersinn durchlaufen.

Viele Naturphanomene zeigen periodisches Verhalten, zum Beispiel mechanischeSchwingungen oder der circadiane Rhythmus in der Physiologie. Wegen der Rotation derErde um ihre Achse und um die Sonne werden in Modellen okologischer Systeme dietages- und jahreszeitlichen Bedingungen periodisch angesetzt. Das Modell ist dannnicht-autonom, x(t) = f(t, x(t)), wobei aber f periodisch in t ist. Auch in solchenModellen kann es periodische Losungen geben.

Attraktoren. Eine Teilmenge des Zustandsraumes, auf die jede Losung zulauft, heißt einAttraktor. Ein global asymptotisch stabiles Gleichgewicht ist ein Attraktor, der nur auseinem Punkt besteht. Aber auch die Kurve einer periodischen Losung konnte einAttraktor sein, in diesem Fall spielt sich jede Losung im Lauf der Zeit mehr und mehrauf einen Zyklus ein, der letztlich durch die periodische Losung wiedergegeben wird. DerAttraktor ist dann eine geschlossene Kurve im Zustandsraum.

Beispiel 9.6. Abbildung 9.1 zeigt eine periodische Losung der Gleichung

d

dtx(t) = (1 −

√

x2 + y2)x(t) − 5y(t),

d

dty(t) = (1 −

√

x2 + y2)y(t) + 5x(t).

Der Zustandsraum ist hier zweidimensional. Die Kurve der periodischen Losung bildeteinen Attraktor fur alle Losungen außer der Gleichgewichtslage bei Null. Zwei andereLosungen, die sich dem Attraktor nahern, sind eingezeichnet. Eine einzige Losung wirdnicht angezogen: Die Gleichgewichtslosung x = y = 0. (Nach dem Satz vonPoincare-Bendixson umlauft jede periodische Losung in einem zweidimensionalenZustandsraum mindestens eine Gleichgewichtslosung.) Abbildung 9.2 zeigt dieselbenLosungen, diesmal werden die beiden Zustandsgroßen gegen die Zeit aufgetragen.


Abbildung 9.2. Losungen aus Abbildung 9.1 in anderer Darstellung

0 0.5 1 1.5 2 2.5 3−1.5

−1

−0.5

0

0.5

1

1.5

t

x

0 0.5 1 1.5 2 2.5 3−2

−1

0

1

2

t

y

Durchgezogen: periodische Losung. Strichliert: zwei weitere Losungen.

Strichpunkt: instabiles Gleichgewicht.

9.5. Bifurkationen.

Die dynamischen Eigenschaften eines Modells hangen von den speziellen Werten seinerParameter ab. Als Beispiel haben wir die verschiedenen Epidemieverlaufe im SIR Modellkennengelernt. Im allgemeinen andern sich die Eigenschaften stetig, wenn ein Parameterstetig variiert wird. Es gibt aber auch die Moglichkeit, dass sich die Eigenschaften einesModells plotzlich und sprunghaft andern, wenn ein Parameter einen kritischen Wertstetig uber- oder unterschreitet. Dieses Phanomen der “Verzweigung” des dynamischenVerhaltens heisst Bifurkation, den kritischen Parameterwert bezeichnet man alsBifurkationspunkt.

Betrachten wir als erstes die Gleichung

x(t) = x(t)2 + r

mit einem Parameter r ∈ R. Auf der Suche nach Gleichgewichtslagen mussen wir dieGleichung f(x) = x2 + r = 0 losen. Diese besitzt fur r > 0 keine relle Losung, aber furr ≤ 0 die beiden Losungen x = ±

√−r. Wegen f ′(x) = 2x (′ bedeutet die Ableitung nach

x), ist die positive Gleichgewichtslage x = +√−r instabil, die negative x = −

√−r ist

stabil (Satz von der linearisierten Stabilitat). Wir fassen dies in einemBifurkationsdiagramm zusammen, in dem in horizontaler Richtung der Parameter raufgetragen ist, in vertikaler Richtung die Lage der Gleichgewichtspunkte, sofern esuberhaupt welche gibt. Die Kurve der stabilen Gleichgewichtspunkte wird durchgezogengezeichnet, jene fur die instabilen ist strichliert, siehe Abbildung 9.4 (a). An Hand derDifferentialgleichung uberlegt man sich, dass, fur negatives r, alle Losungen mitx(0) <

√−r ins negative Gleichgewicht konvergieren. Wenn aber r stetig vergroßert

wird, kommt es bei r = 0 zu einer drastischen Anderung der Dynamik: Fur denkritischen Wert r = 0 ist die Gleichgewichtslage x = 0 instabil, und fur r > 0 gehen alleLosungen gegen +∞.

112

Abbildung 9.3. Graphische qualitative Analyse eines eindimensionalen Modells

x

f(x)

r > r*

x

f(x)

−r* < r < r*

x

f(x)

r < −r*

Als zweites Beispiel betrachten wir die Gleichung

x(t) = rx(t) − x(t)3.

Bei Analyse der Gleichgewichtslagen in Abhangigkeit vom reellen Parameter r erkennenwir, dass die fur r < 0 stabile Lage x = 0 bei Durchgang von r durch den kritischen Wertr = 0 instabil wird, siehe Abbildung 9.4 (b). Gleichzeitig entstehen die stabilenGleichgewichtslagen x = ±√

r, die jede Storung von x = 0 einfangen.

Im dritten Beispiel sei q > 0 eine feste Zahl, und r ∈ R ein variabler Parameter in derDifferentialgleichung

x(t) = −x(t)3 + qx(t) + r.

Die Lage und Stabilitat der Gleichgewichte dieses Modells bestimmen wir qualitativ mitgraphischen Argumenten, wie folgt: In Abbildung 9.3 ist fur drei Werte von r die rechteSeite f(x) = −x(t)3 + qx(t) + r als Funktion von x gezeichnet. Die Variation von rbewirkt die vertikale Parallelverschiebung der Kurve −x3 + qx. Es gibt zwei kritischeWerte r∗ > 0 und −r∗ < 0 . Falls r < −r∗ gibt es genau ein x < 0 mit f(x) = 0. Fallsr > r∗ liegt die einzige Nullstelle von f in x > 0. Fur −r∗ < r < r∗ gibt es beide dieserNullstellen von f und zusatzlich eine weitere dazwischen. Auch die Stabilitat derGleichgewichtslagen lasst sich direkt aus den (x, f(x)) Diagrammen ablesen (dies ist beiModellen mit nur einer Zustandskomponente immer moglich). Nullstellen von f mitpositiver Steigung sind instabil: ruckt man ein kleines Stuck nach links, wird f negativ,daher die Ableitung der Losung negativ, was eine weitere Abnahme der Losung bedeutet.Analog bedeutet eine Abruckung nach rechts, dass f positiv wird, was eine weitereZunahme der zugehorigen Losung bewirkt. Ist aber f an der Nullstelle fallend, dannbewirken Abruckungen nach links bzw. rechts eine Ruckbewegung der Losung in positivebzw. negative x-Richtung, also sind solche Gleichgewichtslagen stabil.

Als Zusammenfassung dieser Analyse erhalten wir das Bifurkationsdiagramm inAbbildung 9.4 (c). (In diesem Beispiel wurde r∗ = 2 gewahlt). Wenn der Parameter vonr < r∗ langsam und stetig auf r > r∗ erhoht wird, wahrend sich das System im negativenGleichgewicht befindet, geht dieses plotzlich verloren und der Zustand konvergiert(zeitlich relativ schnell) in das positive Gleichgewicht. Dies gleicht dem Ubergang einerSubstanz vom festen in den flussigen Aggregatzustand, wenn die Temperatur denSchmelzpunkt uberschreitet. Diesen Typus einer Bifurkation bezeichnet man daher als


Abbildung 9.4. Bifurkationsdiagramme

−4 −2 0 2−2

−1

0

1

2

Parameter r

Lage

der

Gle

ichg

ewic

hte

(a)

−1 0 1 2−2

−1

0

1

2

Parameter r

Lage

der

Gle

ichg

ewic

hte

(b)

−4 −2 0 2 4−4

−2

0

2

4

Parameter r

Lage

der

Gle

ichg

ewic

hte

(c)

5 6 7 80.55

0.6

0.65

Parameter k

Par

amet

er r

(d)

nurRuhe−Ggwstabil

beideGgw stabil

nurAusbruch−Ggwstabil

Phasenubergang. Geht dann der Parameter r wieder langsam zuruck, folgt das Systemder durchgezogenen Kurve der positiven Gleichgewichtslage bis zum Wert r = −r∗ undkonvergiert dann (zeitlich relativ schnell) in das negative Gleichgewicht. Das System liegtalso auf dem Weg r < −r∗ → r∗ in einem anderen Zustand vor als auf dem Ruckwegr > r∗ → −r∗. Dieses Verhalten bezeichnet man als Hysterese: bei gleichemParameterwert r hangt es von der Vergangenheit oder Geschichte der Systementwicklungab, welcher Gleichgewichtszustand vorliegt.

Schließlich betrachten wir ein Beispiel dafur, wie dramatische Entwicklungen in einemOkosystem mit Bifurkationen eines mathematischen Modells erklart werden konnen. InKanada gibt es eine Spezies von Raupen (spruce budworm), die sich von den Nadelneiner Fichtenart (balsam fir) ernahrt. Etwa alle 70 Jahre kommt es zu einerRaupenplage, die fast den ganzen Fichtenbestand vernichtet. Daraufhin geht auch dieRaupenpopulation aus Nahrungsmangel drastisch zuruck. Im Lauf der folgenden 70Jahre koexistieren Raupen und Fichten, beide Populationen wachsen langsam an, bis eswieder zu einer Explosion der Raupenpopulation kommt.Das logistische Modell fur die Raupenpopulation X(t) lautet X = RX(1−X/K) mit derGeburtenrate R und der Kapazitat K, die dem Fichtenbestand entspricht. Ferner wirdmodelliert, dass die Raupen von Vogeln gefressen werden, und zwar mit der Rate

V (X) =BX2

A2 + X2.

114

Diese Funktion von X (sogenannte Hill Funktion zweiter Ordnung) hat einen S-formigenVerlauf zwischen V (0) = 0 und V (∞) = B, ahnlich wie die logistische Funktion. Diesbedeutet eine obere Grenze (Sattigung) fur die Dezimierung der Beute (= Raupen)durch die Rauber (= Vogel).Nach Umformungen y = X/A und x(t) = y(A

Bt) bekommen wir das dimensionslose

Modell

x = rx(1 − x

k) − x2

1 + x2

mit den beiden Parametern r und k. r ist proportional zur Wachstumsrate derRaupen-Population (ohne Rauber) und k ist proportional zum Fichtenbestand(Lebensraum fur die Raupen). Dem langsamen Wachstum der Fichten entsprechend ist kein Parameter der sich langsam andert, wahrend die Raupenpopulation x(t) sich denandernden Gleichgewichtslagen laufend anpasst. Eine qualitative Analyse derGleichgewichtslagen des Modells ist in einem sogenannten Phasendiagramm inAbbildung 9.4 (d) zusammengefasst. Hier wird nicht, wie in (a)-(c), die Lage derGleichgewichte gezeigt, sondern es sind Trennlinien zwischen verschiedenenKombinationen der beiden Parameter in der (r, k) Ebene eingezeichnet, die zugehorigeAnzahl und Stabilitat der Gleichgewichte ist als Text in das Diagramm eingeschrieben.Fur ein und dieselbe Geburtenrate r, z.B. r = 0.6, gibt es zu kleinem k nur einGleichgewicht (Koexistenz). Fur k etwa zwischen 5.9 und 6.7 gibt es zwei Gleichgewichte,beide sind stabil. Fur großere k ist aber die koexistierende Raupenpopulation nicht mehrstabil. Wenn, bei gleichbleibender Geburtenrate r = 0.6, der Fichtenbestand k denkritischen Wert k = 6.7 uberschreitet, verlasst die Raupenpopulation x(t) dasRuhegleichgewicht und vernichtet den Fichtenbestand; dann wird k innerhalb kurzer Zeitso klein, dass wieder nur mehr das Ruhegleichgewicht existiert. In diesem Modell sind dieUnterschiede der beiden Wachstumsraten (Zeitskalen) die Ursache fur abrupteperiodische Dynamik mit sehr großer Periodendauer.

9.6. Chaos.

Auch streng deterministische Systeme konnen ein Verhalten zeigen, das Voraussagenpraktisch unmoglich macht, namlich dann, wenn die Losung so empfindlich von denAnfangsdaten abhangt, daß kleinste Anderungen der Anfangswerte zu totalverschiedenen Losungen fuhren konnen. Die Losungskurven solcher Systeme zeigenbizarre, scheinbar vollig unregelmaßige Formen.In der Mathematik spricht man von deterministischem Chaos, wenn das System dreiEigenschaften zeigt:

(1) Es gibt einen Abstand M > 0, sodaß es zu jeder Losung und jedem noch sokleinen Abstand ǫ > 0 eine zweite Losung gibt, die sich im Lauf der Zeit weiterals M von der ersten Losung entfernt, obwohl die gleichzeitigen Anfangswerteder beiden Losungen kleineren Abstand haben als ǫ.

(2) Zu jedem Punkt x des Zustandsraumes und zu jedem beliebig kleinvorgegebenen Abstand ǫ > 0 gibt es periodische Losungen, die naher als ǫ an xvorbeilaufen. (Es gibt ungeheuer viele periodische Losungen mit denunterschiedlichsten Perioden.)


(3) Zu je zwei Punkten x, y im Zustandsraum und jedem beliebig kleinvorgegebenen Abstand ǫ > 0 gibt es eine Losung, die zuerst naher als ǫ an x,und dann naher als ǫ an y vorbeilauft. (Beinahe, namlich wenn man von beliebigkleinen Toleranzen absieht, konnte man von jedem Zustand des Systems ausjeden anderen erreichen: Daher ist der Verlauf der Losung vollig unvorhersagbar,wenn die Anfangsdaten nicht absolut genau vorgegeben sind: Sie kann uberallhingehen.)

Diskrete chaotische Systeme gibt es schon im eindimensionalen Fall, z.B. die logistischeDifferenzengleichung xk+1 = Rxk(1 − xk) fur R = 3.9 und xk ∈ [0, 1]. Kontinuierlicheautonome Systeme konnen erst ab Dimension 3 chaotisch sein, dh. der Zustand mussmindestens 3 Komponenten haben. Wir betrachten zwei Beispiele, beide aus derMechanik.

Beispiel 9.7. Abbildung 9.5 zeigt eine Vorrichtung, bei der ein Pendel mit einerEisenkugel uber vier gleichstarken Magneten hangt, die in einem Quadrat angeordnetsind. Die Pendelbewegung wird von der Erdanziehungskraft und den Magnetfeldernbeeinflußt. Wir halten das Pendel etwas abseits von der Ruhelage fest und lassen es dannaus. Abbildung 9.6 zeigt die Pendelbewegung, die sich ergibt, von oben gesehen. VerfolgenSie die Bewegung.

Beschreibung der Pendelbewegung: Wir beginnen die Bewegung bei denKoordinaten x = 0.6, y = 2. Das Pendel beginnt zu schwingen. Abgelenkt von denMagneten der rechten Seite, schwingt es aber nicht aufs Zentrum zu, sondern zunachstzwischen diesen beiden Magneten. Im Zuge der dritten Schwingung gerat es aber nahegenug in den Einflußbereich der linken Magneten, daß es anschließend den Magnetenlinks hinten umschwingt. Auch die vierte und funfte Schwingung sind weit nach linksausgelenkt.Wir verstehen leicht, warum dieses System sich chaotisch verhalt. Die Schwerkraft undvier Magneten ziehen in verschiedene Richtungen am Pendel. Wenn sich das Pendel inBereichen bewegt, in denen mehrere Magneten ungefahr gleich stark wirken, entscheidenkleinste Unterschiede der Lage und Bewegung daruber, welcher Magnet die Oberhandgewinnt und fur die nachste Schwingung das Pendel am starksten anzieht.

Abbildung 9.5. Pendel aus Beispiel 9.7

116

Abbildung 9.6. Pendelbewegung aus Beispiel 9.7

−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

2

Start der Bewegung an den Koordinaten x = 0.6, y = 2.

Kreise: Magneten. Stern: Zentrum.

Ubrigens widerspricht dieses System nicht unserer Behauptung, dass Chaos erst abDimension 3 auftritt. Zu den zwei Freiheitsgraden der Lage kommen ja noch dieselbenzwei Freiheitsgrade der Geschwindigkeit. Erst wenn wir Lage und Geschwindigkeit desPendels im Augenblick kennen, konnen wir die weitere Bewegung vorausberechnen. DerZustandsvektor, der Lage- und Geschwindigkeitskoordinaten umfaßt, istvierdimensional. ¤

Beispiel 9.8. Ein periodisch erregter nichtlinearer Oszillator wird durch dieDuffing-Gleichung beschrieben:

d2

dt2x(t) + k

d

dtx(t) + ax(t) + bx(t)3 = B cos t

Fur b = B = 0 und a, k > 0 ist dies die klassische Pendelgleichung fur kleine Auslenkungx(t) aus der Ruhelage x = 0 eines senkrecht aufgehangten Pendels im Schwerefeld: ax(t)ist die Ruckstellkraft und kx(t) ist die Reibungskraft. Wenn die Aufhangung aus einemsteifen aber doch verbiegbaren Stahlband besteht, das am oberen Ende unbeweglich ist,dann werden nichtlineare Ruckstellkrafte wirksam: b > 0. In physikalischenExperimenten wird die periodische Erregung B cos t mit starken Magneten realisiert; dieSimulationsergebnisse stimmen qualitativ gut mit der beobachteten Bewegung uberein.Insbesondere fur starke Erregung, das heisst fur großes B beobachtet man chaotischeDynamik. Hier ist der (x, x)-Zustandsraum zwar nur 2-dimensional, aber das Modell istnicht-autonom, dh. die rechte Seite der Dgl. hangt explizit von der Zeit t ab. DurchEinfuhrung der Zeit als Zustandskomponente, z(t) = t, z = 1 kann jedes nicht-autonomeDgl-System mit n Komponenten in ein autonomes der Dimension n + 1 verwandeltwerden.


Abbildung 9.7. Eine Losung x(t) des erregten nichtlinearen Oszillators

0 50 100 150

−2

0

2

x(t)

x’’ + 0.05x’ + x3 = 7.5 cos t

160 180 200 220 240 260 280 300

−2

0

2

320 340 360 380 400 420 440 460

−2

0

2

480 500 520 540 560 580 600 620

−2

0

2

t

Abblidung 9.7 enthalt die Simulationsergebnisse fur die Losung x(t), t ∈ [0, 200π] beiWahl der Parameter k = 0.05, a = 0, b = 1, B = 7.5. Die Erregung ist periodisch mitPeriode 2π, jede Zeile der Graphik enthalt 25 × 2π Zeiteinheiten. Die x-Komponente derTrajektorie ist jedoch nicht periodisch, auch nicht asymptotisch. Dennoch ist dieTrajektorie nicht ganz unregelmaßig: sie bewegt sich in einem beschrankten Bereich,etwa −3.5 < x(t) < 3.5, und es tauchen immer wieder ahnliche Muster auf, wenn auch inunregelmaßiger Abfolge. Insbesondere befindet sich x(t) zu den Zeitpunktentk = 2kπ, k = 1, 2, ... immer in der Nahe des Maximums und ist niemals negativ. Wiekomplex und doch begrenzt und strukturiert die Bewegung ist, vermittelt dieDarstellung in einem sogenannten Poincare-Schnitt, Abbildung 9.8. Gezeigt sind diewahrend einer Matlab-Simulation gesammelten Punkte (x(tk), x(tk)), k = 1001, ..., 50000in der x, x-Ebene (die ersten 1000 Punkte sind als transiente Einschwingphaseweggelassen). Man sieht also nicht die ganze Trajektorie, sondern nur die auf ihrliegenden Punkte im zeitlichen Abstand von je 2π Zeiteinheiten, der Periode derErregung. Trotz der Unregelmaßigkeit der Bewegung bilden diese Punkte nicht einediffus verteilte Wolke, sondern verdichten sich zu einer interessant strukturiertenPunktmenge. Diese ist der Poincare-Schnitt eines Attraktors. Der Attraktor selbst isttatsachlich global asymptotisch stabil: jede Trajektorie (x(t), x(t)) nahert sich ihm furt → ∞. Sowohl der Attraktor selbst als auch sein Poincare-Schnitt ist eine Mengefraktaler, das heisst nicht-ganzzahliger Dimension.

Definition 9.9. Sei E eine Teilmenge des n-dimensionalen Raumes Rn mit

Durchmesser L = sup{||x − y|| : x, y ∈ E} und N(ℓ) die Mindestanzahl von Mengen mitDurchmesser ℓ, die notig sind, um E zu uberdecken. Dann ist die Dimension von Edefiniert durch den Grenzwert

D(E) = limℓ↓0

log N(ℓ)

log(L/ℓ).

118

Abbildung 9.8. Poincare-Schnitt der Duffing-Gleichung.

1.5 2 2.5 3 3.5−6

−4

−2

0

2

4

6

Poincare Schnitt von x’’(t)+0.05 x’(t)+x(t)3 = 7.5 cos(t), tk=2kπ, k=1001,...,50000

x(tk)

x’(t

k)

Man kann zeigen, dass dies eine Verallgemeinerung des herkommlichen Begriffsraumlicher Dimension als maximale Anzahl linear unabhangiger Vektoren darstellt. EineMenge mit Dimension 1 < D < 2 liegt in der Ebene R

2 dichter als eine (eindimensionale)Kurve, aber weniger dicht als jede (zweidimensionale) kontinuierlich gefullte Flache.Fraktale Mengen konnen sehr viel Struktur besitzen (was der chaotischen Dynamik zuwidersprechen scheint), insbesondere die Eigenschaft der Selbstahnlichkeit, das ist dieUnunterscheidbarkeit der Menge selbst von Vergroßerungen immer kleinerer Ausschnitte.Die Existenz eines fraktalen Attraktors ist so typisch fur chaotische dynamische Systeme,dass in der Literatur die beiden Begriffe Chaos und Fraktale oft in einem Atemzug wieSynonyme verwendet werden, obwohl der Zusammenhang alles Andere als trivial ist.

Wenn nun aber der kleinste, unvermeidbare Fehler bei der numerischen Berechnung vonTrajektorien zu total verschiedenen Ergebnissen fuhrt, welche Bedeutung haben dannnumerische Simulationen chaotischer Dynamik? Zeigen Graphiken wie in Abbildung 9.7,9.8 nur vergroßerte Effekte von Rundungsfehlern? Nein: Zwar weicht die numerischberechnete Losung von der exakten Losung (zur selben Anfangsbedingung) exponentiellzunehmend ab – es kann jedoch gezeigt werden, dass es exakte Losungen (zu anderenAnfangsbedingungen) gibt, die der berechneten sehr nahe sind.

Tatsachlich werden fraktale Attraktoren chaotischer Modelle mit Hilfe numerischerMethoden untersucht: es gibt Algorithmen zur naherungsweisen Berechnung vonLyapunov-Exponenten, die ein Maß fur die Sensitivitat der Abhangigkeit vonAnfangsbedingungen sind und damit fur die Starke der Faltung des Attraktors. Und mitdem ’box counting algorithm’ kann die fraktale Dimension eines Attraktors abgeschatztwerden.


Nicht-mechanische Beispiele fur Modelle mit chaotischer Dynamik sindReaktions/Diffusions Gleichungen (Turbulenz vs. Musterbildung), Klima-Modelle(Lorenz’ Schmetterling-Effekt), Modelle der Okonomie (Aktienkurse) und der Okologie(z.B. Nahrungsketten in aquatischen Okosystemen). Diesen Modellen gemeinsam ist derscheinbare Gegensatz der Einfachheit ihrer mathematischen Formulierung einerseits unddem Reichtum und der Komplexitat ihres dynamischen Verhaltens andererseits. Allediese Modelle sind strikt deterministisch, dh. zu einem Anfangszustand gibt es genaueine eindeutig bestimmte Trajektorie, es sind keine Zufallsprozesse modelliert.

120

10. Verteilte Parameter

In unseren bisherigen Modellen haben wir stets die Einflusse der raumlichenGegebenheiten vernachlassigt. Zum Beispiel wurde im Punktreaktormodell angenommen,daß die Neutronen, instabile Isotopen und das spaltbare Material uber den ganzenReaktor gleichmaßig verteilt sind. Beim Entwurf eines Reaktors ist aber gerade dieraumliche Anordnung der verschiedenen Elemente ein wesentliches Kriterium, daher gibtdas Punktreaktormodell nur eine grobe Bilanz, aber keine Auskunft uber die spezifischenVor- und Nachteile von Konstruktionsdetails.Als ein anderes Beispiel konnten wir die Stickstoffbelastung von Grundwasser durchDungung heranziehen. Naturlich konnen wir wieder eine grobe Bilanz aufstellen, indemwir abschatzen, wie groß die gedungte Flache ist, und welcher Anteil des eingetragenenStickstoffs ins Grundwasser absickert. Wenn aber vorhergesagt werden soll, welcheBrunnen besonders gefahrdet sein werden, mussen wir auch einkalkulieren, an welchenStellen der Boden besonders durchlassig fur Sickerwasser ist, und wohin der Schadstoffdurch die Grundwasserstromungen transportiert wird, die in verschiedenen Schichtenverschieden schnell und in verschiedene Richtungen laufen konnen. Je nachdem, wohindie Stromung den Schadstoff tragt, werden manche Bereiche besonders stark durchStickstoff belastet werden.Wenn wir einen Stahltrager belasten, treten Spannungen auf, die bei Uberlastung zumBruch fuhren. Das Gesamtgewicht, das der Trager aushalten muß, ist nur ein Faktor inder Beurteilung, ob eine Belastung zulassig ist. Ebenso wichtig ist, wie sich die Last aufden Trager verteilt. Auch die Spannungen sind nicht im ganzen Trager gleich. An denStellen mit den starksten Spannungen ist bei Uberlastung der Bruch zu erwarten.

Wir werden uns also jetzt Modellen zuwenden, in denen die Zustandsgroßen nicht nurvon der Zeit, sondern auch vom Ort abhangen. Wir sollten von (raumlich) verteiltenZustanden sprechen, aber es hat sich die Phrase “verteilte Parameter” eingeburgert,selbst fur Modelle, die nur skalare Parameter oder gar keine Parameter enthalten.

10.1. Ein Modell fur die Schadstoffbelastung eines Flusses.

Beispiel 10.1. An einem Flußlauf befindet sich eine punktformige Schadstoffquelle (etwader Abfluß einer Fabrik oder Klaranlage), die zu bestimmten Zeitpunkten organischeStoffe in das Wasser einleitet. Durch Diffusion (zufallige Molekularbewegung) undDispersion (zufallige Bewegung der Wassertropfchen) verteilt sich der Schadstoff imWasser, vor allem aber treibt ihn die Stromung flußabwarts. Im Lauf der Zeit wird derStoff durch Oxidation abgebaut, und dadurch belastet er die Wasserqualitat, denn dieSauerstoffkonzentration im Wasser sinkt. An der Wasseroberflache kann das WasserSauerstoff aus der Luft aufnehmen und sich dadurch regenerieren. (SieheAbbildung 10.1.)Gesucht ist ein Modell, das in der Zeit nach einer Schadstoffausschuttung den Verlaufder Schadstoff- und Sauerstoffkonzentration an verschiedenen Stellen des Flussesvorhersagt. Insbesondere sollen die Stellen erfaßt werden, an denen dieSauerstoffverknappung besonders drastisch verspurt wird.

10. VERTEILTE PARAMETER 121

Abbildung 10.1. Schadstoffbelastung eines Flusses

Wir nehmen gleich vorweg, daß der großte Engpaß fur die Sauerstoffversorgung erst einStuck stromabwarts von der Schadstoffquelle auftritt. Auf dem Weg dahin wird namlichSchadstoff abgebaut, und erst dadurch Sauerstoff aufgebraucht.

Der Zustand unseres Systems besteht aus den Konzentrationen von Schadstoff undSauerstoff, in ihrer raumlichen Verteilung uber ein Stuck des Flußlaufs. DieseKonzentrationen konnen und werden sich mit der Zeit andern, wenn nicht dieSchadquelle zu allen Zeiten dieselbe Menge an Schadstoff ausstoßt.Folgende Faktoren bestimmen die Struktur des Systems:

• Quellen: Das Modell dreht sich um die Auswirkungen einer punktformigenQuelle, den Fabriksabfluß. Aber auch das von weiter oben ankommendeFlußwasser kann mit Schadstoff vorbelastet sein. Entlang der Ufer kann, zumBeispiel als Folge landwirtschaftlicher Aktivitaten, weiteres organisches Materialeindringen, wodurch sich zusatzlich eine verteilte Schadstoffquelle ergibt.

• Advektion: Das ist Transport durch die Stromung. Der treibende Faktor ist alsodie Wasserstromung, und die Richtung ist vorgegeben, namlich flußabwarts.

• Diffusion und Dispersion: Das sind zufallige Bewegungen, die fur sich alleindazu fuhren wurden, das sich das Material gleichformig uber das Wasser verteilt.Die Bewegung des einzelnen Molekuls oder Tropfchens ist rein zufallig, aber inder Gesamtschau bewirken diese Vorgange, daß die Mehrheit der Teilchen vonBereichen großer Konzentration zu Bereichen kleiner Konzentration abwandern.

• Chemische und physikalisch-chemische Reaktionen: An jedem Ort im Wasserfindet die Oxidation des Schadstoffes statt. Treibende Kraft sind dieKonzentration an Schadstoff und Sauerstoff. Die Reaktion bewirkt, daß letztlichaller Schadstoff abgebaut und dabei Sauerstoff aufgebraucht wird. An derWasseroberflache findet ein Austausch zwischen dem gelosten Sauerstoff imWasser und dem Luftsauerstoff statt. Treibende Kraft ist der Partialdruck desLuftsauerstoffs im Vergleich zur Konzentration des gelosten Sauerstoffs. Jeweniger Sauerstoff gelost ist, desto mehr wird durch die Wasseroberflacheaufgenommen.

122

Abbildung 10.2. Einteilung in Zellen

• Nicht berucksichtigt in unserer Systembeschreibung sind zum Beispiel dieWasserzufuhr von Seitenarmen (was zur Verdunnung derSchadstoffkonzentration fuhren wurde), Versickerung und Verdunstung vonWasser, zeitliche Schwankungen des Wasserpegels, Ablagerung des Schadstoffesam Boden und Ufer des Flusses, und vieles mehr.

10.2. Modellierung durch Zellen.

Die Idee dieses einfachen Modellansatzes besteht darin, den Flußlauf (willkurlich) ineinzelne Zellen aufzuteilen, und fur jede Zelle eine Schadstoff- und Sauerstoffbilanzaufzustellen. Wir werden dieses Modell unter vereinfachten Annahmen erstellen, wirnehmen namlich an, daß der Schadstoff uber die Tiefe und Breite des Flusses gleichmaßigverteilt ist, sodaß Konzentrationsunterschiede nur in der Richtung stromauf–stromabauftreten. Außerdem vernachlassigen wir die Effekte der Dispersion und Diffusion.

Wir teilen also den Flußlauf der Lange nach in Zellen, die wir stromabwarts mit denZahlen 1, 2, 3, · · · , N durchnummerieren (siehe Abbildung 10.2). Wir machen dievereinfachende (aber falsche) Annahme, daß in jeder Zelle Schadstoff und Sauerstoffgleichmaßig verteilt sind. Je kleiner die Zellen, desto eher entspricht diese Annahme derWirklichkeit, aber desto umfangreicher wird das Modell.

Die Zellen mit Nummern 1, 2, · · · , N haben also zur Zeit t ihre SchadstoffkonzentrationenL1(t), L2(t), · · · , LN(t) und ihre Sauerstoffkonzentrationen C1(t), C2(t), · · · , CN(t). Dassind die gesuchten Zustandsgroßen. Es bewahrt sich, die Schadstoffmasse nicht in MolSchadstoff, sondern indirekt in Sauerstoffaquivalenten anzugeben, also durch die Mengevon Sauerstoff, die benotigt wird, um den Schadstoff vollig zu oxidieren.Als Parameter hat jede Zelle ihr Wasservolumen V1, · · · , VN und ihre Wasseroberflache(Phasengrenzflache zur Luft) a1, · · · , aN . In jede Zelle sickert pro Zeiteinheit ausverteilten Quellen die Schadstoffmenge S1(t), · · · , SN(t) ein. S1(t) beinhaltet auch den


Schadstoff, der aus dem Fabriksabfluß in den Fluß in Zelle 1 eingeleitet wird. Diegesamte Wassermenge, die den Flußquerschnitt in einer Zeiteinheit durchsetzt, sei F ,und weil wir keine Seitenarme, Versickerung und Verdunstung und auch keine zeitlichenPegelschwankungen berucksichtigen, muß diese Zahl das gesamte betrachtete Flußstuckentlang konstant sein: Dieselbe Menge, die ganz oben zufließt, muß ganz unten wiederabfließen. Die Konzentrationen an Schadstoff und Sauerstoff im Wasser, das auf Zelle 1zufließt, bezeichen wir mit L0(t) und C0(t) und stellen uns auf den Standpunkt, dassbeide Funktionen gegeben sind. Weitere Parameter werden eingefuhrt werden, wenn wirdie einzelnen Prozesse modellieren.

Wir betrachten jetzt eine einzelne Zelle, die Zelle mit Nummer i, und stellen eineMengenbilanz auf. Zum Zeitpunkt t befindet sich in der Zelle die Menge ViLi(t) anSchadstoff und ViCi(t) an Sauerstoff. (Das gilt deshalb, weil in Vi m3 Wasser Sauerstoff

mit einer Konzentration von Ci(t)kmol

m3 gelost ist.)

Modellierung der Advektion: Oberhalb der Zelle i befindet sich die Zelle i − 1 mit einerSchadstoffkonzentration von Li−1(t). Pro Sekunde fließen F Liter des Wassers aus deroberen Zelle in die Zelle i. Das liefert einen Schadstoffeintrag von FLi−1(t) und einenSauerstoffeintrag von FCi−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 nehmenSchadstoff mit sich. Das ergibt einen Schadstoffabgang von FLi(t) kmol pro Sekunde.Ebenso errechnet sich der Sauerstoffabgang FCi(t). Insgesamt ergibt sich als Bilanz:

Schadstoffzuwachs durch Advektion: FLi−1(t) − FLi(t),

Sauerstoffzuwachs durch Advektion: FCi−1(t) − FCi(t).

Beachten Sie, daß die Namen L0(t) und C0(t) fur die Konzentrationen des auf Zelle 1zufließenden Wassers so gewahlt wurden, daß diese Bilanzgleichungen auch in der erstenZelle gelten.

Modellierung der Oxidation: Die treibenden Faktoren der chemischen Reaktion sind dieKonzentrationen von Schadstoff und Sauerstoff. (Wir vernachlassigen die Abhangigkeitvon der Konzentration verschiedener Katalysatoren, oder vom Bestand anMikroorganismen bei Oxidation infolge von biologischen Prozessen.) Die Reaktionsrateeiner Reaktion erster Ordnung ist proportional dem Produkt der Schadstoff- undSauerstoffkonzentration. Wir machen aber die vereinfachende Annahme, daß trotz derSauerstoffverknappung immer noch soviel Sauerstoff vorhanden ist, daß die Oxidationnicht wesentlich beeintrachtigt wird: Der limitierende Faktor ist also dieSchadstoffkonzentration. Auf diese Weise bekommen wir ein lineares Modell. Mit kS

bezeichnen wir die Reaktionsgeschwindigkeit. Pro kmol gelosten Schadstoff (inSauerstoffaquivalenten) werden in der Sekunde kS kmol Schadstoff oxidiert, und dabeidieselbe Menge Sauerstoff verbraucht. kS hangt von der Temperatur ab, wir nehmen

124

aber an, dass die Temperatur im Fluss konstant ist. Weil in der Zelle i insgesamt dieSchadstoffmenge ViLi(t) gelost ist, ergibt sich als Bilanz:

Schadstoffbilanz der Oxidation: − kSViLi(t),

Sauerstoffbilanz der Oxidation: − kSViLi(t).

Modellierung der Sauerstoffaufnahme aus der Luft: Wir nehmen an, daß der Partialdruckvon Sauerstoff in der Luft konstant bleibt, was durchaus realistisch scheint. Nach demHenryschen Gesetz gibt es eine Sattigungskonzentration cS von Sauerstoff im Wasser, diediesem Partialdruck das Gleichgewicht halt. Wird das Gleichgewicht gestort, weicht alsoCi(t) von cS ab, so geht pro m2 Phasengrenzflache und pro Sekunde kO(cS − Ci(t)) kmolSauerstoff aus der Luft in Losung. Ist Ci(t) großer als die Sattigungskonzentration cS, soergibt sich ein negatives Vorzeichen: Sauerstoff entweicht aus der Losung in die Luft. Wirmerken an, daß auch kO temperaturabhangig ist, vernachlassigen dies aber. Die Zelle ihat eine Phasengrenzflache von ai Quadratmetern. Diese Zahl ergibt sich nicht einfachaus einer geodatischen Vermessung des Flusses, es mussen auch dieStromungseigenschaften des Wassers in Rechnung gestellt werden: Die feinen Tropfchenuber einem Wasserfall haben zum Beispiel eine gewaltige Oberflache. Die Bilanz derSauerstoffaufnahme aus der Luft ist daher

Sauerstoffzuwachs aus der Luft: kOai(cS − Ci(t)).

Modellierung der Quellen: Jede Zelle hat eine Schadstoffzufuhr von Si(t) kmol proSekunde (in Sauerstoffaquivalenten), die wir als gegeben ansetzen,

Schadstoffzuwachs aus Quellen: Si(t).

Summieren wir die Bilanzen aller Prozesse, so erhalten wir die Zuwachsrate fur dengesamten Schadstoff und des gesamten Sauerstoff in Zelle i. Wir erinnern uns, das dergesamte Schadstoff durch ViLi(t) und der gesamte Sauerstoff durch ViCi(t) gegeben sind:

d

dt(ViLi(t)) = FLi−1(t) − FLi(t) − kSViLi(t) + Si(t),

d

dt(ViCi(t)) = FCi−1(t) − FCi(t) − kSViLi(t) + kOai(cS − Ci(t)).

Division durch die Konstante Vi liefert die Konzentrationsbilanz der Zelle i:

d

dtLi(t) = −(

F

Vi

+ kS)Li(t) +F

Vi

Li−1(t) +Si(t)

Vi

,

d

dtCi(t) = −F + kOai

Vi

Ci(t) − kSLi(t) +F

Vi

Ci−1(t) +kOaicS

Vi

.

Diese Gleichungen gelten fur jede der Zellen i = 1, · · · , N . Mit 10 Zellen wurden wir einSystem von 20 Differentialgleichungen fur die 20 ZustandsgroßenL1(t), · · · , L10(t), C1(t), · · · , C10(t) erhalten. Mit Simulationssoftware furDifferentialgleichungen laßt sich dann das Zellenmodell auf dem Computer nachbilden.Je kleiner die einzelnen Zellen, desto genauer trifft das Modell die Realitat, aber destogroßer ist auch der Rechenaufwand und damit die Anfalligkeit auf Rundungsfehler.


Tabelle 10.1. Eindimensionales Zellenmodell zu Beispiel 10.1

ModellgroßenGroße Einheit Benennung

t sec ZeitN 1 Anzahl der Zellen

Li(t) kmol/m3 Schadstoffkonzentration in Zelle iCi(t) kmol/m3 Sauerstoffkonzentration in Zelle iL0(t) kmol/m3 Schadstoffkonzentration ZuflußC0(t) kmol/m3 Sauerstoffkonzentration ZuflußSi(t) kmol/sec Schadstoffzufluß zu Zelle i (Sauerstoffaquivalente)Vi m3 Volumen in Zelle iai m2 Wasser-Luft-Oberflache in Zelle iF m3/sec StromungkS 1/sec Reaktionsgeschwindigkeit OxidationkO m/sec Austauschgeschwindigkeit Wasser—LuftcS kmol/m3 Sattigungskonzentration O2

Modellgleichungen

d

dtLi(t) = −(

F

Vi

+ kS)Li(t) +F

Vi

Li−1(t) +Si(t)

Vi

,

d

dtCi(t) = −F + kOai

Vi

Ci(t) − kSLi(t) +F

Vi

Ci−1(t) +kOaicS

Vi

.

(Die Gleichungen gelten fur jede Zellen-Nummer i = 1 · · ·N .)

Abbildung 10.3. Dreidimensionales Zellenmodell

Abbildung 10.3 zeigt eine Moglichkeit, den Fluß in Zellen einzuteilen, wenn dievereinfachende Annahme fallen gelassen werden soll, daß die Schadstoffkonzentrationuber den ganzen Querschnitt des Flußbettes konstant sein soll. Wahrend das obenentwickelte Zellenmodell sich mit einer Raumdimension begnugt, namlich derLangsrichtung des Flusses, bildet Abbildung 10.3 die Grundlage fur ein

126

dreidimensionales Modell. Die Modellierung wurde wieder darin bestehen, fur jede Zelleeine Mengenbilanz aufzustellen. Die Details werden komplizierter, weil jede Zelle mehrals nur 2 Nachbarn hat. Vor allem mußte man die Stromungsverhaltnisse nicht nur inLangsrichtung, sondern auch quer und vertikal, und nach Tiefe und Abstand vom Uferaufgeschlusselt kennen. Der Rechenaufwand fur dreidimensionale Modelle ist sehr vielgroßer als fur eindimensionale. Wenn man zur Erhohung der Genauigkeit dieSeitenlangen der Zellen zehntelt, erhalt man im eindimensionalen Modell zurUberdeckung desselben Flußabschnittes zehnmal so viele Zellen, aber imdreidimensionalen Modell die 1000-fache Anzahl an Zellen.

10.3. Das eindimensionale Zellen-Modell in Matlab.

Dieser Abschnitt soll zeigen, wie man das Zellenmodell recht einfach und kompaktprogrammieren kann. In Matlab kommt uns dabei die Moglichkeit, mit Vektoren zuarbeiten, sehr zugute. Zunachst schreiben wir ein Matlab-Script, genannt zellen.m, dasdie Modell-Parameter festlegt, einen Differentialgleichungsloser (hier ode45) aufruft unddie 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 derParameter nur den file editieren und wieder aufrufen muss. Um z.B. N auf 120 zu setzen,braucht man nur eine Null anhangen, wie man im Programmierer-Jargon sagt. Dann lostdas Programm eben ein System von 240 Differentialgleichungen, was fur dieses Beispielauf einem herkommlichen PC auch nicht langer als eine halbe Sekunde dauert. DieVariablen Breite, Tiefe, V, a, Lanf und Canf sind Vektoren der Lange N. DiePunkt-Operatoren .* und ./ bewirken komponentenweise Multiplikation und Divisionvon Matrizen (hier Spaltenvektoren) derselben Große. Mit [a;b] werden dieSpaltenvektoren a und b zu einem Spaltenvektor aneinandergfugt. Wenn wir dieZellenabmessungen nicht fur alle Zellen gleich wahlen, konnen wir dennoch den Rest desProgramms unverandert verwenden: das Programm ist wie das Modell fur NZellenvolumina Vi formuliert.Die rechte Seite der Modellgleichungen aus Tabelle 10.1 wird von ode45 als functionreseze aufgerufen, die in folgendem file reseze.m enthalten ist:

function y=reseze(t,x)

% rechte Seite des Fluss / Schadstoff Zellenmodells


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


if t >= tein && t <= taus

schadst = L_0;

else

schadst = 0;

end

Der Einfachheit halber haben wir die Terme Si(t) weggelassen - der Schadstoffeintragerfolgt nur in Zelle 1 und zwar uber die function L0(t). Die Schadstoffkonzentration,die in Zelle 1 einstromt, wird zwischen den Zeitpunkten tein und taus konstant auf L 0

gesetzt, sonst auf 0. Wir wahlen L 0=5e-3 um hubsche Graphiken zu erhalten, dieserWert ist aber im Vergleich zu cS = 6.44 × 10−4 unrealistisch hoch (und liefert beigenauerer Rechnung (z.B. N = 120) negative Ci). Den Sauerstoffzustrom in Zelle 1setzen wir konstant auf cS.

128

Abbildung 10.4. Simulationsergebnisse mit 12 Zellen zu den Zeitpunktent=900, 1800, 2700, 3600. * = Schadstoffkonzentration, o = Sauerstoffkon-zentration

1 2 3 4 5 6 7 8 9 10 11 12

2

4

6

8

x 10−4 t = 900:

1 2 3 4 5 6 7 8 9 10 11 12

2

4

6x 10

−4 t = 1800:

1 2 3 4 5 6 7 8 9 10 11 12

2

4

6x 10

−4 t = 2700:

1 2 3 4 5 6 7 8 9 10 11 12

2

4

6x 10

−4 t = 3600:

In zellen.m wahlen wir das Simulationsintervall [0, 3600] und den Belastungszeitraum[0, 360]. Dann geben wir im command window den Befehl zellen und erhalten die vierplots in der Abbildung. Man sieht, wie das Maximum der Schadstoff-Belastungflussabwarts stromt, und der Sauerstoffmangel dem Schadstoffabbau entsprechendmitwandert. Das im Beispiel 10.1 gesuchte Minimum an Sauerstoff liegt in der Abbildungbei t = 2700 und zwar in Zelle 7 und betragt etwa 1.55 × 10−4. Nach Beendigung derBelastung und durch die Aufnahme von Sauerstoff aus der Luft erholt sich dieSauerstoffkonzentration; nach genugend langer Zeit (nicht gezeichnet) wird derAusgangszustand wieder angenahert, namlich Schadstoff Li = 0 und Sauerstoff amSattigungswert Ci = cS, i = 1, . . . , 12.

Es fallt auf, dass die zur Verdeutlichung mit Geraden verbundenen Zellen-Werte ziemlichglatte raumliche Verteilungen suggerieren. Zum Beispiel ist die Schadstoffverteilung beit = 900 schon bis zur Zelle 7 eindeutig positiv, wogegen die Verteilung durch dieStromung physikalisch hochstens bis x = 900v = 1800, also bis Zelle 3 reichen durfte(wegen 1800/Zellenlange = 1800/(10000/12) = 2.16). Dies ist ein charakteristischerFehler in Modellen mit gewohnlichen Differentialgleichungen fur raumlich gemittelteZellen, weil durch die Koppelung der Gleichungen sich Anderungen in einer Zelle


mathematisch sofort auf alle Zellen auswirken. Das Modell verschmiert also dieraumliche Verteilung, was auch eine Verringerung des Maximums von Li und eineVergroßerung des Minimums von Ci bewirkt. Immerhin ist wenigstens die Wanderungdes Maximums von Li selbst mit gut 2 Zellen pro plot (Abstand 900s) gleich derStromungsgeschwindigkeit v =2m/s, wie es sein soll.

10.4. Partielle Ableitungen.

Bevor wir weitere Modelle entwickeln, wiederholen wir kurz, was man unter einerpartiellen Ableitung versteht.

Definition 10.2. Ist w eine Funktion zweier Veranderlicher t und x, so ist die partielleAbleitung nach x an der Stelle (t, x) gegeben durch

∂

∂xw(t, x) = lim

h→0

w(t, x + h) − w(t, x)

h.

Das bedeutet, dass man sich eine der Veranderlichen (hier t) festgehalten denkt, und vonder so entstandenen Funktion w(t, ·) einer Veranderlichen (hier x) die Ableitung wieublich definiert. Der Wert der partiellen Ableitung an der Stelle (t, x) ist der Anstieg derTangentialebene im Punkt

(

t, x, w(t, x))

an die von w(·, ·) gebildete Flache uber 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 abhangt. Diepartielle Ableitung von f nach der Variablen xj errechnet man, in dem man in f alleVariablen außer xj wie Konstante ansieht, und f nach der Variablen xj differenziert.Man schreibt fur die partielle Ableitung

∂

∂xj

f(x1, · · · , xn).

Die partiellen Ableitungen hoherer 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) = 3x4y2 + x2 + 5y.

130

Losung.

∂

∂x(3x4y2 + x2 + 5y) = 12x3y2 + 2x

∂

∂y(3x4y2 + x2 + 5y) = 6x4y + 5

∂2

∂x2(3x4y2 + x2 + 5y) =

∂

∂x(12x3y2 + 2x) = 36x2y2 + 2,

∂2

∂y∂x(3x4y2 + x2 + 5y) =

∂

∂y(12x3y2 + 2x) = 24x3y

∂2

∂x∂y(3x4y2 + x2 + 5y) =

∂

∂x(6x4y + 5) = 24x3y

∂2

∂y2(3x4y2 + x2 + 5y) =

∂

∂y(6x4y + 5) = 6x4

¤

Beachten Sie, daß die Reihenfolge bei den partiellen Ableitungen ∂2

∂x∂yund ∂2

∂y∂xkeine

Rolle spielt. Das gilt fur alle Funktionen, deren partielle Ableitungen 2. Ordnung stetigsind.

10.5. Advektion und Reaktion in einer Dimension.

Die Schwache des Zellenmodells besteht in der Annahme, daß die Konzentrationeninnerhalb einer Zelle, also uber einen gewissen Bereich des Flußlaufes, konstant sind.Erst im Grenzubergang zu “unendlich vielen, unendlich kleinen” Zellen wird das Modellrealistisch. Die mathematische Formulierung eines solchen Grenzuberganges fuhrt zueiner Differentialgleichung, die partielle Ableitungen der Zustandsgroßen enthalt, einersogenannten partiellen Differentialgleichung. Wahrend die Theorie und auch dienumerische Losung von partiellen Differentialgleichungen auf dem Computer ziemlichkompliziert ist, genugt zum Erstellen und Verstandnis der Modelle oft ein bißchennaturwissenschaftlicher Hausverstand.

Wir nehmen wieder an, daß in jedem Querschnitt quer zur Stromungsrichtung alleKonzentrationen konstant sind, sodaß wir uns mit einem eindimensionalen Modellbegnugen konnen. Wir fuhren daher nur eine Raumkoordinate ein: Den Abstand x einesFlußquerschnittes vom Fabriksauslaß. Wir betrachten den Flußabschnitt vom Auslaß bisl Meter stromabwarts, also liegt x im Intervall von 0 bis l.

Zu jedem Zeitpunkt t herrscht am Flußquerschnitt mit Raumkoordinate x eineSchadstoffkonzentration von L(t, x) kmol/m3 Sauerstoffaquivalenten, und eineSauerstoffkonzentration von C(t, x) kmol/m3. Die Zustandsgroßen C und L sind alsoFunktionen, die von zwei Variablen, namlich Zeit und Raum, abhangen.

Die Idee der Modellbildung besteht darin, eine kleine Zelle im Fluß, namlich den Bereichzwischen den Raumkoordinaten x und x + ∆x zu betrachten, wobei ∆x eine Lange ist,


die wir uns sehr klein und zunachst positiv vorstellen. Fur diese Zelle stellen wir eineMengenbilanz auf, und betrachten dann, was geschieht, wenn wir die Lange ∆x gegenNull gehen lassen. (Achtung: Entsprechend einer verbreiteten Konvention haben wir hiereine Große, namlich die Lange der Zelle, mit zwei Buchstaben ∆x bezeichnet. Das isteine einzige Große, und nicht etwa ein Produkt einer Zahl ∆ mit x.)

Geometrie und Hydrodynamik des Flusses: Mit A(x) bezeichnen wir dieQuerschnittsflache des Flusses an der Stelle mit Abstand x stromabwarts vom Auslaß.Die Stromungsgeschwindigkeit des Wassers wird naturlich nahe am Ufer und am Grunddes Flusses langsamer sein als in der Mitte. Wir bilden den Mittelwert derStromungsgeschwindigkeit uber den ganzen Querschnitt an der Stelle x und benennenihn mit u(x). Die Wassermenge, die dann pro Sekunde diesen Querschnitt durchsetzt, istA(x)u(x). Wir nehmen an, daß A und u nicht von der Zeit abhangen, also daß der Flußim beobachteten Zeitraum immer die gleiche Menge Wasser fuhrt. Zur Modellierung derSauerstoffaufnahme benotigen wir wieder die Wasser-Luft-Grenzflache. Wir nehmen an,daß an der Stelle x auf ein kurzes Langenstuck ∆x die Oberflache b(x)∆x kommt. DieGroße b(x) hat die Dimension einer Lange. Sie hangt von der Breite des Flusses und demStromungsverhalten ab: Wellen und schwebende Tropfchen vergroßern diePhasengrenzflache.

Schadstoff- und Sauerstoffbestand in der Zelle: Wir gehen davon aus, daß die Lange ∆xder Zelle ziemlich klein ist. Wenn sich der Querschnitt A entlang dieses kurzen Stuckesnicht stark andert, enthalt die Zelle ungefahr ein Wasservolumen von A(x)∆x. In diesemWasservolumen sind Sauerstoff und Schadstoff mit Konzentrationen C(t, x) und L(t, x)gelost. 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 stromaufwarts dringen pro Sekunde A(x)u(x) m3

Wasser durch den Querschnitt am Zellanfang x in die Zelle ein. Dieses Wasser hat eineSauerstoffkonzentration von C(t, x) und eine Schadstoffkonzentration von L(t, x), dahertreten mit dem Wasser in jeder Sekunde A(x)u(x)C(t, x) kmol Sauerstoff undA(x)u(x)L(t, x) kmol Schadstoff (in Sauerstoffaquivalenten) ein. Am unteren Ende derZelle, an der Stelle x + ∆x, werden nach derselben Schlußrechnung jede SekundeA(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 folgendeBilanz

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

132

Modellierung der verteilten Quellen: Wir nehmen an, daß an der Stelle x auf ein kurzesLangenstuck der Lange ∆x pro Sekunde die Schadstoffmenge s(t, x)∆x vom Ufer herzufließt. Wir erhalten daher fur unsere Zelle:

Schadstoffzufluß aus verteilten Quellen: s(t, x)∆x.

Auf die Modellierung der punktformigen Quelle bei x = 0 kommen wir spater zusprechen. Wenn 0 nicht gerade innerhalb unserer sehr kleinen Zelle zu liegen kommt,spielt diese Quelle fur die Zelle ohnehin nur indirekt uber die Advektion eine Rolle.

Modellierung der Oxidation: Wie im Zellenmodell beschreiben wir die Oxidation als eineReaktion erster Ordnung in Abhangigkeit von der Schadstoffkonzentration, das heißt, wirgehen davon aus, daß das Sauerstoffangebot ausreichend ist. Mit kS bezeichnen wirwieder die Reaktionsgeschwindigkeit der Oxidation. In der Zelle ist die SchadstoffmengeA(x)∆xL(t, x) zur Oxidation vorhanden, davon oxidiert in der SekundekSA(x)∆xL(t, x). Daraus ergibt sich die folgende Bilanz:

Schadstoffbilanz der Oxidation: − kSA(x)∆xL(t, x),

Sauerstoffbilanz der Oxidation: − kSA(x)∆xL(t, x).

Modellierung der Sauerstoffaufnahme aus der Luft: Als Phasengrenzflache stehen b(x)∆xQuadratmeter zur Verfugung, das Konzentrationsgefalle von Sattigungskonzentration cS

zur tatsachlichen Sauerstoffkonzentration ist cS − C(t, x). Mit einerAustauschgeschwindigkeit kO ergibt sich die Bilanz

Sauerstoffzuwachs aus der Luft: kOb(x)∆x(cS − C(t, x)).

Mengenbilanz: Die Mengenbilanz fur die Zelle gibt nun die Anderungsraten derSchadstoff- und Sauerstoffmengen in der Zelle an. Dazu werden die Beitrage vonAdvektion, Quellen, Oxidation und Reaeration aus der Luft summiert. DieAnderungsraten sind, wie gewohnt, die Ableitungen der Mengen nach der Zeit. Da jetztaußer der Zeit noch die Raumvariable zur Verfugung steht, schreiben wir die Ableitungenals partielle Ableitungen:

∂

∂t(A(x)L(t, x)∆x) = A(x)u(x)L(t, x) − A(x + ∆x)u(x + ∆x)L(t, x + ∆x)

+ s(t, x)∆x − kSA(x)∆xL(t, x),

∂

∂t(A(x)C(t, x)∆x) = A(x)u(x)C(t, x) − A(x + ∆x)u(x + ∆x)C(t, x + ∆x)

− kSA(x)∆xL(t, x) + kOb(x)∆x(cS − C(t, x)).


Wir kurzen beide Gleichungen durch ∆x und erhalten

∂

∂t(A(x)L(t, x)) = −A(x + ∆x)u(x + ∆x)L(t, x + ∆x) − A(x)u(x)L(t, x)

∆x− kSA(x)L(t, x) + s(t, x),

∂

∂t(A(x)C(t, x)) = −A(x + ∆x)u(x + ∆x)C(t, x + ∆x) − A(x)u(x)C(t, x)

∆x− kSA(x)L(t, x) + kOb(x)(cS − C(t, x)).

Fur ∆x < 0 drehen sich die Vorzeichen der Advektionsbilanz um, alle anderen Termeenthalten die Lange −∆x als Faktor. Kurzen durch −∆x fuhrt dann also auf genaudieselben Gesamtbilanzen fur die Zelle zwischen x + ∆x und x. Schließlich lassen wir ∆xgegen Null gehen und wenden die Definition 10.2 der partiellen Ableitung an (mit h = ∆xund w(t, x) = A(x)u(x)L(t, x) bzw. A(x)u(x)C(t, x)). Wir erhalten im Grenzwert:

A(x)∂

∂tL(t, x) = − ∂

∂x(A(x)u(x)L(t, x)) + s(t, x) − kSA(x)L(t, x),

A(x)∂

∂tC(t, x) = − ∂

∂x(A(x)u(x)C(t, x)) − kSA(x)L(t, x) + kOb(x)(cS − C(t, x)).

Die Modellbildung fuhrt also auf ein System von zwei partiellen Differentialgleichungen,in denen sowohl die Zeit- als auch die Raumableitungen der Zustandsgroßen L(t, x) undC(t, x) auftreten. Zu jedem Zeitpunkt t ist der Zustand des Modells ein Paar von zwei“verteilten” Funktionen (L(t, ·), C(t, ·)).

Merksatz 10.5. Zur Modellierung zeitlich und raumlich kontinuierlich veranderlicherZustande betrachtet man die dynamische Mengenbilanz fur ein kleines Raumelement. Furdie Anderung in der Zeit t verwendet man die partielle Ableitung nach t und schreibt dieZu- und Abflusssraten in Form von Differenzenquotienten in den Raumkoordinaten.Lasst man das Volumen (bzw. die Flache bzw. die Lange) gegen Null gehen, ergeben sichpartielle Ableitungen nach den Raumkoordinaten, also partielle Differentialgleichungen.

Bemerkung. Eine alternative (hier nicht vorgefuhrte) Vorgangsweise ist es, die zeitlicheAbleitung des Integrals des Zustandes uber ein Raumelement dem Oberflachenintegralseiner Anderungsrate gleichzusetzen. Dies fuhrt - fur hinreichend differenzierbareZustande - mit dem Satz von Gauß zu partiellen Differentialgleichungen fur dieIntegranden und zu denselben Modellen wie im Merksatz.

Modellierung der Zuflusse und der punktformigen Quelle: Wie schon im Zellenmodell,bezeichnen wir mit C0(t) die Sauerstoffkonzentration des zufließenden Wassersunmittelbar oberhalb des Fabriksauslasses, und mit L0(t) die Schadstoffkonzentrationdes zufließenden Wassers. Mit S(t) bezeichnen wir den Schadstoffausstoß der Fabrik, inkmol Sauerstoffaquivalenten pro Sekunde. Wir nehmen an, daß wir diese Großen kennen.

134

An der Stelle x = 0, unmittelbar nach dem Fabriksauslaß, betragt also dieSauerstoffkonzentration

C(t, 0) = C0(t).

Unmittelbar unter dem Auslaß fließen pro Sekunde u(0)A(0) m3 Wasser. Wenn wirannehmen, daß der Beitrag des Auslasses zur Wassermenge verschwindend klein ist,fließt fast dieselbe Menge von oberhalb des Auslasses zu. Sie tragt pro Sekunde vonflußaufwarts L0(t)u(0)A(0) kmol Schadstoff (Sauerstoffaquivalente) heran. Dazu kommenjede Sekunde S(t) kmol Schadstoff von der Fabrik. Daher fließen unmittelbar unter demAuslaß jede Sekunde S(t) + u(0)A(0)L0(t) kmol Schadstoff, gelost in u(0)A(0) m3

Wasser. Das ergibt die Konzentration (die wir der Kurze halber mit L1(t) bezeichnen):

L(t, 0) = L0(t) +S(t)

u(0)A(0)= L1(t).

Die Information uber die Zuflußkonzentrationen beschreibt die beiden Zustandsgroßen Lund C an einem Randpunkt des betrachteten Gebietes, namlich bei x = 0. Sie sindRandbedingungen.

Anfangsbedingungen: Um die Entwicklung der Schadstoffkonzentration im Flußvorhersagen zu konnen, brauchen wir letztlich noch Informationen uber denAnfangszustand. Wir nehmen an, daß wir die Schadstoffkonzentration La und dieSauerstoffkonzentration Ca zum Zeitpunkt t = 0 im ganzen untersuchten Flußgebietkennen.

L(0, x) = La(x),

C(0, x) = Ca(x).

Diese Bedingungen legen die Zustandsgroßen zu Anfang des betrachteten Zeitraumesfest, sie sind Anfangsbedingungen.

Insgesamt erhalten wir ein Modell aus zwei partiellen Differentialgleichungen, mit zweiAnfangsbedingungen, und zwei Randbedingungen. Das vollstandige Modell ist inTabelle 10.2 zusammengefaßt.

Merksatz 10.6. Eine partielle Differentialgleichung beschreibt eine unbekannteFunktion, die von mehreren Variablen abhangt, durch eine Beziehung zwischen derFunktion und ihren partiellen Ableitungen. In der Modellierung verteilter Zustande tretenpartielle Differentialgleichungen zusammen mit Anfangsbedingungen undRandbedingungen auf. Anfangsbedingungen geben die Werte der unbekanntenFunktion(en) zu Anfang des Beobachtungszeitraumes an. Randbedingungencharakterisieren die unbekannten Funktionen am Rand des beobachteten Gebietes.


Tabelle 10.2. Reaktions-Advektionsgleichung zu Beispiel 10.1

ModellgroßenGroße Einheit Benennung

t sec Zeitx m Distanz vom Auslaß

L(t, x) kmol/m3 Schadstoffkonzentration (O2-Aquivalente)C(t, x) kmol/m3 SauerstoffkonzentrationL0(t) kmol/m3 Schadstoffkonzentration ZuflußC0(t) kmol/m3 Sauerstoffkonzentration ZuflußS(t) kmol/sec Schadstoffzufluß von Auslaß

s(t, x) kmol/(m sec) Schadstoffzufuhr verteilte Quellen pro LangeLa(x) kmol/m3 Schadstoffkonzentration zu AnfangCa(x) kmol/m3 Sauertoffkonzentration zu Anfang

kS 1/sec Reaktionsgeschwindigkeit OxidationkO m/sec Austauschgeschwindigkeit Wasser—LuftcS kmol/m3 Sattigungskonzentration O2

A(x) m2 Querschnittsflacheu(x) m/sec Stromungsgeschwindigkeitb(x) m Phasengrenzflache pro Lange

Partielle Differentialgleichungen

(10.1)

A(x)∂

∂tL(t, x) = − ∂

∂x(A(x)u(x)L(t, x))

+ s(t, x) − kSA(x)L(t, x),

A(x)∂

∂tC(t, x) = − ∂

∂x(A(x)u(x)C(t, x))

− kSA(x)L(t, x) + kOb(x)(cS − C(t, x)).

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

Bemerkung. In der soeben durchgefuhrten Modellbildung haben wir alle betrachtetenTeil-Bilanzen (fur die Stromung, die Oxidation, den Sauerstoffeintrag und denSchadstoffzutritt aus den Ufern) zusammengefasst. Falls man nur die Advektion(=Stromung) berucksichtigt und alle anderen Effekte einmal weglasst, dann steht in derdurch ∆x dividierten Mengenbilanz fur 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. Aus der Inkompressibilitat

136

des Wassers folgt A(x)u(x) = F mit einer von x und t unabhangigen Konstanten F .Diese kann also, nach Grenzubergang ∆x → 0, aus der partiellen Ableitung nach xherausgezogen werden. Kurzt man die gesamte Gleichung dann durch A(x), so erhaltman eine partielle Differentialgleichung der Gestalt

∂

∂tL(t, x) + u(x)

∂

∂xL(t, x) = 0, t > 0, x ∈ (0, xL).

Dies ist die klassische Advektionsgleichung fur die Bewegung einer inkrompressiblenFlussigkeit mit der raumlich und zeitlich veranderlichen Eigenschaft L(t, x). DieStromungsgeschwindigkeit u(x) ist umgekehrt proportional zur Breite desStromungskanals A(x), namlich u(x) = F/A(x). Fur u > 0 benotigt manRandbedingungen am Anfang des Kanals bei x = 0, fur u < 0 kann hingegen nur eineRandbedingung am Ende bei x = xL vorgegeben sein. Eine Geschwindigkeitsverteilungu(x) mit wechselndem Vorzeichen ist wegen der Inkompressibilitat unmoglich (beiKompressibilitat verwendet man verallgemeinerte Konzepte fur zulassige ’Losungen’).

10.6. Charakteristische Kurven.

Die Hauptstoßrichtung dieses Kapitels ist die Modellbildung, und die ist nunabgeschlossen, sofern wir Diffusion und Dispersion nicht berucksichtigen wollen. AufDiffusion und Dispersion kommen wir spater zuruck. Wir werden jetzt die vorliegendenModellgleichungen losen. Das Modell in Tabelle 10.2 sieht zwar abschreckend aus, aber esist ein sehr einfaches Beispiel einer partiellen Differentialgleichung. PartielleDifferentialgleichungen konnen außerordentlich schwierig sein, sowohl bei dertheoretischen Behandlung als auch bei der numerischen Losung auf dem Computer.

Um die Darstellung der Berechnungen einfacher zu gestalten, machen wir die folgendenAnnahmen: Querschnitt A(x) = A, Stromungsgeschwindigkeit u(x) = u, und diePhasenoberflache pro Langeneinheit b(x) = b sind konstant, also unabhangig von x.(Statt von einem Fluß sprechen wir wohl eher von einem regulierten Kanal.) Die einzigeSchadstoffquelle ist der Fabrikauslaß: s(t, x) = 0.

Der Beobachterstandpunkt bei der Modellbildung war fest: Wir haben aus dem Flußlaufeine kurze Zelle herausgegriffen, und berechnet, wieviel Schadstoff und Sauerstoff dasWasser im Lauf der Zeit durch diese Zelle hindurchtragt. Wir konnten uns vorstellen,daß ein Beobachter an einer festen Stelle x am Ufer steht, und an einer Angel seineMeßsonden ausgehangt hat. Er mißt damit den Verlauf der Konzentrationen L(t, x) undC(t, x) im Lauf der Zeit. Wenn sich diese Konzentrationen im Lauf der Zeit andern, kanndas zwei Grunde haben: Einerseits verschieben die chemischen Reaktionen dieKonzentrationen. Andererseits kann der Schadstoffausstoß der Quellen zu verschiedenenZeiten unterschiedlich sein, und dadurch kann Wasser, das erst auf die Meßstelle zufließt,eine andere Vorgeschichte haben und andere Schadstoffwerte liefern als Wasser, das ebenvorbeigeflossen ist. Weil der Beobachter stets am selben Standort x steht, ist x fur ihneine Konstante. Die Anderungsraten von L und C, die er wahrnimmt, sind dieAbleitungen der Konzentrationen, die man erhalt, wenn man nach der Zeit differenziertund x konstant laßt: die partiellen Ableitungen ∂

∂tL(t, x) und ∂

∂tC(t, x).


Die Idee des folgenden Losungswegs besteht in einer geschickten Veranderung desBeobachterstandpunktes: Wir stellen uns vor, daß der Beobachter zu einem gegebenenZeitpunkt t0 an der Stelle x = 0 in ein Boot steigt und sich mit der Stromung treibenlaßt. Vom Boot aus hangt er seine Meßsonden in das Wasser, sodaß sie mit der Stromungmittreiben. Weil die Stromung vollig von Zufalligkeiten frei ist (laut Modellannahme),sind die Sonden immer vom selben Wasser umgeben, und alle Anderungen derKonzentrationen, die dieser Beobachter wahrnimmt, sind auf die chemischen Reaktionenzuruckzufuhren.

Weil der Beobachter zum Zeitpunkt t = t0 an der Stelle x = 0 ablegt, und mit derStromungsgeschwindigkeit u fortgetragen wird, befindet er sich zur Zeit t > 0 an derStelle 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 konnen sich diese Werte andern, und dieAnderungsraten bezeichnen wir mit d

dtL(t, (t − t0)u) und d

dtC(t, (t − t0)u). Das sind nicht

mehr die partiellen Ableitungen nach der Zeit, denn die Raumkoordinate x verandertsich mit der Zeit: das Boot treibt weiter. Wir sprechen von totalen Ableitungen: Sieresultieren aus der Veranderung von Raum und Zeit zugleich.

Um die totale Ableitung mathematisch zu beschreiben, halten wir uns die Bedeutung derpartiellen Ableitungen vor Augen: Die partielle Ableitung ∂

∂xL(t, x) ist die

Anderungsrate bezogen auf Lange, die sich ergibt, wenn man die Konzentration Lgleichzeitig an der Stelle x und ein “unendlich kleines” Stuck ∆x weiter stromabwartsmißt. Die partielle Ableitung ∂

∂tL(t, x) ist die Anderungsrate bezogen auf Zeit, die sich

ergibt, wenn man die Konzentration an derselben Stelle zur Zeit t und eine “unendlichkleine” Zeitspanne ∆t spater mißt:

L(t, x + ∆x) − L(t, x) ≈ ∆x∂

∂xL(t, x),

L(t + ∆t, x) − L(t, x) ≈ ∆t∂

∂tL(t, x).

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 wurde sich eine Anderung derMeßwerte von

∆t∂

∂tL(t, x)

ergeben. Außerdem erfolgt aber eine Anderung des Standortes um den Betrag ∆x = u∆tund schlagt sich in einer zusatzlichen Verschiebung der Meßwerte von

u∆t∂

∂xL(t, x)

nieder. Daher ist die gesamte Anderung

∆t∂

∂tL(t, x) + u∆t

∂

∂xL(t, x).

138

Die Anderungsrate, bezogen auf das Zeitintervall ∆t, ist dann die totale Ableitung:

(10.2)

d

dtL(t, u(t − t0)) =

∂

∂tL(t, u(t − t0)) + u

∂

∂xL(t, u(t − t0)),

d

dtC(t, u(t − t0)) =

∂

∂tC(t, u(t − t0)) + u

∂

∂xC(t, u(t − t0)).

Merksatz 10.7. Ein Beobachter durchlauft den Raum entlang einer Kurve(

x1(t), · · · , xn(t))

. Die totale Ableitung einer Funktion f(x1(t), · · · , xn(t)) entlang dieserKurve errechnet sich durch

d

dtf =

∂

∂x1

f · d

dtx1 + · · · + ∂

∂xn

f · d

dtxn.

Die totale Ableitung beschreibt die Anderungsrate einer raumlich verteilten Eigenschaftf , die ein Beobachter feststellt, der sich zu den Zeitpunkten t jeweils an der Stelle mitden Raumkoordinaten x1(t), · · · , xn(t) befindet.

Wir ordnen die partiellen Differentialgleichungen (10.1) aus dem Modell in Tabelle 10.2etwas um und heben den konstanten Faktor A(x) = A heraus. Unter Beachtung, daßauch u konstant und s = 0 ist, erhalten wir

∂

∂tL(t, x) + u

∂

∂xL(t, x) = −kSL(t, x),

∂

∂tC(t, x) + u

∂

∂xC(t, x) = −kSL(t, x) +

kOb

A(cS − C(t, x)).

Wir setzen jetzt in die totalen Ableitungen (10.2) der Schadstoff- undSauerstoffkonzentration die partiellen Differentialgleichungen (10.1) aus dem Modell inTabelle 10.2 ein:

d

dtL(t, (t − t0)u) =

∂

∂tL(t, (t − t0)u) + u

∂

∂xL(t, (t − t0)u)

= −kSL(t, (t − t0)u),

d

dtC(t, (t − t0)u) =

∂

∂tC(t, (t − t0)u) + u

∂

∂xC(t, u(t − t0))

= −kSL(t, (t − t0)u) +kOb

A(cS − C(t, (t − t0)u)

Alle partiellen Ableitungen sind verschwunden, ubrig bleiben nur mehr die totalenAbleitungen, die Anderungsraten, die der Beobachter im Boot feststellt. Und diesetotalen Ableitungen resultieren aus den Konzentrationsanderungen durch Oxidation undSauerstoffaufnahme aus der Luft. Die Effekte der Stromung werden vom Boot aus nichtwahrgenommen, weil es mitschwimmt. Statt zwei partiellen Differentialgleichungen hatder Beobachter im Boot nur zwei gewohnliche Differentialgleichungen zu losen, um seineMeßwerte zu erklaren:


(10.3)

d

dtL(t, (t − t0)u) = −kSL(t, (t − t0)u),

d

dtC(t, (t − t0)u) = −kSL(t, (t − t0)u) +

kOb

A(cS − C(t, (t − t0)u))

Das Boot bewegt sich auf dem Weg x(t) = (t − t0)u durch den Raum. Diese Kurve heißtin der Mathematik eine charakteristische Kurve der partiellenDifferentialgleichungen (10.1): Eine charakteristische Kurve verlauft so im Raum, daß dietotale Ableitung der Losung entlang einer charakteristische Kurve gerade alle partiellenAbleitungen aus der Differentialgleichung zusammenfaßt. Die Losung wird entlang dercharakteristischen Kurve also durch eine Differentialgleichung beschrieben, die nur eineAbleitung enthalt, namlich die totale Ableitung nach der Zeit.

Wir wiederholen: Die Konstruktion von charakteristische Kurven, und damit dieUmwandlung der partiellen Differentialgleichungen in gewohnlicheDifferentialgleichungen, war nur moglich, weil die Stromungsverhaltnisse a priori klarund frei von Zufalligkeiten waren. Die Methode scheitert, wenn Diffusions- undDispersionseffekte zum Tragen kommen.

10.7. Losung der Reaktions-Advektions-Gleichungen.

Die Losung der gewohnlichen Differentialgleichungen (10.3) ist nun eine Standardaufgabefur eine Vorlesung aus Differentialgleichungen. Zur Abkurzung schreiben wir

k1 =kOb

A.

Wir erhalten die Losungen

(10.4)

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

− kSL1(t0)

k1 − kS

[e−kS(t−t0) − e−k1(t−t0)].

was man durch Einsetzen uberprufen kann.

Der Vollstandigkeit halber geben wir den Losungsweg an:Wir losen zuerst die Differentialgleichung fur L, denn sie ist von C entkoppelt: DieSauerstoffkonzentration kommt in der Gleichung gar nicht vor:

d

dtL(t, (t − t0)u) = −kSL(t, (t − t0)u)

Das ist eine homogene lineare Differentialgleichung

d

dty(t) = −kSy(t)

mit y(t) = L(t, (t − t0)u). Die Losung ergibt die Exponentialfunktion

y(t) = e−kS(t−t0)y(t0),

140

alsoL(t, (t − t0)u) = e−kS(t−t0)L(t0, 0) = e−kS(t−t0)L1(t0).

Wir setzen in die Gleichung fur die Sauerstoffkonzentration ein:

d

dtC(t, (t − t0)u) = −kSL(t, (t − t0)u) +

kOb

A(cS − C(t, (t − t0)u)

= −k1C(t, (t − t0)u) + k1cS − kSL1(t0)e−kS(t−t0).

Diese Gleichung ist eine inhomogene lineare Differentialgleichung der Form

d

dty(t) = −k1y(t) + g(t),

mit y(t) = C(t, (t − t0)u) und g(t) = k1cS − kSL1(t0)e−kS(t−t0). Ihre Losung ist nach der “Variation der

Konstanten”-Formel

y(t) = e−k1(t−t0)y(0) +

∫ t

t0

e−k1(t−s)g(s) ds,

das heißt also

C(t, (t − t0)u) = e−k1(t−t0)C(t0, 0) +

∫ t

t0

e−k1(t−s)[k1cS − kSL1(t0)e−kS(s−t0)] ds

= e−k1(t−t0)C0(t0) + cS [1 − e−k1(t−t0)]

− kSL1(t0)

k1 − kS

[e−kS(t−t0) − e−k1(t−t0)].

Um die Losung L(t, x) zu bestimmen, mussen wir jetzt nur den Zeitpunkt t0 so wahlen,daß (t − t0)u = x gilt, also

t0 = t − x

u.

Dann ist

L(t, x) = L(t, (t − t0)u) = e−kS(t−t0)L1(t0)

= e−kSx/uL1(t −x

u).

Ebenso errechnet man

(10.5) C(t, x) = e−k1x/uC0(t −x

u) + cS

[

1 − e−k1x/u]

− kSL1(t − xu)

k1 − kS

[

e−kSx/u − e−k1x/u]

.

Das sind die Streeter-Phelps-Gleichungen. Sie gelten fur solche (t, x), fur die t − xu≥ t0,

also x ≤ u(t − t0), andernfalls konnten die Funktionen der Randbedingungen L1 und C0

nicht ausgewertet werden. Fur x ≥ u(t − t0) muß man die Anfangsbedingungen statt derRandbedingungen verwenden. Die Notwendigkeit dieser Fallunterscheidung istanschaulich verstandlich: Fur grosse x und kleine t sind die Anfangsbedingungenmaßgeblich, die Information aus dem linken Rand kommt erst nach genugend langer Zeitnach x. Fur grosses t und kleines x sind die Anfangsbedingungen schon zur Ganzevorbeigestromt, nur mehr die Zustromung aus dem linken Rand ist maßgeblich.

Dies wird in Abbildung 10.5 illustriert. Wir sehen die beiden Teilbereiche der (t, x)Ebene, die von der Randbedingung (oberes Bild) bzw. von der Anfangsbedingungbeeinflusst werden (unteres Bild). Beachten Sie, dass hier keine Losungen oder Werte vonL,C gezeigt sind, sondern nur eine Draufsicht auf die Ebene der unabhangigen Variablen


Abbildung 10.5. Einflussgebiete der Rand- und Anfangsbedingung

-

6

©©©©©©©©©©©©©©©©©©©©©©©©

t

x

0

xL

t0

©©©©©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©

©©©*©©©©©©©©©©

©©©*©©©©©©©

©©©*©©©©©©©*

©©

t

• (t, u(t − t1))

t1

bei x = 0 ist die Randbedingung L1 gegeben, insbesondere L1(t1)

-

6

©©©©©©©©©©©©©©©©©©©©©©©©

t

x

0

xL

t0

©©©©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©©©©

©©©*©©©©©©©©©©©©©

©©©*©©©©©©©©©©

©©©*©©©©©©©©

©©©*©©©©©

©©©*©©

t

• (t, x1 + u(t − t0))

x1

bei t = t0 ist die Anfangsbedingung La gegeben, insbesondere La(x1)

– Funktionen dieser Variablen mussten senkrecht zu dieser Ebene aufgetragen werden,oder am Papier in einer perspektivischen ’dreidimensionalen’ Darstellung. DieInformation wandert entlang der charakteristischen Kurven, die Richtung ist durchPfleile gekennzeichnet. Im Bereich oberhalb von x = 0 geht die Information von derRandbedingung aus, im Dreieck rechts von t = t0 stammt die Information aus derAnfangsbedingung. Die charakteristischen Kurven sind hier wegen der konstantenStromungsgeschwindikeit u Geraden mit Anstieg u. Fur variableStromungsgeschwindigkeiten (bei variablem Flussquerschnitt) waren diecharakteristischen Kurven entsprechend gekrummte Kurven in der (t, x)-Ebene.

Als Anwendung der berechneten Losung in (10.5) konnte man das Minimum von C(Nullsetzen der Ableitung nach x), also fur festes t jenen Punkt im Fluß bestimmen, andem die Sauerstoffkonzentration minimal ist, und damit das maximale Sauerstoffdefizitvorhersagen. Daraus kann man dann umgekehrt berechnen, wie groß L1 hochstens sein

142

Abbildung 10.6. Zufallsbewegung

−5 0 5 10−5

0

5

10

−5 0 5 10−4

−2

0

2

4

6

8

10

Linkes Bild: Verlauf der Bewegung. Rechtes Bild: Anfangs- und Endpositionen.

Kreise: Anfangsposition, Kreuze: Endposition

darf, damit die Sauerstoffkonzentration einen gewissen vorgegebenen Standard nichtunterschreitet.

Beobachten Sie auch einige qualitative Aussagen der Formeln: Die Voraussagen furC(t, x) und L(t, x) hangen nur vom gesamten Schadstoffzufluß L1 (dadurch indirekt vomSchadstoffausstoß L0) zur Zeit t − x/u ab. Die Zeitspanne x/u ist genau die Zeit, die derSchadstoff braucht, um mit der Stromung vom Auslaß zur Stelle x abzutreiben. Folgtman dem Schadstoff auf seinem Weg stromabwarts, so ergibt sich als Grenzwert

limt→∞

L(t, (t − t0)u) = 0,

limt→∞

C(t, (t − t0)u) = cS.

Weit stromabwarts regeneriert sich das Wasser auf die volle Sauerstoffsattigung, undaller Schadstoff wird abgebaut.

10.8. Modellierung von Diffusion und Dispersion.

Diffusion und Dispersion sind Zufallsbewegungen. Diffusion ist Bewegung der einzelnenMolekule in der Losung, Dispersion ist Bewegung der Tropfchen. Obwohl dasphysikalisch zwei vollig verschiedene Vorgange sind, bewirken beide eine Durchmischungder Losung. Im Lauf der Zeit verteilt sich jeder geloste Stoff gleichmaßig in der Losung.Wir werden im Folgenden nur von Diffusion reden, aber Dispersion wird genausomodelliert. Nur die Zahlenwerte der Parameter sind verschieden.

Abbildung 10.6 ist das Ergebnis einer Simulation einer Zufallsbewegung. 20 Teilchenstehen anfangs in einem Rechteck in der Mitte der gezeigten Flache angeordnet undunterliegen zufalligen Bewegungen. Die Bewegung des einzelnen Teilchens ist reinzufallig, aber insgesamt ist das Ergebnis, daß sich die Teilchen mehr oder mindergleichmaßig uber die Flache ausbreiten. Von den Regionen, wo viele Teilchen sitzen,


Abbildung 10.7. Die Wirkung der Diffusion

Durchgezogen: Konzentration zu Beginn. Strichliert: Konzentration nach einiger Zeit.

Die Pfeile symbolisieren Starke und Richtung des Teilchenflusses.

konnen sich naturgemaß 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 naturlich nicht die einzelnen Molekule oderTropfchen, sondern sie geht davon aus, daß die Gesamtwirkung der Diffusion vielerTeilchen einen Teilchenstrom (Fluß) von den Gebieten hoher Konzentration zu denGebieten niederer Konzentration ergibt. Mit Φ(t, x) bezeichnen wir die Teilchenmenge,die pro Flacheneinheit und Sekunde einen Querschnitt durchsetzt. Der einfachsteModellansatz ist das Ficksche Gesetz. Wir formulieren es nur fur den eindimensionalenFall.

Merksatz 10.8 (Ficksches Gesetz). Sei L(t, x) die Konzentration eines diffundierendenStoffes und Φ der Teilchenfluß auf Grund einer reinen Zufallsbewegung. Positives Φbedeutet Fluß in Richtung auf großere x, negatives Φ bedeutet Fluß in Gegenrichtung.Dann ist

(10.6) Φ(t, x) = −D∂

∂xL(t, x).

Dabei ist D eine Konstante, der Diffusionskoeffizient. Unter Umstanden kann D(x) auchvom Ort abhangen.

Das Ficksche Gesetz gibt zumindest qualitativ das typische Verhalten der Diffusionwieder: Wird die Konzentration in positive Richtung, also fur steigende x, großer, so istdie partielle Ableitung ∂

∂xL positiv. Damit wird Φ negativ, was einem Fluß in negative

Richtung, also in Richtung auf kleinere x, entspricht. Je großer dasKonzentrationsgefalle, also je großer die Ableitung der Konzentration, desto starker istauch der Fluß. Abbildung 10.7 gibt eine schematische Darstellung dieses Sachverhaltes.Mit Methoden der Wahrscheinlichkeitsrechnung laßt sich das Ficksche Gesetz fundierter

144

begrunden. Ubrigens gilt dasselbe Gesetz auch fur die Warmeausbreitung durchWarmeleitung in einem Korper.

Wir komplettieren jetzt unser Schadstoffmodell durch Einbau der Diffusion oderDispersion. Seien DS(x) und DO(x) die Diffusions- beziehungsweiseDispersionskonstanten fur Schadstoff und Sauerstoff. Die Diffusionskonstanten konnenauf Grund der verschiedenen Beweglichkeit der Molekule verschieden sein. DieDispersionskonstante ware in jedem Falle dieselbe, denn beide Stoffe werden mitdenselben zufallig schwirrenden Tropfchen transportiert. Zusatzlich zum Transport durchdie Stromung erhalten wir auf Grund der Zufallsbewegung fur beide geloste Stoffe diefolgenden Flusse:

ΦS(t, x) = −DS(x)∂

∂xL(t, x),

ΦO(t, x) = −DO(x)∂

∂xC(t, x).

Positives Φ zeigt nach unserer Konvention in Richtung auf großere x. Wir betrachtenwieder die “unendlich kleine” Zelle zwischen x und x + ∆x. Am oberen Zellrand, also ander Stelle x, durchsetzt der Fluß ΦS(t, x) den Querschnitt A(x) in Richtung auf hoherex, also ins Innere der Zelle. Das ergibt einen Schadstoffzufluß von A(x)ΦS(x) kmolSchadstoff 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 ergibteinen 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)∂

∂xL(t, x) + A(x + ∆x)DS(x + ∆x)

∂

∂xL(t, x + ∆x),

Sauerstoffbilanz durch Diffusion:

A(x)ΦO(x) − A(x + ∆x)ΦO(x + ∆x)

= −A(x)DO(x)∂

∂xC(t, x) + A(x + ∆x)DO(x + ∆x)

∂

∂xC(t, x + ∆x).

Dieser Term muß in der Modellbildung zu den anderen Mengenbilanzen der Zelle addiertwerden. Im Grenzubergang zu unendlich kleinen Zellen haben wir die Mengenbilanzdurch ∆x dividiert und dann den Grenzwert fur ∆x → 0 ausgewertet. Der Beitrag desDiffusionsterms ist dabei

lim∆x→0

A(x + ∆x)DS(x + ∆x) ∂∂x

L(t, x + ∆x) − A(x)DS(x) ∂∂x

L(t, x)

∆x

=∂

∂x

(

A(x)DS(x)∂

∂xL(t, x)

)

.

Dieser Term wird den anderen Termen der Advektions-Reaktions-Gleichung hinzugefugt.Die Diffusion des Sauerstoffs wird ebenso modelliert.


Das Schadstoffmodell unter Einbeziehung der Diffusion besteht also aus den beidenpartiellen Differentialgleichungen

A(x)∂

∂tL(t, x) =

∂

∂x

(

A(x)DS(x)∂

∂xL(t, x)

)

− ∂

∂x

(

A(x)u(x)L(t, x))

−kSA(x)L(t, x) + s(t, x),

A(x)∂

∂tC(t, x) =

∂

∂x

(

A(x)DO(x)∂

∂xC(t, x)

)

− ∂

∂x

(

A(x)u(x)C(t, x))

−kSA(x)L(t, x) + kOb(x)(cS − C(t, x)).

In den neuen Gleichungen werden L und C zweimal nach x partiell differenziert:Diffusion fuhrt auf partielle Ableitungen zweiter Ordnung in den Raumvariablen. Um dieGleichung zu behandeln, wurden wir auch je eine weitere Randbedingung fur Schadstoffund Sauerstoff brauchen, zum Beispiel die Konzentrationen am unteren Ende desbetrachteten Flußstuckes. Durch die Diffusion kann ja ein kleiner Teil des Schadstoffeswieder stromaufwarts wandern, sodaß auch vom unteren Ende her theoretischgeringfugige Einflusse auf die Schadstoffkonzentration im betrachteten Flußstuck zuerwarten sind. Weil man diese Konzentrationen ja nicht weiß, kann man sich in derPraxis behelfen, indem man ein ausreichend langes Flußstuck berechnet und am Endedie Schadstoffkonzentration mit 0 und die Sauerstoffkonzentration mit cS ansetzt. Dabeigehen wir davon aus, daß weit stromabwarts von der Schadstoffquelle die Wasserqualitatweitgehend regeneriert ist. Das Festlegen der Werte der gesuchten Funktionen am Randnennt man Dirichlet Randbedingungen. Ein anderer haufiger Typ sind sogenannteNeumann Randbedingungen, die die Ableitung der gesuchten Funktionen am Randvorschreiben. Wir konnten zum Beispiel verlangen, dass die Ableitung derKonzentrationen am unteren Ende null sind,

(10.7)∂

∂xL(t, xL) =

∂

∂xC(t, xL) = 0,

was bedeutet, dass sich die Konzentrationen an dieser Stelle raumlich nicht andern sollen.

Tabelle 10.3 zeigt zur Wiederholung noch einmal die Gestalt einer derReaktions-Diffusions-Advektionsgleichungen.

Durch die Einfuhrung der Diffusion andert sich das qualitative Verhalten deutlich.Stellen wir uns vor, daß die Fabrik zum Zeitpunkt 0 einen Augenblick lang eine gewisseSchadstoffmenge in den Fluß laßt, und dann ihren Abfluß wieder schließt. An der Stelle xlassen wir einen Beobachter mit seinen Meßsonden am Ufer warten. DasAdvektionsmodell ohne Diffusion sagt voraus, daß der Schadstoff zur Zeit x/u(Wegstrecke durch Stromungsgeschwindigkeit) die Stelle x erreicht. Bis zu diesemZeitpunkt nimmt der Beobachter nichts wahr. Zur Zeit x/u schlagen seine Instrumentekurz aus, einen Augenblick spater ist der Schadstoff vorbeigetrieben und der Beobachternimmt nichts mehr wahr.Diffusion und Dispersion bewirken, daß manche Schadstoffteilchen der großen Massevorauseilen, manche nachhinken. Das mathematische Modell wurde voraussagen, daß der

146

Tabelle 10.3. Reaktions-Diffusions-Advektionsgleichung

∂∂t

A(x)L(t, x) Zeitliche Anderungpartielle Ableitung nach t

=∂∂x

(

DS(x)A(x) ∂∂x

L(t, x))

Diffusionpartielle Ableitungen zweiter Ordnung nach x

− ∂∂x

A(x)u(x)L(t, x) Advektionpartielle Ableitung nach x von

Konzentration × Querschnitt

× Stromungsgeschwindigkeit

+s(t, x) Beitrag der Quellen−kSA(x)L(t, x) Beitrag der chemischen Reaktion

Gestalt hangt von der Art der Reaktion ab

Beobachter unmittelbar nach der Zeit 0 bereits einen minimalen Anstieg des Schadstoffeswahrnehmen wurde. Die Konzentration wurde weiter ansteigen, bis die große Massevorbeigetrieben ist, und wieder langsam gegen Null absinken. Je weiter stromabwarts derBeobachter steht, desto langsamer gehen Anstieg und Abfall der Konzentration vor sich.

Die Methode der charakteristische Kurven scheitert am Diffusionsmodell. Die Grundideedieser Methode war, daß ein Beobachter, der mit der Stromung mittreibt, dieStromungseffekte nicht wahrnimmt: Er ist immer vom selben Wasser und denselbenSchadstoffteilchen umgeben. Die Durchmischung durch Diffusion und Dispersion machtaber gerade diesen Effekt zunichte. Wenn man versucht, die totalen Ableitungen derbeobachteten Konzentrationen zu berechnen, heben sich nicht mehr alle partiellenAbleitungen weg.

Bemerkung. Wenn man bei der Modellbildung nur die Diffusion (oder Dispersion)berucksichtigt und alle anderen Effekte weglasst, kommt man (fur konstantenQuerschnitt) zur klassischen Diffusionsgleichung der Gestalt

∂

∂tL(t, x) =

∂

∂x

(

D(x)∂

∂xL(t, x)

)

, t > 0, x ∈ (0, xL),

fur eine Eigenschaft L, die sich durch kleine Zufallsbewegungen andert. Dies ist zumBeispiel auch fur die Warmeverteilung der Fall, man nennt diese Gleichung daher auchWarmeleitungsgleichung, wobei dann L die Temperatur ist. Um unter den Losungen derpartiellen Differentialgleichung eine eindeutig auszuwahlen, mussen Anfangsbedingungenfur t = t0 und Randbedingungen bei x = 0, x = xL gegeben sein.


Abbildung 10.8. Finite Elemente in einer Dimension

10.9. Numerische Losung von partiellen Differentialgleichungen.

Partielle Differentialgleichungen sind schwierig zu losen, und benotigen vor allem beizwei- und dreidimensionalen Problemen viel Speicherplatz und Rechenzeit. Wir habengesehen, daß sie ganz verschiedene Phanomene beschreiben konnen und ihre Losungenqualitativ ganz unterschiedliches Verhalten zeigen. Jede partielle Differentialgleichunghat sozusagen ihren ganz eigenen Charakter. Daher gibt es auch ganz verschiedeneAnsatze zur numerischen Losung, und Programmpakete werden oft fur spezielleAnwendungen erstellt. So gibt es zum Beispiel Programme, die auf die Berechnung vonStromungs- und Diffusionsvorgangen im Grundwasser spezialisiert sind, Programme zurBerechnung von Spannungen und Dehnungen in elastischen Strukturen, Programme zurSimulation der Luftstromung um eine Tragflache, und viele andere.

Eines der bekanntesten und erfolgreichsten Verfahren ist die Methode der finitenElemente. Sie baut darauf auf, die partielle Differentialgleichung durch eine ArtZellenmodell naherungsweise darzustellen. Ein Element besteht aus einem Raumstuck,sozusagen einer Zelle, und einer Funktion auf dieser Zelle, die noch durch geeigneteKoeffizienten modifiziert werden kann. Das gesamte untersuchte Gebiet wird in Elementezerlegt. Als Naherungslosungen werden nur Funktionen zugelassen, die auf jedereinzelnen Zelle durch ein Element dargestellt werden konnen. Auf diese Weise wird jedezulassige Naherungslosung durch endlich viele Koeffizienten beschrieben. Gegebenenfallskonnen noch zusatzliche Bedingungen die Auswahl der Naherungslosung einschranken.

Wir erklaren das an einem Beispiel fur ein eindimensionales Problem: EineDifferentialgleichung auf dem Intervall [0, l] soll gelost werden. Als Element betrachtenwir ein Intervall mit einer Geradengleichung y = kx + d als Funktion. Zwei Koeffizienten,namlich k und d, stehen zur Anpassung zur Verfugung. Wir teilen das gesamte Gebiet,also das Intervall [0, l], zum Beispiel in 20 kleine Teilintervalle. Als Naherungslosungenwerden nur Funktionen zugelassen, die auf jedem der Teilintervalle Geradenstucke sind.Daher wird jede Naherungslosung durch 40 Koeffizienten, namlich ki und di fur jedeeinzelne Zelle (mit Nummer i = 1, · · · , 20), beschrieben. Als Zusatzbedingung fordernwir, daß die gesamte Naherungslosung stetig ist, also die einzelnen Geradenstucke sichstetig treffen, wo die Zellen zusammenstoßen. Abbildung 10.8 zeigt ein einzelnes Elementund eine zulassige Naherungslosung. Statt Geradenstucken konnte man auch zumBeispiel quadratische oder vor allem kubische Parabeln heranziehen (sogenanntekubische Splines). Zweidimensionale Gebiete konnen zum Beispiel in Rechtecke oder

148

Dreiecke zerlegt werden, dreidimensionale Gebiete in Quader, Tetraeder oder Prismen(wie in Abbildung 10.3).

Die exakte Losung der partiellen Differentialgleichung ist naturlich normalerweise nichtgenau aus Elementen aufgebaut, aber wenn die Zellen ausreichend klein gewahlt sind,gibt es zulassige Naherungslosungen, die sich der exakten Losung recht genauanschmiegen. Unter den zulassigen Naherungslosungen sucht man eine, die die partielleDifferentialgleichung wenigstens bis auf einen moglichst kleinen Fehler genau erfullt. Dasfuhrt rechnerisch auf die Bestimmung von endlich vielen Koeffizienten als Losung einesMinimumproblems. Je kleiner die Zellen, desto genauer wird die Naherung, aber destomehr Koeffizienten sind zu bestimmen, und desto großer wird auch der Rechenaufwand.Die mathematische Durchfuhrung dieses Konzeptes im Detail ist naturlich ziemlichkompliziert. Ausgefeilte Programmpakete suchen nicht nur die Koeffizienten derNaherungslosung, sie bestimmen auch automatisch eine Einteilung des Gebietes in Zellen,bieten geeignete ODE-solver zur Auswahl an und stellen die Losung graphisch dar.

Die Methode der finiten Elemente ist – fur zeitabhangige Probleme – ein sogenanntessemidiskretes Verfahren, weil nur die raumliche Verteilung mit endlich vielen Elementendiskretisiert wird; dadurch entsteht fur eine partielle Differentialgleichung, in der auchzeitliche Anderungen modelliert sind, ein großes System von gewohnlichenDifferentialgleichungen, in dem nur mehr Ableitungen nach der Zeit vorkommen. Wiedieses System dann numerisch gelost wird, ist im Prinzip ein unabhangiges Problem,man sucht sich je nach Anforderungen einen passenden ODE-solver.

Anders als bei semidiskreten Verfahren erfolgt die Diskretisierung bei Verfahren derfiniten Differenzen vollstandig, indem auch die Zeit (in Abstimmung mit denRaumkoordinaten) explizit diskretisiert wird. Solche Verfahren kommen vor allem dannzur Anwendung, wenn es auf genaue Berechnung beinahe unstetiger Anderungen derZustande ankommt (z.B. Schockwellen). Durch die Anpassung der Diskretisierung anspezifische Eigenheiten des Modells (etwa entlang von charakteristische Kurven) kannnicht nur die Genauigkeit der Naherungslosungen verbessert werden, sondern auch dieEffizienz ihrer numerischen Berechnung.

10.10. Losungen des Reaktions-Advektions-Diffusions-Systems in Femlab.

Femlab ist ein Softwarepaket, das zur bequemen Behandlung von partiellenDifferentialgleichungen in Matlab verwendet wird (FEM=finite element method). Eskann Systeme von PDGen mit bis zu 3 Raumkoordinaten bearbeiten. Es enthalt Modulefur spezielle Probleme, etwa Warmeleitung, Stomungen, Elastizitat, Elektromagnetismus,etc. Wir verwenden fur unser Schadstoff-Modell ein allgemeines Modul fur eine lineareGleichung mit Zeitableitungen erster Ordnung und Raumableitungen erster und zweiterOrdnung.

Wenn wir in Matlab den Befehl femlab eingeben, erhalten wir die aus etlichen Fensternund klickbaren Menus bestehende graphische Arbeitsumgebung von Femlab. Hier konnen


wir den Modell-Typ wahlen und die raumliche Geometrie auf recht anschauliche Weiseeingeben (fur zwei oder drei Raumdimensionen mit Hilfe vorbereiteter Bausteine). EinGitter-Generator diskretisiert das Raumgebiet mit wahlbar Feinheit. Die Koeffizientender PDG und die gewunschten Anfangs- und Randbedingungen werden durch Eintragenvon Matlab-Formeln in Formulare programmiert, die je nach PDG-Typ gestaltet sind. Essteht eine Auswahl von finiten Elementen und ODE-solvern zur Verfugung. NachdemAlles vollstandig ausgewahlt und konsistent eingegeben wurde, startet die Simulation aufKnopfdruck. Zur Visualisierung der Ergebnisse werden verschiedene variierbareDarstellungs-Moglichkeiten angeboten.

Fur unserer Schadstoffmodell haben wir nur eine raumliche Dimension, wir klicken 1D.Dann legen wir fest, dass es 2 Zustandsgroßen gibt, die hier als Spaltenvektor u mit denKomponenten u1 und u2 bezeichnet werden. Wir wahlen als Geometrie einzusammenhangendes Intervall, 0 < x < xL = 10000. Als Einheiten verwenden wir Meter,Kilogramm, Sekunde, also ist der betrachtete Flußabschnitt zehn Kilometer lang. BeimAusfullen der Formulare fur die (immer passend angezeigten) Koeffizienten, legen wiruns darauf fest, dass mit u1 durchgehend die Schadstoffkonzentration und mit u2durchgehend die Sauerstoffkonzentration bezeichnet wird. Dann lassen sich die in denFormularen vorgegebenen Bezeichnungen eindeutig auf unsere Bezeichnungen abbilden(fur die Stromungsgeschwindigkeit weichen wir auf die Bezeichnung v aus). Es werdenauch Terme angeboten, die in unserem Modell garnicht vorkommen, deren Koeffizientensetzen wir einfach 0.

Zum Vergleich mit den Zellen-Simulationen (Abschnitt 10.3) wahlen wir die Parameterwie dort, insbesondere zeit- und ortsunabhangig. Es ist zweckmaßig, die Werte derbenotigten Konstanten in die Add/Edit Constants Tabelle einzutragen, sodass ihreNamen dann zur Bildung der Koeffizienten in den Formularen verwendet werden konnen.Die Schadstoff-Zufuhr modellieren wir der Einfachheit halber hier nicht alsRandbedingung, sondern als kurzzeitigen Eintrag vom Ufer am oberen Ende. Diesgeschieht mit der im Modell mit s(t, x) bezeichneten Inhomogenitat: Wir wahlen eineKonstante L0 und setzen s(t, x) = L0/(Av) fur 0 < x < 100, 0 < t < 360, und s(t, x) = 0sonst. Um orts- und zeitabhangige Koeffizienten zu programmieren, mussen wir in denFormeln die Variablen x und t verwenden. In die erste Zeile fur die Inhomogenitatschreiben wir also die Matlab Formel fur dieses s(t, x), namlichL 0*(0<x&x<100)*max(0,sign(360-t))/A/v. In der zweiten Zeile steht kO*b*cS (diesist die Inhomogenitat auf der rechten Seite der C-Gleichung). Durch Probieren wahlenwir die Paramter L0 und DS, DO so, dass hubsche, mit Abschnitt 10.3 vergleichbareSimulationsergebnisse entstehen. Dies fuhrt 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. AlsRandbedingungen nehmen wir L(t, 0) = 0, C(t, 0) = cS bei x = 0 und (10.7) bei x = xL.Wir legen den Simulationszeitraum 0 ≤ t ≤ 3600 fest, alles Andere lassen wir zunachstwie voreingestellt (z.B. ode15s als ODE-solver). Dann drucken wir gespannt auf dieTaste Solve und hoffen, dass nur wenige Fehlermeldungen kommen. Nach derenBeseitigung lauft die Simulation dann endlich und zeigt uns die Ergebnisse zu den vonuns gewahlten Zeitpunkten an, sagen wir t = 900, 1800, 2700, 3600, wie beim

150

Zellenmodell. Wir konnen die Zeitpunkte aber auch dichter wahlen, sodass wir dann mitder Taste Animate eine bewegte Graphik sehen, die die zeitliche Entwicklung je einer derKomponenten der Losung zeigt. Weitere graphische Moglichkeiten sind perspektivischeplots der Flachen L(t, x), C(t, x) uber der (t, x)-Ebene, oder die Farbung einerKomponente je nach der Große der anderen an derselben Stelle (t, x).

Die Abbildung zeigt die Ergebnisse fur die gerade geschilderten Parameter. Um glatteKurven zu erhalten, haben wir 120 Elemente (fur jede der beiden Komponenten)gewahlt. Ausserdem haben wir ausprobiert, mit welchem der wahlbaren Verfahren fur dieLosung der linearen Gleichungssysteme (fur die Koeffizienten der Naherungslosung) diebesten Rechenzeiten erzielt werden: Diese betragen auf einem Standard PC insgesamtetwa 6 Sekunden.

Die Kurven mit den Maxima sind die Schadstoffkonzentrationen L(ti, x), jene mit denMinima die zugehorigen Sauerstoffkonzentrationen C(ti, x). Die negativen Werte vonC(900, x) sind naturlich chemisch unsinnig und stammen daher, dass fur so grosseSchadstoffmengen das lineare Reaktionsmodell nicht zutreffend ist. Im Ubrigen schauendie Simulationsergebnisse aber vernunftig aus und entsprechen unseren Erwartungen: Dader Schadstoffzutritt nur fur kurze Zeit am linken Ende stattfindet, bildet sich einzunachst gut lokalisiertes L-Paket am linken Ende. Ohne Diffusion und Dispersion (kurzD&D) ware die Lange des Pakets 360 ∗ v Meter, bei Stromungsgeschwindigkeit v = 2 also720m. Durch D&D verbreitert sich das Paket, seine Hohe nimmt dadurch ab.Gleichzeitig verschiebt sich das Paket wegen der Stromung nach rechts. Zusatzlich nimmtdie L-Konzentration wegen der Oxidation ab, die auch C verbraucht. Dasdurchhangende C-Paket verbreitert sich durch seine Kopplung an L und mit seinereigenen D&D Konstanten DO. Wegen DS = 60 > 20 = DO bleibt das Maximum derL-Kurve starker hinter der Stromung zuruck als das Minimum der C-Kurve.

Die Verbreiterung der Pakete kann im Modell mit D&D im Gegensatz zu demZellen-Modell (das nur Advektion modellierte) also plausibel gemacht werden. Trotzdemruhrt die Verbreiterung der Pakete in den Simulationsergebnissen nicht nur aus denD-Termen, sondern sie beruht auch auf systematischen Problemen mit finitenElementen: Wie im Zellenmodell kommt es zu einer Verschmierung von lokalkonzentrierten Anderungen. Dies konnen wir zeigen, indem wir im selben Femlab-Modelldie Diffusionskoeffizienten einfach 0 setzen, DS = DO = 0. Dann mussten die Fronten desSchadstoff-Pakets genau mit der Stromungsgeschwindigkeit wandern, also wegen derspeziellen Wahl von s(t, x) die vordere zum Zeitpunkt t genau die Stelle x = 100 + v ∗ tpassieren, und die hintere Front x = v ∗ (t − 360) (vgl. die exakte Losung derAdvektionsgleichung im Abschnitt uber charakteristische Kurven). Bei genauererUntersuchung der zugehorigen Simulationsergebnisse stellen wir jedoch fest, dass diePaket-Langen grosser als 720 sind und mit der Zeit zunehmen, selbst bei DS = DO = 0.

Es ware naiv, sich ein Simulations-Programm zu wunschen, das alle PDGen zur vollenZufriedenheit im Griff hat, zu unterschiedlich sind die Eigenschaften der verschiedenenGleichungstypen. Legt man Wert auf beinahe unstetige Wellenfronten, sind angepasstefinite Differenzen Verfahren am geeignetsten. Andererseits ist software wie Femlab fur


Abbildung 10.9. Simulationsergebnisse des Reaktions-Advektions-Diffusions-Modells in Femlab

0

5

10

15x 10

−4

t=90

0

0

2

4

6x 10

−4

t=18

00

0

2

4

6x 10

−4

t=27

00

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 100000

2

4

6x 10

−4

t=36

00

Anwender, die sich nicht mit der Numerik partieller Differentialgleichung auseinandersetzen wollen, und die moglichst 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 systematischenumerische Tests verringert werden.

152

11. Diskrete deterministische Modelle

Angenommen, der Zustand eines Systems wird durch n Zustandsgroßen beschrieben, diein einem Zustandsvektor x = (x1, . . . , xn) zusammengefasst werden. Falls sich das Systemmit der Zeit t kontinuierlich verandert, und wir dies mathematisch beschreiben wollen,brauchen wir Funktionen der Zeit x1(t), . . . , xn(t) fur t ≥ t0, mit einem Anfangszeitpunktt0. Bisher haben wir nur Modelle betrachtet, deren Zustandsvektor x(t) koninuierlich vonder Zeit abhangt (im Falle raumlich verteilter Zustande sind fur jedes t dieKomponenten von x(t) selbst Funktionen, die den Raumkoordinaten Werte zuordnen).

Es kann aber vorteilhaft und in vielen Situationen durchaus angebracht sein, denZustand des Systems nur zu diskreten Zeitpunkten t0 < t1 < t2 < . . . zu modellieren. Einmathematisches Modell, das dies leistet, konnte zum Beispiel von der folgenden Gestaltsein:

(11.1)x0 gegebenxk+1 = F (xk), k = 0, 1, 2, . . .

Dabei ist y = F (x) eine Funktion, die einem n-Vektor x einen n-Vektor y zuordnet undxk ist der Zustandsvektor zum Zeitpunkt tk. Es ist denkbar (und in manchen Fallenunabdingbar), dass auf der rechten Seite von (11.1) die Zeit explizit vorkommt und dasszur Berechnung von xk+1 auch fruhere Zustande xk−1, xk−2, . . . benotigt werden. In (11.1)haben wir uns aber auf Modelle beschrankt, in denen xk+1 aus dem Vorganger-Zustandallein berechenbar ist. Insbesondere in der Form xk+1 = xk + G(xk) nennt man so einModell eine Differenzengleichung erster Ordnung. Die Funktion G beschreibt dieAnderung des Zustands zwischen diskreten Zeitpunkten. Ausgehend vom Anfangszustandx0 bestimmt das Modell sukzessive die Zustande x1, x2, . . . zu den diskreten Zeitpunktent1, t2, . . . . Die Werte der Zustandsgroßen sind im Allgemeinen nicht diskret, dh. dieKomponenten von xk sind im Allgemeinen keine ganzen Zahlen, sondern beliebig reell.

In diesem Kapitel beschaftigen 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 zuganglich, man braucht nur F zuprogrammieren und wiederholt anzuwenden. Es ist aber auch moglich, eine qualitativeAnalyse der durch x0 und F erzeugten Folge x1, x2, . . . durchzufuhren; es geht dabei umGleichgewichtslagen, deren Stabilitat, periodische Losungen, asymptotisches Verhaltenbei k → ∞, etc. Im Besonderen ist es chaotische Dynamik, die an Hand von einfachendiskreten deterministischen Modellen studiert werden kann: das beruhmteste Beispieldafur ist die Differenzengleichung xk+1 = αxk(1 − xk) mit der logistischen FunktionF (x) = αx(1 − x) fur α > 1 und x ∈ R.

11.1. Ein diskretes deterministisches Epidemie-Modell.

Ein fruhes 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, nichtansteckenden Immunen (wie im SIR-Modell des ersten Semesters). Wir betrachten diePopulation zu diskreten Zeitpunkten tk in gleichbleibenden Abstanden von zwei Wochen,sagen wir. Mit xk bezeichnen wir die Anzahl der Ansteckenden zum Zeitpunkt tk, fur die

11. DISKRETE DETERMINISTISCHE MODELLE 153

Anzahl der Individuen in den Klassen y, z schreiben wir yk, zk. Wir nehmen an, dass dieGesamtgroß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 ganzenZahlen sind. Dies ist insbesondere bei großen Klassen zulassig, weil es dann nicht aufeinzelne Individuen ankommt. Vielmehr konnen wir uns in diesem Fall die Werte vonxk, yk, zk als Anteile an der Gesamtbevolkerung vorstellen, die in geeigneten Einheitenangegeben sind.

Sei 0 < p < 1 die gleichbleibende Wahrscheinlichkeit, dass ein ansteckbares Individuumeffektiven Kontakt mit einem ansteckenden Individuum hat, dh. dass es wahrend einesBeobachtungsintervalls angesteckt wird. Dann ist die Wahrscheinlichkeit, dass es keineneffektiven Kontakt mit einem der ansteckenden Individuen hat, gleich 1 − p, und dieWahrscheinlichkeit, dass es den effektiven Kontakt mit allen yk Ansteckenden vermeidetist gleich (1 − p)yk . Wir setzen daher

xk+1 = (1 − p)ykxk.

Dies ist fur sich allein als Modell fur xk+1 nicht brauchbar, weil ja yk benotigt wird. ZurAbkurzung fuhren wir einen Parameter a ein, und zwar a = −ln(1 − p) > 0. Dann iste−a = 1 − p. Der Anteil der wahrend des Intervalls neu Infizierten, namlich (1 − e−ayk)xk

ist zum Zeitpunkt tk+1 in der Klasse der Angesteckten y. Aufgrund der Heilung verlassenein Teil der Angesteckten y und kommen zur Klasse z. Wahrend [tk, tk+1] betreffe dies(1 − b)yk Individuen, wobei wir einen konstanten Parameter 0 < b < 1 ansetzen. Damiterhalten wir

(11.3)xk+1 = e−aykxk

yk+1 = (1 − e−ayk)xk + byk

zk+1 = zk + (1 − b)yk

Die Klasse z ist strikt an die Klasse y gekoppelt und spielt fur die weitere Entwicklungder Epidemie keine Rolle mehr. Zur Kontrolle uberprufen wir die Eigenschaft (11.2)indem wir die drei Zeilen summieren: xk+1 + yk+1 + zk+1 = xk + yk + zk, die Summebleibt im Modell unverandert.

Das Modell enthalt zwei positive reelle Zahlen, a und b, als Parameter. Obwohl bei derErstellung des Modells mit Wahrscheinlichkeiten argumentiert wurde, ist es keinprobabilistisches Modell. Der Parameter a ist fest vorgegeben, das Modell bestimmtxk+1, yk+1, zk+1 eindeutig aus xk, yk, zk und ist daher deterministisch.

In einer interaktiven Programmiersprache wie Matlab kann, ausgehend von Startwertenfur 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 ausgefuhrt werden. Wir wollen dies hier nicht weiterverfolgen, stellen aber eine qualitative Frage: Wird die Zahl der Angestecktenzunehmen ? Mit anderen Worten: Sagt das Modell den Ausbruch einer Epidemie voraus ?

154

Man kann die Frage auch so formulieren: wenn y0 = 1 (dies ist durch Wahl der Einheitenimmer erreichbar), wird dann y1 großer als 1 sein ? Laut Modell gilt fur die Differenz

y1 − y0 = (1 − e−a)x0 + b − 1.

Dies ist großer als 0 genau dann, wenn

x0 > (1 − b)/(1 − e−a) = x∗.

Der kritische Wert x∗ ist ein Schwellenwert fur die Große der Klasse der Ansteckbaren.Die Signifikanz eines solchen Schwellenwerts (threshold) ist unabhangig vonmathematischen Modellen aus Felddaten bekannt, und wird als Test fur die Gute desModells verwendet. Wir sehen, dass der Zusammenhang zwischen den Parametern unddem Schwellenwert wenigstens plausibel ist: je kleiner (1 − b), desto mehr Ansteckendeuberdauern den Beobachtungszeitraum und desto kleiner ist x∗. Je großer (1 − e−a),desto großer ist die Wahrscheinlichkeit der Ansteckung p und desto kleiner ist x∗. Oballerdings die Großenordnungen passen und mit realistischen, moglichst unabhangig zueruierenden Parametern erreichbar sind, ist eine andere Frage. Ein weiteres Kriterium,das vom Modell erfullt wird, ist die Forderung nach einer wohlbestimmten Bedeutungder Parameter: a und b sind anschaulich interpretierbar und konnen als Systemgroßenvon der Sanitatsverwaltung gemessen und beeinflusst werden. Sind die Parameter einmalfestgelegt, kann das Modell im Prinzip dazu verwendet werden, die Starke der Epidemievorherzusagen (bei allem Vorbehalt, der solch einfachen Modellen entgegenzubringen ist).

Zusammenhang mit dem kontinuierlichen SIR Modell. Bei der Modellbildung zuKermack-McKendrick und zum SIR Modell werden dieselben Argumente verwendet, alsosollte es Ubergange 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, dassxk = S(kh), yk = I(kh) und zk = R(kh), zumindest in guter Naherung fur kleineZeitintervalle h = tk+1 − tk.

Im Modell (11.3) machen wir den folgenden linearen Ansatz fur die Parameter a und b:

a = λh, b = 1 − γh.

Dies ist plausibel, denn fur h → 0 sollten die Populationen gleich bleiben. Mit t = khsieht (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),

alsoS(t + h) − S(t) = (e−λhI(t) − 1)S(t) = −λhI(t)S(t) + O(h2),I(t + h) − I(t) = −γhI(t) + (1 − e−λhI(t))S(t)

= −γhI(t) + λhI(t)S(t) + O(h2),R(t + h) − R(t) = γhI(t).

Dabei haben wir die Reihenentwicklung fur die Exponentialfunktion verwendet und furderen Terme hoherer Ordnung das O-symbol geschrieben (insbesondere gilt O(h2)/h → 0


fur h → 0). Wenn wir durch h dividieren erhalten wir links Differenzenquotienten. Lassenwir dann h gegen 0 gehen, verschwinden die Terme hoherer Ordnung und es bleibt:

ddt

S(t) = −λI(t)S(t),ddt

I(t) = λI(t)S(t) − γI(t),ddt

R(t) = γI(t).

Dies ist das von fuher bekannte SIR-Modell ohne Zuwanderung (β = 0) und ohneSterblichkeit (µ = 0). Zur Kontrolle uberprufen wir die Eigenschaft (11.2) indem wir diedrei Zeilen summieren: d

dt(S(t) + I(t) + R(t)) = 0, die Summe bleibt im Modell

unverandert.

Wegen des engen Zusammenhangs von gewohnlichen Differentialgleichungen undDifferenzengleichungen, hat der/die Modellierer/in unter gegebenen Umstanden dieMoglichkeit, den Modell-Typ zu wahlen, je nach Fragestellung, Vorlieben undVerfugbarkeit von Methoden.

11.2. Zellular-Automaten.

Wenn wir einen Vorgang in Zeit und Raum diskret beschreiben mochten, dann mussenwir auch das betrachtete Raumgebiet in endlich viele Teile zerlegen, oder, was aufdasselbe hinauslauft, nur endlich viele seiner Punkte berucksichtigen. Betrachten wirzunachst den eindimensionalen Fall, also ein Raumintervall [0, L] der Lange L. In diesemIntervall wahlen wir N + 1 aquidistante Stutzstellen, namlich, mit ∆x = L/N ,xj = j∆x, j = 1, . . . , N − 1 im Inneren und die beiden Rander x0 = 0, xN = L.

Viele dynamische Prozesse sind durch lokale Regeln bestimmt, dh. dass nur raumlichbenachbarte Teile des Systems einander beeinflussen (dies ist eine enorm generelleAussage und sollte uberdacht werden). Denken wir zum Beispiel an die Stromung einerFlussigkeit in einem Kanal [0, L] mit Stromungsgeschwindigkeit v von 0 nach L. DieFlussigkeit besitze eine zeitlich und raumlich veranderliche Eigenschaft C(t, x) (z.B. dieKonzentration einer gelosten Substanz, wir lassen hier aber chemische Reaktionen oderDiffusion weg und betrachten nur die reine Advektion). Wir geben eine diskreteraum-zeitliche Beschreibung der Entwicklung von C, indem wir den Wert Ck

j von C zujedem Zeitpunkt tk = k∆t an jeder Stelle xj = j∆x angeben. Stimmen wir ∆t und ∆xaufeinander ab, und zwar so, dass ∆x = v∆t. Dann ist der Wert von C zur Zeit tk+1 ander Stelle xj genau der Wert von C zur Zeit tk an der Stelle xj−1. Die Entwicklung derdiskreten Werte von C folgt also folgender lokaler Regel:

Ck+1j = Ck

j−1, j = 1, . . . , N, k = 0, 1, . . .

Die Anfangswerte C0j zur Zeit t0 sowie die Randwerte Ck

0 bei x0 mussen vorgegeben sein.Dann ist durch die lokale Regel der Wert von C fur alle Zeiten tk an allen Stellen xj

eindeutig bestimmt. Das Zeit-Raum-Gitter (tk, xj) zusammen mit den Randbedingungenund den Regeln zur Berechnung der Zustande Ck+1

j (Update) aus den Zustanden Ckj

bezeichnet man als Zellular-Automat (engl. cellular automaton, dt. auch zellularerAutomat). Die raumliche Diskretisierung (die wir gerade mit Hilfe von Stutzstellen

156

beschrieben haben) denkt man sich dabei als Unterteilung des Raumgebiets in Zellenund Ck

j als Mittelwerte von C wahrend [tk, tk+1) in [xj−1, xj).

Der soeben beschriebene Zellular-Automat lasst sich als eine finite DifferenzenDiskretisierung der Advektionsgleichung entlang der Charakteristiken x = vtinterpretieren. Auch Diffusion und Reaktion kann mit Zellular-Automaten simuliertwerden. Fur viele wichtige Zellular-Automaten besteht jedoch kein Bezug zu einerpartiellen Differentialgleichung.

Insbesondere gibt es Zellular-Automaten, deren Zellen nur eine endliche Anzahl vonZustanden annehmen; wenn es nur zwei mogliche Zustande pro Zelle gibt, spricht manvon boolschen Automaten. Ein Zellular-Automat ist gepragt durch die Regeln fur dasUpdate, wobei die Anderung des Zustands einer Zelle sich typischerweise aus denZustanden ihrer Nachbarzellen ergibt (und zwar eindeutig im deterministischen Fall).Dazu wird insbesondere auch festgelegt, welche Zellen zur Nachbarschaft gehoren. Inunserem Stromungsbeispiel benotigten wir nur die links liegende, westliche Nachbarzelle.Andere eindimensionale Prozesse, etwa die Diffusion, benotigen die Berucksichtigungbeider Nachbarzellen (westlich und ostlich). Denkbar und durchfuhrbar sind auchgrossere Nachbarschaften, etwa je zwei westliche und ostliche Zellen. Im Fall vonzwei-oder dreidimensionalen raumlichen Gittern gibt es eine Vielzahl geometrischerMoglichkeiten fur Nachbarschaften.

Die lokalen Regeln zum Update konnen oft eindeutig und leicht verstandlich formuliertwerden, auch ohne mathematische Formalisierung. Die Ubersetzung in eine globaleFunktion F wie in (11.1) ist vergleichsweise viel komplizierter, bringt wenig Erkenntnisund ist fur eine numerisch-graphische Implementierung am Computer gar nichtnotwendig. Ein deterministischer Zellular-Automat gilt als vollstandig gegeben, wenn zubeliebigen Anfangszustanden gehorende Simulationsergebnisse eindeutig festgelegt sind.

Der bekannteste Zellular-Automat ist zweifellos Game of Life, der intensivnumerisch-experimentell aber auch analytisch untersucht wurde. Es handelt sich um eimstilisiertes raumlich zweidimensional verteiltes Populations-Modell, in dem der Tod auflokale Unter- und Uberbevolkerung folgt. Uberleben und Geburt setzen sehreingeschrankte lokale Bevolkerungsdichten voraus. Betrachten wir ein Rechteck, das inviele kleine gleich große rechteckige Zellen unterteilt ist. Jede der Zellen kann zu jedemZeitpunkt tk, k = 0, 1, . . . entweder tot oder lebendig sein. In einer graphischenDarstellung des Zustands zu einem festen Zeitpunkt werden tote Zellen leer gelassen,lebendige eingefarbt oder mit einem Symbol besetzt. Ausgehend von einer gegebenenAusgangsverteilung wird der Zustand einer Zelle zur Zeit tk+1 aus ihrem und ihrerNachbarn Zustand zur Zeit tk bestimmt. Die Regel lautet:

eine Zelle ist zur Zeit tk+1 genau dann lebendig, wenn zur Zeit tksie selbst und zwei ihrer Nachbarn, oder drei ihrer Nachbarzellen lebendig waren.

Zu den Nachbarzellen einer Zelle zahlen die 8 nachstliegenden Zellen. Also: Hat einelebendige Zelle zu wenige oder zu viele lebendige Nachbarn, stirbt sie. Eine lebendigeZelle mit 2 oder 3 lebendigen Nachbarn uberlebt. Eine tote Zelle mit 3 lebendigenNachbarn wird lebendig.


Abbildung 11.1. Ein glider bewegt sich in vier Zyklen um einen Schrittnach Sud-West

Zeichnen Sie auf kariertem Papier einige Beispiele oder eine Sequenz von Schritten. EinBeispiel ist in Abbildung 11.2 zu sehen, der sogenannte glider, eine Konfiguration, diesich in 4 Schritten um je eine Zelle nach Suden und Westen verschiebt (uberprufen 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 verandernde) Konfigurationen.

Da wir mit einem endlich großen Rechteck arbeiten, gibt es Randzellen, die zunachstweniger als 8 Nachbarzellen haben. Eine mogliche Vervollstandigung besteht inperiodischen Randbedingungen: fur das Update der Zellen in der obersten Zeileverwendet man die unterste Zeile als nordlichen Nachbarn, fur die rechteste Spalte wirddie linkeste Spalte als ostlicher Nachbar verwendet, und analog fur die unterste Zeile unddie linkeste Spalte (durch dieses Zusammenkleben der Rander macht man aus demRechteck einen Schlauchring). Ein glider, der den sudlichen Rand uberschreitet kommtdann im Norden wieder herein, usw.

Wenn zum Beispiel ein glider in eine statische/periodische Konfiguration eindringtkonnen die erstaunlichsten Entwicklungen ausgelost werden.... Viel besser als durch eineBeschreibung lernen Sie Game of Life durch Experimentieren kennen: siehe Matlab>>demo. Naturlich finden Sie mit einer Suchmaschine im Internet jede Menge sites mitGame of Life, etwa www.bitstorm.org/gameoflife. Nur als Andeutung zeigen wir inAbbildung 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 dielebendigen Zellen nach 100 Zyklen.

Die Faszination, die Game of Life auf viele seiner Betrachter und Erforscher ausubt,beruht wohl darauf, dass eine einfache, in einem Satz formulierbare Regel, derartkomplexe und reichhaltige Entwicklungsszenarien in sich birgt. So gibt es z.B. shooter,das sind periodische, ortsfeste Konfigurationen, die alle 30 Zyklen einen glider ausstossen.Darauf aufbauend wurde gezeigt, dass es sich reproduzierende Konfigurationen gibt: istdas Zellenfeld genugend groß, dann entstehen aus so einer Konfiguration zwei identischesolche. Naturlich wurden als Varianten zu Game of Life andere Regeln untersucht undRegeln systematisch klassifiziert oder mit Parametern stufenlos geandert. Bei der ersten,legendaren Tagung zu diesen Themen wurde der Begriff “artificial life” eingefuhrt.

158

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

Bemerkung. In diesem Abschnitt des Skriptums hat sich - unangekundigt - eineErweiterung des Begriffs “Mathematisches Modell” ereignet. Wahrend wir bisher untereinem mathematischen Modell eine vereinfachende, idealisierende Abbildung einesTeilbereichs der Realitat verstanden, erheben die Zellular Automaten des artificial lifenicht den Anspruch, in der Natur existierende Systeme abzubilden oder moglichst gutnachzuahmen. Ausserdem haben wir das Modell gar nicht vollstandig mathematischausformuliert: fur Simulationen genugt die programmierbare Beschreibung einesAlgorithmus und die graphische Darstellung der Simulation. Vielleicht sollte man solcheModelle einfach Simulations-Modelle nennen.

Die Moglichkeit, mit solchen Modellen zu arbeiten, eroffnete sich der Wissenschaftinsbesonders seit es Computer gibt. Zusatzlich zu Theorie und Experiment mit naturlichexistierenden Objekten gibt es nun Computer-Simulationen in Welten unserersprachlichen Phantasie. Dies ist auch fur sehr grundlegende Fragen von großerwissenschaftlicher Bedeutung. So wird an Hand von Simulationen untersucht, ob (undwie) in einer stilisierten Ursuppe chemischer Bestandteile Makromolekule als Vorlauferdes Lebens entstehen konnen. Prozesse, die der biologischen Evolution ahnlich sind,finden auch in recht einfachen virtuellen, digitalen Universen statt. Soziologische undkulturelle Entwicklungen werden mit vereinfachten Gesellschaften simuliert (multi agentmodels). Scheinbar (?) komplexe und komplizierte Phanomene wie Lernen, Anpassungoder Selbstorganisation konnen mit Simulations-Modellen imitiert werden, und dies kannwesentlich zum Verstandnis der Grundlagen solcher Prozesse beitragen.

12. DISKRETE PROBABILISTISCHE MODELLE 159

12. Diskrete probabilistische Modelle

Probabilistische Modelle - immer haufiger auch stochastische Modelle genannt -beinhalten Elemente des Zufalls, in der einen oder anderen Form. Dies kann aus vielerleiGrunden angebracht sein. Zum Einen kann es sich um die Modellierung von Systemenhandeln, in denen der Zufall explizit eine grosse Rolle spielt. Dies ist etwa beiGlucksspielen der Fall, der Ubertragung von Krankheitserregern oder der Mutationgenetischer Information. Zum Anderen kann es sich um Modelle handeln, die dieAnderung von Modellgroßen nicht bis ins letzte (deterministische) Detail beschreiben,sondern wie vom Zufall gesteuert behandeln. Dies ist etwa der Fall bei der Modellierungunregelmaßiger Schwankungen von Umweltgroßen wie Niederschlag, Temperatur etc.

In probabilistischen Modellen ist die Entwicklung der Zustandsgroßen nur gemaßgewisser Wahrscheinlichkeitsverteilungen bestimmt. Bei Simulationen werden dieseGroßen mittels Zufallsgeneratoren festgelegt, sodass auch die Simulationsergebnisse vomZufall abhangen und von Simulation zu Simulation verschieden sind. Im Rahmen vonsogenannten Monte-Carlo Simulationen fuhrt man eine Serie von Einzelsimulationendurch und macht statistische Aussagen uber die Ergebnisse.

Es gibt stochastische Modelle, in denen zeitlich kontinuierliche Modellgroßen vom Zufallabhangen und ein Kontinuum von Werten annehmen konnen. In dieser Vorlesung werdenaber keine kontinuierlichen Zufallsprozesse behandelt: Unsere Zustandsgroßen werdennur endlich viele diskrete Werte annehmen, und wir werden nur zeitlich diskrete Modellebehandeln. Der Ubergang zwischen den Zustanden wird mit Wahrscheinlichkeitenmodelliert. Das Paradebeispiel fur ein System, das so modelliert werden kann, ist einSpielwurfel, der vor bzw. nach jedem Wurf in einem von 6 moglichen Zustanden ist. Fureinen perfekten Wurfel ist die Ubergangswahrscheinlichkeit zwischen allen Zustandengleich, namlich 1

6.

12.1. Markov-Ketten.

Eine Markov-Kette oder Markov-Prozess n-ter Ordnung ist ein zeitlich diskretesdynamisches Modell bestehend aus einer Menge {S1, . . . , Sn} von Zustanden und n2

Ubergangswahrscheinlichkeiten pij, i = 1, . . . , n, j = 1, . . . , n. Der Prozess kann zu jedemZeitpunkt k nur in einem der n Zustande sein. Ist der Prozess zum Zeitpunkt k imZustand Si, dann ist pij die Wahrscheinlichkeit, dass der Prozess zum Zeitpunkt k + 1 imZustand Sj ist. Im einfachsten Fall hangen die Ubergangswahrscheinlichkeiten nicht vonder Zeit k oder vom Zustand zur Zeit k ab (lineare MK mit konstantenUbergangswahrscheinlichkeiten). Der Prozess startet in einem gegebenenAnfangszustand.

Beispiel. Wir wollen das Wetter in einer bestimmtem Stadt mit Hilfe dreier moglicherZustande S, W und R, namlich Sonne, Wolken, Regen beschreiben. Je einer dieserZustande soll fur jeden Tag k zutreffen. In Abbildung 12.1 sind die (fur diese Stadttypischen) Ubergangswahrscheinlichkeiten mit Hilfe von Pfeilen und dazugeschriebenenWahrscheinlichkeiten angegeben. Dh. zB. falls es heute wolkig ist, dann gibt es eineChance von 25%, dass es morgen wieder wolkig ist, eine Chance von 25%, dass es morgen

160

Abbildung 12.1. Wettermodell mit Ubergangswahrscheinlichkeiten

regnet und eine Chance von 50%, dass es morgen sonnig ist. Wir konnen dieUbergangswahrscheinlichkeiten zu einer Tabelle zusammenfassen:

S W RS 1/2 1/2 0W 1/2 1/4 1/4R 0 1/2 1/2

Zu einem heutigen Wetter finden wir in der entsprechenden Zeile dieWahrscheinlichkeiten fur 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 dasModell 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 jederRunde k des Spiels gewinnt A mit der Wkt p einen Jeton von B, oder aber A verlierteinen Jeton an B, namlich mit der Wkt q = 1 − p (dh. in jeder Runde wechselt genau einJeton den Besitzer). Sobald einer der beiden Spieler alle Jetons hat, kann sich nichtsmehr andern. Betrachten wir ein Spiel mit 4 Jetons: es gibt 5 Zustande, dir wir einfachmit der Anzahl der Jetons bezeichnen, die in Besitz von A sind. Es gibt also dieZustande 0,1,2,3 und 4. Das Pfeile-Diagramm sieht wie in Abbildung 12.1 aus, diezugehorige Tabelle ist

0 1 2 3 40 1 0 0 0 01 q 0 p 0 02 0 q 0 p 03 0 0 q 0 p4 0 0 0 0 1

Da es aus jedem der Zustande heraus immer weitergehen soll, sollte die Summe der auseinem Knoten des Diagramms herausfuhrenden Wahrscheinlichkeiten uberall 1 sein. Esist klar, wie das gleiche Modell fur eine andere Anzahl von Jetons ausssieht.


Abbildung 12.2. Spielermodell mit Ubergangswahrscheinlichkeiten

Allgemein gehort zu einer Markov Kette n-ter Ordnung eine n × n-Matrix P vonUbergangswahrscheinlichkeiten

P =

p11 p12 · · · p1n

p21 p22 · · · p2n...

......

pn1 pn2 · · · pnn

Alle Elemente der Matrix sind reelle Zahlen zwischen 0 und 1, 0 ≤ pij ≤ 1. Ferner sindalle Zeilensummen gleich 1, denn jeder Zustand muss in irgendeinen der moglichenZustande ubergehen. Eine quadratische Matrix mit diesen beiden Eigenschaften heissteine stochastische Matrix. Und ein Vektor x ∈ R

n heisst Wahrscheinlichkeitsvektor, wennalle seine Komponenten nichtnegativ kleiner gleich 1 sind und ihre Summe 1 ist. Wennder Zeilenvektor x ein Wahrscheinlichkeitsvektor ist und P eine stochastische Matrix,dann ist das Vektor mal Matrix Produkt xP wieder ein Wahrscheinlichkeitsvektor:

n∑

j=1

(xP )j =n

∑

j=1

n∑

i=1

xi(P )ij =n

∑

i=1

n∑

j=1

xi(P )ij =n

∑

i=1

xi

(

n∑

j=1

pij

)

=n

∑

i=1

xi = 1.

Daher sind stochastische Matrizen die naturlichen linearen Transformationen zwischenWahrscheinlichkeitsvektoren. Wir haben, wie bei Markov Ketten ublich, die Vektoren alsZeilen verwendet und multiplizieren diese von rechts mit quadratischen Matrizen(genausogut konnten wir die fur MK unubliche Schreibweise mit Spaltenvektoren undden transponierten Matrizen verwenden). Fur das Element einer Matrix A in der i-tenZeile und in der j-ten Spalte schreiben wir wie ublich (A)ij.

Wenn ein Markov Prozess mit Ubergangsmatrix (P )ij = pij im Zustand Si startet, dannist 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 einWahrscheinlichkeitsvektor, die i-te Zeile von P , und entsteht aus P durch Multiplikationmit dem Vektor

ei = (0 · · · 010 · · · 0)

von links; in ei steht der Einser an der i-ten Stelle. Wie gross ist die Wkt p(2)ij , dass Si in

zwei Schritten nach Sj ubergeht ? Falls dazwischen sicher Sk lage, ware die Wkt

p(2)ij = pkj. Aber die Wkt, dass Sk dazwischen liegt ist nur pik, und wir mussen die

162

Summe uber alle moglichen Wege bilden,

p(2)ij =

n∑

k=1

pikpkj = (P 2)ij.

Hier bezeichnet P 2 das Quadrat der Matrix P , das dadurch entsteht, dass man die n × nMatrix P mit sich selbst multipliziert, P 2 := P · P (dies ist nicht die Matrix, die durchQuadrieren der einzelnen Elemente entsteht). Also ist P 2 die 2-Schritt Ubergangsmatrix,deren Element (P 2)ij die Wkt angibt, dass Si in zwei Schritten nach Sj ubergeht. Wennwir diese Uberlegung induktiv fortsetzen, sehen wir dass die Matrix Pm, die dadurchentsteht, dass man P m mal mit sich selbst multipliziert, die m-Schritt Ubergangsmatrixist: (Pm)ij ist die Wkt, in m Schritten von Si nach Sj zu kommen.

Anstatt die m-Schritt Ubergangsmatrix auszurechnen, konnen wir auch einen iterativenProzess als diskretes dynamisches System formulieren: sei x(k) ein Zeilenvektor, dessenj-te Komponente die Wkt angibt, dass im kten Schritt der Zustand Sj vorliegt. Dannergibt 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) legtdieses diskrete dynamische System die Wahrscheinlichkeitsverteilung x(k) fur alle Zeitenk fest. Das Resultat fur die Wktverteilung x(m) ausgehend von einer Anfangsverteilungx(0) ist naturlich das gleiche wie mit der m-Schritt Ubergangsmatrix Pm, denn

x(m) = x(m − 1)P = x(m − 2)P · P = · · · = x(1)Pm−1 = x(0)Pm.

Achtung, x(k) ist nicht der k-te Zustand des Markov Prozesses, mogliche Zustande sindnur S1, . . . , Sn. x(k) enthalt aber die Wahrscheinlichkeiten, dass die MK zur Zeit k einender Zustande annimmt. Zur Simulation einer Zustandsentwicklung musste man in jedemSchritt einen Zustand gemaß den Wkten in x(k) auswahlen (siehe Monte CarloSimulation im nachsten Abschnitt).

Fur eine andere Interpretation von x(k), die vielleicht etwas anschaulicher ist, stellen wiruns eine große Anzahl N gleichartiger Systeme vor, die alle mit dem gleichen Modellbeschrieben werden. Sei zB. N eine Anzahl von Stadten, fur die wir das gleicheWetter-Modell anwenden. Jede Stadt wird mit der gleichen MK modelliert, wobei derMarkov Prozess in jeder Stadt unabhangig von den anderen Stadten ablauft. Auch wennalle Stadte im gleichen Anfangszustand starten, werden sie sich im Lauf der Zeitunterschiedlich entwickeln. Das Modell liefert Erwartungswerte fur eine Verteilung derZustande zu spateren Zeitpunkten: zum Zeitpunkt k sind x(k)jN Stadte im Zustand Sj

zu erwarten. Obwohl das Modell keine verlassliche Aussage uber den Zustand in einereinzelnen Stadt macht, sagt es fur die ganze Menge von Stadten eine Verteilung voraus.

Abgesehen von der Berechnung der m-Schritt Ubergangswahrscheinlichkeiten in Pm oderder Iteration (12.1) ist vor Allem die qualitative Analyse von Markov Ketten vonInteresse. Folgende Fragen liegen auf der Hand und konnen fur lineare Markov Kettenmit konstanten Ubergangswahrscheinlichkeiten beantwortet werden: Langzeitverhalten,Grenzzustand, mittlere Anzahl von Schritten, um spezielle Zustande zu erreichen, Wkt,


dass ein spezieller Zustand je erreicht wird, etc. Ein wichtiges Resultat in diese Richtungformulieren wir als Satz:

Satz 12.1. Sei P die Ubergangsmatrix einer regularen Markov Kette, das heisst, dass furein m > 0 alle Elemente der Matrix Pm positiv sind. Dann gilt

a) Es gibt einen eindeutig bestimmten Wahrscheinlichkeitsvektor p mit derEigenschaft

pP = p.

b) Fur beliebigen Anfangszustand x(0) = ei gilt

limm→∞

eiPm = p,

insbesondere ist die Grenzverteilung p unabhangig vom Anfangszustand.c)

limm→∞

Pm = P :=

pp...p

.

Der Beweis dieses Satzes ist eine Anwendung der Vektorrechnung, der Vektor p ist einlinker Eigenvektor von P zum Eigenwert 1. Die Aussagen b) und c) sind Aussagen uberdas Langzeitverhalten der Kette. Wenn wir uns wieder vorstellen, dass N Kettengestartet werden, dann konvergiert nach vielen Schritten die Anzahl der Ketten, die sichin Zustand Sj befinden, nach pjN , unabhangig von der Anfangsverteilung der Zustande.

Beispiel. Fur die obige Wetter Matrix P sind alle Elemente von P 2 positiv, wie Sieleicht nachrechnen konnen. Also handelt es sich um eine regulare MK. Durch einfacheinteraktive Iterationen zB. in Matlab kommt man schnell zu der Vermutung, dass

P =

0.4 0.4 0.20.4 0.4 0.20.4 0.4 0.2

.

Uberzeugen Sie sich, dass p = (0.4 0.4 0.2) tatsachlich ein linker Eigenvektor von Pzum Eigenwert 1 ist. Die Aussage des Satzes bedeutet folgende Langzeit-Statistik lautModell: 40% der Tage sonnig, 40% wolkig, 20% regnerisch.

Wir gehen hier nicht auf andere qualitative Fragestellungen in Bezug auf Markov Kettenein, ebensowenig auf Markov Ketten mit unendlich vielen diskreten Zustanden.

Die beiden einfachen Beispiele dieses Abschnitts machen vielleicht den Eindruck, dassMarkov Ketten triviale mathematische Modelle fur Entwicklungen sind, die auch ohnemathematische Modellierung vorhergesagt werden konnen. Doch gibt es sehr wohl MKModelle mit konstanten Ubergangswahrscheinlichkeiten, die wesentliche Beitrage zumVerstandnis von Zufallsprozessen geleistet haben: Jahrelang spekulierten Experten derGentechnologie, warum die einem Bakterium zusatzlich zum Typ A implantiertenPlasmide vom Typ B im Lauf der Generationen wieder separiert werden. Plasmidetragen zusatzliche Erbinformation ausserhalb des Chromosoms des Bakteriums. Bei dernaturlichen Zellteilung werden die N vorhandenen Plasmide zunachst kopiert und dann

164

zufallig auf die beiden Tochterzellen verteilt, und zwar erhalt jede Tochter N Plasmide.Im Lauf von vielen Zellteilungen entstehen aus einer Population, in der jedes derBakterien gleichviele Plasmide vom Typ A wie B besitzt, zwei Unterpopulationen vonBakterien, die jeweils nur einen Typ Plasmide besitzen. Als man schließlich dieZellteilungen als Markov Kette modellierte, wobei die Wahrscheinlichkeiten strengkombinatorisch hergeleitet wurden, stellte sich heraus, dass die Plasmid Separation einProzess ist, der sich von selbst einstellt, allein auf Grund binomialerWahrscheinlichkeitsverteilungen (random genetic drift). Irgendwelche Spekulationenbzgl. Plasmid-Inkompatibilitaten u. dgl. haben sich mit diesem Modell erubrigt. DasModell wurde gut in diese Vorlesung passen und ist nicht besonders aufwandig, wirbetrachten nun aber eine andere, nichtlineare Markov Kette mit zustandsabhangigenUbergangswahrscheinlichkeiten.

12.2. Infektion einer Familie.

Wir betrachten nun wieder eine Population ansteckbarer (S), ansteckender (I) undgeheilter (R) Individuen, wie in Abschnitt 11.1. Die Reihenfolge von Infektion undHeilung ist durch S → I → R gegeben. Hier wollen wir nun eine kleine Population mitganzzahligen Klassen modellieren und die Zufallsprozesse bis ins Detail verfolgen. Seienk = 0, 1, . . . die Zeitpunkte, zu denen das System beobachtet wird, der Abstand zwischenzwei 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. Wiedernehmen wir an, dass die Große der Gesamtpopulation konstant N bleibt:

Sk + Ik + Rk = N.

Die moglichen Zustande des Modells sind also alle Tripel mit nichtnegativenganzzahligen Komponenten, deren Summe N ist. Zum Beispiel fur N = 4 gibt es 15solche Zustande. Die Modellbildung besteht darin, Formeln fur die Berechnung derUbergangswahrscheinlichkeiten zwischen diesen Zustanden zu entwickeln.

Sei p die Wkt des effektiven Kontakts eines der S mit einem der I wahrend einesBeobachtungszeitraums. Dann ist 1 − p die Wkt, dass ein S den effektiven Kontakt miteinem der I vermeidet. Wenn es Ik Infizierte gibt, dann ist also die Wkt, dass ein S denKontakt mit allen Infizierten vermeidet gleich

qk = (1 − p)Ik ,

und qnk 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 dieWahrscheinlichkeit, dass beim Ubergang k → k + 1 gerade n der Ansteckbaren jedenKontakt mit Infizierten vermeiden und gerade Sk − n aber effektiven Kontakt haben alsdas Produkt

qnk (1 − qk)

Sk−n


berechnet werden. Aus der Kombinatorik ist bekannt, wieviele Moglichkeiten es gibt, ausSk Indiviuen jeweils n auszuwahlen. Diese Anzahl ist gegeben durch denBinomialkoeffizienten

(

Sk

n

)

=Sk!

n!(Sk − n)!,

wobei fur jede nichtnegative ganze Zahl K, K Fakultat gegeben ist durchK! = K · (K − 1) · (K − 2) . . . 3 · 2 · 1. Somit ist die Wkt fur Sk+1 = n gleich

(

Sk

n

)

qnk (1 − qk)

Sk−n.

Der Einfachheit halber nehmen wir nun an, dass die Infektionsdauer in jedem Fall mitSicherheit gleich dem Beobachtungszeitraum ist. Das heisst, jedes infizierte Individuumist genau 2 Wochen ansteckend und wechselt dann in die Klasse der nicht mehransteckenden und nicht mehr ansteckbaren Geheilten. Zusammenfassend erhalten wir

Pr(Sk+1 = n|Sk und Ik) =

(

Sk

n

)

qnk (1 − qk)

Sk−n(12.2)

Falls Sk+1 = n, dann Ik+1 = Sk − n und Rk+1 = Rk + Ik.

Hier haben wir eine Schreibweise aus der Wahrscheinlichkeitsrechnung verwendet: mitder Probabilitat Pr(A|B) bezeichnet man die Wahrscheinlichkeit von A unter derBedingung B, dh. die Wahrscheinlichkeit fur das Eintreten von A unter derVoraussetzung, dass B erfullt ist. Die Ausdrucke auf der rechten Seite von (12.2) sindnur sinnvoll, wenn n ≤ Sk. Andere Falle brauchen auch nicht modelliert zu werden, weildie Anzahl an Ansteckbaren nur abnehmen oder hochstens gleichbleiben kann.

Die Gleichungen (12.2) heissen das Reed-Frost Modell und geben dieUbergangswahrscheinlichkeiten fur eine Markov Kette. Zu jedem Zeitpunkt k ist derZustand ein Vektor mit drei nichtnegativen ganzzahligen Komponenten, deren Summe Nist. Die Ubergangswahrscheinlichkeiten sind hier aber keine Konstanten, sondern hangenstark nichtlinear vom Zustand ab, und zwar von Sk und, uber qk, von Ik. Eine qualitativeAnalyse des Reed-Frost Modells ist daher nicht moglich. Wohl aber konnen Beipiele derEntwicklung simuliert werden und die Auswertung vieler Beispiele fuhrt zu statistischenAussagen uber das bzw. mit Hilfe des Modells.

Monte Carlo Simulationen. Wegen der Annahme, dass Infizierte nach zwei Wochengeheilt sind, ist es sicher, dass die Infektion der Familie nach spatestens 2N Wochenvoruber ist. Ist namlich einmal Ik = 0, dann kann im Modell kein Krankheitsfall mehrauftauchen. Wie lange wird die Infektion aber wirklich dauern, wievieleFamilienmitglieder werden krank werden und wieviele werden die Sache gesunduberstehen? Schauen wir uns ein Beispiel an, das mit Reed-Frost Wahrscheinlichkeitenarbeitet.

Die Familie hat N = 5 Mitglieder und zum Zeitpunkt k = 0 ist genau eines davoninfiziert, also S0 = 4, I0 = 1, R0 = 0. Die Krankheit sei wenig ansteckend, sagen wirp = 0.05. Mit Hilfe der Formel (12.2) berechnen wir die Wahrscheinlichkeiten fur allemoglichen Nachfolgezustande (unmoglich ist zB. S1 = 5 oder R1 6= 1) und tragen diese ineine Tabelle ein:

166

Beobachtung WahrscheinlichkeitS1 I1 R1 dieser Beobachtung

a 4 0 1 0.814b 3 1 1 0.172c 2 2 1 0.0135d 1 3 1 0.000495e 0 4 1 0.000005

Da ja irgendeiner der moglichen Nachfolgezustande angenommen wird, muss die Summeder Wahrscheinlichkeiten 1 sein, was tatsachlich der Fall ist. Zur Veranschaulichungzeichnen wir eine (nicht ganz maßstabgetreue) den Wahrscheinlichkeiten entsprechendeUnterteilung des Intervalls [0,1], und benennen die Teilintervalle wie die Zeilen derTabelle: Nun ziehen wir eine Zufallszahl z ∈ [0, 1], das heisst wir verwenden einen

a 0.814

1 0.999995

e d

0.9995

c

0.986

b

0

Mechanismus oder eine Methode, zum Beispiel einen Zufallszahlengenerator, der unterden Zahlen in [0,1] eine auswahlt, wobei alle der Zahlen die gleiche Chance haben,ausgewahlt zu werden. Lage z zB. im Intervall c, dann ware S1 = 2 ausgewahlt. Durchdiese Vorgangsweise wird erreicht, dass genau eine der Zeilen ausgewahlt wird und zwarmit der Wahrscheinlichkeit, die der letzten Spalte der Tabelle entspricht. (Wir konntendie Teilintervalle auch in anderer Reihenfolge in der Zeichnung anordnen, wichtig ist nur,dass ihre Langen gleich den Wahrscheinlichkeiten sind, dass alle vorkommen, und dass[0,1] luckenlos uberdeckt ist.) Wegen der geringen Ansteckungswahrscheinlichkeit p ist esbei Weitem am wahrscheinlichsten, dass I1 = 0, womit die Infektion schon nach dem demersten Beobachtungszeitraum beendet ist.Angenommen wir ziehen z = 0.96, unser z liegt also im Teilintervall b, und das bedeutetdie Auswahl

S1 = 3 und I1 = 1.

Dies war der erste Simulationsschritt und die Simulation kann nun analog fortgesetztwerden: wir berechnen die Wahrscheinlichkeiten Pr(S2 = n|S1 = 3 und I1 = 1) fur n ≤ 3nach (12.2):


a 3 0 2 0.857b 2 1 2 0.135c 1 2 2 0.0075d 0 3 2 0.0005

Die dazupassende Unterteilung von [0,1] in Teilintervalle der Lange a, b, c, d hat dieTrennungspunkte bei 0.875, 0.992, 0.9995. Wir ziehen eine Zufallszahl, sagen wir z = 0.9.Damit ist b ausgewahlt, also

S2 = 2 und I2 = 1.


Von hier ausgehend finden wir folgende WahrscheinlichkeitenPr(S3 = n|S2 = 2 und I2 = 1) fur n ≤ 2:


a 2 0 3 0.902b 1 1 3 0.095c 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, 3Mitglieder der Familie sind erkrankt, 2 haben die Ansteckungsgefahr gesunduberstanden.

Dies war eine Monte Carlo Simulation. Klarerweise arbeitet man am Computer undschreibt sich ein passendes Programm zur Berechnung der Wahrscheinlichkeiten. Es isteffizienter, in jedem Schritt die Zufallszahl z zuerst zu bestimmen (wenigRechenaufwand) und die Wahrscheinlichkeiten nur soweit zu berechnen, bis erstmals zkleiner ist als die Summe der Wahrscheinlichkeiten.

Man kann nun viele Beispiele rechnen – gleichwertig damit ist die Denkweise, dass manviele Familien mit dem gleichen Modell modelliert – und dann eine Statistik derErgebnisse ziehen. Auf diese Weise konnen wir etwa die durchschnittliche Dauer derKrankheit, oder die zu erwartende Anzahl an Infizierten bestimmen. Zur Illustrationbetrachten wir die Krankheit in Familien mit 5 Mitgliedern unter drei verschiedenenAnnahmen uber die Ansteckungswahrscheinlichkeit: p = 0.1, p = 0.3 und p = 0.5. Injedem der drei Falle fuhren wir 40 Simulationen durch, die alle mit S0 = 4, I0 = 1starten. Wir bewerten die Starke der Infektion einer Familie mit Hilfe von S∞, das ist dieAnzahl von (noch) Ansteckbaren zu dem Zeitpunkt k, zu dem die Krankheit aus derFamilie verschwindet, also Ik = 0 wird (dies ist spatestens nach 5Beobachtungszeitraumen der Fall). In obigem ausfuhrlichen Beispiel mit p = 0.05 hattenwir S∞ = 2. Ergebnisse solcher Simulationen sind in folgender Tabelle aufgelistet, wobeidie Anzahl von Beipielen angegeben ist, die mit S∞ endeten.

S∞ p=0.1 p=0.3 p=0.50 1 12 231 1 8 122 5 3 13 11 3 14 22 14 3

Die Spaltensumme ist jeweils 40, das ist die Anzahl der gerechneten Beispiele (Familien)fur jedes der p. Fur p = 0.1 endete nur eine der 40 Simulationen ohne verbliebeneAnsteckbare, aber in 22 Beispielen wurde kein weiters Familienmitglied angesteckt. Furp = 0.5 wurden in 23 der Beispiele alle Familienmitglieder angesteckt und in nur 3Beispielen war die Gefahr nach dem ersten Beobachtungszeitraum vorbei. Interessant istdie Anderung der Verteilung auf die funf Zeilen, wenn die Ansteckungswahrscheinlichkeitzunimmt: das Maximum ist fur p = 0.1 eindeutig bei S∞ = 4, das heisst, dass keine

168

weitere Erkrankung in einer Familie zu erwarten ist, die mit einem Infizierten startet.Fur p = 0.5 ist das Maximum bei S∞ = 0, das heisst, dass am wahrscheinlichsten dieganze Familie erkrankt. Aber fur p = 0.3 macht das Modell keine eindeutige Aussage: esgibt zwei Maxima und die Anzahl der 12 Falle, in denen alle erkranken ist etwa gleichgroß wie die Anzahl der 14 Falle, in denen nur ein Mitglied erkrankt. Gibt es in der Nahevon p = 0.3 vielleicht so etwas wie einen Schwellwert ? Dies konnten wir imdeterministischen Modell in 11.1 analytisch beantworten, im probabilistischen Modellmussen wir uns mit Simulationen und Statistik begnugen.

Implementierung in Matlab. Es ist nicht schwierig, das Reed-Frost Modell in Matlabumzusetzen und damit Simulations-Reihen auszufuhren. Der Befehl Y=rand(m,n) erzeugteine m×n Matrix Y deren Elemente im Intervall (0,1) gleichverteilt sind (rand steht furrandom, es gibt auch die Funktion randn fur Normalverteilung). Fur unsere Zweckegenugt es rand immer wieder ohne Argument aufzurufen, jedesmal wird ein Skalarerzeugt und die Folge der so erzeugten Skalare ist de facto gleichverteilt im Intervall (0,1).Die folgende function im file ReedFrost.m fuhrt fur die Argumente p, S0, I0,

num of fam in einer Schleife num of fam Monte Carlo Simulationen durch undretourniert einen Zeilen-Vektor S 00, dessen jte Spalte die Anzahl der Familien enthalt,in denen S∞ = j − 1 Familienmitglieder nie angesteckt wurden.

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

S_00(S+1)=S_00(S+1)+1;

end


Fur die Große der Familie wird also S0+I0 genommen (R0 = 0). DieBinomialkoeffizienten werden mit Hilfe der Matlab function nchoosek berechnet, dieWahrscheinlichkeiten im Vektor a gespeichert und die Unterteilungspunkte von [0,1]stehen in b. Die Zufallszahlen werden in der Zeile z=rand generiert, erst nach derBerechnung aller Wahrscheinlichkeiten: da die 120 Simulationen fur die S∞ Tabelle nur0.3 Sekunden Rechenzeit verbrauchen, kummern wir uns hier nicht um Effizienz. Wennman dieses Programm im command window je drei mal fur p=.1, p=.3, p=.5 aufruft,erhalt man zB. folgende Ergebnisse:

>> S_00 = ReedFrost(.1,4,1,40)

S_00 =

1 0 3 9 27

>> S_00 = ReedFrost(.1,4,1,40)

S_00 =

0 1 2 11 26

>> S_00 = ReedFrost(.1,4,1,40)

S_00 =

0 2 3 8 27

>> S_00 = ReedFrost(.3,4,1,40)

S_00 =

10 12 5 3 10

>> S_00 = ReedFrost(.3,4,1,40)

S_00 =

10 9 0 6 15

>> S_00 = ReedFrost(.3,4,1,40)

S_00 =

12 7 5 7 9

>> S_00 = ReedFrost(.5,4,1,40)

S_00 =

30 3 1 1 5

>> S_00 = ReedFrost(.5,4,1,40)

S_00 =

30 5 2 2 1

>> S_00 = ReedFrost(.5,4,1,40)

S_00 =

27 9 1 1 2

>>

Das heisst z.B. fur die letzte Simulation mit p = 0.5, und 40 Familien, deren jede mitvier S und einem I starten: in 27 Familien wurden 0 Familienmitglieder nicht angesteckt,in 9 Familien wurde 1 Mitglied nicht angesteckt, in je 1 Familie wurden 2 bzw. 3Mitglieder nicht angesteckt, und in 2 Familien wurden 4 Mitglieder nicht angesteckt.

Bemerkung. Da der output eines numerischen Programms eindeutig vorbestimmt ist,ist es prinzipiell unmoglich, einen echten Zufallsgenerator mit Hilfe deterministischermathematischer Operationen zu programmieren. Es gibt zwar software, die aufphysikalisches Rauschen im Computer zugreift, fur viele Simulationszwecke genugen aber(im Prinzip vorhersagbare) Generatoren wie rand in Matlab. Um immer wieder andere

170

Zahlenfolgen zu bekommen, kann man den Zufallszahlengenerator am Beginn einer jedenFolge immer wieder anders initialisieren, etwa mit Hilfe von Datum und Uhrzeit.

12.3. Ein probabilistischer Zellular-Automat.

Als Beispiel fur ein raumlich verteiltes System, dessen Elemente zufalligenEntwicklungen ausgesetzt sind, denken wir an einen Waldbrand. Ein einfaches aberqualitativ reichhaltiges Modell dazu ist der ZA Forest Fire. Das Rechteck aller Zellenstellt einen Wald dar und jede der Zellen kann zu den diskreten Zeitpunkten in einemvon drei moglichen Zustanden sein:

T : Baum (tree)F : brennender Baum (fire)E: leerer Platz (empty site)

Diese Zustande konnen nur in der Reihenfolge T → F → E → T → . . . durchlaufenwerden, wobei der Zustand F immer nur einen Zeitschritt lang anhalt. Die Verweildauereiner Zelle in E und T hangt uber vorgegebene Wahrscheinlichkeiten p und f vom Zufallab. 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 mussen wir noch festlegen, was wir unter einer Nachbarschaft verstehen. Anders alsbei Game of Life zahlen wir hier nur jene 4 Zellen zur Nachbarschaft einer Zelle, die andie Seiten angrenzen; jene mit nur einem gemeinsamen Eckpunkt zahlen wir nicht dazu.Als Randbedingung betrachten wir die Situation, dass das Rechteck von leeren Platzenumgeben ist, auf denen keine Baume wachsen konnen.

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 = [zeros(m,1) F zeros(m,1)];

FE = [zeros(1,n+2); FE; zeros(1,n+2)];

FE = FE(no,:)+FE(:,e)+FE(s,:)+FE(:,w); % +FE(no,e)+FE(s,e)+FE(s,w)+FE(no,w);

FN = FE(2:m+1,2:n+1);

NFF = (T & FN);


NFN = ~FN;

NFL = (T & NFN & (rand(m,n) <= f));

E = (~T & ~F);

NT = (E & rand(m,n) <= p);

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 erhalt als Argumente die Ausgangsmatrizen T und F, die ausser Nullen nurEinser an den entsprechenden Platzen enthalt, ferner die Wahrscheinlichkeiten p und f ,sowie die Anzahl der zu simulierenden Zyklen. Alle Operationen werden alsMatrix-Gleichungen geschrieben, dies ist wesentlich effizienter als Schleifen fur dieIndizes. Die Zeilen- und Spaltenzahl der Matrizen sind m,n. In den Vektoren no, e, s,

w werden die Indizes der nordlichen und sudlichen Zeilen jeder Zeile und der ostlichenund westlichen Spalten jeder Spalte einer m + 2, n + 2 Matrix gespeichert, und zwarperiodischen Randbedingungen entsprechend. Dies wird dann zur Berechnung derAnzahl brennender Nachbarn eines jeden Platzes in FE verwendet, wobei FE dieVergroßerung von F mit leeren Platzen am Rand ist. Die mit einem Kommentarzeichen% unschadlich gemachte Summe +FE(no,e)+FE(s,e)+FE(s,w)+FE(no,w) wurde auchdie Ecknachbarn dazuzahlen. FN enthalt also die Anzahl brennender Baume in derNachbarschaft eines jeden der m × n Platze. Mit Hilfe des logischen und-Operators &erhalt NFF Einser an den Platzen mit brennenden Nachbarn wo ausserdem ein Baumsteht, dieser beginnt nun zu brennen. Mit Hilfe der Negation˜enthalt NFN jene Platze, inderen Nachbarschaft kein brennender Baum ist. Falls dort ein Baum steht wird diesermit Wahrscheinlichkeit f zu brennen beginnen, was in NFL mit Hilfe einer in jedemZeitschritt neu gebildeten Zufallsmatrix rand(m,n) und der Kleiner-Gleich-Relation <=

realisiert ist. Leere Platze E befinden sich genau dort wo kein Baum und kein brennenderBaum steht. Mit Hilfe einer zweiten Zufallsmatrix und der Wahrscheinlichkeit p wirdfestgelegt, ob dort ein neuer Baum NT wachst. Alte uberlebende Baume OT sind jene, inderen Nachbarschaft kein brennender Baum ist und die nicht gerade vom Blitz getroffenwerden. Insgesamt gibt es zum nachsten Zeitpunkt die Baume OT+NT und die brennendenBaume NFF+NFL. Der Rest ist Graphik.

Am besten lernen Sie die Dynamik dieses ZA durch eigene Experimente kennen, von derSeite teaching propst.html konnen Sie eine Matlab function forestfire.m

herunterladen. Als matten Abglanz der dynamischen Entwicklungen zeigt Abbildung12.3 drei Momentaufnahmen, die mit diesem Programm erzeugt wurden. Zu sehen istjeweils ein Wald mit 100 × 100 Zellen, wobei die weissen Zellen leere Platze sind, dunkle(rote) Zellen brennende Baume und hellgraue (grune) Zellen nicht brennende Baumedarstellen. Fig. a zeigt einen Ausgangszustand gleichverteilter T und F , wobei 10% derPlatze mit T besetzt sind, T=rand(m,n) <= 0.1, und 10% der Platze brennen, sofern sienicht schon von einem T besetzt sind, F=(rand(m,n) <= 0.1) & ˜ T. Wir wahlen

172

Abbildung 12.3. Zustande des Zellular-Automaten Waldbrand

a b

c d

p = 0.05 und f = 0, dh. es gibt zunachst keinen Blitzschlag. In diesem Fall sind dreiverschiedene Szenarien moglich: Erstens, die Anfangsverteilung ist so, dass die F schonnach wenigen Schritten verschwinden, ab dann fullt sich der Wald nur mehr mit grunenZellen. Zweitens, die Anfangsverteilung ist so, dass sich eine oder mehrere Feuerfrontenausbilden, und zwar dann und dort, wo der Baumbestand schon dicht genug dafur ist. Indiesem Szenario kommt es typischerweise zu großen Ringwellen, das sind Karo-formigeKetten dunkler Zellen, die sich in Richtung Rand ausdehnen und eine Zone leerer Platzehinter sich lassen, die sich nur langsam wieder mit grunen Baumen fullt. Fig. b zeigt dieSituation nach etwa 100 Zyklen, ausgehend von a. Zwei großer werdende Ringwellen sindin der Mitte zusammengewachsen und im Begriff das Rechteck uber den Rand zuverlassen. Wenn die zuruckbleibenden kleinen Brande verloschen, wachst nach derFeuersbrunst wieder einen dichter werdender Baumbestand heran. Fur p = 0.05 und der10%igen Anfangsverteilung ist es aber, drittens, haufig der Fall, dass sich nach dergroßen Anfangswelle viele kleinere Brande bilden und nie mehr ein feuerfreier Zustandeintritt. Es entstehen wandernde Muster aus Baumen, Feuern und leeren Platzen,ahnlich wie in Fig. c. Da sich diese Muster (insbesondere auch in b) aus einemzufallsverteilten Zustand entwickeln, liegen Beispiele fur die Selbstorganisation einesSystems vor. In ahnlich gebauten deterministischen Automaten kommt es zur


Ausbildung periodischer Raum-Zeit-Strukturen (Spiralwellen), die hier aber durchZufallsprozesse gestort sind. Zu Spiralwellen, also sich nicht geradlinig oder konzentrischfortpflanzenden Strukturen, kommt es dann, wenn jede Zelle angeregt und damit dieNachbarn anregend sein kann (F ), oder ruhend (T ), wobei dazwischen ein refraktarerErholungs-Zustand (E) angenommen werden muss, der nicht angeregt werden kann.

Fig. c zeigt eine Momentaufnahme fur p = 0.05 und f = 0.0001, die am Ende einerlangeren Simulation gemacht wurde. In diesem Stadium entstehen und vergehenStrukturen verschiedener Großenordnungen. Obwohl sich die Verteilungen dauerndandern spricht man von einem quasistationaren Zustand, weil das Gesamt-Musterqualitativ unverandert bleibt. Im Unterschied zum Fall f = 0 drittens, scheint diesesStadium fur f = 0.0001 aber etwas unruhiger und unregelmaßiger, und ausserdem kanndas Feuer wegen anhaltender Blitzgefahr nie endgultig verschwinden. Erhohen wirschließlich die Blitzschlag-Haufigkeit auf f = 0.01, dann stellen sich im Lauf der Zeitbewegte Verteilungen wie in Fig. d ein (wieder p = 0.05 ). Die Strukturen verschwinden,das Nachwachsen der Baume ist zu sehr gestort, die Gleichverteilung von leeren Platzenverhindert die Ausbildung von Feuerfronten.

Wie in Abschnitt 11.2, beruht die Faszination an ZA des Typs Forest Fire nicht auf derAnwendbarkeit auf reale Waldbrande, sondern auf der Moglichkeit, Mechanismen derMusterbildung zu erkunden, insbesondere ihre Abhangigkeit von stufenlos einstellbarenParametern. So spricht man bei der Ausbildung quasistationarer Strukturen (c und,unter Umstanden, b) von selbstorganisierter Kritizitat (SOC = self-organized criticality).Und bei einer grundsatzlichen Anderung des qualitativen Verhaltens spricht man vonPhasenubergangen: Erstarrung im Fall b, wenn alle F verschwinden, oder Verdampfungim Fall d, wenn alle Strukturen aufgelost werden. Lebendige strukturelle Vielfalt gibt esnur in einem kleinen Bereich aufeinander abgestimmter Parameter, dort entstehen sieallerdings von selbst, auch aus strukturlosen Ausgangsverteilungen.

9. qualitative analyse mathematischer modelle 103...• das gleichgewicht x0 heißt stabil, wenn...

Documents