Grundlagen der mathematischen Modellierung · Grundlagen der Mathematischen Modellierung Skript zur...

Grundlagen

der Mathematischen Modellierung

Skript zur Vorlesung fur Studierende der Masterstudiengange

Marine Umweltwissenschaften und Umweltmodellierung

an der Carl von Ossietzky Universitat Oldenburg

im Wintersemester 2009/2010

Cora Kohlmeier

25. Januar 2010

GMM WS2009/2010

Inhaltsverzeichnis

1 Einfuhrung 1

1.1 Was ist Modellbildung – mathematische Modellierung . . . . . . 1

1.2 Was ist ein Modell ? . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2.1 Merkmale von Modellen . . . . . . . . . . . . . . . . . . 2

1.2.2 Arten von Modellen . . . . . . . . . . . . . . . . . . . . . 2

1.2.3 Modell-Simulation . . . . . . . . . . . . . . . . . . . . . 2

2 Empirische Modelle 3

2.1 Parameteranpassung . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 Lineare Regression . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2.1 Methode der kleinsten Quadrate . . . . . . . . . . . . . . 7

2.3 Logarithmischer Zusammenhang am Beispiel der Karzinogenese 9

2.4 Verallgemeinerung der linearen Regression ersten Grades . . . . 11

2.4.1 Fourieranalyse . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Diskrete Modelle 13

3.1 Bakterienwachstum . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.2 Fibonacci und die Kaninchen . . . . . . . . . . . . . . . . . . . 20

3.2.1 Schritte beim Aufstellen eines Modells . . . . . . . . . . 21

3.2.2 Formalisierung des Kaninchenmodells . . . . . . . . . . . 21

3.2.3 Losung/Simulation des Systems . . . . . . . . . . . . . . 23

3.3 Bienen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.4 Zusammenfassung der Wachstumstypen . . . . . . . . . . . . . . 25

3.5 Diskrete 2D-Modelle . . . . . . . . . . . . . . . . . . . . . . . . 26

3.5.1 Fibonacci und der goldene Schnitt . . . . . . . . . . . . . 28

3.5.2 Uberlagerung von Losungen – Anfangswertproblem . . . 28

4 Aufstellen eines Modells 31

4.1 Verbale Beschreibung, Ziele . . . . . . . . . . . . . . . . . . . . 35

4.2 Zustandsvariable, Prozesse . . . . . . . . . . . . . . . . . . . . . 35

4.3 Randbedingungen, Antriebe . . . . . . . . . . . . . . . . . . . . 39

4.4 Prozessbeschreibungen . . . . . . . . . . . . . . . . . . . . . . . 41

4.5 Parameterliste und Parameterwerte . . . . . . . . . . . . . . . . 43

4.6 Programmierung, Simulation . . . . . . . . . . . . . . . . . . . . 47

3

5 Modellierung mit Differentialgleichungen 49

5.1 Die Wachstumsrate . . . . . . . . . . . . . . . . . . . . . . . . . 50

6 Systeme mit einer Variablen 51

6.1 Aufstellen einer DGL am Beispiel auslaufender Gefaße . . . . . 51

6.1.1 Der auslaufende Zylinder . . . . . . . . . . . . . . . . . . 51

6.1.2 Der auslaufende Trichter . . . . . . . . . . . . . . . . . . 52

6.2 Das logistische Wachstum . . . . . . . . . . . . . . . . . . . . . 55

6.3 Stationare Zustande . . . . . . . . . . . . . . . . . . . . . . . . 56

6.3.1 Stabilitat stationarer Zustande . . . . . . . . . . . . . . 56

6.4 Skalierung der logistischen Differentialgleichung . . . . . . . . . 59

6.5 Wachstums und Sterbeprozesse . . . . . . . . . . . . . . . . . . 60

7 Systeme mit 2 Variablen 63

7.1 Rauber-Beute-Modelle . . . . . . . . . . . . . . . . . . . . . . . 63

7.1.1 Stationare Zustande . . . . . . . . . . . . . . . . . . . . 65

7.1.2 Richtungsfeld . . . . . . . . . . . . . . . . . . . . . . . . 65

7.1.3 Stabilitat . . . . . . . . . . . . . . . . . . . . . . . . . . 66

7.1.4 Modellverbesserungen . . . . . . . . . . . . . . . . . . . 67

7.1.5 Das skalierte allgemeine Rauber-Beute-Modell . . . . . . 70

7.2 Verhalten 2-dimensionaler autonomer Systeme . . . . . . . . . . 71

7.2.1 Stationare Zustande . . . . . . . . . . . . . . . . . . . . 71

7.2.2 Geschlossene Bahnen . . . . . . . . . . . . . . . . . . . . 72

7.2.3 Satz von Poincare . . . . . . . . . . . . . . . . . . . . . . 72

7.2.4 Teilgleichgewichte (Isoklinen) . . . . . . . . . . . . . . . 72

7.2.5 Stabilitat stationarer Zustande . . . . . . . . . . . . . . 73

7.2.6 Stabilitat im Beispielsystem . . . . . . . . . . . . . . . . 75

8 Raumliche Systeme 79

8.1 Zellulare Automaten . . . . . . . . . . . . . . . . . . . . . . . . 79

8.1.1 Conway’s Life . . . . . . . . . . . . . . . . . . . . . . . . 79

8.1.2 Charakteristika von zelluaren Automaten . . . . . . . . . 81

8.1.3 Zirkularer Raum . . . . . . . . . . . . . . . . . . . . . . 81

8.1.4 Erweiterungen . . . . . . . . . . . . . . . . . . . . . . . . 82

8.2 Transportprozesse . . . . . . . . . . . . . . . . . . . . . . . . . . 82

8.2.1 Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

8.2.2 Advektion (Transportgleichung) . . . . . . . . . . . . . . 85

8.2.3 2D Advektion . . . . . . . . . . . . . . . . . . . . . . . . 86

9 Ausbreitung von Borkenkafern (W. Ebenhoh) 87

9.1 Borkenkafer-Dynamik . . . . . . . . . . . . . . . . . . . . . . . . 87

9.2 Wald-Dynamik . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

9.3 Raumliche Ausdehnung des Modells . . . . . . . . . . . . . . . . 93

9.3.1 Bekampfung . . . . . . . . . . . . . . . . . . . . . . . . . 93

10 Mathematische Epidemiemodelle 95

10.1 Das klassische Epidemiemodell nach Kermack & McKendrick . . 96

10.1.1 Stabilitatsanalyse des SIR-Modells . . . . . . . . . . . . 97

10.2 Masern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

10.2.1 Simulationsergebnis . . . . . . . . . . . . . . . . . . . . . 101

A Funktionen 103

A.1 Polynome . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

A.2 Periodische Funktionen – Winkelfunktionen . . . . . . . . . . . 104

B Differential- und Integralrechnung 106

B.1 Kettenregel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

B.2 Substitutionsregel . . . . . . . . . . . . . . . . . . . . . . . . . . 107

B.3 Partielle Integration . . . . . . . . . . . . . . . . . . . . . . . . 108

C Differentialgleichungen 109

C.1 Numerisches Losen von Differentialgleichungen . . . . . . . . . . 113

C.1.1 Verbesserungen . . . . . . . . . . . . . . . . . . . . . . . 115

D Zweidimensionale lineare Abbildungen 115

E Stabilitatsanalyse autonomer Systeme (2D) 119

E.1 System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

E.2 Stationare Zustande . . . . . . . . . . . . . . . . . . . . . . . . 119

E.3 Geschlossene Bahnen . . . . . . . . . . . . . . . . . . . . . . . . 119

E.4 Satz von Poincare . . . . . . . . . . . . . . . . . . . . . . . . . . 119

E.5 Teilgleichgewichte . . . . . . . . . . . . . . . . . . . . . . . . . . 120

E.5.1 Stabilitat stationarer Zustande . . . . . . . . . . . . . . 120

E.6 Veranschaulichung des Systemverhaltens . . . . . . . . . . . . . 122

GMM WS2009/2010 Kapitel 1 – Einfuhrung

1 Einfuhrung

1.1 Was ist Modellbildung – mathematische Modellierung

Die mathematische Modellbildung oder mathematische Modellierung

• bezeichnet eine Methode

• ist nicht an eine spezielle Wissenschaft gebunden und wird in Naturwis-

senschaften und Technik und in der Okonomie angewendet

• versucht Teile der Realitat mathematisch begreifbar zu machen

(Natur)-Wissenschaft ist Modellierung: In der Wissenschaft werden Modelle

aufgestellt, um eine vorgebenen Fragestellung zu beantworten.

Jedes Modelliervorhaben braucht eine Leitfrage oder ein Ziel!

Dies ist wichtig, da die Art und die Komplexitat eines Modells von dieser

Zielvorgabe abhangt. Ein Modell soll einen Teilaspekt der Realitat so darstel-

len, dass es die Information liefern kann, die zur Beantwortung einer Leitfrage

notwendig ist.

Beispiel:

Ein Klimamodell taugt nicht zur Wettervorhersage, ein Regentropfenmodell

auch nicht.

1.2 Was ist ein Modell ?

Hierfur gibt es keine eindeutige Definition. Man konnte es wie folgt beschrei-

ben:

Ein Modell ist ein Objekt oder Konzept, das benutzt wird, um einen realen

Aspekt so darzustellen, dass er in Hinblick auf eine Zielfrage verstanden werden

kann. Beispiele fur Modelle sind:

• ein Landschaftsbild

• Landkarten

• Wettermodelle

1

1.2. Was ist ein Modell ? GMM WS2009/2010

1.2.1 Merkmale von Modellen

• Vereinfachung. Es werden nur die wesentlichen Aspekte durch das Modell

beschrieben. Will man z.B. das Volumen eines Wurfels berechnen, so ist

dessen Farbe unwichtig.

• Skalierung in Raum und Zeit. Ein Atommodell stellt das Atom in einer

Große dar, die wir sehen konnen, ein Globus ist so klein, dass er bequem

auf dem Schreibtisch stehen kann. Modelle zur Wettervorhersage mussen

schneller sein als das Wetter selbst.

• Modelle haben einen begrenzten Gultigkeitsbereich, z.B. gelten die New-

ton’schen Gesetze nicht nahe der Lichtgeschwindigkeit

1.2.2 Arten von Modellen

• einfache Modelle, die Prinzipien erklaren, z.B. Beschreibung des radio-

aktiven Zerfalls durch die Exponentialfunktion

• komplexe Modelle, wie z.B. Klimamodelle, die die Atmosphare und die

Ozeanstromungen berucksichtigen

• empirische Modelle, die einen funktionalen Zusammenhang zu Messdaten

liefern

• prozessorientierte Modelle

• deterministische Modelle

• stochastische Modelle

Diese Einteilung ist nicht eindeutig! Z.B. Konnen einfache Modelle sowohl

deterministisch als auch stochastisch gebildet werden.

1.2.3 Modell-Simulation

Ergebnisse von Modellen erhalt man durch eine Simulation. Die Simulation

liefert eine Realisierung des Modells. Das Modell gibt die Einsicht in die Zu-

sammenhange, die Simulation liefert ein Ergebnis, das z.B. mit Messdaten

verglichen werden kann. Wird dieser Vergleich fur ausreichend gut befunden,

so kann das Modell den beschriebenen Sachverhalt reproduzieren. Zeigt der

2

GMM WS2009/2010 Kapitel 2 – Empirische Modelle

Vergleich von Simulationsergebnis mit den Daten Abweichungen, so muss dass

Modell verbessert werden.

2 Empirische Modelle

Bei der empirischen Modellierung wird ausgehend von einem oder mehreren

Datensatzen versucht, einen funktionalen Zusammenhang zu bestimmen, der

diese Daten ausreichend gut reproduziert.

2.1 Parameteranpassung

Bei der Parameteranpassung wird eine Funktion vorgegeben und die Parameter

der Funktion anhand der Messdaten bestimmt. Es seien z.B. folgende Werte

fur die Wassertemperatur an der Nordseekuste gegeben (Abbildung 2.1):.

J F M A M J J A S O N D1 Jahr

0

5

10

15

20

25

30

Tem

pera

tur

in °

C

Abbildung 2.1: Jahresgang der Wassertemperatur an der Nordseekuste (Bei-spiel)

Diese Werte sollen nun durch eine Kosinusfunktion approximiert werden. Aus

der Grafik erkennt man, dass der kalteste Zeitraum Mitte Februar bis Mitte

Marz ist. Wir nehmen als kaltesten Tag den Tag 60 an. Die kalteste gemessene

Temperatur betragt ca. 4 ◦C die warmste ca. 22 ◦C. Die Temperaturdifferenz

betragt somit 18 ◦C, die mittlere Temperatur ca. 13 ◦C.

Die Kosinusfunktion soll nun so verandert werden, dass sie zu den Daten passt:

3

2.1. Parameteranpassung GMM WS2009/2010

Die Kosinusfunktion: T = cos(t)

π/2 (90°) π (180°) 3/2 π (270°) 2 π (360°)

-1

0

1

Skalierung auf [0, 1] T = cos(2π t)

0.25 0.5 0.75 1

-1

0

1

Skalierung auf 1 Jahr T = cos(2π t365 )

91 182 274 365

Tage im Jahr

-1

0

1

Im Winter soll es kalt sein T = − cos(2π t365 )

91 182 274 365

Tage im Jahr

-1

0

1

Der kalteste Tag ist Tag 60 T = − cos(2π(t − 60)

365 )

91 182 274 365

Tage im Jahr

-1

0

1

Die Amplitude betragt 9 ◦C T = −9 cos(2π(t − 60)

365 )

91 182 274 365

Tage im Jahr-10

-5

0

5

10

Der Mittelwert betragt 13 ◦C T = 13 − 9 cos(2π(t − 60)

365 ) 91 182 274 365

Tage im Jahr

0

5

10

15

20

25

4


Man erhalt so eine Kurve, die recht gut zu den Daten passt (Abbildung 2.2).

J F M A M J J A S O N D1 Jahr

0

5

10

15

20

25

30

Tem

pera

tur

in °

C

Abbildung 2.2: Jahresgang der Wassertemperatur an der Nordseekuste unddaran angepasste Kosinusfunktion.

2.2 Lineare Regression

Im vorherigen Beispiel haben wir die Parameter der Kosinuskurve anhand von

wenigen Charaketristika, wie die Lage der Extrema bestimmt. Im folgenden

wird ein Verfahren beschrieben, dass die Anpassung einer Geraden an gegebene

Messwerte unter Berucksichtigung aller Messwerte erreicht.

Fur das Wachstum einer Bohne seien folgende Messwerte fur die Zeit in Tagen

seit der Pflanzung und die zu dieser Zeit gemessene Hohe der Bohne gegeben:

Zeit in Tagen 3 5 10 15 20 30 40 50 60 70 80 100Hohe in cm 0,5 1 2 7 15 30 70 130 170 230 248 252

In Abbildung 2.4 erkennt man, dass im Bereich zwischen dem 20. Tag und dem

70. Tag wahrscheinlich ein linearer Zusammenhang zwischen der Wachstums-

zeit und der Wachstumshohe besteht. Diesen Zusammenhang kann man durch

eine Geradengleichung beschreiben:

y = m · x + b

Wir wollen nun m und b bestimmen. Dabei interessiert uns insbesondere m,

die Steigung der Geraden. Sie ist ein Maß fur die Wachstumsgeschwindigkeit.

Eine solche Gerade kann man z.B. erhalten, wenn man zwei beliebige Punkte

5

2.2. Lineare Regression GMM WS2009/2010

0 20 40 60 80 100Zeit in Tagen

0

50

100

150

200

250

300

Höh

e in

cm

Abbildung 2.3: Wachstum einer Stangenbohne.

in diesem Bereich miteinander verbindet. Handelt es sich tatsachlich um einen

linearen Zusammenhang und sind die Messungen nicht fehlerbehaftet, so kann

man so die Steigung und damit die Wachstumsgeschwindigkeit erhalten. Liegen

aber nicht alle Punkte auf dieser Geraden, so wird die Lage der Geraden von

der Wahl der Punkte abhangen.

In der Realitat wird aber weder das Wachstum exakt einem linearen Zusam-

menhang folgen (Bohnen wissen gar nicht, was das ist) noch wird man die Hohe

derart exakt messen. In Abschnitt 2.2.1 wird daher ein Verfahren vorgestellt,

eine Gerade zu bestimmen, die allen Punkten in diesem Bereich Rechnung

tragt.

Aber nehmen wir mal beispielsweise die Gerade, die durch die Punkte bei t=20

und t=70 lauft, dann gilt:

15 = m · 20 + b

230 = m · 70 + b

Auflosen ergibt m = 4.3 und b = −71.0 .

Man kann die Steigung m auch erhalten, indem man ein Steigungsdreieck

betrachtet:

m =∆y

∆x=

230 cm − 15 cm

70 Tage− 20 Tage

Im betrachteten Zeitintervall wachst die Bohne also 4.3 cm/Tag.

Die Methode, nur zwei Punkte zu berucksichtigen, ist naturlich falsch, da sich

Messfehler und Abweichungen in diesen Punkten stark auf die Lage der Gera-

den auswirken. Besser ist es alle Punkte zu berucksichtigen, die im Bereich des

6


0 20 40 60 80 100Zeit in Tagen

0

50

100

150

200

250

300

Höh

e in

cm

Abbildung 2.4: Wachstum einer Stangenbohne. Dargestellt sind die Messwerteund die Gerade, die durch die Punkte (20,15) und (70,230) geht.

linearen Wachstums liegen. Hierzu gibt es die Methode der kleinsten Quadrate

auch lineare Regression genannt.

2.2.1 Methode der kleinsten Quadrate

Seien t1 < t2 < · · · < tn die n Zeitpunkte, an denen die n Messwerte y1, y2, · · · , yn

gemessen wurden. Der Datenbereich des Beispiels aus 2.2, in dem lineares

Wachstum angenommen wird, ist durch

Zeit ti 20 30 40 50 60 70

Hohe yi 15 30 70 130 170 230

gegeben.

Es soll nun eine Gerade y(t) = m · t + b so bestimmt werden, dass die Summe

der Quadrate der Abweichungen von Messdaten und Geraden minimal wird.

Q(m, b) :=n∑

i=1

(yi − y(ti))2 =

n∑

i=1

(yi − (m · ti + b))2

soll also minimal werden. Im Beispiel muss also

Q(m, b) := (15 − (m · 20 + b))2 + (30 − (m · 30 + b))2 + (70 − (m · 40 + b))2

+ (130 − (m · 50 + b))2 + (170 − (m · 60 + b))2 + (230 − (m · 70 + b))2

minimal werden.

Um das Minimum zu bestimmen, werden die partiellen Ableitungen von Q

nach m und b bestimmt und gleich null gesetzt:

7

2.2. Lineare Regression GMM WS2009/2010

∂Q(m, b)

∂m= 2

n∑

i=1

(yi − (m · ti + b)) · (−ti) =: 0

∂Q(m, b)

∂b= 2

n∑

i=1

(yi − (m · ti + b)) · (−1) =: 0

Man erhalt also zwei Gleichungen mit den zwei Unbekannten m und b. Formal

muss man naturlich zeigen, dass es sich tatsachlich um ein Minimum handelt!

Satz 2.2.1 Methode der keinsten Quadrate (lineare Regression)

Es seien n Messwerte yi, i = 1 · · ·n zu den Zeitpunkten ti, i = 1 · · ·n gegeben.

Die Summe der Abstandsquadrate der Messwerte von der Geraden

y(t) = m · t + b wird durch

m =

n ·n∑

i=1

ti · yi − T · Y

n ·n∑

i=1

t2i − T · Tb =

1

n(Y − mT )

minimiert, wobei

T :=

n∑

i=1

ti = t1 + · · ·+ tn Y :=

n∑

i=1

yi = y1 + · · ·+ yn

gilt.

Mit Hilfe der Großen Mittelwert, Varianz und Kovarianz erhalt man

m =Cov(t, y)

V (t)b = y − mt

mit

t =1

n

n∑

i=1

ti y =1

n

n∑

i=1

yi V (t) =1

n

n∑

i=1

(ti − t)2 Cov(t, y) =1

n

n∑

i=1

(ti − t) · (yi − y)

Im Beispiel erhalt man die Gerade

y = 4.44 · t − 92.43

8


Das so bestimmte Wachstum betragt also 4.44 cm/Tag. In Abbildung 2.5 sieht

man, dass die Gerade zwar keinen der Messpunkte genau trifft, aber trotzdem

besser ist als die Gerade in Abbildung 2.4. Dies liegt daran, dass Ausreisser

weniger stark einfliessen.

0 20 40 60 80 100Zeit in Tagen

0

50

100

150

200

250

300

Höh

e in

cm

Abbildung 2.5: Wachstum einer Stangenbohne. Dargestellt sind die Messwerte,fur die ein lineares Wachstum angenommen wird, und die Ausgleichgerade, diesich nach der Methode der kleinsten Quadrate ergibt.

Das Verfahren der linearen Regression ist in vielen wissenschaftlichen Taschen-

rechnern und Tabellenkalulationsprogrammen implementiert. Diese geben mei-

stens zusatzlich zu den Parametern das sogenannte Bestimmtheitsmaß r2 an.

Hierbei ist r durch

r =Cov(t, y)

√

V (t) · V (y)

gegeben. Das Bestimmtheitsmaß r2 gibt den Anteil der durch die unabhangige

Große t erklarten Varianz an der Gesamtvarianz von y an. Liegen die Werte

exakt auf einer Geraden, so ist r2 = 1.

2.3 Logarithmischer Zusammenhang am Beispiel der Karzinogenese

Haufig folgen naturliche Prozesse keinem linearen Zusammenhang. Am folgen-

den Beispiel soll gezeigt werden, wie die Methode der linearen Regression auf

Prozesse angewendet werden kann, die einem logarithmischen Zusammenhang

folgen.

9

2.3. Logarithmischer Zusammenhang am Beispiel der Karzinogenese GMM WS2009/2010

Beispiel: Karzinogenese bei Ratten.

Durch Gabe von Diethylnitrosamin, kann bei Ratten Krebs ausgelost werden.

Je hoher die tagliche Dosis ist, desto schneller entwickeln die Ratten einen Tu-

mor. Die Zeit von der Diethylnitrosamingabe bis zur Ausbildung des Tumors

heisst Latenzzeit. Es besteht folgender Zusammenhang

Dosis in mg/kg/Tag 0.075 0.15 0.3 0.6 1.25 2.5 5.0 10.0

Latenzzeit in Tagen 1020 645 480 360 265 225 150 120

Auf den ersten Blick kann man nicht entscheiden, welcher Zusammenhang

besteht (Abbildung 2.6).

0 1 2 3 4 5 6 7 8 9 10 11Dosis in mg/kg/Tag

0

200

400

600

800

1000

Late

nzze

it in

Tag

en

Abbildung 2.6: Karzinogenese bei Ratten. Aufgetragen sind die Dosis an Die-thylnitrosamin und die zugehorige Zeit bis zur Ausbildung eines Tumors (La-tenzzeit)

Tragt man die Datenpunkte doppelt logarithmisch auf so erkennt man den

Zusammenhang (Abbildung 2.7). Doppelt logarithmisch bedeutet, dass sowohl

die x- als auch die y-Achse logarithmisch sind. Hierbei ist es unerheblich, ob

man die Achsenskalierungen logarithmisch wahlt oder die Werte in der Tabelle

logarithmiert und die Achsen linear skaliert (siehe Achsen in Abbildung 2.7).

Es scheint ein linearer Zusammenhang zwischen dem Logarithmus der Dosis

und dem Logarithmus der Latenzzeit zu bestehen:

log(Latenzzeit) = m · log(Dosis) + b

Die Parameter m und b kann man nun mit der Methode der kleinsten Quadrate

10


0,01 0,1 1 10Dosis in mg/kg/Tag

1

10

100

1000

Late

nzze

it in

Tag

en

-2 -1,5 -1 -0,5 0 0,5 1log(Dosis)

0

0,5

1

1,5

2

2,5

3

log(

Late

nzze

it)

Abbildung 2.7: Karzinogenese bei Ratten. Aufgetragen sind die Dosis an Die-thylnitrosamin und die zugehorige Zeit bis zur Ausbildung eines Tumors (La-tenzzeit) in doppelt logarithmischer Darstellung

bestimmen. Hierzu macht man sich zuerst die Wertetabelle der logarithmierten

Werte. Man beachte, dass diese Werte keine Einheit haben!!:

log(Dosis) -1.12 -0.82 -0.52 -0.22 0.1 0.4 0.7 1

log(Latenzzeit) 3.009 2.810 2.681 2.556 2.423 2.352 2.176 2.079

Man erhalt m = −0.42 und b = 2.48.

Aus log(Latenzzeit) = m · log(Dosis) + b folgt

Latenzzeit = 10m·log(Dosis)+b = 10b · Dosism = eb·ln(10)+m·ln(Dosis)

Dieser funktionale Zusammenhang ist in Abbildung 2.8 dargestellt. Man hat

nun eine Funktion, mit der man zu einer gegebenen Dosis die Latenzzeit be-

stimmen kann. In der Praxis wird die tatsachliche Latenzzeit davon abweichen

(Ratten sind keine Computer). Auch weiss man nicht, welchen Gultigkeitsbe-

reich das Gesetz hat (bei sehr hohen Dosen werden die Ratten vermutlich viel

schneller Tumore entwickeln, bei sehr niedrigen eventuell gar keine).

2.4 Verallgemeinerung der linearen Regression ersten Grades

Weiss man, dass der Zusammenhang zwischen abhangiger und unabhangiger

Variablen einem anderen Gesetz folgt, so kann man auch hierauf optimieren.

Viel Softwarepakete und Taschenrechner haben hierzu eine Reihe von Funk-

tionen implementiert:

11

2.4. Verallgemeinerung der linearen Regression ersten Grades GMM WS2009/2010

0 1 2 3 4 5 6 7 8 9 10 11Dosis in mg/kg/Tag

0

200

400

600

800

1000

Late

nzze

it in

Tag

en

Abbildung 2.8: Karzinogenese bei Ratten. Aufgetragen sind die Dosis an Die-thylnitrosamin und die zugehorige Zeit bis zur Ausbildung eines Tumors (La-tenzzeit) sowie der funktionale Zusammenhang, der sich nach der Methode derkleinsten Quadrate ergibt.

f(x) = a1x + a0 linear

f(x) = a2x2 + a1x + a0 quadratisch

f(x) = a3x3 + a2x

2 + a1x + a0 kubisch

... ...

Bei der Optimierung mit Hilfe der angegebenen Funktionen handelt es sich

auch um lineare Regressionen, da die Parameter ai nicht in hoherer Potenz

vorkommen. Allgemein handelt es sich um eine lineare Uberlagerung von ein-

fachen Funktionen:

f(x) =

J∑

j=1

aj · hj(x) ,

wobei hier hj(x) = xj gilt.

Die Losung eines Optimierungsproblems lauft analog wie im Falle der Geraden.

Man leitet die Funktion der Summe der Abstandquadrate nach den Parametern

ab, setzt diese gleich null und erhalt J Gleichungen mit J Unbekannten.

2.4.1 Fourieranalyse

Nun ist es nicht zwingend, dass die Funktionen hj Potenzen von x sind. Hat

man Messwerte, von denen man weiss, dass sie einer Periodizitat unterliegen,

wie z.B. der Temperaturverlauf eines Jahres, so kann man auch periodische

Funktionen als Basis wahlen. Das bekannteste Beispiel hierfur ist die Fourier-

12

GMM WS2009/2010 Kapitel 3 – Diskrete Modelle

Analyse. Man wahlt als Funktion eine Uberlagerung aus Sinus- und Kosinus-

funktionen

f(t, a0, a1, . . . , aJ , b1, . . . , bJ ) = a0 +J∑

k=1

(

ak · cos(2πkt

T) + bk · sin(2πk

t

T)

)

wobei T die Periodenlange ist. Sind die Stutzwerte ti der Messwerte yi, i =

1, . . . , N aquidistant uber die Periode T verteilt, so ergibt sich fur die optimalen

Parameter

a∗0 =

1

N

N∑

i=1

yi a∗k =

2

N

N∑

i=1

yi · cos(2πktiT

) b∗k =2

N

N∑

i=1

yi · sin(2πktiT

)

3 Diskrete Modelle

Im vorherigen Kapitel haben wir ausgehend von Messwerten einen funktiona-

len Zusammenhang hergestellt. Manchmal hat man aber das Problem, dass es

keine Messungen gibt. In diesem Fall fragt man die Experten z.B. Biologen

oder Physiker wie sich das System, das man durch ein Modell beschreiben

mochte, verhalt1.Modelle, die auf den dem System zu Grunde liegenden Pro-

zesse basieren, nennt man auch prozessorientierte Modelle.

3.1 Bakterienwachstum

Aus der Mikrobiologie sei bekannt, dass sich Bakterien einer bestimmten Art

ca. alle 20 Minuten teilen. Wir nehmen vereinfachend an, dass sich dabei auch

die Bakterienbiomasse alle 20 Minuten verdoppelt. Ein Bakterium hat in etwa

einen Durchmesser von 0.5 µm.

Zielfrage: Wieviele Bakterien entstehen aus einer vorgegebenen Anzahl in ei-

ner vorgegebenen Zeit?

1Das Wissen der Experten basiert dabei im allgemeinen naturlich auch auf Messungen undderen Interpretation. Die Messdaten sind aber vielleicht im Laufe der Zeit verloren gegangen.

13

3.1. Bakterienwachstum GMM WS2009/2010

Anhand der uns zur Verfugung stehenden Information stellen wir folgendes

Wachstumsmodell auf:

Bakterienzahl zur Zeit t0 = 0 min x0

t1 = 20 min x1 = 2 · x0

t2 = 40 min x2 = 2 · x1 = 2 · 2 · x0

......

tn = n · 20 min xn = 2 · xn−1 = 2 · · · ·2 · x0 = 2n · x0

Man erhalt das Wachstumsmodell des exponentiellen Wachstums in expliziter

Darstellung 2

xn = 2n · x0 ,

oder in impliziter (rekursiver) Darstellung

xn = 2 · xn−1 ,

wobei jeweils x0 als Anfangswert vorgegeben wird.

Uberprufung der Plausibilitat des Modells

Wieviele Bakterien gibt es nach einem Tag und welches Volumen nehmen sie

ein ?

Ein Tag hat 24 · 3 · 20 Minuten, enspricht also 72 Zeitschritten. Unter der An-

nahme, dass zu Anfang ein Bakterium existiert, hat man nach 72 Zeitschritten

x72 = 272 = 1072 log 2 ≈ 1021.674 ≈ 4.7 · 1021

Bakterien.

Nun, solch eine Zahl sagt uns anschaulich nicht mehr viel. Schauen wir uns

das Volumen an. Wir nehmen der Einfachheit halber an, dass die Bakterien

2Beweis der expliziten Formel durch vollstandige InduktionInduktionsanfang: n = 0 : x0 = 20 · x0

Induktionsschritt n − 1 → n :Es gelte xn−1 = 2n−1 · x0 , dann ist xn = 2 · xn−1 = 2 · 2n−1 · x0 = 2n · x0 .

14


kugelformig sind und trotzdem dicht an dicht (ohne Lucken) gepackt liegen.

Damit unterschatzen wir das tatsachliche Volumen.

Das Volumen VB eines Bakteriums mit einem Durchmesser von

d = 0.5 µm = 0.5 · 10−6 m betragt

VB =4

3π

(d

2

)3

=4

3π(0.25 · 10−6

)3m3 =

1

48π · 10−18m3 .

Das Gesamtvolumen V betragt also nach einem Tag

V = x72 · VB = 47 · 1020 · 1

48· π · 10−18m3 ≈ 3 · 102 m3 .

Dieses Volumen entspricht einem Quader der Kantenlange 3 m× 10 m× 10 m,

also in etwa einem großen Seminarraum.

Vielleicht ist die von uns angenommene Wachstumsrate zu groß.

Bisher ist die Zahl der Bakterien, die pro Zeitschritt hinzukommt, genauso

groß wie die Zahl bereits existierender Bakterien. Die Anderung der Zellzahl

ist also zu jeder Zeit die Zellzahl selbst:

xn+1 = 2 · xn = xn + xn = xn + 1 · xn

Die Wachstumsrate betragt 1.

Vielleicht mussen wir diese Wachstumsrate verringern. Setzen wir anstelle von

1 nun den Kontrollparameter r ein, so erhalten wir folgendes Modell:

xn+1 = xn + r · xn = (1 + r) · xn r ∈ R+

Wir erhalten nun folgende explizite Darstellung:

xn+1 = (1 + r) · xn = (1 + r) · (1 + r) · xn−1 = . . . = (1 + r)n+1 · x0

oder anders geschrieben

xn = x0 · (1 + r)n = x0 · eln(1+r)n

= x0 · en·ln(1+r)

Es handelt sich hierbei um das Modell des exponentiellen Wachstums. Fur jede

positive Wachstumsrate r wachst die Zahl der Bakterien schliesslich uber alle

Grenzen. Also wachst die Bakterienzahl auch dann uber alle Grenzen, wenn

wir die Wachstumsrate verringeren (Abbildung 3.1).

15


0 5 10 15 20 25Zeitschritte

0

5×105

1×106

Bak

terie

nzah

l

Wachstumsrate r=1,0Wachstumsrate r=0,8

Exponentielles Wachstum

Abbildung 3.1: Vergleich des exponentiellem Wachstums fur die Wachstums-rate r = 1 (*) und r = 0, 8 (+).

Dies ist typisch fur das exponentielle Wachstum der Form

xn = x0 · λn = x0 · en·ln(λ)

wobei nun λ der Kontrollparameter ist. Fur λ > 1 wachst das System uber

alle Grenzen, fur 0 < λ < 1 bricht das System zusammen.

Nun, irgendetwas ist also im Bakterienmodell immer noch falsch. Wir haben

bisher eine sehr grobe Zielfrage im Auge gehabt. Anfangs wurde nicht spezi-

fiziert, fur welche Zeitraume das Modell Gultigkeit haben soll. Wenn wir die

Zahl der Bakterien nach wenigen Stunden bestimmen, so scheint das Modell

plausible Ergebnisse zu liefern.

Wir mussen also entweder den Gultigkeitsbereich einschranken oder das Modell

an den gewunschten Gultigkeitsbereich anpassen. Zuerst stellt sich die Frage,

warum der Gultigkeitsbereich eingeschrankt ist.

Bakterien brauchen zum Wachstum Nahrstoffe und i.a. Sauerstoff. Lasst man

eine Bakterienkultur in einer Petrischale auf Nahrlosung wachsen, so wird die

Nahrlosung nach und nach verbraucht. Je dicker der Bakterienrasen wird, de-

sto schlechter werden innen liegende Zellen mit Sauerstoff und Nahrstoffen

versorgt. Es gibt also eine Nahrstofflimitierung und eine Limitierung durch

den zur Verfugung stehenden Platz. Bakterien werden also mit abnehmendem

Nahrstoff- und Raumangebot immer langsamer wachsen und sich dementspre-

chend auch immer seltener teilen.

Wir mussen die Wachstumsrate mit zunehmender Bakterienzahl abnehmen

16


0 5 10 15 20 25Zeitschritte

0

5×105

1×106

Bak

terie

nzah

l

Exponentielles WachstumLogistisches Wachstum

Abbildung 3.2: Vergleich zwischen exponentiellem Wachstum und logistischemWachstum. Startwert x0 = 1.

lassen , also

xn+1 = xn + R(xn) · xn

wobei R nun eine Funktion der Bakterienzahl ist, die monoton fallend sein soll,

also z.B.

R(x) = 1 − x

K

mit z.B. K = 106. Fur sehr kleine Bakterienzahlen betragt R nahezu 1, und

wird immer kleiner, je mehr sich die Zahl der Bakterien einer Million nahert.

Wir erhalten also das verbesserte Modell:

xn+1 = xn +(

1 − xn

K

)

· xn

K heisst Kapazitat der Bakterienpopulation.

xn+1 = xn +(

1 − xn

K

)

· xn

Dieses ist das Modell des logistischen Wachstums. Fur kleine Bakterienzahlen

wachst die Population nahezu ungebremst. Je großer die Population wird, desto

langsamer wachst sie (Abbildung 3.2).

Auch das logistische Wachstum kann man verallgemeinern:

xn+1 = xn + r ·(

1 − xn

K

)

· xn r ∈ R+ (3.1)

17


0 K

Bakterienzahl zur Zeit n

K

Bak

terie

nzah

l zur

Zei

t n+

1

g(x)=xf(x) = x + r*x(1-x/K)

r=1

x0 x1 x2 x3

x1

x2

x3

Abbildung 3.3: ”Rechte Seite” des Modells zum logistischen Wachstum mitgrafischer Losung

Ist die Zellzahl sehr klein gegenuber der Maximalkapzitat, betragt die Wachs-

tumsrate nun in etwa r.

Die rechte Seite kann man nun als Funktion der Zellzahl x schreiben. Sie gibt

an, wie groß die Zellzahl im jeweils nachsten Zeitschritt ist:

f(x) = x + r ·(

1 − x

K

)

· x = − r

Kx2 + (1 + r) · x

Die Zellzahl ist also als Funktion der Zellzahl im vorherigen Zeitschritt zu

verstehen. In Abbildung 3.3 ist der Graph der Funktion f fur K = 106 und

r = 1 gegeben. Es handelt sich um eine nach unten geoffnete, gestauchte und

nach rechts verschobene Parabel.

Es gilt:

• f(0)=0, d.h wenn keine Bakterien da sind, entstehen auch keine,

• f(K)=K, wenn die Bakterienzahl gerade gleich der Maximalkapazitat ist,

verandert sie sich nichts mehr,

• fur 0 < x < K, gilt f(x) > x, die Zellzahl nimmt zu,

• fur x > K , gilt f(x) < x, die Zellzahl nimmt ab.

Vergrossert man die Wachstumsrate r passieren merkwurdige Dinge. Das Sy-

stem erreicht nicht mehr zwingend den Fixpunkt bei K sondern osszilliert

oder macht sogar Chaos (Abbildung 3.4). Dies ist eine Folge des diskreten

18


Zeitschritts. Befindet sich das System zur Zeit n−1 unterhalb der Maximalka-

pazitat und ist die Wachstumsrate gross, schiesst es mit dem nachsten Schritt

n uber die Maximalkapazitat hinaus. Im Schritt n + 1 fallt es dann wieder

unter K.

0 5 10 15 20

Zeit

0

0.5

1

1.5

2

Bak

terie

n in

Mio

.

r=1

0 5 10 15 20

Zeit

0

0.5

1

1.5

2

Bak

terie

n in

Mio

.

r=2.1

0 5 10 15 20

Zeit

0

0.5

1

1.5

2

Bak

terie

n in

Mio

.

r=2.7

Abbildung 3.4: Logistisches Wachstum fur die Wachstumsraten r = 1.0 (links),r = 2.1 (Mitte) und r = 2.7 (rechts). Startwert x0 = 0.1 Mio., K= 1 Mio.

Um dieses Verhalten zu vermeiden, muss man den Zeitschritt verkleinern, in-

finitisimal werden lassen. Dies wird in Kapitel 5 behandelt.

Tragt man die Losung des Systems (ohne die ersten Schritte, in den sich das

System noch einschwingt) uber der Wachstumsrate r auf so erhalt man ein

sogenanntes Bifurkationsdiagramm (Abbildung 3.5)..

Abbildung 3.5: Bifurkationsdiagramm des logistischen Wachstums.

19

3.2. Fibonacci und die Kaninchen GMM WS2009/2010

3.2 Fibonacci und die Kaninchen

Leonardo da Pisa, auch Fibonacci (figlio di Bonacci, Sohn des Bonacci, gebo-

ren ca. 1180) genannt, war italienischer Mathematiker und gilt als der bedeu-

tendste Mathematiker des Mittelalters. Er brachte das arabische Zahlensystem

nach Europa. Bekannt sind heute vor allem die nach ihm benannten Fibonacci-

Zahlen.

Die Fibonacci-Zahlen sind eine Zahlenfolge, die sich daraus ergibt, dass eine

Zahl jeweils die Summe ihrer beiden Vorganger ist. Begonnen wird hierbei mit

1 1:

1 1 2 3 5 8 13 21 35 55 ...

Das Bildungsgesetz kann man wie folgt formalisieren:

xn+1 = xn−1 + xn

wobei x0 = 1 und x1 = 1 die sogenannten Anfangswerte sind.

In seinem ”Buch der Rechenkunst” (ca. 1227) beschreibt er die Zahlenfolge

anhand des Beispiels eines Kaninchenzuchters der herausfinden will, wie viele

Kaninchenpaare innerhalb eines Jahres aus einem einzigen Paar entstehen,

wenn jedes Paar ab dem zweiten Lebensmonat ein weiteres Paar pro Monat zur

Welt bringt. Er kann somit als Begrunder der mathematischen Modellierung

angesehen werden.

Fibonacci veranschaulichte diese Folge durch die mathematischen Modellie-

rung des Wachstums einer Kaninchenpopulation :

1. Zu Beginn gibt es ein Paar geschlechtsreifer (reproduktiver) Kaninchen.

2. Jedes neugeborene Paar wird im zweiten Lebensmonat geschlechtsreif.

3. Jedes geschlechtsreife Paar wirft pro Monat ein weiteres Paar.

4. Es kommen keine Tiere von außen in die Population oder verlassen diese

Das erste Paar erzeugt seinen Nachwuchs im ersten Monat. Jeden Folgemonat

kommt dann zu der Anzahl der Paare, die im letzten Monat gelebt haben, eine

Anzahl von neugeborenen Paaren hinzu, die gleich der Anzahl der Paare ist,

die bereits im vorletzten Monat gelebt haben, da genau diese geschlechtsreif

sind und sich nun vermehren.

Dieses Modell wollen wir nun formalisieren:

20


3.2.1 Schritte beim Aufstellen eines Modells

1. Die Zielsetzung des Modells definieren (haufig noch vage). Eine verbale

Beschreibung des Modells erstellen. Dies ist schwierig, da die Zielset-

zung noch nicht konkretisiert Eine prazise Formulierung setzt bereits ein

starkes Systemverstandnis voraus. Die erste Formulierung legt meistens

schon die Struktur des Modells fest und last sich spater schwer andern.

Man hat hier meistens auch schon die eigene Denkstruktur festgelegt.

2. Man muss das Wesentliche (z.B. Vermehrungszyklus) vom Unwesentli-

chen(z.B. Farbe der Kaninchen) trennen.

3. Festlegung der Zustandsvariabeln und der Prozesse. Man nimmt Verein-

fachungen an:

• Kaninchen sterben nicht

• Jedes reproduktive Paar produziert exakt ein neues Paar pro Monat

• Jedes neue Paar wird exakt nach einem Monat reproduktiv

Solche Vereinfachungen sind unumganglich, da nicht alle Fakten des

betrachteten Systems bekannt sind. Haufig werden Modelle durch die

Berucksichtigung vieler Details nicht unbedingt realistischer sondern nur

komplizierter.

Die Kunst des Modellierens besteht darin, genau die fur die Zielfrage

wesentlichen Faktoren zu berucksichtigen.

4. Losung/Simulation des Systems

3.2.2 Formalisierung des Kaninchenmodells

Die Zielfrage und die verbale Beschreibung sind bereits durch Fibonacci gege-

ben worden.

Festlegung der Zustandsvariabeln

Die Kaninchenpaare (im folgenden kurz Kaninchen genannt) unterscheiden

sich durch die Eigenschaft der reproduktivitatsfahigkeit. Es gibt Neugeborene,

die sich noch nicht vermehren konnen und bereits Reproduktive. Da wir Aus-

sagen uber die Anzahl der Kaninchen von Monat zu Monat bestimmen wollen,

liegt es Nahe die Zeit in Monatsschritten zu betrachten. Somit liegt es nahe,

21

3.2. Fibonacci und die Kaninchen GMM WS2009/2010

als Zustandsvariablen die Neugeborenen zur Zeit (Nt) und die Reproduktiven

zur Zeit t (Rt) zu beschreiben.

Wir erhalten also die Zustandsvariablen

Nt, n ∈ N0 die Anzahl der Neugeborenen(paare)

Rt, n ∈ N0 die Anzahl der Reproduktiven(paare)

Die Gesamtzahl der Kaninchenpaare Ft ergibt sich dann zu

Ft = Nt + Rt

Ft ist keine Zustandsvariable, sondern eine abgeleitete Große.

Prozessbeschreibungen

Gesucht ist eine Funktion f (hier eine Folge), die zu jeder Zeit t die Werte Nt

und Rt angibt:

f : N −→ N × N (3.2)

t 7→ (Nt, Rt)

Damit weiss man, wie sich der Zustand des Systems in der Zeit andert. Man

hat die Dynamik des Systems beschrieben.

Die Prozesse, mit denen wir es zu tun haben sind:

Reproduktion: Reproduktive Kaninchen bringen Neugeborene hervor

Altern: Neugeborene Kaninchen werden reproduktiv

also erhalt man folgende Gleichungen

Nt+1 = Nt + Rt − Nt = Rt

Rt+1 = Rt + Nt = Rt + Nt

Diese Gleichungen sind nun die Vorschrift, wie man aus den Kaninchenzahlen

zur Zeit t die Zahl zur Zeit t+1 bestimmt. Zusammen mit der Anfangsbedin-

gung N0 = 1, R0 = 0 hat man das Systen nun vollstandig beschrieben.

Mathematisch handelt es sich hierbei um ein lineares Differenzenmodell:

22


(Nt+1

Rt+1

)

=

(

0 1

1 1

)

·(

Nt

Rt

)

3.2.3 Losung/Simulation des Systems

Viele einfache Modelle kann man noch analytisch losen. Dies bedeutet, dass

man eine Funktion findet, die den Systemzustand fur jeden Zeitpunkt explizit

angibt. Bei großeren Modellen ist dies meistens nicht mehr moglich. man lost

sie iterativ, d.h. man berechnet den Verlauf von Schritt zu Schritt. Wenn man

dies fur das Kaninchenmodell tut, stellt man fest, dass die abgeleitete Große F,

die die Gesamtzahl der Kaninchen angibt, gerade die Fibonacci-Zahlen sind.

N: 1 0 1 1 2 3 ...

R: 0 1 1 2 3 5 ...

F: 1 1 2 3 5 8 ...

0 1 2 3 4 5 6 7 8 9 10Zeit in Monaten

0

20

40

60

80

100

Kan

inch

enpa

are

0 1 2 3 4 5 6 7 8 9 10Zeit in Monaten

0

1

10

100

Kan

inch

enpa

are

Abbildung 3.6: Das Kaninchenmodell nach Fibonacci. Die ersten zehn Schritte,links linear, rechts logarithmisch.

Wie Abbildung 3.6 zeigt, besteht auch hier vermutlich ein linearer Zusammen-

hang, wenn man den Logarithmus der Ergebnisse betrachtet:

logFt ∼ t

Um eine analytische Losung zu finden, macht also der Ansatz

Ft = λt

23

3.3. Bienen GMM WS2009/2010

Sinn. Wir haben es also wieder mit dem exponentiellen Wachstum zu tun. Wie

man λ bestimmt, wird in der Ubung gezeigt.

3.3 Bienen

Bei den Bakterien und den Kaninchen wachst die Population schliesslich ex-

ponentiell (wenn man mal die Modellverbesserung bei den Bacs weglasst). Im

folgenden soll das Wachstum einer Bienenkolonie beschrieben werden. Hierzu

wird zwischen Koniginnen und Arbeiterinnen unterschieden.

Es wird angenommen, dass die Zahl der Koniginnen konstant ist, d.h man

muss sie nicht explizit modellieren. Die Zahl der Arbeiterinnen hangt von der

bisherigen Zahl der Arbeiterinnen und der Zahl der gerade Geschlupften ab.

Die Zahl der gerade Geschlupften ist proportional zur Zahl der Koniginnen K

mulipliziert mit der Zahl der Nachkommen pro Konigin r, die ja konstant ist.

Fur die Arbeiterinnen erhalt man also das Modell

At+1 = At + r · K = At + α = A0 + (t + 1) · α

man hat es hier mit linearem (arithmetischem) Wachstum zu tun (Abbil-

dung 3.7)

0 1 2 3 4 5 6 7 8 9 10Zeit

A

Abbildung 3.7: Bienenwachstum.

24


3.4 Zusammenfassung der Wachstumstypen

Exponentielles Wachstum

rekursive Darstellung: xn+1 = λ · xn

explizite Darstellung: xn = x0 · λn = x0 · eln(λ)·n = x0 · 10log(λ)·n

n

xn

λ > 1

n

log(xn)

λ > 1

n

xn

λ <1

n

log(xn)

λ < 1

Arithmetisches Wachstum

rekursive Darstellung: xn+1 = xn + α

explizite Darstellung: xn = x0 + α · n

n

xn

α > 0

n

log(xn)

α > 0

n

xn

α < 0

n

log(xn)

α < 0

25

3.5. Diskrete 2D-Modelle GMM WS2009/2010

Logistisches Wachstum

rekursive Darstellung: xn+1 = xn + α · (1 − xn)

n

xn

r=1

n

log(xn)

r > 0

In der logarithmischen Darstellung sieht man, dass das Wachstum zu Beginn

nahezu exponentiell ist.

3.5 Diskrete 2D-Modelle

Tragt man Nt und Rt des Fibonacci-Modells im sogenannten Phasenraum auf

Abbildung 3.8, so sieht man, dass sich das System nach kurzer Einschwingphase

nur noch in eine Richtung bewegt. Die Frage ist, wie man diese Richtung

bstimmt.

0 2 4 6 8N

0

2

4

6

8

10

12

14

R

Abbildung 3.8: Phasendiagramm des Fibonacci-Modells

Iterativ bestimmt man die Losung des Fibonacci-Modells, indem man die

Fibonacci-Matrix

A =

(

0 1

1 1

)

immer wieder auf den Vektor x =

(

N

R

)

anwendet, wobei mit x0 =

(

1

0

)

gestartet wird:

26


x1 = A · x0 =

(

0 1

1 1

)

·(

1

0

)

=

(

0

1

)

x2 = A · x1 =

(

0 1

1 1

)

·(

0

1

)

=

(

1

1

)

Man erhalt so die Trajektorie

[(

1

0

)

,

(

0

1

)

,

(

1

1

)

,

(

1

2

)

, · · ·]

Wir haben es hier also mit einem linearen System der Form

xn+1 = A · xn

zu tun.

Eine Losung eines solchen linearen Systems lasst sich leicht finden, indem man

die Eigenwerte und Eigenvektoren der Fibonacci-Matrix bestimmt.

Man sucht eine Zahl λ und einen Vektor v fur die gilt

A · v = λ · v

Dies bedeutet, dass die Matrix A den (Eigen-) Vektor v gerade um den (Eigen-)

Wert λ streckt. Hierzu lost man die charakteristische Gleichung

det(A − λ · I) = det

((

0 − λ 1

1 1 − λ

))

= −λ · (1 − λ) − 1 = 0

Zur Bestimmung der Eigenwerte und Eigenvektoren siehe Anhang D.

Im Falle der Fibonacci-Matrix erhalt man 2 verschiedene reelle Eigenwerte mit

ihren zugehorigen Eigenvektoren:

λ1 =1 +

√5

2= 1.618 · · · =: Φ v1 =

(

1

Φ

)

27


In die Richtung von v1 findet eine Streckung um den Faktor Φ statt.

λ2 =1 −

√5

2= −0.618 · · · =: φ v2 =

(

1

φ

)

In die Richtung von v2 findet eine alternierende Stauchung um den Faktor φ

statt.

Nach einer gewissen Zeit ”wirkt” nur noch Φ: xn+1 = Φ · xn.

damit wachsen Rn,Nn und Fn in jedem Zeitschritt um den Faktor Φ.

Jetzt wissen wir auch, in welche Richtung das System wachst: Rt

Nt

= Nt+1

Nt

≈ Φ

3.5.1 Fibonacci und der goldene Schnitt

Zwei aufeinanderfolgende Fibonacci-Zahlen stehen also im Verhaltnis Φ zuein-

ander: Fn+1 = Φ ·Fn. Die Zahl Φ = 1+√

52

≈ 1.618 heisst auch goldener Schnitt.

Sie ist eine beruhmte Zahl und beschreibt in Kunst und Architektur ein be-

sonderes Verhaltnis. Ein Rechteck mit den Seitenlangen a und b (a>b) genugt

dem goldenen Schnitt wenn, gilt:

a + b : a = a : b

Dieses Verhaltnis wird als besonders harmonisch empfunden und findet sich

z.B. bei Bilderrahmen.

3.5.2 Uberlagerung von Losungen – Anfangswertproblem

Wir haben in unseren Uberlegungen bisher nicht berucksichtigt, dass unser

Modell naturlich auch den Anfangswerten Nt = 1und Rt = 0, d.h. x0 =

(

1

0

)

genugen soll. Wir wissen aber, dass beide Ansatze xn = λn1 · v1 und auch

xn = λn2 · v2 die Dynamik des Systems beschreiben.

Da es sich bei unserem Modell um ein lineares Modell handelt, ist aber auch

jede Linearkombination von Losungen wiederum eine Losung, also

xn = α1 · λn1 · v1 + α2λ

n2 · v2

fur beliebige α1, α2 ∈ R.

28


Man muss nun also α1 und α2 so bestimmen, dass x0 =

(

1

0

)

gilt:

x0 = α1 · λ01 · v1 + α2λ

02 · v2 = α1 · v1 + α2 · v2

29


30

GMM WS2009/2010 Kapitel 4 – Aufstellen eines Modells

4 Aufstellen eines Modells

Der Wasserstand eines Sees (W. Ebenhoh)

Wir erfinden einen (Modell-)See. Der Zufluss unterliege starken jahreszeitli-

chen Schwankungen, insbesondere tritt ein Fruhjahrshochwasser auf, und im

Hochsommer ist der Zufluss sehr gering oder verschwindet. Der jahrliche Zu-

fluss konnte den See etwa zweimal fullen. Der Abfluss hangt vom Wasserstand

ab. Sinkt dieser unter eine Schwelle, versiegt er ganz. Ein Teil des zugeflos-

senen Wassers verdunstet, besonders im warmen Sommer. Die Verdunstung

steigt mit der Temperatur stark an. Im Jahresmittel moge etwa die Halfte des

zugeflossenen Wassers verdunsten. Die Oberflache des Sees betrage bei mittle-

rem Wasserstand 2 km2, er habe flache Ufer. Seine maximale Tiefe betragt 6 m.

Gesucht sind die jahreszeitlichen Wasserstandsschwankungen des Sees!

Da die Angaben qualitativ und unvollstandig sind, mussen sie durch sinnvolle

Annahmen erganzt werden. Auch die Gleichungen fur Zufluss, Abfluss und

Verdunstung mussen erfunden werden! Die Aufgabe lasst damit viel Freiheit.

Es werden aber fur die Modellierung zwei technische Forderung gestellt:

• Die Wasserstandsanderungen werden von Tag zu Tag berechnet (die Zeit

schreitet in Tagesschritten fort, diskretes Modell).

• Der Zufluss zum See soll durch zwolf Monatsmittel beschrieben werden.

Der Modellzufluss”springt“ dann jeweils zum Monatswechsel (12 Monate

zu je 30 Tagen).

Bei diesem Beispiel geht es nicht nur um die mathematische Technik, vielmehr

mussen alle Arbeitsschritte eines Modellierers durchgefuhrt werden. Es werden

acht Modellierungsschritte vorgestellt und kommentiert und fur das einfache

Modell”See“ ausgefuhrt.

Acht Schritte bei der Aufstellung eines Modells

1. Verbale Beschreibung des Systems, Fragen, Ziele des Modells

2. Zustandsvariable auswahlen, Prozesse benennen, Modellgleichungen (top

level)

31

Kapitel 4 – Aufstellen eines Modells GMM WS2009/2010

3. Randbedingungen, Antriebe

4. Prozessbeschreibungen, Modellgleichungen (lower levels)

5. Parameterliste, Parameterwerte beschaffen

6. Programmierung und Simulation

7. Modellverbesserung: Iteration der vorangegangenen Schritte

8. Anwendungen des Modells, Dokumentation der Ergebnisse und Erkennt-

nisse

1. Eine verbale Beschreibung des Modells ist dann unerlasslich, wenn man

nicht allein am Modell arbeitet, oder wenn der Modellierer einen Auftrag-

geber hat. Man arbeitet sonst aneinander vorbei. Aber auch sonst deckt

der Versuch, das aufzuschreiben, was man im Kopf zu haben glaubt,

erst die Lucken der Vorstellungen auf. Man schafft durch die verbale

Beschreibung eine hilfreiche Vorstufe des mathematischen Modells: ein

konzeptionelles Modell. Insbesondere ist es wichtig, die Ziele und die Fra-

gen, die mit dem Modell beantwortet werden sollen, gut zu formulieren.

Man modelliert ja nicht einfach”ein System“, sondern modelliert es fur

einen Zweck, und der bestimmt die Art des Modells.

2. Die Zustandsvariablen (auch Systemvariablen) sind die Komponenten

eines Zustandsvektors. Sie werden durch die Fragestellung bestimmt,

hangen also vom Zweck des Modells ab. Die Auswahl kann sich spater im

Modellierungsprozess noch andern (Schritt”Modellverbesserung“), aber

man muss mit einer Entscheidung fur einen Satz von Zustandsvariablen

beginnen, um die nachsten Schritte ausfuhren zu konnen. Bei einem dy-

namischen System andert sich der Zustand mit der Zeit. Deshalb mussen

als nachstes die Prozesse aufgelistet werden, die die Zustandsvariablen

andern. Beide (Zustandsvariable und Prozesse) werden nun verbunden,

indem die Modellgleichungen in einer wenig spezifizierten Form aufge-

stellt werden. Das Ergebnis nennen wir”top level“-Gleichungen, da sie

die Ubersicht uber das gesamte Modell enthalten und noch alle Einzel-

heiten offen lassen. Zum Beispiel konnte man schreiben:

Anderung der Bevolkerungszahl = Geburten - Todesfalle + Zuwanderung

32


Die Bevolkerungszahl ist die Zustandsvariable, Geburten, Todesfalle und

Wanderungsbewegungen sind die Prozesse, die sie andern.

3. Bei einem”System“ muss man immer klar trennen, was dazu gehort und

was als extern anzusehen ist. Die Umgebung beeinflusst die Entwicklung

des Modellsystems, aber – und das ist eine Folge der (kunstlichen) Mo-

dellabgrenzung – das Modellsystem kann das”Externe“ nicht andern,

auch wenn es”in Wirklichkeit“ eine Ruckwirkung uber die Systemgren-

zen hinaus gibt. Soll diese modelliert werden, sind die Systemgrenzen

weiter zu fassen. Die externen Bedingungen, die fur den Fortgang der

internen Entwicklung wesentlich sind, werden”Antriebskrafte“ genannt

(”driving forces“). Bei Okosystemen sind die wichtigsten Antriebskrafte

Licht, Temperatur und Niederschlag. Sie andern sich mit der Jahreszeit

und der Tageszeit. Aber auch andere”Randbedingungen“ gehoren zu den

Antriebskraften. Das sind z.B. Konzentrationen von Stoffen außerhalb

der raumlichen Systemgrenzen, wenn sie diese durch Transportprozesse

uberschreiten und dadurch im System die Zustandsvariablen andern.

4. Die”top-level“-Gleichungen mussen durch jede Menge Einzelbetrachtun-

gen prazisiert werden, und zwar so weit, dass man sie auf dem Rech-

ner programmieren kann. Es mussen Formeln fur die Prozesse gefunden

oder erfunden werden (”lower level“). Dazu ist notig, dass man sich klar

macht, wovon die Prozessraten abhangen, und wie die Einflussfaktoren

zu gestalten sind. Dabei werden Parameter in den neu entwickelten Glei-

chungen auftauchen, von denen einige noch unbekannte Werte haben.

5. Den Parametern, die in den spezifizierten Gleichungen von Schritt (4)

enthalten sind, mussen Werte zugeordnet werden. Viele Werte kann man

aus der Literatur, aus Lehrbuchern oder durch Expertenbefragung gewin-

nen. Allerdings muss der Modellierer meistens die Werte einiger Parame-

ter schatzen, da sie nicht bekannt sind. Diesen Prozess, das”Erraten“

von Parameterwerten, nennt man”educated guessing“. Es handelt sich

nicht um ein blindes Herumraten, sondern um den Einsatz des gesam-

ten Systemverstandnisses, das sich der Modellierer auf den bisherigen

Schritten und durch Literaturstudium erworben hat. Regelmaßig ist ein

Teil der Parameter in den Modellgleichungen erst durch Aufstellung der

Gleichungen definiert worden, also durch die Einsicht des Modellierers in

das Problem. Fur diese Parameter konnen gar keine Werte vorher durch

33

Kapitel 4 – Aufstellen eines Modells GMM WS2009/2010

Messungen bestimmt worden sein! Durch die Modellierung liegt nun eine

theoretische Vorstellung des Systems vor, die es erlaubt, Experimente zu

ersinnen, mit denen diese neuen Parameter gemessen werden konnten.

Zunachst aber muss man diese Werte schatzen.

6. Die Programmierung des Modells ist eine ganz andere Kunst, andere

Fertigkeiten sind gefordert, die man im Prinzip auch zu anderen Zei-

ten (vorher) lernen kann. Die Durchfuhrung von Simulationslaufen und

der Vergleich mit Beobachtungsergebnissen oder experimentellen Resul-

taten erfordert dann wieder den Modellierer im Wissenschaftler. Man

kann davon ausgehen, dass die Modellresultate der ersten Versuche nicht

das gewunschte Bild zeigen. Nun kommt es darauf an, herauszufinden,

woran das liegt. Ist es ein Programmierfehler, ein Fehler in den Parame-

terabschatzungen oder ein Denkfehler beim Aufstellen der Modellglei-

chungen? Wenn es sich um Programmierfehler handelt, so erfordert es

Geduld, Ubersicht und detektivische Kleinarbeit, sie in einem langeren

Code aufzuspuren. Fehlersuche ist anspruchsvoll und wird durch Erfah-

rung erleichtert.

7. Spatestens an dieser Stelle erweist sich der Modellierungsprozess als Ite-

rationsschleife. Man muss wieder und wieder die Modellgleichungen in

Frage stellen (nicht nur den Code oder die Parameterwerte). Man muss

sogar die Ziele des Modells hinterfragen. Der Vergleich zwischen Simu-

lation und Messdaten erzeugt Einsicht in das Problem, auch und gerade

dann, wenn die Ubereinstimmung schlecht ist.

8. Wenn das Modell vertrauenswurdig”lauft“, kann man es in verschie-

denen Anwendungen nutzen, um weitere Einsichten in das System zu

erhalten. Durch Sensitivitatsanalysen sind die Parameter und Modelltei-

le zu identifizieren, die das Simulationsergebnis empfindlich beeinflussen.

Man kann Szenarien aller Art berechnen, und insbesondere das Modell

unter Bedingungen laufen lassen, die man dem realen System nicht zu-

muten mag. Die Dokumentation des Modells und die Publikation der

Ergebnisse sind unabdingbar.

34


4.1 Verbale Beschreibung, Ziele

Die Wasserstandsanderungen eines Sees sollen modelliert werden. Es handelt

sich nicht um einen realen See, sondern um einen erdachten See. Wir sind seine

”Schopfer“ und wollen ihn nicht zu langweilig gestalten. Die Einleitungsaufga-

be enthalt eine sehr knappe Beschreibung des Sees und einige Zahlenangaben

(Tiefe, Flache, Zufluss, Verdunstungsanteil). Er soll einen mit der Jahreszeit

stark wechselnden Wasserstand haben. Der Zufluss ist in der kalteren regen-

reichen Zeit hoch, insbesondere gibt es ein Hochwasser (z.B. nach der Schnee-

schmelze). Im Sommer kann eine Periode ohne Zufluss existieren, und in der

Sommerhitze werden die großten Wasserverluste durch Verdunstung hervorge-

rufen. Der Seespiegel sinkt dann so stark, dass der Abfluss versiegt. Im Ex-

tremfall konnte der See auch austrocknen. Der uber das Jahr integrierte Zufluss

soll etwa zweimal das mittlere Seevolumen betragen, davon verdunstet etwa die

Halfte. Das Modell dient zu Lehrzwecken, deshalb wird ein sogenanntes”kli-

matologisches Jahr“ zugrunde gelegt, also mittlere Verhaltnisse im Zeitraum

der letzten 30 Jahre (mindestens). Dieses klimatologische Jahr wiederholt sich

im Modell immer wieder, und die Losung wird entsprechend periodische jahres-

zeitliche Schwankungen aufzeigen. Die Freiheit eines”erdachten“ Sees ist nicht

nur angenehm, sie erzwingt auch, uber Große, Querschnitt und alles andere

erfinderisch nachzudenken.

Das Wasserstandsmodell konnten wir als erste Stufe einer See-Modellierung an-

sehen. In einer Modellerweiterung wurde ein Betrieb am Ufer einen Schadstoff

uber Abwasser einleiten. Dann interessiert die Zeitentwicklung der Schadstoff-

konzentration im See. Offenbar sind die Zeiten mit niedrigem Wasserstand und

ohne Abfluss besonders kritisch. Man konnte das erweiterte Modell benutzen,

um herauszufinden, wie die Schadstoffeinleitung zu steuern und zu begrenzen

ware, wenn ein vorgegebener Grenzwert im See nicht uberschritten werden

darf.

In weiteren Fortsetzungen konnten dann die Auswirkung des Schadstoffes (oder

einer Eutrophierung) auf das Okosystem des Sees modelliert werden.

4.2 Zustandsvariable, Prozesse

Dynamische Systeme haben einen zeitabhangigen”Zustand“. Der Zustands-

raum ist die Menge aller moglichen Zustande, die das Modellsystem im Prin-

zip annehmen kann. Der aktuelle Zustand des Systems zu einem ausgewahlten

35

4.2. Zustandsvariable, Prozesse GMM WS2009/2010

Zeitpunkt t ist ein Punkt im Zustandsraum. Mit Ablauf der Zeit t verschiebt

er sich darin. Das Entwicklungsgesetz gibt die Regeln an, nach denen diese

Verschiebung erfolgt.

Wenn die Zustandsvariablen definiert sind, werden die Modellgleichungen auf-

gestellt, die es ermoglichen, ihre Zeitentwicklung zu berechnen. Im Falle des

oben beschriebenen Sees bietet sich auf den ersten Blick als einzige Zustandsva-

riable der Wasserstand an. Spater kann als zweite Zustandsvariable die Schad-

stoffkonzentration hinzukommen. Als”Wasserstand“ H bezeichnen wir einen

Pegel, der sich nur durch eine Nullpunktverschiebung von der aktuellen Tiefe

des Sees unterscheidet.

Zustandsvariable: Wasserstand H

Der Name H fur die Zustandsvariable ist willkurlich (H wie Hohe des Was-

serstandes). T fur Tiefe geht nicht, weil T fur Temperatur benutzt werden

soll, und P fur Pegel kollidiert mit dem in aquatischen Okosystemmodellen

ublichen P fur Phytoplankton. Derartige Entscheidungen haben oft langanhal-

tende Auswirkungen. Modellierer halten sich moglichst an die Bezeichnungen,

die ihre Vordenker eingefuhrt haben. Sie erhohen damit die Akzeptanz ihrer

Arbeit.

Mit dem Wasserstand verknupft sind Seeoberflache F und Wasservolumen V,

die beide fur die Modellierung gebraucht werden:

Oberflache: F (H) Volumen: V (H) =H∫

−L

F (h) dh

Der Pegelstand, bei dem der See leer ist, wird hier mit -L (leer) bezeichnet,

er tritt als untere Integrationsgrenze auf und ist negativ, wenn der Pegel-

Nullpunkt daruber liegt. Wird die Seeoberflache F (H) als Funktion des Was-

serstandes H vorgegeben, so ist alles bekannt, was das Modell uber die Geo-

metrie des Sees wissen muss. Die horizontale Form spielt keine Rolle, aber der

vertikale Querschnitt, das Tiefenprofil. Als einfaches Beispiel kann man sich

einen See mit einem parabelformigen vertikalen Querschnitt vorstellen (4.1),

dann ist die Flache proportional zur Tiefe, und die Tiefe ist H + L.

F (H) = FNH + L

HN + L(4.1)

36


In diesem Kapitel soll der Zusammenhang 4.1 gelten. Dabei ist HN ein Referenz-

Wasserstand (etwa der mittlere Wasserstand) und FN die zugehorige Flache.

Die Prozesse sind die Vorgange, die die Zustandsvariablen andern. Der Was-

serstand H wird durch Zufluss, Abfluss und Verdunstung geandert:

Prozesse: Zufluss zu

Abfluss ab

Verdunstung vd

Alle drei Großen haben die Dimension Volumen pro Zeiteinheit, also z.B. die

Einheit m3d−1 (Kubikmeter pro Tag). Es ist von allergroßter Bedeutung, sich

beim Modellieren immer uber Dimension und Einheiten aller auftretender

Großen im Klaren zu sein und bei jedem Term, bei jeder Gleichung sofort

die Einheiten auf Konsistenz zu uberprufen! Das ist ein wichtiges Hilfsmittel,

um Fehler beim Aufstellen der Gleichungen fruhzeitig aufzuspuren.

Die Modellgleichungen (top-level) verbinden die Zustandsvariablen mit den

Prozessen in Form von”dynamischen“ Gleichungen, das sind Gleichungen,

die die Anderungen der Zustandsvariablen beschreiben. In der Regel sind das

Differentialgleichungen oder Differenzengleichungen. Die Prozesse andern den

Wasserstand, indem sie das Volumen andern:

Volumenanderung = Zufluss - Abfluss - Verdunstung

Hier soll die Zeit in festen Tagesschritten ∆t fortschreiten. Statt mit der Ande-

rungsgeschwindigkeit V = dV/dt rechnen wir mit der Naherung”Volumenande-

rung pro Tag“

∆V = (zu − ab − vd) · ∆t (4.2)

Die Multiplikation mit ∆t auf der rechten Seite ist aus zwei Grunden bemer-

kenswert: Zum einen werden die Prozesse mit Dimension Volumen pro Zeit

durch die Multiplikation mit einer Zeit in Volumen verwandelt, wie es die linke

Seite erfordert. Zum anderen erkennt man eine”Linearitat“: Die Volumenande-

rung ist proportional zur Lange des Zeitschrittes. Genau hierin besteht der

Approximationscharakter der Differenzengleichung, denn offensichtlich ist die-

37

4.2. Zustandsvariable, Prozesse GMM WS2009/2010

se Linearitat nur fur kleine Zeitschritte ∆t akzeptabel, namlich solange sich

zu, ab und vd in ∆t nur unwesentlich andern.

In Abbildung 4.1 ist die dynamische Gleichung in Form eines Diagramms dar-

gestellt.

Abbildung 4.1: Schematische Darstellung des Modellsystems. Die Pfeile stellenhier materielle Flusse dar (Wassermengen pro Zeit).

Sprachlich unscharf und sogar falsch ware die Ausdrucksweise:”Das Volumen

hangt von Zufluss, Abfluss und Verdunstung ab“. Richtig ist aber:”Das Volu-

men wird durch Zufluss, Abfluss und Verdunstung geandert“. Das Volumen V

ist keinesfalls eine”Funktion“ von zu, ab und vd! Dann musste sich namlich

aus den augenblicklichen Werten von zu, ab und vd auf das augenblickliche

Volumen schließen lassen, was nicht der Fall ist. Man sieht es schon daran,

dass das Volumen unterschiedlich sein kann, wenn z.B. alle drei Flusse ver-

schwinden. Nicht V selbst sondern die Zeitableitung von V ist eine”Funktion“

der Prozesse!

Als Zustandsvariable hatten wir eigentlich den Wasserstand H und nicht das

Volumen gewahlt! Der Integral-Zusammenhang zwischen V und H uber die

Oberflache F hat die differentielle Form V = F · H. Damit lasst sich (??a)

entsprechend umschreiben:

H =1

F (H)(zu − ab − vd) (4.3a)

Die Dimensionsprufung dieser Gleichung befriedigt, denn (Volumen/Zeit)/Flache

ergibt Lange/Zeit, also z.B. m d−1 fur die Wasserstandsanderung H. Mit dem

festen, endlichen Zeitschritt ∆t erhalt man die Differenzengleichung (4.3b) als

gesuchte Modellgleichung. Sie ist”top-level“, alle notigen Prazisierungen der

38


Prozessbeschreibungen fehlen noch.

∆H =1

F (H)(zu − ab − vd) · ∆t (4.3b)

Abbildung 4.2: Im Zeitintervall ∆t, z.B. einen Tag, andert sich das Volumenum ∆V = (zu − ab − vd) · ∆t. Bei gegebener Oberflache F des Sees sind Vo-lumenanderung und Wasserstandsanderung durch ∆V = F · ∆H verknupft.Damit ergibt sich (4.3b).“

4.3 Randbedingungen, Antriebe

Als”Antriebe“ werden alle Variablen bezeichnet, die das System von

”außen“

beeinflussen. Im Falle des Sees sind das die Niederschlage, die den Zufluss be-

stimmen, und die Temperatur, die die Verdunstung reguliert. Beide andern sich

mit der Jahreszeit. Der Abfluss hangt im Gegensatz zum Zufluss vom System-

zustand (Wasserstand) ab, und die Verdunstung ist neben der Temperatur-

Empfindlichkeit proportional zur Oberflache. In Wirklichkeit wirkt ein großer

See durch Verdunstung auf das Regionalklima und somit auf Wolken, Nieder-

schlag, Zufluss, Luftfeuchtigkeit und Verdunstung zuruck. Durch Festlegung

der Systemgrenzen wird diese Ruckwirkung im Modell aber gerade ausgeschlos-

sen (unterdruckt).

Zufluss: Statt des Niederschlages geben wir den Zufluss selbst als Antrieb vor.

Das vereinfacht die Sache erheblich. Wir ersparen uns, uber die Verzogerung

zwischen Niederschlag und Zufluss nachzudenken. Da sich der Niederschlag und

entsprechend verzogert der Zufluss von Tag zu Tag und von Jahr zu Jahr stark

andern, findet man in der Fachliteratur als Angaben fur die Wasserfuhrung

von Flussen oft Monatsmittel uber mehrere Jahre (klimatologische Mittel).

Sie glatten diese Variabilitat. Auch hier soll so vorgegangen werden. Fur einen

realen See erhalten wir die Angaben aus dem Internet oder durch Anfrage bei

einem Amt fur Hydrologie. Wir dagegen mussen uns fur den hypothetischen See

39

4.3. Randbedingungen, Antriebe GMM WS2009/2010

die 12 Monatsmittel MM(m), m = 1 . . . 12 ausdenken. Dabei ist es sinnvoll,

den absoluten Wert von den relativen Monatsmitteln relMM(m) zu trennen:

MM(m) = relMM(m) · vzu , wobei

12∑

m=1

relMM(m) = 12 (4.4)

Das Jahresmittel vzu des Zuflusses ist in der Aufgabenstellung gegeben:

vzu =2VN

360 d(4.5)

wobei VN = V (HN) das Volumen beim Referenz-Wasserstand HN ist und mit

12 Monaten zu je 30 Tagen gerechnet wird. In einem halben Jahr konnte der

Zufluss also den See fullen, gabe es keine Verluste. Die 12 relativen Monats-

mittel addieren sich zu 12 (im Mittel 1). Im Marz soll ein erhohter Zufluss zu

beobachten sein und im Sommer durch Trockenheit eine zuflusslose Periode.

Das wird durch die in Abbildung 4.3 dargestellten Annahmen erfullt.

Abbildung 4.3: Relative Monatsmittel des Zuflusses zum hypothetischen See(Summe ist 12).

Temperatur: Auch die Temperatur setzen wir klimatologisch an. Mehr noch,

sie soll sich wie eine verschobene Sinusfunktion verhalten. Dann benotigt man

nur wenige Angaben, um die Temperatur an jedem Tag t im Jahr auszurechnen,

die jahrliche Mitteltemperatur Tm, die Amplitude der Temperaturschwankung

Ta und den kaltesten Tag tk im Jahr:

T (t) = Tm − Ta · cos

(

2πt − tk365 d

)

(4.6)

Am kaltesten Tag, das ist etwa der 1. Marz (Wassertemperatur!), ist T (tk) = Tm − Ta.

Die Formel 4.6 liefert eine recht grobe Naherung, die aber fur dieses Modell

vollkommen ausreicht (Abbildung 4.4).

40


Abbildung 4.4: Temperaturverlauf im Jahresgang. Kaltester Tag: 1. Marz.

4.4 Prozessbeschreibungen

Zufluss: Der Zufluss ist schon als”Antrieb“ im vorangegangenen Modellie-

rungsschritt behandelt worden (ersatzweise fur die Niederschlage). Mit Glei-

chung (4.4) ist fast alles gesagt. Es fehlt nur noch, wie man zu einem gegebenen

Tag t den Monat m ausrechnet. Das ist gar nicht so trivial und ware bei un-

gleich langen Monaten noch umstandlicher:

m =

[1

30(t mod 360)

]

+ 1

Hier ist (t mod 360) der Rest bei der Division der Zeit t (in Tagen) durch

360, und [..] symbolisiert das Abrunden auf die nachste ganze Zahl. Der Tag

t = 400 gehort somit zum Monat 2.

Abfluss: Wenn der Wasserstand unter einen Grenzwert sinkt, liegt der Abfluss

trocken. Um Parameter zu sparen, setzen wir den Pegel-Nullpunkt auf diese

kritische Grenze. Steigt der Wasserstand uber den Pegelwert 0 an, nimmt der

Abfluss erst langsam, dann schneller zu. Es gilt, fur dieses Verhalten eine For-

mel zu finden (4.7):

ab (H) = vab ·(

max

(

0,H

HN

))2

(4.7)

Diese Formulierung enthalt eine Konstante mit den”richtigen“ Einheiten und

einen dimensionslosen Ausdruck. Hinter diesem Ansatz steckt ein extrem wich-

tiges Prinzip, das man beim mathematischen Modellieren dringend beachten

sollte: Die Formeln werden so strukturiert, dass die einzelnen Teile eines Pro-

duktes Dimensionen mit sinnvoller Bedeutung haben. Hier ist vab eine”Ab-

flusskonstante“ mit der Dimension Volumen/Zeit und den Einheiten m3d−1.

Ihre Bedeutung ist glasklar: Wenn der Wasserstand H auf dem Referenzwert

41

4.4. Prozessbeschreibungen GMM WS2009/2010

HN > 0 steht, gibt sie gerade den Abfluss an, da dann der dimensionslose

Faktor 1 ist. Nach der Aufgabenstellung in der Einleitung ist sogar ungefahr

ihr Wert bekannt. Da namlich etwa die Halfte des Zuflusses verdunsten und

die andere Halfte abfließen soll, wird vab ≈ 1/2vzu gelten. vzu ist in 4.4 definiert

und hat die gleichen Einheiten wie vab.

Verdunstung: Die Verdunstung vd ist proportional zur aktuellen Oberflache

F des Sees. Die auf die Flacheneinheit bezogene Verdunstung ist damit vd/F .

Sie steigt mit der Temperatur stark an und hangt in Wahrheit außerdem

vom Wind und von der Luftfeuchtigkeit ab. Wir ignorieren diese zusatzlichen

Abhangigkeiten und setzen eine exponentielle Zunahme mit der Temperatur

an (4.8).

vd (t, H) = avd · F (H) · exp

(T (t) − 20 ◦C

Tvd

)

(4.8)

Wie schon beim Abfluss haben auch in (4.8) die Faktoren des Produktes klare

Dimensionen. Die eigentliche Temperaturabhangigkeit steckt in einem dimen-

sionslosen Faktor, der bei T = 20 ◦C den Wert 1 annimmt und sich bei einer

Temperatursteigerung um Tvd ver-e-facht. Damit ist die Bedeutung des Pa-

rameters Tvd beschrieben. Die Verdunstungskonstante avd hat die Dimension

Lange/Zeit mit den Einheiten md−1. Sie gibt an, um wie viel der Seespie-

gel bei 20 ◦C taglich allein durch Verdunstung sinken wurde. Das Produkt

vvd = avd · F hat wieder die Dimension Volumen/Zeit und ist mit vzu und vab

vergleichbar. Weil sich F mit H verandert, ist vvd allerdings keine Konstante.

Es kann gar nicht genug betont und wiederholt werden, wie wichtig es beim

Erfinden von Gleichungen ist, den neu auftretenden Parametern”Bedeutun-

gen“ und sinnvolle Dimensionen zu geben, um sie abschatzbar und messbar,

auch erklarbar und verstandlich zu machen! Ein grundsatzlicher”Trick“ dabei

sind dimensionslose Ausdrucke, wie sie in (4.7) und (4.8) auftreten, wo sie die

wesentlichen Abhangigkeiten vom Wasserstand bzw. von der Temperatur ent-

halten.

Entwicklungsgesetz

Nachdem die drei Prozesse durch (lower-level) Gleichungen prazisiert wurden,

lasst sich das Modell als Ganzes durch Zusammenfassung von 4.1 bis (4.8) in

42


Gleichungsform darstellen (4.9). Das System 4.9 ist die mathematische For-

mulierung des Modells, das”Entwicklungsgesetz“, da es den Zustand zur Zeit

t+∆t aus dem Zustand zur Zeit t zu berechnen erlaubt. Aus gegebenen “alten“

Werten werden”neue“ berechnet, und das wird iterativ wiederholt.

Bei einer Iteration muss ein”Iterationsanfang“ vorgegeben werden. Zu Beginn

muss fur einen Anfangszeitpunkt t0 der zugehorige Wasserstand H0 vorgegeben

werden. Dann kann 4.9 immer wiederholt werden: Fur gegebene Werte t und

H(t) werden zuerst die Hilfsgroßen Flache F und Temperatur T und Monat

m berechnet. Diese Hilfsgroßen werden fur die Berechnung der Prozesse zu, ab

und vd benotigt. Zum Schluss wird mit ihnen die Anderung des Wasserstandes

∆H und der neue Wasserstand H (t + ∆t) = H + ∆H bestimmt. Die Zeit

schreitet um ∆t fort (in 4.9 nicht aufgefuhrt). Nach dem Durchlaufen kann

diese Abfolge mit den”neuen“ Werten t und H erneut gestartet und immer

wiederholt werden, bis ein Endzeitpunkt erreicht ist. Der Iterationsanfang bei

t0 mit dem Startwert H0 legt so eine ganze Iterationsfolge fest.

Es seien t und H = H (t) gegeben:

F = FNH + L

HN + L

T = Tm − Ta · cos

(

2πt − tk365 d

)

m =

[1

30(t mod 360)

]

+ 1

zu = vzu · relMM (m) (4.9)

ab = vab ·(

max

(

0,H

HN

))2

vd = avd · F · exp

(T (t) − 20 ◦C

Tvd

)

∆H =1

F(zu − ab − vd) · ∆t

H(t + ∆t) = H(t) + ∆H

4.5 Parameterliste und Parameterwerte

Wie man aus dem System 4.9 entnehmen kann, hat das Modell die Parameter

HN , FN , L, vab, avd, Tvd, vzu, relMM, Tm, Ta, tk

HN , FN beschreiben die Seegroße, L die Tiefe bei Pegel 0, vzu, relMM, Tm,Ta,

tk sind Angaben uber die externen Antriebe. Dabei hangt vzu durch die Vor-

43

4.5. Parameterliste und Parameterwerte GMM WS2009/2010

gaben von der Seegroße ab. Nur die drei Parameter vab,avd, Tvd beeinflussen

unmittelbar die Prozesse, ohne von externen Daten abhangig zu sein. Wir

nennen sie die dynamischen Parameter (genauer: Dynamik-Parameter). Sie

mussen nicht geandert werden, wenn beispielsweise das klimatologische Jahr

durch eine Abfolge unterschiedlicher”realer“ Jahre ersetzt wird. Wenn ein rea-

ler See modelliert werden sollte, mussten die Parameterwerte durch Befragung

von Experten, durch Literaturstudium und zum Teil durch eigene Schatzungen

bestimmt werden. Das ist eine schwierige und verantwortungsvolle Aufgabe.

Hier aber denken wir uns einen See aus. Auch dies erfordert tiefes Nachdenken

uber die Parameterwerte, wenn sie einigermaßen konsistent und realistisch sein

sollen.

Seegroße: Unser See soll bei Normalfullung 2 km2 Oberflache haben und 6 m

Wassertiefe, das Tiefenprofil wird als parabelformig angesetzt. Den Pegel-

Nullpunkt (wo der Abfluss gerade versiegt) wollen wir so festlegen, dass der

Normalpegel HN bei 1m liegt.

FN = 2 km HN + L = 6 m HN = 1 m L = 5 m VN = 6 · 106 m3 (4.10)

Antriebe: Der Zufluss soll so hoch sein, dass das Wasser im Mittel in einem

Jahr etwa zweimal erneuert wird. Dadurch wird der Parameter vzu festgelegt.

vzu =2VN

360 d≈ 30 · 103 m3 d−1 (4.11)

Fur die relativen jahreszeitlichen Schwankungen relMM wurden bereits in Bild

2 Annahmen gemacht. Fur die Parameter der Temperaturschwankung nach

Gleichung (4.6) wurden in Bild 3 auch schon Werte festgelegt, die fur einen

eher warmen Ort mit ausgepragtem Landklima gelten konnten:

Tm = 10 ◦C Ta = 12 ◦C tk = 60 d (4.12)

Abfluss: Vom Parameter vab hangt ab, um wieviel der Wasserstand uber den

kritischen Pegel 0 steigen kann. Die Formel (4.7) fur den Abfluss wurde so

konstruiert, dass bei einem Pegel von HN = 1 m (1 m in der Abflussrinne)

der Abfluss vab betragt. Da sich bei mittlerem Zufluss vzu bei diesem Pegel ein

44


Gleichgewicht einstellen soll, muss man (4.13) ansetzen, denn entsprechend der

verbalen Beschreibung verdunstet die Halfte des zufließenden Wassers.

vab ≈1

2vzu (4.13)

Wenn vab um einen Faktor 2 vergroßert oder verkleinert wird, so hat das wegen

der quadratischen Abhangigkeit in (4.7) die Auswirkung, dass sich bei gleichem

vzu ein Gleichgewichtspegel von 1.4 m bzw. 0.70 m einstellen wurde. Das sind

relativ geringe Verschiebungen, vab ist also kaum als kritischer Parameter zu

bezeichnen.

Verdunstung: Wir beginnen mit einer groben Schatzung unter Benutzung

des gesunden Menschenverstandes. Es ist plausibel, dass bei 20 ◦C (Wasser-

temperatur!) der Wasserspiegel durch Verdunstung taglich um etwa 2 cm sinkt

– ohne Zufluss naturlich.

avd = 0.02 m d−1 (4.14)

Der wahre Wert hangt sicher von Wind und Luftfeuchtigkeit ab und ist mogli-

cherweise um einen Faktor 2 falsch, aber wohl nicht viel mehr. Auch die

Abhangigkeit der Verdunstung von der Temperatur muss erraten werden, wenn

wir nicht tief in Physiklehrbuchern graben wollen. Die funktionelle Form die-

ser Abhangigkeit wurde bereits in (4.8) als Exponentialfunktion angesetzt. Bei

welcher Temperaturerhohung wird sich wohl die Verdunstung verdoppeln? Es

scheint nicht unvernunftig zu sein, dafur 5 ◦C anzunehmen. Das bedeutet

2 ≈ exp

(5 ◦C

Tvd

)

⇒ Tvd ≈ 5 ◦C

ln 2≈ 7.2 ◦C (4.15)

Prufung auf Konsistenz: Da die Schatzungen fur Zufluss und Verdunstung ganz

unabhangig voneinander gemacht wurden, ist keineswegs klar, ob die Vorga-

be aus der verbalen Beschreibung auch nur annahernd erfullt ist, nach der die

Verdunstung etwa die Halfte des Zuflusses ausmacht. Wir schatzen deshalb die

jahrliche Verdunstung ab. Sie konzentriert sich auf die vier Sommermonate, in

denen die Temperatur bis zu 27oC erreicht, also mit bis zu 2 · avd ≈ 4 cm d−1

erfolgt. In vier Monaten (∼ 124 d) macht das etwa 5 m aus. Der See konnte sich

also weitgehend leeren. Der Zufluss vzu wurde in (4.11) festgelegt. Er fullt den

See im Jahr zweimal. Damit ist das Verhaltnis von Zufluss und Verdunstung

45

4.5. Parameterliste und Parameterwerte GMM WS2009/2010

im Sinn der Vorgabe akzeptabel. Diese Abschatzung wurde hier eingeschoben,

weil derartige Uberlegungen einen wichtigen Teil der Modellierungsarbeit aus-

machen. Es geht eben nicht nur um das Aufstellen von Gleichungen, Nachschla-

gen von Parameterwerten, Programmierung und Darstellung der Simulations-

ergebnisse, sondern auch um den kreativen Umgang mit Konsistenzprufungen,

Vermutungen, Hypothesen, fachubergreifenden Uberlegungen u.a..

In der Tabelle 4.1 werden alle Parameter mit Werten, Einheiten und Erklarun-

gen zusammengefasst.

Symbole Einheiten Werte ErklarungGeometrie:

HN m 1 Normalpegel uber derAbflussschwelle

FN m2 2·106 Seeoberflache beiNormalpegel

VN m3 6·106 Volumen beiNormalpegel

L m 5 Tiefe unter derAbflussschwelleDynamischeParameter:

vab m3·d−1 15 · 103 Abfluss beimNormalpegel

avd m · d−1 0.02 Verdunstung bei 20 ◦CTvd

◦C 7.2 Ver-e-fachung derVerdunstungAntriebskrafte:

vzu m3 · d−1 30·103 Jahresmittel desZuflusses

relMM 1 1.0 - 0.7 -4.3 - 1.9 -0.8 - 0.1 -0.0 - 0.0 -0.0 - 0.6 -1.4 - 1.2

Relative Monatsmitteldes Zuflusses (Summeergibt 12)

Tm◦C 15 Mittlere

JahrestemperaturTa

◦C 12 Amplitude derTemperaturschwankung

tk d 60 Kaltester Tag (1.3.)

Tabelle 4.1: Zusammenfassung der Parameter mit Einheiten, Werten undErlauterungen

46


4.6 Programmierung, Simulation

In Abbildung 4.5 sind die Simulationsergebnisse fur verschiedene Startwerte

(t0, H0) (Zeit und Wasserstand) dargestellt.

Abbildung 4.5: Bild 5: Wasserstand eines Sees im Jahresgang. UnterschiedlicheAnfangspegel

”vergisst“ das System nach wenigen Monaten. Es stellt sich bei

periodischen Antrieben (klimatologisches Jahr) nach einiger Zeit eine periodi-sche Losung ein.

Man gewinnt aus Abbildung 4.5 eine wesentliche Erkenntnis: Das System”ver-

gisst“ den Anfangswert nach einiger Zeit, hier bald nach dem ersten Hoch-

wasser. Solche Systeme nennt man”dissipativ“. Nach einiger Zeit stellt sich

in jedem Fall die gleiche periodische Losung ein, da die Antriebe (Tempera-

tur und Zufluss) als klimatologische Mittel periodisch angesetzt wurden. Es

gibt offenbar eine ausgezeichnete periodische Losung, die zum Anfangswert

H0 = H(1080 d) ≈ −1.074 m gehort, und zu der alle anderen Losungen hin-

streben.

47

4.6. Programmierung, Simulation GMM WS2009/2010

48

GMM WS2009/2010 Kapitel 5 – Modellierung mit Differentialgleichungen

5 Modellierung mit Differentialgleichungen

Gehen wir noch einmal zuruck zu dem Bakterienmodell des exponentiellen

Wachstums:

xn+1 = xn + r · xn ,

wobei x0 als Anfangswert vorgegeben wird. Dieses Modell berechnet die Bak-

terienzahl von Zeitschritt zu Zeitschritt (n → n + 1).

Betrachtet man nun die Bakterienzahl zu einer Zeit und wahlt den Zeitschritt

∆t, so bleibt mit

x(t + ∆t) = x(t) + r · x(t) · ∆t ,

fur ∆t = 1 alles beim alten. Nun stellen wir die Gleichung um und lassen ∆t

immer kleiner werden. Der Differenzenquotient geht uber in den Differential-

quotienten

x(t+∆t)−x(t)∆t

= r · x(t)

↓dx(t)

dt= r · x(t)

Man erhalt die Differentialgleichung (DGL) des exponentiellen Wachstums:

dx(t)

dt= r · x(t) kurz x′ = r · x oder x = r · x , (5.1)

sie hat die Losung

x(t) = C · er·t .

C ist eine freie Konstante, die durch den Anfangswert bestimmt wird. Der

Parameter r heisst Wachstumsrate.

Fur alle, die noch keine Vorlesung zu Differentialgleichungen gehort haben, sei

an dieser Stelle der Anhang C empfohlen.

49

5.1. Die Wachstumsrate GMM WS2009/2010

5.1 Die Wachstumsrate

Betrachten wir noch einmal den diskreten Fall

xn+1 = (1 + r) · xn = x0 · (1 + r)n+1 = x0 · e(n+1)·ln(1+r)

Fur n=0 ergibt sich x1 = x0 · eln(1+r). Im stetigen Fall gilt x(1) = x0 · eα.

Um also mit dem diskreten und dem stetigen Modell dasselbe Wachstum zu

erreichen, muss α = ln(1 + r) gesetzt werden.

x(t) = x0 · e(ln(1+r))·t = x0 · eln(1+r)t

= x0 · (1 + r)t ,

Will man ganz speziell nach einem Zeitschritt eine Verdopplung (r=1) errei-

chen, so gilt α = ln 2:

x(t) = x0 · e(ln 2)·t = x0 · eln 2t

= x0 · 2t ,

50

GMM WS2009/2010 Kapitel 6 – Systeme mit einer Variablen

6 Systeme mit einer Variablen

6.1 Aufstellen einer DGL am Beispiel auslaufender Gefaße

6.1.1 Der auslaufende Zylinder

Ein Gefaß mit festem Querschnitt F (Zylinder) habe ein kleines Loch am Bo-

den. Anfangs sei es bis zur Hohe H0 gefullt. Die zeitlich veranderliche Fullhohe

H(t) ist die Systemvariable, d.h. es ist die Funktion H(t) gesucht, die das Sy-

stem beschreibt. Diese erhalt man, indem man die Anderung der Hohe mit der

Zeit beschreibt. Wenn Reibung und Zahigkeit dominieren (das Loch ist klein),

ist die Auslaufgeschwindigkeit (Anderung des Fullvolumens pro Zeiteinheit)

proportional zum Druck am Loch, und dieser ist proportional zur Hohe H.

Die Proportionalitatskonstante hangt von Zahigkeit, Lochform und Lochgroße

ab. Das sind Großen, die sich wahrend des Auslaufprozesses nicht andern. Die

Auslaufgeschwindigkeit hat die Dimension Volumen/Zeit, man muss deshalb

zur Herleitung der Differentialgleichung fur H den Umweg uber die Anderung

des Fullvolumen ∆V im Zeitintervall ∆t gehen:

Es gilt, da die Zylinderflache F konstant ist ∆V = F · ∆H , also

∆V = F · ∆H = −k · H · ∆t also ∆H = − k

F· H · ∆t . (6.1)

wobei k die Proportionalitatskonstante ist. Der Grenzubergang liefert

H = − k

F· H , (6.2)

Gleichung 6.2 hat dieselbe Form wie Gleichung 5.1. Sie hat fur die Anfangsbe-

dingung H(0) = H0 die Losung

H(t) = H0 · exp

(

−k · tF

)

. (6.3)

Die Fullhohe des Zylinders nimmt also mit der Zeit exponentiell ab.

Um den Wert von k zu bestimmen und um die Bedeutung von k besser zu

erkennen, kann man in Gedanken experimentieren:

Uberlegen wir zuerst, welche Einheit k hat: V hat die Einheit Volumen pro

51

6.1. Aufstellen einer DGL am Beispiel auslaufender Gefaße GMM WS2009/2010

Zeit, z.B. Kubikmeter pro Sekunde , H hat die Einheit Lange z.B. Meter, dann

muss nach Gleichung 6.2 der Parameter k die Einheit Flache pro Zeit also z.B.

Quadratmeter pro Sekunde haben.

Angenommen, ein Zylinder mit einem Radius von 5 cm sei anfangs bis

H0 = 30 cm mit Wasser gefullt, und wir beobachten, dass der Wasserstand

in der Zeit ∆t = 1 min um ∆H = −1 cm gesunken ist. Dann kann man k

abschatzen. Das in ∆t ausgelaufene Volumen ∆V ist naherungsweise gleich

dem Produkt aus Anfangs-Auslaufgeschwindigkeit kH0 und Zeitintervall :

∆V = F · ∆H ≈ −kH0 · ∆t

Damit gilt

k ≈ −∆H · FH0 · ∆t

=−(−1 cm) · (π · (5 cm)2)

30 cm · 1 min≈ 78, 54 cm2

30 min≈ 2, 62 cm2/min

6.1.2 Der auslaufende Trichter

Im Falle des Zylinders anderte sich die Flache der Wasseroberflache nicht.

Anders beim Trichter. Hier verandert sich die Flache F in Abhangigkeit von

der Fullhohe H . Fur das Wasservolumen im Trichter bei der Fullhohe H gilt

V (H) =

∫ H

0

F (h) dh

Nimmt man wieder an, dass die Auslaufoffnung sehr klein ist, so ist die Aus-

laufgeschwindigkeit proportional zum Druck, als proportional zur Fullhohe H .

Es gilt also wieder

∆V = −k · H · ∆t

also

V = −k · H

Nach der Kettenregel gilt

V ≡ dV (H(t))

dt=

dV (H)

dH· dH

dt= F (H) · H

52


Damit gilt

H = −kH

F (H)(6.4)

Der Radius des Trichters sei bei der anfanglichen Fullhohe H0 durch den

Wert R0 gegeben. Der Trichter besitzt gerade Seiten, damit ist der Radius

in Abhangigkeit von der Fullhohe durch eine Geradengleichung gegeben:

R(H) = m · H + b

Da nun R(H0) = H0 gilt und R(0) = 0, folgt b = 0 und a = R0

H0(Steigung),

also

R(H) =R0

H0· H

Fur die Flache in Anhangigkeit von der Hohe gilt also:

F (H) = π

(R0

H0· H)2

= F0 ·(

H

H0

)2

Man erhalt die Differentialgleichung fur die Fullhohe:

H = −k · H20

F0

· 1

H(6.5)

Aus Gleichung 6.5 folgt H · H = −A, mit A =kH2

0

F0. Integration auf beiden

Seiten fuhrt zu

∫ t

0

H(s) · H(s) ds = −A · t

Partielle Integration der linken Seite fuhrt zu

H(s) · H(s)|t0 −∫ t

0

H(s) · H(s) ds

︸︷︷︸

−A·t

= −A · t

H(t)2 − H20 = −2A · t

H(t) =√

H20 − 2A · t = H0 ·

√

1 − 2k

F0· t mit F0 = F (H0)

53

6.1. Aufstellen einer DGL am Beispiel auslaufender Gefaße GMM WS2009/2010

Man erhalt als Losung

H(t) = H0 ·√

max

(

0, 1 − 2k

F0· t)

(6.6)

54


6.2 Das logistische Wachstum

In Abschnitt haben wird das das logistische Wachstum kennen gelernt. Im

diskreten Fall hatte das logistische Wachstum folgende Form

Yn+1 = Yn + r ·(

1 − Yn

K

)

· Yn .

Diese Gleichung geht nun in

dY (t)

dt= r ·

(

1 − Y

K

)

· Y

uber.

Die Losungen der logistischen DGL weisen im Unterschied zur Differenzenglei-

chung unabhangig von der Parametrisierung nun kein Chaos mehr auf.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

0.2

0.4

0.6

0.8

1

1.2Richtungsfeld der logistischen DGL

Abbildung 6.1: Richtungsfeld des logistischen Wachstums mit einer moglichenTrajektorie.

Am Richtungsfeld (Abbildung 6.1) kann man drei verschiedene Losungstypen

erkennen:

1. Losungen, die zwischen 0 und K bleiben und monoton wachsen

2. Losungen, die sich monoton fallend K annaheren

3. Losungen, die sich von 0 nach unten entfernend gegen −∞ streben

55

6.3. Stationare Zustande GMM WS2009/2010

Das logistische Wachstum weist fur kleine Anfangswerte Y0 zunachst ein ex-

ponentielles Wachstum, denn in diesem Fall gilt

dY (t)

dt= r ·

(

1 − Y

K

)

· Y = r · Y − r

K· Y 2 ≈ r · Y .

Fur Startwerte in der Nahe von K/2 ist das Wachstum zunachst linear.

Die sigmoiden Losungen (0 < Y0 < K/2) der logistischen Differentialgleichung

sind punktsymmetrisch am Wendepunkt. Daher nahern sich Losungen in der

Nahe von K exponentiell dem Wert K.

Die logistische Differentialgleichung kann man noch analytisch losen, die be-

schrankten Losungen lauten

Y (t) = K · er(t−m)

1 + er(t−m)−∞ < m < ∞

Dies ist eine Funktionenschar, der Scharparameter m kann alle reellen Werte

annehmen und hangt fur das Anfangswertproblem Y (0) = Y0 von Y0 ab.

6.3 Stationare Zustande

Autonome Systeme der Form Y = f(Y ) haben gewohnlich spezielle Losungen

der Form

Y (t) = Y ∗ = konstant ,

d.h. dass die Anderung von Y gleich Null ist, also Y = 0 Die Konstante Y ∗

erfullt also die Gleichung

f(Y ∗) = 0 .

Diese zeitlich unveranderlichen Losungen heißen stationare Zustande, sta-

tionare Losungen, Gleichgewichtszustande oder kurz Fixpunkte.

Im Fall des logistischen Wachstums gibt es zwei stationare Zustande:

Y ∗ = 0 und Y ∗ = K.

6.3.1 Stabilitat stationarer Zustande

Es gibt zwei Arten von stationaren Zustanden:

56


1. stabil: kleine Abweichungen wachsen nicht, genauer:

asymptotisch stabil: kleine Abweichungen verschwinden

marginal stabil: kleine Abweichungen wachsen nicht

2. instabil: es gibt kleine Abweichungen, die wachsen

Fur das logistische Wachstum gilt, dass Y ∗ = 0 instabil ist. Sobald man vom

stationaren Zustand auch nur ein wenig abweicht, fangt die Losung an zu wach-

sen. Y ∗ = K ist stabil. Sobald man ein wenig von K abweicht, lauft das System

zu K zuruck.

Bestimmung der Stabilitat eines stationaren Zustands

Bei autonomen Systemen mit einer Zustandsvariablen kann schnell festgestellt

werden, ob ein stationarer Zustand stabil ist oder nicht. man schaut sich hierzu

die Steigung an (Abbildung 6.2): Betrachtet man den Fipunkt Y ∗ = 0, so sieht

man, dass die Steigung positiv wird, sobald man nur ein wenig von Y ∗ nach

rechts abweicht. Damit wird Y immer weiter erhoht. Man sagt der Fixpunkt

Y ∗ = 0 ist instabil, weil kleine Storungen dazu fuhren, dass die Losung vom

Fixpunkt weg lauft. Anders ist es bei Y ∗ = K. Verringert man Y ein wenig, so

ist die Steigung positiv, Y nimmt also wieder zu. Damit kommt das System von

Y ∗ = K nicht weg, der Fixpunkt ist stabil. Analytisch kann man die Stabilitat

y

y’

=y*

(1-y

)

Y*=0 Y

*=K

Abbildung 6.2: ”Steigung” des logistischen Wachstums. Hierbei tragt man y’uber y auf!

57

6.3. Stationare Zustande GMM WS2009/2010

also wie folgt untersuchen: Man bildet die Ableitung der rechten Seite der

Differentialgleichung nach der Zustandsvariablen im stationaren Zustand:

γ =df(Y ∗)

dY

Es gilt :

γ > 0: Y ∗ ist instabil

γ < 0: Y ∗ ist stabil

γ = 0: weitere Analysen sind erforderlich

Fur das logistische Wachstum

Y = f(Y ) = r ·(

1 − Y

K

)

· Y r > 0

gilt

df(Y )

dY= r − 2r

K· Y

Fur Y ∗ = 0 gilt: γ = r > 0, der Zustand ist instabil

Fur Y ∗ = K gilt: γ = r − 2r < 0, der Zustand ist stabil

Was im Detail passiert

Man betrachte Y als Abweichung vom Fixpunkt Y ∗ mit der Storung h: Y =

Y ∗ + h. Fur Y gilt

Y =dY

dt=

d(Y ∗ + h)

dt=

d(Y ∗)

dt+

d(h)

dt=

d(h)

dt= h

Andererseits gilt die Taylorentwicklung

Y = f(Y ) = f(Y ∗ + h) = f(Y ∗) +df(Y ∗)

dY· h + · · ·

Betrachtet man nun die Taylorentwickung nur bis zum ersten Glied1, so erhalt

1

f(x0 + h) ≈ f(x0) + f ′(x0) · h

58


man eine neue, lineare DGL fur die Storung h = Y ∗ − Y .

h = f(Y ∗) +df(Y ∗)

dY· h =

df(Y ∗)

dY· h = γ · h (6.7)

mit der Losung

h(t) = h0 · eγ·t

Damit wachst die”Storung“ fur γ > 0 exponentiell (instabil) und zerfallt ex-

ponentiell fur γ < 0 (stabil). Dieses Verfahren nennt man auch”Linearisierung

am stationaren Zustand“. Man sagt, man habe die ursprungliche DGL lineari-

siert. Nun ist auch klar, warum man fur γ = 0 keine Aussage uber die Stabilitat

des Fixpunkts der ursprunglichen Gleichung treffen kann. Da die Storung un-

verandert bleibt, hangt die Stabilitat nun womoglich von den hoheren Termen

der Taylorentwicklung ab, die beim linearisierten System vernachlassigt wur-

den.

6.4 Skalierung der logistischen Differentialgleichung

Wie schon im diskreten Fall kann man auch bei der logistischen Differential-

gleichung zwei Parameter loswerden:

Setze y = YK

, dann gilt

y =1

KY = r · Y

K· (1 − Y

K) = r · y · (1 − y)

Damit liegt der stabile stationare Zustand nun bei y∗ = 1. Man kann dies

so interpretieren, dass die dimensionslose Zustandsgroße y den Anteil an der

Maximalkapazitat K angibt. Jetzt spielt es keine Rolle mehr, wie groß K ist.

Um auch noch r loszuwerden, muss man die Zeit skalieren. Fur die”neue“ Zeit

s = r · t gilt nach der Kettenregel

dy

dt=

dy(s(t))

dt=

dy(s)

ds· ds(t)

dt

Mit ds(t)dt

= r folgt

dy(s)

ds=

1

r· dy

dt= y · (1 − y)

59

6.5. Wachstums und Sterbeprozesse GMM WS2009/2010

Betragt die maximale Wachstumsrate z.B. 5 pro Minute, und lauft die Zeit t

in Minuten, so lauft die neue Zeit nun in 5 Minuten Abschnitten, so dass die

neue maximale Wachstumsrate gerade eins wird.

6.5 Wachstums und Sterbeprozesse

Betrachtet man die Veranderung einer Population, so kann man fur die Popu-

lationsgroße ganz allgemein die Differentialgleichung

N = Wachstum − Sterben (6.8)

aufstellen. Im allgemeinen ist das Wachstum von der aktuellen Populations-

große abhangig, also

Wachstum = r · N (6.9)

Es findet kein Wachstum aus dem nichts statt.

Auch der Anteil der stirbt hangt von der aktuellen Populationsgroße ab

Sterben = µ · N (6.10)

Sind r und µ Konstanten, kann man sie zu einer Nettowachstumsrate d zu-

sammenfassen

N =d·N (6.11)

Anders ausgedruckt, bedeutet es, dass

N

N= konstant (6.12)

gilt.

Im Fall des logistischen Wachstums gilt:

N = R(N) · N R(N) = r · (1 − N

K) (6.13)

also

N

N= R(N) (6.14)

60


Dieses ist das Wachstumsmodell nach Verhulst (1838). Es handelt sich um

ein dichtereguliertes Wachstum. Fur R(N) > 0 nimmt die Population zu, fur

R(N) < 0 nimmt sie ab. Unter idealen Bedingungen betragt die Wachstums-

rate r.

61

6.5. Wachstums und Sterbeprozesse GMM WS2009/2010

62

GMM WS2009/2010 Kapitel 7 – Systeme mit 2 Variablen

7 Systeme mit 2 Variablen

7.1 Rauber-Beute-Modelle

Beim Beispiel des Bakterienwachstums sind wir davon ausgegangen, dass die

Abnahme der Wachstumsrate von der Bakterienbiomasse selbst abhangt. Je

mehr Bakterien da sind, umso schlechter kann die Population wachsen. Dies

nennt man auch intraspezifische Konkurrenz. Im folgenden wollen wir uber-

legen was passiert, wenn eine andere Art die Population dezimiert und de-

ren Biomasse also von der ersten Art abhangt. Solche Modelle werden auch

Rauber-Beute-Modelle genannt. Traditionell spricht man von Hasen als Beu-

te und Fuchsen als Rauber. In Wirklichkeit fressen Fuchse wohl keine Hasen

und es ware sinniger, Luchse als Rauber zu nehmen, aber das interessiert den

Mathematiker nicht wirklich In diesem Abschnitt soll also das Zusammenspiel

zwischen einer Beutepopulation, den Hasen und einer Rauberpopulation, den

Fuchsen simuliert werden (Abbildung 7.1). Hierzu machen wir folgende An-

nahmen.

Hasen(Beute)

Füchse(Räuber)

wachsen fressen sterben

Abbildung 7.1: Diagramm des Rauber-Beute-Systems. Die Pfeile geben denBiomassenfluss an.

1. Die Hasen wachsen in Abhangigkeit der Anzahl der aktuell vorhandenen

Hasen.

2. Die Fuchse fressen einen Teil der Hasen, umso mehr, je mehr Hasen da

sind.

3. Ein Teil der Fuchse stirbt

Sei B(t) (Beute) die Anzahl der Hasen und R(t) (Rauber) die Anzahl der

Fuchse zur Zeit t.

Zur Zeit t0 = 0 sei die Zahl der Hasen B0 = 200, die Zahl der Fuchse R0 = 20.

Wir nehmen in Kauf, dass auch nicht-ganzzahlige Werte vorkommen konnen.

Wen das stort, der kann als Einheit auch kg Hasen und kg Fuchse wahlen.

63

7.1. Rauber-Beute-Modelle GMM WS2009/2010

Das Wachstum der Hasen wird wie in Abschnitt 6.2 durch

B = r · B

beschrieben. Hierbei ist r die Wachstumsrate der Hasen. Wir nehmen an, dass

die Wachstumsrate der Hasen r = 0.01 pro Tag betragt. Das bedeutet, dass die

Hasenpopulation jeden Tag um ca. 1 %1 wachst, wenn keine Hasen gefressen

werden. Die Abnahme der Rauber ist durch

R = −s · R

beschrieben.Hierbei ist s die Sterberate der Rauber. Wir nehmen an, dass

die Sterberate der Fuchse s = 0.05 pro Tag betragt . Das bedeutet, dass die

Fuchspopulation jeden Tag um 5 % abnimmt, wenn nichts gefressen wird.

Nun mussen wir noch berucksichtigen, wie viel ein Fuchs frisst. Klar ist, dass

er nichts zu fressen findet, wenn es keine Hasen gibt, und dass er umso mehr

frisst, je mehr Hasen da sind (er muss dann nicht so lange suchen). Die Wachs-

tumsrate der Fuchse w wird also von der Hasenzahl abhangen, z.B.

w(B) = b · B

Setzen wir fur den Parameter b = 0.0004, dann gilt bei 100 Hasen gerade

w = 0.04. Stehen den Fuchsen 100 Hasen als Nahrung zur Verfugung, kann

ihre Population am Tag 4% wachsen (dies ist unrealistisch, aber es ist ein sehr

einfaches Modell). Die Differentialgleichungen lauten nun:

B = r · B − b · B · R (7.1)

R = b · B · R − s · R

Simuliert man nun die Hasen- und Fuchszahl (Abbildung 7.2), so fallt folgendes

auf: Beide Populationen oszillieren. Zuerst steigt die Zahl der Hasen. Dadurch

steht den Fuchsen mehr Nahrung zur Verfugung und ihre Zahl steigt kurz

darauf ebenfalls. Dies reduziert die Zahl der Hasen wieder und den Fuchsen

steht wieder weniger Nahrung zur Verfugung usw..

1genauer: R(1) = R0 · (1 + 0.01 + ...)Verdopplungszeit T: R(T ) = 2R0 = R0 · er·T , also T = ln2

r= ln2 · 100.

64


0 500 1000 1500 2000Zeit in Tagen

0

50

100

150

200

Anz

ahl

Hasen (Beute)Füchse (Räuber)

Abbildung 7.2: Losungskurven des Rauber-Beute Systems.

7.1.1 Stationare Zustande

Ahnlich wie beim eindimensionalen Modell kann man die stationaren Zustande

des Modells bestimmen. Hierzu werden beide Gleichungen gleich Null gesetzt.

Man erhalt ein Gleichungssystem aus zwei Gleichungen mit zwei Unbekannten.

B = r · B − b · B · R = (r − b · R) · B = 0 (7.2)

R = b · B · R − s · R = (b · B − s) · R = 0

Eine Losung (trivial) ist durch (B∗, R∗) = (0, 0) gegeben. Eine weitere Losung

ist durch

r − b · R = 0

b · B − s = 0

bestimmt, sie lautet (B∗, R∗) = ( sb, r

b). Wenn das System genau diesen Zustand

hat, wird es sich nichts mehr andern. Die Frage ist, wie sich kleine Storungen

auswirken, d.h. wie stabil dieser Zustand ist.

7.1.2 Richtungsfeld

Fur ein autonomes 2-dimensionales System kann man ebenfalls ein Richtungs-

feld zeichnen. Man betrachet hierzu

65


dR

dB=

dRdtdBdt

d.h. die Anderung des Raubers nach der Beute. Im Falle des Lotka-Volterra-

Modells gilt

dR

dB=

b · B · R − s · Rr · B − b · B · R

Abbildung 7.3: Richtungsfeld des Lotka-Volterra-Modells mit Losungskurvenzu verschiedenen Anfangswerten.

Die R-B-Ebene nennt man auch Zustandsraum oder Phasenraum des Systems,

eine Losungskurve auch Bahn oder Trajektorie2

7.1.3 Stabilitat

Das Modell ist ein singulares Modell, da der nicht-trivial stationare Zustandes

stabil aber nicht asymptotisch stabil ist. Bei kleinen Veranderungen findet dass

Modell nicht in den ursprunglichen Zustand zuruck.

Anmerkung: Will man das Modell mittels Euler-Verfahren losen, ist Vorsicht

geboten. In diesem Fall wird die Zahl der Hasen und Fuchse immer starker

2Diese Bezeichnung ist etwas ungenau. Streng genommen ist die Trajektorie eine Abbildungder Zeit in den Zustandsraum, eine Losung des Anfangswertproblem. Die Bahn ist das Bild,enthalt also die Zeitinformation nicht mehr.

66


oszillieren. Dies ist ein numerischer Artefakt.

Wie man die Stabilitat nachweisen kann wird spater gezeigt.

7.1.4 Modellverbesserungen

Dieses Modell kann naturlich verbessert werden. Folgende Verbesserungen fal-

len einem sofort ein:

• Die Hasen durfen nicht unbegrenzt wachsen, wenn keine Fuchse da sind.

Es muss eine Maximalkapazitat wie in Abschnitt 6.2 geben.

• Die Fuchse mussen satt werden. Die Zunahme der Wachstumsrate der

Fuchse darf nicht linear von der Zahl der Hasen abhangen.

• Nur ein Teil der Hasen kann als Nahrung verwertet werden (Fell und Kno-

chen werden ausgeschieden).

Maximalkapazitat

Verbessern wir das Modell nun dahingehend, dass wir eine Maximalkapazitat

fur Hasen berucksichtigen. In diesem Fall wird das Wachstum der Hasen durch

Gleichung 3.1 beschrieben, und das vollstandige Modell lautet:

B = r · B ·(

1 − B

K

)

− b · B · R (7.3)

R = b · B · R − s · R

Die Modellergebnisse sind in Abbildung 7.3 fur K = 250 und K = 500 dar-

gestellt (andere Parameter wie bisher). Die Oszillationen sind nun mehr oder

weniger stark gedampft. Die Zahl der Hasen und Fuchse strebt letztlich auf

konstante Werte zu.

Sattigung

Will man nun zusatzlich berucksichtigen, dass die Fuchse satt werden, so muss

man anstelle der monoton steigenden Fressrate b ·B eine Fressrate annehmen,

die durch einen maximalen Wert begrenzt ist. Diese maximale Fressrate gibt an

wie viel ein Fuchs pro Zeiteinheit maximal fressen kann. Dies ist u.a. von der so

genannten”handling time“ abhangig, der Zeitspanne, die ein Fuchs benotigt,

die Nahrung aufzunehmen (den Hasen zu zerlegen). Ein weiterer Faktor ist

die Stoffwechselgeschwindigkeit, d.h. wie schnell Biomasse aus der Nahrung

67


0 500 1000 1500 2000Zeit in Tagen

0

50

100

150

200A

nzah

l


0 500 1000 1500 2000Zeit in Tagen

0

50

100

150

200

Anz

ahl


Abbildung 7.4: Losungskurven des Rauber-Beute Systems mit einer Maximal-kapazitat fur die Hasen (Beute) von K = 250 (links) und K = 500 (rechts).

Abbildung 7.5: Richtungsfeld des verbesserten Lotka-Volterra-Modells mit ei-ner Maximalkapazitat fur die Hasen(K = 250) mit Losungskurven zu verschie-denen Anfangswerten.

aufgebaut werden kann. Die tatsachliche Fressrate soll vom Nahrungsangebot

abhangen. Hierbei kann man davon ausgehen, dass die Rate anfangs mehr oder

weniger linear ansteigt und sich dann dem Maximalwert nahert.

Hierzu geeignet ist eine Monod-Funktion oder auch Michaelis-Menten-Kinetik:

68


w(B) = bMB

B + M

oder allgemeiner

w(B) = µB

B + M

Die erste Gleichung geht fur M → ∞ in die alte Form uber. In der zweiten

Form ist der Bruch dimensionslos und µ ist die maximale Wachstumsrate in

1/Zeiteinheit. Der Bruch geht fur B → ∞ gegen 1.

Die Michaelis-Menten-Kinetik wird in vielen Modellen verwendet, z.B.

• Algenwachstum als Funktion der Nahrstoffkonzentration (Monod-Modell)

• Alkoholabbau als Funktion der Blutalkoholkonzentration

• Aktenbearbeitung als Funktion der Aktenmenge

Eigenschaften

• nimmt bei M die Halfte des Maximalwerts an. M heißt deshalb auch

Halbsattigungskonstante

• die Tangente im Ursprung schneidet den Grenzwert µ bei M

• die Annaherung an den Grenzwert ist sehr langsam

0 M 2M 3M 4M 5M

µ/2

µ

Abbildung 7.6: Graph der Monod-Funktion (Michaelis-Menten-Funktion) mitder Halbsattigungskonstanten M und dem Grenzwert µ.

Man erhalt hiermit folgendes allgemeinere Rauber-Beute-Modell nach Rosen-

zweig und MacArthur:

69


B = r · B ·(

1 − B

K

)

− bMB

B + M· R (7.4)

R = η · b MB

B + M· R − s · R

Der Effizienz-Parameter η gibt an, wie viel der aufgenommenen Nahrung tatsachlich

in Biomasse umgesetzt werden (Knochen werden ausgespuckt). Ein Simulati-

onsergebnis ist in Abbildung 7.7 gegeben. Man sieht, dass das System bei der

gegebenen Parameterkonstellation einen Grenzkreis erreicht.

0 5000 10000Zeit in Tagen

0

50

100

150

200

250

300

Anz

ahl


0 100 200 300 400 500

Beute

0

20

40

60

Räu

ber

Abbildung 7.7: Simulationsergebnis fur das allgemeine Rauber-Beute-Modellnach 7.4 mit den Parametern r = 0.05, K = 250, M = 30, s = 0.05, η = 0.7,b = 0.003. Links ist das Zeitdiagramm, rechts das Phasendiagramm dargestellt.Fur andere Parameterwerte kann sich das System vollig anders verhalten.

Das das System nun sehr viele Parameter enthalt, ist es einfacher, mit der

skalierten Form zu arbeiten (Herleitung siehe Skript von Wolfgang Ebenhoh,

http://eagle.icbm.de/~mathmod/mm/skript_mm.html).

7.1.5 Das skalierte allgemeine Rauber-Beute-Modell

In der skalierten Form erhalt man folgende Gleichungen

B = (1 − εB) · B − R

1 + κB· B (7.5)

R = δB

1 + κB· R − δR

Das System hat nun die Parameter ε, κ und δ, fur die angenommen wird, dass

sie alle großer oder gleich Null sind. Fur ε = 0 und κ = 0 hat man wieder ein

normales Lotka-Volterra-Modell.

Das System hat 3 stationare Zustande

70


1. B∗=0 und R∗ = 0 (trivial)

2. B∗ = 1ε und R∗ = 0 (Beute ohne Rauber, bei ε > 0)

3. B∗ = 11 − κ und R∗ = 1 − κ − ε

(1 − κ)2 (positiver stationarer Zustand)

Der positive stationare Zustand existiert im Zustandsraum (B, R ≥ 0), wenn

κ + ε < 1

erfullt ist. Fur κ + ε = 1 fallt der Zustand (3) mit Zustand (2) zusammen.

Das Verhalten des Systems fur δ = 1, κ = 0.5 und verschiedene ε ist in

Abbildung 7.8 dargestellt.

B

86420

R

3

2,5

2

1,5

1

0,5

0

B

86420

R

3

2,5

2

1,5

1

0,5

0

B

86420

R

3

2,5

2

1,5

1

0,5

0

Abbildung 7.8: Phasendiagramm fur das skalierte allgemeine Rauber-Beute-Modell nach 7.5 mit den Parametern δ = 1, κ = 0.5 und ε = 0.15 (links),ε = 0.3 (Mitte) und ε = 0.6 (rechts) und den Startwerten B0 = 6 und R0 = 1.

Fur δ = 1, κ = 0.5 hat das System die stationaren Zustande

ε 0.15 0.3 0.6

(B∗1 , R

∗1) (0,0) (0,0) (0,0)

(B∗2 , R

∗2) (6.6, 0) (3.3, 0) (1.6, 0)

(B∗3 , R

∗3) (2,1.4) (2,0.8) neg.

7.2 Verhalten 2-dimensionaler autonomer Systeme

7.2.1 Stationare Zustande

Wie bereits gezeigt, erhalt man im 2-dimensionalen Fall die stationaren Zustande

des (autonomen) Systems

X = f(X, Y ) (7.6)

Y = g(X, Y )

71

7.2. Verhalten 2-dimensionaler autonomer Systeme GMM WS2009/2010

in dem

f(x, y) = 0 (7.7)

g(x, y) = 0

gesetzt wird, und die Losungen (X∗, Y ∗) dieses Gleichungssystems bestimmt

werden.

7.2.2 Geschlossene Bahnen

Im Gegensatz zu eindimensionalen Systemen konnen im zweidimensionalen

Fall geschlossene Losungsbahnen auftreten. Im Normalfall ist die geschlossene

Bahn das Bild einer periodischen Losung. Man nennt eine solche periodische

Bahn stabil, wenn alle Bahnen, die in ihrer Nahe starten, zu dieser hin laufen.

Man nennt solche stabilen Bahnen auch Grenzkreise. Periodische Bahnen

konnen auch instabil sein.

7.2.3 Satz von Poincare

Bei autonomen Systemen mit zwei Zustandsvariablen lauft jede Losung ent-

weder

• auf einen stationaren Zustand zu, oder

• nahert sich einer geschlossenen Kurve (Grenzkreis) an, oder

• lauft nach unendlich

7.2.4 Teilgleichgewichte (Isoklinen)

Wir betrachten das System aus Gleichung 7.6. Zuerst wollen wir einen Uber-

blick uber das Losungsverhalten gewinnen. Hierzu betrachten wird die Menge

aller Punkte fur die X = 0 bzw. Y = 0 gilt:

Teilgleichgewicht von X = {X = 0} = {(X, Y )|f(X, Y ) = 0}

Das Teilgleichgewicht von X trennt die Bereiche des Zustandsraums, in denen

die Bewegung nach rechts bzw. links erfolgt. Analog gilt

Teilgleichgewicht von Y = {Y = 0} = {(X, Y )| g(X, Y ) = 0}


die Bewegung nach oben bzw. unten erfolgt.

72


Die Schnittpunkte der Teilgleichgewichte sind gerade die stationaren Zustande.

B

86420

R

3

2,5

2

1,5

1

0,5

0

B

86420

R

3

2,5

2

1,5

1

0,5

0

B

86420

R

3

2,5

2

1,5

1

0,5

0

Abbildung 7.9: Teilgleichgewichte fur das skalierte allgemeine Rauber-Beute-Modell nach 7.5 mit den Parametern δ = 1, κ = 0.5 und ε = 0.15 (links),ε = 0.3 (Mitte) und ε = 0.6 (rechts). Die senkrechte Linie zeigt das Teilgleich-gewicht {R = 0}, die Parabel zeigt das Teilgleichgewicht {B = 0}.

7.2.5 Stabilitat stationarer Zustande

Die Stabilitat eines stationaren Zustands kann wie im eindimensionalen Fall

untersucht werden. Man betrachtet hierzu das System als Vektor

(X

Y

)

=

(f(X, Y )

g(X, Y )

)

(7.8)

Analog zum eindimensionalen Fall entwickelt man die rechte Seite des Differen-

tialgleichungssystems am stationaren Zustand (X∗, Y ∗) nach Taylor bis zum

ersten Glied und erhalt so das linearisierte System (ist fur jeden stationaren

Zustand unterschiedlich!!):

(x

y

)

= J(X∗, Y ∗) ·(

x

y

)

(7.9)

Die Matrix J ist die Jacobi-Matrix:

J(X∗, Y ∗) =

∂f(X∗, Y ∗)∂X

∂f(X∗, Y ∗)∂Y

∂g(X∗, Y ∗)∂X

∂g(X∗, Y ∗)∂Y

(7.10)

73


Das linearisierte System ist am stationaren Zustand (X∗, Y ∗) entwickelt und

gibt die Abweichung (Storung) von diesem an. Zu linearen Abbildungen, Ei-

genwerten etc. siehe Anhang D.

Hat die Matrix J zwei verschiedene Eigenwerte λ1 und λ2, so kann man eine

allgemeine Losung des linearen Systems E.4 angeben:

~z(t) = C1 · eλ1t · ~v1 + C2 · eλ2t · ~v2

Hierbei ist ~v1 Eigenvektor zum Eigenwert λ1 und ~v2 Eigenvektor zum Eigenwert

λ2, C1 und C2 sind freie Konstanten.

Sind beide Eigenwerte reell so handelt es sich bei den Summanden entweder

um exponentielles Wachstum oder exponentiellen Zerfall. Ist also einer der

Eigenwerte großer Null, wachst die Storung, der Zustand ist instabil. Sind beide

Eigenwerte negativ, so verschwindet eine anfangliche Storung, der Zustand ist

stabil.

Fur einen komplexen Eigenwert λ = a + ib gilt

eλt = e(a+ib)t = eat · eibt = eat · (cos bt + i sin bt)

Das System vollzieht also Schwingungen. Ist der Realteil a kleiner Null, so ver-

schwindet die Storung, ist der Realteil a großer Null, so explodieren Storungen.

Bei a = 0 bleiben Storungen erhalten. Uber die Stabilitat des Ursprungssy-

stems kann dann nicht entschieden werden.

Bezuglich der Stabilitat des Systems 7.6 gilt der folgende Satz:

Satz 7.2.1 Stabilitat von stationaren Zustanden autonomer Systeme

1. Der stationare Zustand (X∗, Y ∗) von 7.6 ist asymptotisch stabil, falls

alle Eigenwerte von 7.10 negative Realteile besitzen.

2. Der stationare Zustand (X∗, Y ∗) ist instabil, falls mindestens ein Ei-

genwert von 7.10 einen positiven Realteil besitzt.

3. Sind alle Realteile aller Eigenwerte von 7.10 kleiner oder gleich Null und

mindestens ein Realteil gleich Null, so kann nicht uber die Stabilitat

des stationaren Zustands (X∗, Y ∗) entschieden werden.

74


Die Eigenwerte der Jacobimatrix J := J(X∗, Y ∗) sind durch

λ1,2 =spur J

2±√

(spur J)2

4− det J

gegeben. Fur die Stabilitat des Zustands (X∗, Y ∗) gilt:

spur J det J Stabilitat

spur J < 0 det J > 0 stabil

spur J > 0 beliebig instabil

beliebig det J < 0 instabil

spur J = 0 det J ≥ 0 keine Entscheidung moglich

spur < 0 det J = 0 keine Entscheidung moglich

7.2.6 Stabilitat im Beispielsystem

Das skalierte allgemeine Rauber-Beute-System (7.5)

B = (1 − εB) · B − R

1 + κB· B

R = δB

1 + κB· R − δR

hat die drei stationaren Zustande

1. (B∗, R∗) = (0, 0) (trivial)

2. (B∗, R∗) = (1ε , 0) (Beute ohne Rauber, bei ε > 0)

3. (B∗, R∗) = ( 11 − κ,

1 − (κ + ε)(1 − κ)2 ) (positiver stationarer Zustand)

Der positive stationare Zustand existiert im Zustandsraum (B, R ≥ 0), wenn

κ + ε < 1

erfullt ist. Fur κ + ε = 1 fallt der Zustand (3) mit Zustand (2) zusammen.

Die Jacobi-Matrix des System ist durch

J(B, R) =

1 − 2 ε B − R

(1 + κ B)2 − B

1 + κ B

δR

(1 + κ B)2 δ

(B

1 + κ B− 1

)

(7.11)

75


Stabilitat von (B∗1 , R

∗1) = (0, 0)

Fur (B∗1 , R

∗1) = (0, 0) erhalt man

J(B∗1 , R

∗1) =

1 0

0 −δ

(7.12)

Damit gilt det J < 0, der Zustand ist instabil.


∗2) = (1

ε, 0) (Beute ohne Rauber)

Fur (B∗2 , R

∗2) = (1

ε, 0) erhalt man

J(B∗2 , R

∗2) =

−1 − 1

κ + ε

0 δ

(1

κ + ε− 1

)

(7.13)

Es gilt

det J = δ

(

1 − 1

κ + ε

)

1.Fall κ + ε < 1 (der dritte Zustand existiert):

det J < 0, der Zustand 2 ist instabil.

2. Fall κ + ε > 1 (der dritte Zustand existiert nicht ):

det J > 0 und spur J = −1 + δ

(1

κ + ε− 1

)

< 0, der Zustand 2 ist stabil.

3. Fall κ + ε = 1 (der zweite und dritte Zustand fallen zusammen): det J = 0

keine Entscheidung uber die Stabilitat des Zustands 2 moglich.


∗3) =

(

11 − κ,

1 − (κ + ε)(1 − κ)2

)

Der Zustand 3 existiert fur κ + ε < 1. Fur die Jacobimatrix erhalt man

J(B∗3 , R

∗3) =

−ε1 + κ

1 − κ+ κ −1

δ (1 − (κ + ε)) 0

(7.14)

76


Es gilt

det J = δ (1 − (κ + ε)) > 0 (κ + ε < 1)

und

spur J = −ε1 + κ

1 − κ+ κ

Der Zustand 3 ist also stabil, wenn die Spur negativ ist, also fur

−ε1 + κ

1 − κ+ κ < 0

Dies gilt fur

ε > κ1 − κ

1 + κ=: εkrit(κ)

Das Systemverhalten in Abhangigkeit von ε und κ ist in Abbildung 7.10 dar-

gestellt.

0 1 κ

0

1

ε

Verhalten des Zustands 3

stabil

stabil (Oszillationen)

instabil (Osszilationen)

instabil

Abbildung 7.10: Verhalten des stationaren Zustands 3 (Rauber und Beute exi-stieren) des verallgemeinerten skalierten Rauber-Beute-Systems in Abhangig-keit von κ und ε fur δ = 0.2.

77


78

GMM WS2009/2010 Kapitel 8 – Raumliche Systeme

8 Raumliche Systeme

8.1 Zellulare Automaten

Bei der Beschreibung von Prozessen in der Natur spielt nicht nur die Verande-

rung der Zustandsgroßen im Laufe der Zeit eine Rolle. Haufig mochte man

zusatzlich auch Informationen uber die raumliche Verteilung haben. Stellt man

sich eine Rasenflache vor, und mochte man die Grasmenge auf dieser Flache

modellieren, so kann man naturlich ein Modell entwickeln, dass die Grasmenge

auf dieser Flache im Laufe der Jahreszeiten beschreibt. Dann hat man aber

keinerlei Informationen daruber, ob auf dieser Rasenflache an einigen Stellen

braune Stellen oder dicke Grasbuschel existieren. Hierzu benotigt man ein Mo-

dell, das die Grasflache raumlich beschreibt. Hierzu kann man die Rasenflache

zum Beispiel in Quadrate einteilen. Man kann nun Wachstumsregeln aufstellen,

die das Graswachstum beschreiben. Ist z.B. ein Quadrat von braunen Stellen

umgeben, d.h., dass auf den Nachbarquadraten sehr wenig Gras ist, so wird

das Gras auch in der Mitte schlecht wachsen, da der Schutz vor Wind fehlt.

Sind um ein Quadrat dicke Grasbuschel, so werden die Nahrstoffe knapp und

auch dann kann das Gras schlecht wachsen.

8.1.1 Conway’s Life

Ein erstes solches Modell, genannt LIFE, wurde von Conway im Jahre 1970

zur Beschreibung einer fiktiven Bevolkerung beschrieben. Jedes Quadrat, auch

Zelle genannt, kann hierbei zwei Zustande annehmen, tot (0) oder lebendig

(1). Conway hat nun folgende Regeln aufgestellt:

1. Eine Zelle wird geboren, wenn sie drei lebende Nachbarn hat.

2. Eine Zelle bleibt am Leben, wenn sie zwei oder drei lebende Nachbarn

hat.

3. Eine Zelle stirbt sie an Vereinsamung oder Uberbevolkerung, wenn sie

weniger als zwei oder mehr als drei lebende Nachbarn hat.

4. Nachbarn sind die 8 Zellen um eine Zelle herum.

Mit diesen Regeln ergeben sich interessante Ergebnisse. Es gibt Konstellatio-

nen, die sich uber die Zeit nicht andern, andere die Zyklen durchlaufen, welche

79

8.1. Zellulare Automaten GMM WS2009/2010

die sich uber das Feld bewegen (Gleiter) und solche, die in gewissen Abstanden

Gleiter aussenden.

Abbildung 8.1: Gleiter in Life

Abbildung 8.2: Zyklus in Life

80


8.1.2 Charakteristika von zelluaren Automaten

Allgemein kann man die Charakteristika von zelluaren Automaten wie folgt

beschreiben:

• Entwicklung in Raum und Zeit

• diskrete Anzahl an Zellen

• endliche Anzahl an Zustanden fur jede Zelle (kann man aufweichen)

• Anderung der Zustande in diskreten Zeitschritten

• Zustand hangt von den Nachbarzellen ( und evtl. der Zelle selbst) ab

• alle Zellen Verhalten sich nach den gleichen Regeln

Nachbarschaftsbeziehungen

Der Zustand einer Zelle hangt von den Zustanden seiner Nachbarn ab. Man

kann verschiedene Nachbarschaftsbeziehungen betracheten.

• 4 Nachbarn (oben,unten, rechts und links)

• 8 Nachbarn (alle direkt angrenzenden Zellen)

• die Zelle selbst kann zu den Nachbarn gezaht werden

• entferntere Zellen werden hinzugenommen

Randbedingungen

Man kann die Rander offen lassen, so dass die Randzellen einfach weniger

Nachbarn haben. In diesem Fall wird der Zustand der Randzellen je nach

Automat anders ausfallen als der der Zellen in der Mitte. Will man dies ver-

meiden schliesst man das Feld zu einem Torus. Dann haben alle Zellen gleich

viele Nachbarn.

8.1.3 Zirkularer Raum

Der zirkulare Raum ist ein beruhmtes Beispiel fur Selbstorganisation. Nach

einer geringen Anzahl von Schritten bilden sich Strukturen heraus.

• Es gibt N mogliche Zustande

• Der Zustand der Zelle wird um Eins erhoht, wenn mindestens ein Nachbar

den Zustand der Zelle plus Eins hat.

• Jede Zelle hat 4 Nachbarn (links, rechts, oben und unten)

• Der Zustand 0 wird mit dem Zustand N gleichgesetzt

81

8.2. Transportprozesse GMM WS2009/2010

Abbildung 8.3: Selbstorganisation im zirkularen Raum mit N=6 Zustanden.Zwischen den Bildern liegen jeweils 6 Generationen

8.1.4 Erweiterungen

Man kann neben diskreten Zustanden auch kontinuierliche Werte fur die ein-

zelnen Zellen zulassen. Damit lassen sich Prozesse wie Diffusion und Advektion

beschreiben (siehe Abschnitt 8.2).

8.2 Transportprozesse

Lasst man in zellularen Automaten auch kontinuierliche Werte zu, und betrach-

tet die einzelnen Zellen als Diskretisierung eines kontinuierlichen Bereichs (z.B

eines Sees), so kann man mit der ”Automatenstruktur” auch Stromungspro-

zesse beschreiben.

Von besonderer Bedeutung sind hierbei die Advektion (von lat. advehi = her-

anbewegen; bezeichnet allgemein das Heranfuhren von Dingen) und die Diffu-

sion.

Mit Advektion bezeichnet man z.B. den Transport von einer in Wasser gelosten

Substanz.

Beispiel: Man nimmt einen durchsichtigen Gartenschlauch und iniiziert an ei-

ner Stelle etwas Tinte in den Schlauch.

Macht man dieses Experiment in der Realitat, wird man feststellen, das die

Bande aus Tinte mit der Zeit breiter und heller wird. Dies liegt daran, dass

die Tintenteilchen im realen System das Bestreben haben, die Konzentration

im Schlauch auszugleichen. Dieses Bestreben nennt man Diffusion.

Dreht man nun das Wasser auf, so stromt die Tinte entsprechend der angeleg-

ten Geschwindigkeit durch den Schlauch. Dieses bezeichnet man als Advektion.

82


8.2.1 Diffusion

Die Diffusion kann man wie folgt formalisieren: Angenommen wir betrachten

den Schlauch als eindimensional (x-Richtung) und geben die Tinte als Kon-

zentration c an. Dann ist c(t, x) die Konzentration zur Zeit t am Ort x.

Nach dem 1. Fickschen Gesetz ist die Teilchenstromdichte (der Fluss) J pro-

portional zum Konzentrationsgradienten entgegen der Diffusionsrichtung:

J = −D∂ c

∂x

Hierbei ist D die Proportionalitatskonstante, die Diffusionskonstante. Das Mi-

nuszeichen bewirkt, dass der Fluss entgegen dem Konzentrationsgradienten

wirkt, also von hoherer zu niedrigerer Konzentration. Die ”runden” ∂’s be-

zeichnen die partielle Ableitung.

Einheiten: Ist die Konzentration c z.B. in mol/m3 angegeben und die Zeit in s,

so ist die Einheit des Flusses J mol/m2s. Die Einheit der Diffusionskonstanten

ist m2/s.

Die zeitliche Anderung der Konzentration ergibt sich aus der Kontinuitatsglei-

chung, die besagt, dass keine Masse verloren geht oder entsteht:

∂c

∂t+

∂J

∂x= 0

∂c

∂t= −∂J

∂x= − ∂

∂x

[

−D · ∂ c

∂x

]

Fur eine Diffusionskonstante unabhangig von x ergibt sich

∂c

∂t= D · ∂2c

∂x2

Um die Diffusionsgleichung fur ein System zu losen, muss man eine Anfangs-

bedingung vorgeben. In unserem Fall die Konzentration an jeder Stelle x zur

Zeit t=0 c(0,x).

Bisher haben wir noch keine Angabe uber die Lange des Schlauches gemacht.

Wir sind von einem unendlich langen Schlauch ausgegangen. Will man die

Verteilung in einem endlichen Schlauchstuck beschreiben, so muss man auch

Randbedingungen definieren.

Hier gibt es zwei Moglichkeiten: Entweder man setzt die Konzentration am

Anfang (x=0) und am Ende des Schlauchstucks (x=L) fest, also

83


c(t, 0) und c(t, L) f.a t >= 0 Dirichlet-Problem

oder man gibt die Ableitungen am Rand vor

∂c

∂tc(t, 0) und

∂c

∂tc(t, L) f.a t >= 0 Neumann-Problem

Der Unterschied ist wie folgt:

Sind die Konzentrationen am Rand vorgegeben, so bestimmen diese letztlich

nach einer gewissen Zeit die Konzentration im Schlauchstuck. Vom Rand her

diffundiert Substanz (Tinte) in den Schlauch hinein.

Halt man den Gradienten am Rand fest, z.B setzt in gleich null, dann bedeutet

dies, dass keine Substanz aus dem Schlauch oder in den Schlauch diffundieren

kann. Die Menge im Schlauch bleibt konstant und verteilt sich mit der Zeit

gleichmassig.

Im allgemeinen ist es sehr schwierig, Anfangs- und Randwertprobleme analy-

tisch zu losen. Daher ist es meistens sinnvoll solche Probleme numerisch zu

behandeln. Hierzu gibt es hervorragende Verfahren. Das Grundprinzip soll im

folgenden mit einem sehr elementaren Verfahren erlautert werden:

Numerische Losung

Der Schlauch wird durch ein eindimensionales Gitter simuliert. Hierzu wird

die die x-Achse in Abschnitte der Lange ∆x unterteilt. Jeder Abschnitt be-

schreibt ein kurzes Stuckchen des Schlauches (Gitterzelle), n dem die Substanz

als homogen verteilt angenommen wird. Es wird nun die zeitliche Anderung

der Konzentration in der i-ten Gitterzelle bestimmt. Hierzu werden zuerst die

Gradienten zu den benachbarten Gitterzellen diskretisiert, indem man die je-

weiligen Differenzenquotienten betrachtet:

84


∂c(t, x)

∂t≈ D ·

∂c(t,x+∆x)∂x

− ∂c(t,x)∂x

∆x

≈ D ·c(t,x+∆x)−c(t,x)

∆x− c(t,x)−c(x−∆x)

∆x

∆x

= D · c(t, x + ∆x) − 2 · c(t, x) + c(x − ∆x)

∆x2

Die zeitliche Diskretisierung (Euler) fuhrt zu:

c(t + ∆t, x) = c(t, x) + D · c(t, x + ∆x) − 2 · c(t, x) + c(x − ∆x)

∆x2· ∆t

Betrachtet man nun Gitterzellen mit Abstand ∆x so erhalt man folgende For-

mel:

c(t + ∆t, i) = c(t, i) + D · c(t, i + 1) − 2 · c(t, i) + c(i − 1)

∆x2· ∆t

Diese Art der Diskretisierung nennt man forward in time und centered in space.

Man muss nun ∆x und ∆t so wahlen, dass tatsachlich ein Konzentrationsaus-

gleich geschaffen wird, und dass das diskretisierte System nicht osszilliert. Dies

ist fur

D <1

2· ∆x2

∆t

der Fall.

8.2.2 Advektion (Transportgleichung)

Die Beschreibung der Advektion, also einer Stromung, die durch eine konstante

Stromungsgeschwindigkeit v > 0 bestimmt ist, lasst sich im eindimensionalen

fall wie folgt beschreiben

∂c(t, x)

∂t+ v · ∂c(t, x)

∂x= 0

gegeben.

85


Eine sehr einfache Diskretisierung ist durch

c(t + ∆t, i) = c(t, i) + v · c(t, i − 1) − c(t, i)

∆x· ∆t

gegeben. Hierbei ist die raumliche Ableitung nicht symmetrisch sondern strom-

aufwarts verschoben. Man nennt dies daher auch Upstream-Schema. Es ist

stabil, wenn

v · ∆t

∆x<= 1

gilt.

Das Problem bei der Diskretisierung der Advektion ist, dass ∆x und ∆t genau

zueinander ”passen” mussen. Dies kann man aber nur bei konstanter Geschwin-

digkeit v erreichen. Ansonsten ”verschmiert” das Signal, d.h. die Bande aus

Tinte wird immer breiter. Dieses sieht ahnlich wie eine Diffusion aus. Daher

wird dieses Phanomen, dass sich aufgrund der Diskretisierung ergibt, auch

numerische Diffusion genannt.

8.2.3 2D Advektion

Eine sehr einfache Moglichkeit, diese Transportprozesse im zweidimensionalen

Fall zu beschreiben, funktioniert ahnlich wie bei denAutomaten.

Man betrachtet jetzt nicht nur die Nachbarzellen rechts und links, sondern

auch oben und unten. Man definiert also fur jede Zelle vier Austauschflusse.

Diese Transportprozesse konnen mit Differentialgleichungen gekoppelt werden.

86

GMM WS2009/2010 Kapitel 9 – Ausbreitung von Borkenkafern (W. Ebenhoh)

9 Ausbreitung von Borkenkafern (W. Ebenhoh)

Naturlicherweise stirbt ein Teil eines Waldes gelegentlich, das ist Teil seiner re-

gularen Erneuerung. Er stirbt an Waldbranden, an Sturmschaden, an Trocken-

heit, an Schadlingsbefall, und stets folgt dem Sterben eine Erneuerung. Jeder

Teil des naturlichen Waldes erneuert sich etwa alle 30 bis 200 Jahre. Alter Wald

verhindert das Aufkommen von Samlingen. Viel zu effektiv sammeln die Kro-

nen das Licht und die Wurzeln Nahrstoffe und Wasser. Samlinge haben nur in

Lichtungen und an den Randern eine Chance – und eben nach einem kollekti-

ven Absterben eines großeren Bestandes von altem Wald. Alternder Wald wird

gegen die oben genannten Gefahren immer anfalliger, irgendwann setzt dann

das von außen induzierte Sterben ein, fur jeden Standort je nach Anfalligkeit

zu einem anderen Zeitpunkt, aber infolge der oft weiten raumlichen Ausdeh-

nung der Ursachen (z.B. ein besonders trockener Sommer in Mitteleuropa) an

vielen Orten gleichzeitig.

In diesem Kapitel wird ein sehr einfaches Modell einer Walderneuerung kon-

struiert. Im ersten Abschnitt wird die Entwicklung der Borkenkaferpopulation

isoliert betrachtet. Das heißt, zunachst wird die Ruckwirkung der Borkenkafer

auf den Wald vergessen. Naturlich ist das der wichtigste Prozess im Gesamt-

system”Wald“, jedenfalls im Modell. Aber um das

”Subsystem“ Borkenkafer

zu verstehen, wird zunachst die Borkenkaferpopulation in einen Wald mit fest-

gehaltener Altersstruktur gesetzt. Wir lernen dann, in welcher Art Wald sie

sich schnell oder langsam vermehrt. Im zweiten Abschnitt beschreiben wir die

Wirkung einer hohen Borkenkaferdichte auf den Wald als”Verjungung“. Die

Annahme dabei ist, dass die alten Baume sterben und Platz fur Samlinge ma-

chen. Das Modell ist naturlich ubertrieben einfach, da es nur zwei Variable hat,

namlich Borkenkaferdichte und Alter des Waldes. Zuletzt wird das Modell als

raumliches Modell auf einer Automatenstruktur betrachtet.

9.1 Borkenkafer-Dynamik

Borkenkafer vermehren sich und werden von Vogeln gefressen oder sterben an

eigenen Krankheiten.

Als erste (und in diesem Abschnitt einzige) Systemvariable fuhren wir die

Borkenkaferdichte B ein. Als leicht vorstellbare Maßeinheit fur B konnte man

etwa”Zahl der Kafer pro Ar“ verwenden. Unser Modell wird aber nicht an die

Wirklichkeit heranreichen, in der diese Zahl wohl um viele Großenordnungen

87

9.1. Borkenkafer-Dynamik GMM WS2009/2010

schwankt. Die Prozesse, die die Anderung der Borkenkaferdichte B beschrei-

ben, sind Wachstum und Mortalitat. Beide Prozesse sind dichteabhangig, und

außerdem noch abhangig vom Alter A des Waldes. Das Modell besteht in dieser

Entwicklungsstufe aus Gleichung (9.1).

B = Wachstum(B, A) − Mortalitat(B, A) (9.1)

Fur das Wachstum bietet sich als fast einfachste Moglichkeit eine Art logisti-

scher Ansatz an, wobei die Wachstumsrate bei niedrigen Dichten vom Alter A

des Waldes abhangt. Alte Baume sind anfalliger, deshalb wird r(A) in Ansatz

(9.2) mit dem Alter A steigen.

Wachstum = r(A) · B ·(

1 − B

K

)

(9.2)

Hier ist K die Kapzitat der Kafer, eine Wachstumsgrenze die nie erreicht wird,

da ja eine Mortalitat in (9.1) vorgesehen ist. Nur dieser Teil der Gesamtmorta-

litat wird explizit durch eine Formel beschrieben (9.3). Die Wachstumsgrenze

ist dagegen Ausdruck einer”impliziten“ Mortalitat, die nicht weiter spezifiziert

ist und zusatzlich zur explizit angegebenen Mortalitat (9.3) wirkt.

Mortalitat = µB2

B2 + M2(9.3)

Der Ansatz (9.3) ist eine sigmoide Kurve. Ohne die Quadrate ware (9.3) ein

Monod-Ansatz, also ein Predationsterm mit Sattigung. Der Sattigungswert µ

gibt an, wie viel Borkenkafer maximal von seinen Fressfeinden, den Vogeln, ge-

fressen werden konnen (pro Tag und Ar). Naturlich ist die Sache komplizierter.

Einerseits zieht ein reichhaltiges Nahrungsangebot mehr Vogel an, andererseits

fuhrt Uberfluss und Einseitigkeit der Nahrung auch zur Abwanderung von

Vogeln. Im jungen Wald sind wahrscheinlich mehr Vogel als im alten. Deshalb

konnte man vermuten, dass die (maximale) Sterberate µ mit A abnimmt. Wir

ignorieren jedoch der Einfachheit halber eine mogliche Abhangigkeit µ(A). Die

Michaeliskonstante M gibt an, bei welcher Borkenkaferdichte diese mit halb-

maximaler Intensitat gefressen werden.

Die quadratische Abweichung bewirkt, dass bei sehr niedrigen Dichten die Mor-

talitat der Borkenkafer sehr klein ist. Die Schadlinge konnen so von den Vogeln

nicht ausgerottet werden. Man kann vermuten, dass die Kafer umso besser ver-

steckt sind, je weniger von ihnen da sind. Außerdem konzentrieren Rauber ihre

88


Wahrnehmung auf den haufigsten Beutetyp und neigen dazu”seltene“ Beute

zu ignorieren.

Das Alter des Waldes A wird spater die zweite Systemvariable. Hier ist A

eine”Konstante“. Wir normieren das Alter auf den Bereich 0 bis 1. Damit ist

A eine dimensionslose Große. Man kann sich etwa vorstellen, dass mit A =

1 ein Wald aus hundertjahrigen oder alteren Baumen beschrieben wird. Da

unser Modell so einfach wie moglich sein soll, wird die Wachstumsrate der

Borkenkafer proportional zu A angesetzt, (9.4).

r(A) = r · A (9.4)

Die Doppelfunktion des Buchstaben r als Wachstumsfunktion r(A) und als

Ratenkonstante r wird keine Schwierigkeiten machen, da im folgenden nur noch

letztere auftritt. Wird (9.2) bis (9.4) in die Systemgleichung (9.1) eingesetzt,

ergibt sich Gleichung (9.5).

B = r · A · B ·(

1 − B

K

)

− µB2

B2 + M2(9.5)

Die Fixpunkte des Systems erhalt man als Schnittpunkte des Wachstumsterms

und des Mortalitatsterms (Abbildung 9.1). Je nach Alter des Waldes gibt es

zwei, drei oder vier Schnittpunkte. Die erste ist die triviale Losung B∗1 =

0, die anderen sind die Losung einer kubische Gleichung. Da wir analytisch

den Zusammenhang B∗(A) nur schwer berechnen konnen (wer lost schon gern

kubische Gleichungen), betrachten wir die Umkehrfunktion A (B∗) und plotten

die Kurve mit vertauschten Achsen:

A(B∗) =µ · B∗

r ·(1 − B∗

K

)· (B∗2 + M2)

B∗ > 0 (9.6)

In Abbildung 9.2 ist die Losung fur K = 5000 Kafer pro Ar, M = 500 Kafer

pro Ar, r = 0.15 pro Tag und µ = 100 pro Ar und Tag dargestellt. Es gibt

ein breites Intervall auf der A-Achse, in dem zu jedem Alter drei positive

stationare Zustande B∗(A) koexistieren. (Fur großere M verkleinert sich das

Intervall, fur kleinere M dehnt es sich bis uber A = 1 aus). Das Modell kann

zu einer sogenannten Hysterese fuhren. Um die Hysterese zu verstehen, drehen

wir im Gedanken am Parameter A. Bei A = 0.4 existiert (neben B∗ = 0)

zunachst nur einen Gleichgewichtswert mit sehr kleiner Borkenkaferdichte B∗.

Bei Erhohung von A steigt B∗ nur langsam an und erreicht bei A ≈ 0.75 etwa

89

9.2. Wald-Dynamik GMM WS2009/2010

0 1000 2000 3000 4000 5000Anzahl Borkenkäfer pro Ar

0

1×104

2×104

3×104

4×104

5×104

6×104

7×104

Anz

ahl p

ro A

r un

d Ja

hr

MortalitätKäferwachstum, Waldalter = 1.0Käferwachstum, Waldalter = 0.7Käferwachstum, Waldalter = 0.5

Wachstum und Mortalität

Abbildung 9.1: Vergleich von Wachstum und Mortalitat fur verschiedene AlterA des Waldes. Bei A = 0.7 haben die Graphen von Mortalitat (blau) undWachstum vier Schnittpunkte, in den anderen Fallen zwei, jeweils einschließlichdes Nullpunktes. Die Schnittpunkte entsprechen den stationaren Zustanden,sie sind abwechselnd stabil und instabil , wobei der Nullpunkt in jedem Falleinstabiler stationarer Zustand ist.

den Wert B∗ = 550. Bei noch hoherem A hort dieser Ast von B∗(A) auf zu

existieren. Das System muss sich einen anderen stabilen stationaren Zustand

suchen. Faktisch bedeutet das, die Borkenkaferdichte explodiert und pendelt

sich auf dem hohen Wert B∗ ≈ 4200 ein. Eine nachfolgende Erniedrigung

von A muss bis unter den Wert von etwa 0.5 fuhren, ehe die hohe Dichte

zusammenbricht.

9.2 Wald-Dynamik

Im diesem Abschnitt wird das Alter A des Waldes dynamisch modelliert. Die

Hysterese nach Abbildung 9.2 legt nahe, dass man mit Borkenkaferausbruchen

rechnen kann, wenn die Parameter entsprechend gewahlt werden. Neben Glei-

chung (9.1) bzw. (9.5) benotigt man eine weitere Systemgleichung, mit der das

”Alter des Waldes“ A durch die Prozesse Altern und Verjungung verandert

wird. Das Alter eines Baumes ist eine einfache Sache, aber nicht das Alter

des Waldes. Wir wollen als Systemvariable A ein irgendwie gemitteltes Alter

der Baume verstehen, wobei uns klar ist, dass es nicht dasselbe ist, ob der

Wald aus gleich alten 50-jahrigen oder zur halben Flache aus 10- und 90-jahri-

90


0 0.2 0.4 0.6 0.8 1A

0

1000

2000

3000

4000

5000

B*

Lage der Fixpunkte B*

Abbildung 9.2: Stationare Borkenkaferdichten B∗ als Funktion des Alters Ades Waldes. Im Bereich 0.55 < A < 0.75 koexistieren vier stationare Zustande(einschließlich des trivialen Falles B∗ = 0). Die Zustande sind abwechselndinstabil (gepunktet) und stabil, beginnend mit B∗ = 0 instabil.

gen Baumen besteht. A ist eine integrierte Waldeigenschaft, die wegen ihrer

Globalitat nicht leicht messbar ist, aber ein intuitives Verstandnis ermoglicht.

Die Normierung von A auf das Intervall 0 bis 1 soll beibehalten werden. Das

Waldsystem mit B und A wird nun durch (9.7) beschrieben:

B = Wachstum(B, A) − Mortalitat(B) (9.7)

A = Altern(A) − Verjungung(B, A)

Die Spezifizierung der Prozesse Altern und Verjungung erfordert wie immer

das Erfinden von Formeln. Die Alterungsrate wird als linear angenommen. Je

Alter der Wald desto geringer die Alterungsrate:

Altern(A) = α · (1 − A) (9.8)

Der Wachstumsparameter α wird etwa bei 0.01 pro Jahr liegen, was dem Baum

hundert Jahre Zeit gibt, seine volle Große zu erreichen. Anders als das Al-

tern ist die Verjungung biologisch ein schwieriger Prozess. Zuerst mussen die

alten Baume durch Borkenkafernbefall absterben. Im Modell treten Sturme,

Trockenheit, Luftverschmutzung und andere Krankheiten nicht auf. Dann mussen

die Baumleichen umfallen und zersetzt und die Lucken im Wald durch Samlinge

91

9.2. Wald-Dynamik GMM WS2009/2010

neu besiedelt werden. Der Wald enthalt junge und alte Baume nebeneinander.

Die Formel fur den Prozess Verjungung muss eine extreme Zusammenfassung

der dazu notigen Teilprozesse sein. Das Tatsache, dass A eine integrierte, qua-

litative Waldeigenschaft ist, erfordert auch fur die Verjungung eine ebenso

integrierte Betrachtung. Wir konnen annehmen, dass sie nur bei hohen Bor-

kenkaferdichten stattfindet, aber dann sehr schnell erfolgt:

Verjungung(B, A) = b · A ·(

B

K

)3

(9.9)

Durch den Exponent 3 in (9.9) schaffen geringe Borkenkaferdichten B kaum

eine Verjungung.

Die Modellgleichungen sind in (9.10) zusammengefasst.

B = r · A · B ·(

1 − B

K

)

− sB2

B2 + M2(9.10)

A = a · (A < 1) − b · A ·(

B

K

)2

Mit diesen Gleichungen ist es moglich, periodische Ausbruche von Borkenkafer-

plagen zu simulieren (Abbildung 9.3).

0 10 20 30 40 50Jahre

0

1000

2000

3000

4000

5000

Anz

ahl p

ro A

r

Borkenkäferdichte

0 10 20 30 40 50Jahre

0.0

0.2

0.4

0.6

0.8

1.0 Alter des Waldes

Abbildung 9.3: Simulation von Borkenkaferausbruchen.

Ganz befriedigend ist das in Abbildung 9.3 dargestellte Simulationsergebnis

noch nicht. Die Lange der Borkenkaferplagen ist recht groß gegenuber den”nor-

malen“ Intervallen. Es ware realistischer, wenn nach die Plage noch schneller

wieder voruber ware. Ob das mit diesem einfachen Modell uberhaupt moglich

ist, kann nur mit einer Sensitivitatsanalyse herausgefunden werden.

92


9.3 Raumliche Ausdehnung des Modells

Mit Hilfe der Struktur der zellularen Automaten kann das Modell erweitert

werden. Hierzu teilen den Wald in Quadrate der Große 10 m x 10 m (1Ar2) ein.

Wir nehmen an, dass jeweils ein Teil der Kafer in die benachbarten Zellen (4

Nachbarn) abwandert. Zusatzlich kann man annehmen, dass diese Wanderung

durch eine vorgegebene Windrichtung nicht symmetrisch ist.

9.3.1 Bekampfung

Die Bekampfung der Kafer erfolgt durch das Schlagen von Schneisen. Im Mo-

dell kann man das simulieren, indem man die Kaferzahl in einigen benachbar-

ten Spalten stark dezimiert oder auf Null setzt und das Waldalter herabsetzt

(Schlagen alter Baume).

Man darf bei Betrachtung dieser Modellergebnisse naturlich keinen Moment

vergessen, dass sie mit einem sehr einfachen Modell erzielt wurden, das nur

einige prinzipielle Mechanismen erklaren und nicht die Realitat quantitativ

widerspiegeln soll. Seine Parameter sind nicht an die Realitat angepasst und

vielleicht nicht einmal anpassbar, da wesentliche Komponenten des naturlichen

Systems im Modell nicht erfasst wurden.

93

9.3. Raumliche Ausdehnung des Modells GMM WS2009/2010

94

GMM WS2009/2010 Kapitel 10 – Mathematische Epidemiemodelle

10 Mathematische Epidemiemodelle

Unter einer Epidemie (aus dem Griechischen fur”im Volk verbreitend

”) ver-

steht man ein stark gehauftes ortlich und zeitlich begrenztes Vorkommen einer

Erkrankung, insbesondere von Infektionskrankheiten. Bei einer Ausbreitung

uber meherer Lander und Kontinente spricht man von einer Pandemie.

Mathematische Epidemiemodelle beschreiben die Dynamik und Ausbreitung

von Epidemien. Sie liefern ein Verstandnis der Auswirkungen krankheitssep-

zifischer Eigenschaften auf den Epidemieverlauf, wie z.B. die Auswirkung der

Lange der Inkubationsszeit. Die Inkubationszeit, ist die Zeitspanne zwischen

der Infektion (Eindringen des Erregers in den Korper) bis zum Auftreten der

ersten Symptome (bei Masern 8-14 Tage). Die Inkubationszeit ist zu unter-

scheiden von der Latenzzeit, Zeitspanne zwischem den Eindringen des Erregers

und dem Beginn der Infektiositat.

Historischer Uberblick (sehr unvollstandig)

• Hippokrates 459 - ca. 370 v. Chr.: Begrunder der wissenschaftlichen Medi-

zin, erste Aufzeichnungen von Krankheitsverlaufen

• Louis Pasteur 1822-1895: Begrunder der Mikrobiologie, Entwicklung der

Schutzimpfungen gegen Huhnercholera, Milzbrand und Tollwut

• Robert Koch 1843 - 1910: Entdecker des Tuberkullose-Bazillus

• 19. Jahrhundert: Aufzeichnung von Krankheitsverlaufen und Erklarungsver-

suche von Epidemien (Cholera, Typhus, Pocken)

• Ende des 19 Jahrhundert: erste Modelle

• Kermack & McKendrick (1927): Klassisches Epidemiemodell (SIR), Schwel-

lenwertsatz der Epidemiologie

Wichtige Faktoren bei der Modellierung von Epidemien

• Ubertragungswege Zwischenwirte, Blutkontakt, Tropfcheninfektion

• Immunitat durch Erkrankung, durch Impfung

• Latenzzeit und Inkubationszeit

• Krankheitsverlauf Dauer, Letalitat, bei Kindern/Erwachsenen

• Populationsstruktur Dichte, Verteilung

• Reaktion auf Fallzahlen Schulschließungen, Impfungen, Quaratane, Verhal-

tensanderungen

95

10.1. Das klassische Epidemiemodell nach Kermack & McKendrick GMM WS2009/2010

Statistische Kenngroßen von Epidemien

• Pravalenz

relative Pravalenz: Anteil der Infizierten/Erkrankten an der Gesamtbevolke-

rung

absolute Pravalenz: Anzahl der Infizierten/Erkrankten

• Inzidenz

relative Inzidenz: Anteil der Neuinfizierten/Neuerkrankten an der Gesamt-

bevolkerung in einem bestimmten Zeitraum

absolute Inzidenz: Anzahl der Neuinfizierten/Neuerkrankten an der Ge-

samtbevolkerung in einem bestimmten Zeitraum

kumulative Inzidenz: Gesamtzahl der Neuinfizierten/Neuerkrankten bis zu

einem bestimmten Zeitpunkt

10.1 Das klassische Epidemiemodell nach Kermack & McKendrick

Es wird eine Bevolkerung mit einer konstanten Große N zugrunde gelegt.

In dieser Bevolkerung seien alle Personen fur die modellierte Krankeit gleich

anfallig. Die Ubertragung der Krankheit findet durch den Kontakt einer infi-

zierten Person I (engl. infectious) mit einer anfalligen Person S ( engl. suszep-

tible) statt. Eine Person, die infiziert wurde, gesundet nach einer gewissen Zeit

und wird eine immune Personen R (engl. removal). Es wird hierbei nicht zwi-

schen infizierten und erkrankten Personen unterschieden. Aus den englischen

Bezeichnungen leitet sich der Name SIR-Modell ab:

S(t) Anzahl der anfalligen Personen zur Zeit t,

I(t) Anzahl der infizierten Personen zur Zeit t und

R(t) Anzahl der gesundeten Personen zur Zeit t

Man nimmt an, dass zu Beginn einer Epidemie eine sehr kleine Anzahl von

Personen infiziert ist: I(0) << N und die restlichen Personen anfallig sind

S(0) = N − I(0), R(0) = 0.

Wird eine Infektion durch den Kontakt eines Anfalligen mit einem Infizierten

hervorgerufen, so ist das Infektionsrisiko eines Anfalligen umso großer, je mehr

Infizierte in der Bevolkerung sind. Die Ansteckungsrate B wird mit der Zahl

der Infizierten steigen:

B = β I

96


Hierbei ist β ein Parameter, der die Infektiositat angibt. Geht man von einer

Population mit konstanter Populationsgroße N aus, so betragt der Anteil der

Infizierten IN

. Betrachtet man eine anfallige Person, so hat diese innerhalb eines

Zeitintervalls k Kontakte zu anderen Personen. Die Anzahl der Kontakte zu

infizierten betragt somit k · IN

. Von diesen potentiell gefahrlichen Kontakten

fuhrt nur ein Teil p zur Infektion (Ansteckungswahrscheinlichkeit oder An-

steckungsrisiko). Die Zahl der fatalen Kontakte (also derer, die zur Infektion

fuhren) betragt somit fur eine anfallige Person im betrachteten Zeitintervall

p · k · IN

. Dies ist die Infektionsrate B

B = p · k · I

N= β I mit β =

p · kN

(10.1)

Die Konvertierung von infiziert zu immun hangt von der Gesundungsrate γ

und der Anzahl der aktuell Infizierten ab. Man erhalt somit folgendes Modell:

S = −β I S

I = β I S − γ I (10.2)

R = γ I

Fuhrt man zusatzlich den Parameter ρ = γβ

ein, so erhalt man das Modell

S = −β S I (10.3)

I = β (S − ρ) I

R = β ρ I

Fur eine konstante Bevolkerungsgroße gilt zu jeder Zeit S + I + R = N . Das

System lasst sich also auf zwei Gleichungen reduzieren.

Diesem System sieht man den Epidemieverlauf bereits an. Fur S > ρ gilt I > 0,

fur S < ρ gilt I < 0. Dass bedeutet das die Epidemie ausbrechen kann, wenn

die Anzahl der Anfalligen oberhalb des Schwellenwerts ρ liegt. Anderenfalls

ebbt sie ab. Das Systemverhalten ist in Abbildung 10.1 dargestellt.

10.1.1 Stabilitatsanalyse des SIR-Modells

Es reicht, das zweidimensionale Modell zu untersuchen:

S = −β S I (10.4)

I = β (S − ρ) I

97

10.1. Das klassische Epidemiemodell nach Kermack & McKendrick GMM WS2009/2010

800Suszeptible S

0

10

20

30

40

Infiz

iert

e I

Abbildung 10.1: Losungen im S-I-Phasenraum des klassischen Epidemiemo-dells nach Kermack & McKendrick fur β = 0.005 und ρ = 800.

Teilgleichgewichte

{

S = 0}

= {(S, I) | I = 0 oder S = 0}{

I = 0}

= {(S, I) | I = 0 oder S = ρ}

Stationare Zustande

Die stationaren Zustande des Modells sind alle Punkte (S∗, I∗) = (S, 0), S ≥ 0.

Damit besteht die gesamte S-Achse aus Fixpunkten.

Stabilitat

Da die gesamte S-Achse aus Fixpunkten besteht, kann keine asymptotische

Stabilitat erwartet werden.

Die Linearisierung des Systems 10.4 ergibt die Jacobimatrix

J(S∗, I∗) =

−β I∗ −β S∗

β I∗ β (S∗ − ρ)

(10.5)

Fur (S∗, I∗) = (S, 0) gilt spur J = β (S − ρ) und det J = 0, und man erhalt

die Eigenwerte

λ1 = β (S − ρ)

λ2 = 0

Fur S > ρ gilt λ1 > 0. In diesem Fall ist der Fixpunkt (S, 0) instabil, die

98


Epidemie bricht aus.

Schwellenwertsatz der Epidemiologie (ohne Beweis)

Fur das Epidemiemodell 10.3 gilt: Sei die Anzahl der Anfalligen zu Beginn der

Epidemie S(0) = ρ + ε fur ein ε > 0 und die Anzahl der Infizierten sehr viel

kleiner als die Anzahl der Anfalligen I(0) << S(0), so werden schliesslich ca.

2 ε Personen erkranken.

Hierbei handelt es sich um eine Approximation. Diese unterschatzt den tatsachli-

chen Epidemieverlauf etwas.

10.2 Masern

Masern sind hochansteckende fieberhafte, exanthematische Viruserkrankun-

gen, die nur beim Menschen vorkommen. Schwere Krankheitsverlaufe mit Kom-

plikationen in Form von Pneumonien, Mittelohrentzundungen, Bronchitiden

sowie der lebensbedrohenden akuten postinfektiosen Enzephalitis (0,1 % al-

ler Erkrankungsfalle) sind moglich. Daruber hinaus kommt es bei etwa 1 pro

100.000 Erkrankungen zum Auftreten der sog. subakuten sklerosierenden Pa-

nenzephalitis (SSPE), die immer zum Tod fuhrt. Masernerkrankungen konnen

durch eine Impfung effektiv verhindert werden. Die Eliminierung der Masern

bis zum Jahr 2010 ist ein erklartes Gesundheitsziel der WHO. Um dies zu er-

reichen, sollten 95 % der Bevolkerung durch Impfung bereits im Kindesalter

geschutzt sein (s. a. RKI, Epid. Bull. 10/2004). Masern sind weltweit verbrei-

tet. Aufgrund der hohen Infektiositat treten Masern meist als Kinderkrankheit

auf und hinterlassen eine lebenslange Immunitat. Masern werden durch Tropf-

cheninfektion ubertragen, also z.B. durch Husten, Niesen oder Sprechen. Die

Inkubationszeit betragt zwischen ca. 9 und 14 Tagen. Weltweit sind die Ma-

sern mit jahrlich 31 Millionen Erkrankungen und 614.000 Todesfallen (2002)

weiterhin eine Hauptursache fur Todesfalle im Kindesalter, die durch Impfung

vermeidbar waren (Quelle: RKI).

Epidemiologisch gesehen ist die Bevolkerungsdichte ein wesentlicher Parame-

ter fur den Verlauf der Epidemie (siehe Spektrum der Wissenschaft, 1984.

Insulare Epidemien). Die Ansteckungsrate variiert im Laufe des Jahres, sie ist

z.B wahrend der Schulferien geringer.

Um die Auswirkung der Bevolkerungsdichte zu berucksichtigen, muss im klas-

99

10.2. Masern GMM WS2009/2010

sischen Epidemiemodell die Ansteckungsrate B modifiziert werden. Desweitern

wird eine Zu- und Abwanderung von Teilen der Bevolkerung berucksichtigt.

Masern-Modell

Es wird angenommen, dass die Bevolkerungsgroße N konstant ist. Es findet

eine Migration statt, wobei pro Tag µ · N gesunde Personen pro Quadratkilo-

meter einwandern und dieselbe Personenzahl das Gebiet verlasst. Die auswan-

derenden Personen setzen sich anteilig aus Gesunden, infizierten und immunen

zusammen. Es wird im Modell nicht zwischen infizierten und erkrankten Per-

sonen unterschieden. Erkrankte Personen werden mit einer Gesundungsrate ρ

immun.

S = µ · N − B · S − µ · SI = B · S − ρ · I − µ · IR = ρ · I − µ · R

Zustandsvariablen:

S: Dichte der Suszeptiblen in Einwohner/km2

I: Dichte der Infizierten in Einwohner/km2

R: Dichte der Immunen in Einwohner/km2

Parameter:

N: Bevolkerungsdichte in Einwohner/km2

µ: Zuwanderungsrate in in 1/Tag

B: Infektionsrate in 1/Tag

ρ: Gesundungssrate von infiziert nach immun in 1/Tag

Geht man im klassischen Epidemiemodell von einer konstanten Populations-

große N aus, so kann man die Infektionsrate B entsprechend Gleichung 10.2

interpretieren:

B = p · k · I

N

Hierbei ist p das Ansteckungsrisiko bei einem fatalen Kontakt und k die Anzahl

der Kontakte einer anfalligen Person in einem Zeitintervall.

Bei hoher Bevolkerungsdichte (z.B. in Großstadten) ist die Zahl der Kontak-

te pro Zeiteinheit großer als bei niedriger Bevolkerungsdichte (z.B. auf dem

100


Land). Da Masern uber Tropfcheninfektion ubetragen werden, kann man sich

vor einer Infektion nicht wirklich schutzen, so dass die Infektionsrate bei hoher

Dichte zunehmen wird.

Die Kontaktzahl k soll also abhangig von der Bevolkerungsdichte N sein. Ist

κ die Anzahl der Kontakte bei einer Bevolkerungsdichte Nκ, so ist folgende

Annahme sinnvoll:

k(N) =κN · N

Nκ

Dann erhalt man die Infektionsrate B:

B = p · κN · NNκ

· I

N= p · κN

Nκ· I

Man kann nun uber Nκ oder κ die Kontaktzahl steuern. Angenommen die

Kontaktzahl k betragt bei einer Bevolkerungsdichte von Nκ = 200 Personen

pro Quadratkilometer κ = 5 Kontakte pro Tag, dann wird sie fur N=1000

Personen pro Quadratkilometer das Funffache betragen.

Die Anzahl der Kontakte κ bei der Bevolkerungsdichte Nκ hangt von der

Jahreszeit ab. Wahrend der Sommermonate, wenn Schulferien sind, ist die

Zahl der Kontakte geringer. Hierzu modifiziert man κ so, dass es im Sommer

niedriger als im Winter ist. Dies kann durch eine Rechteckkurve oder auch

durch eine Cosinuskurve geschehen:

κ(t) = κ ·(

1 + c

2+

1 − c

2· cos

2π t

365

)

Der Parameter c ∈ [0, 1] gibt an, auf wieviel Prozent des ursprunglichen Wert

κ gesenkt wird.

10.2.1 Simulationsergebnis

Die Simulationsverlaufe in Abbildung 10.2 wurden mit folgender Parametri-

sierung berechnet:

101

10.2. Masern GMM WS2009/2010

N = 200, 500, 900 Bevolkerungsdichte in Einwohner/km2

µ = 0.001 Migrationsrate in 1/Tag

p = 0.05 Ansteckungsrate in 1/Tag

ρ = 0.1 Gesundungsrate in 1/Tag

κ = 5 Kontaktzahl bei der Bevolkerungsdichte Nκ

Nκ = 200 Referenzbevolkerungsdichte in Einwohner/km2

I(0) = 1.0 anfangliche Infiziertendichte in Einwohner/km2

R(0) = 0.8 · N anfangliche Immunendichte in Einwohner/km2

S(0) = N − R(0) − I(0) anfangliche Suszeptiblendichte in Einwohner/km2

0 5 10 15 20Zeit in Jahren

0

2

4

6

Infiz

iert

enan

teil

in %

200 Einw. pro km2

500 Einw. pro km2

900 Einw. pro km2

Konstante Kontaktzahl

0 5 10 15 20Zeit in Jahren

0

2

4

6

Infiz

iert

enan

teil

in %

Sinusförmige Kontaktzahl (c=0.7)

Abbildung 10.2: Simulationsergebnisse fur die Masernepidemie mit konstanterKontaktzahl (oben) und sinusformiger Kontaktzahl mit einem Minimum imSommer (c=0.7, unten).

Das Model beschreibt ausschliesslich die Abhangigkeit der Infiziertenzahlen

von der Kontaktzahl. In einem realistischeren Modell muss zusatzlich der Impf-

status der Bevolkerung berucksichtigt werden.

102

GMM WS2009/2010 Anhang A – Funktionen

A Funktionen

Eine Funktion f ist eine Vorschrift, die jedem Element einer Menge D genau

ein Element der Menge W zuordnet. Die Menge D heisst Definitionsbereich,

die Menge W Wertebereich:

f : D −→ W (A.1)

x 7→ f(x)

x heisst unabhangige Variable oder Argument. Eine Funktion wird auch Ab-

bildung genannt.

Die Menge aller Punkte G := {(x, f(x))|x ∈ D} der Funktion (A.1) heisst

Graph der Funktion.

Beispiel A.0.1

f : [0, 5] −→ R (A.2)

x → x2

Das Argument x ist anschaulich gesehen ein Platzhalter und ist beliebig durch

einen anderen Buchstaben austauschbar. Durch

f : [0, 5] −→ R

t → t2

ist dieselbe Funktion wie in (A.2) definiert.

In der Schule schreibt man statt f(x) haufig einfach y. Wenn klar ist, welches

in der Vorschrift die unabhangige Variable ist, so ist das ok:

y = 5 · x2 (A.3)

Aber was ist,wenn die Vorschrift y = a · b2 geben ist. Man kann so nicht

erkennen, ob a, b oder sogar beide als Argumente dienen sollen. Moglicherweise

ist a ein Parameter, der einen festen Wert hat, und b das Argument oder

umgekehrt. Daher wird das Argument in der Vorschrift mit angegeben:

y(b) = a · b2

103

A.1. Polynome GMM WS2009/2010

Ist nun a = 5, so haben wie dieselbe Vorschrift wie in (A.3).

Nun gibt es einige Konventionen, die aber nicht strikt eingehalten werden:

So weiss man im Allgemeinen, dass bei y = mx + b eine Geradengleichung

gemeint ist, mit dem Argument x, der Steigung m und dem Abzissenwert b.

Meistens ist, wenn nicht anders angegeben x die unabhangige Variable.

In der Physik gibt es haufig Funktionen, in denen t als Argument auftritt. In

den meisten Fallen ist damit die Zeit gemeint.

Bildet man die Ableitung einer Funktion, so ist mit der Schreibweise y′ die

Ableitung nach dem (einzigen) Argument gemeint, dies ist meistens x. Ins-

besondere beschreibt f ′(x) die Ableitung der Funktion f nach x. Will man

dies besonders hervorheben so schreibt man ddx

f(x). In der Physik wird haufig

die Ableitung nach der Zeit, die sogenannte Zeitableitung, mit einem Punkt

dargestellt:

f(t) :=d

dtf(t) .

Ist der Definitionsbereich einer Abbildung N, so nennt man die Abbildung eine

Folge und schreibt

a : N −→ W

n 7→ an .

oder kurz (an)n∈N.

A.1 Polynome

Abbildungen der Form p(x) = am · xm + · · · + a1 · x + a0 , m ∈ N heissen

Polynome. Beruhmtestes Beispiel ist die Normalparabel p(x) = x2 mit m = 2,

a2 = 1 und a1 = a0 = 0.

A.2 Periodische Funktionen – Winkelfunktionen

Eine Funktion heisst periodisch mit der Periodenlange p , wenn fur alle x ∈ D

f(x + p) = f(x) gilt. Beispiele fur periodische Funktionen sind die Winkel-

funktionen (trigonometrische Funktionen) Sinus und Kosinus.

104

GMM WS2009/2010 Anhang A – Funktionen

1 2 3 4 5 6

-1

0

1

sin(x)cos(x)

α

sin(

α)

cos(α)

Abbildung A.1: Der Einheitskreis und die Winkelfunktionen

Betrachtet man im Einkeitskreis, das in Abbildung A.1 dargestellte Dreieck, so

ist die Lange der Hypotenuse 1. Der resultierende Abschnitt auf der x-Achse ist

also der Kosinus, des Winkels α, der Abschnitt auf der y-Achse ist der Sinus.

Die Steigung der Geraden betragt sin(α)cos(α)

= tan(α).

Das Argument der Winkelfunktionen kann sowohl in Grad als auch in Bo-

genmaß (Radiant) angegeben werden. Ein Vollkreis von 360 ◦ entspricht dabei

2π rad.

Wichtige Werte der Sinus und Kosinusfunktion:

Grad 0 30 45 60 90 180 270 360

Bogenmaß 0 π/6 π/4 π/3 π/2 π 3π/2 2π

Sinus 0 1/2√

2/2√

3/2 1 0 -1 0

Kosinus 1√

3/2√

2/2 1/2 0 -1 0 1

Eigenschaften der Sinus und Kosinusfunktion:

Periodenlange 2 π sin(x) = sin(x + 2π)

cos(x) = cos(x + 2π)

Sinus ist symmetrisch zum Ursprung sin(−x) = − sin(x)

Kosinus ist symmetrisch zur y-Achse cos(x) = cos(−x)

sind gegeneinander um π/2 verschoben cos(x) = sin(x + π/2)

Schnittpunkte bei x = π/4 + n · π, n ∈ Z sin(π/4) = cos(π/4)

Bemerkung: Berechnet man Sinus und Kosinus mit dem Taschenrechner so muss

105

Anhang B – Differential- und Integralrechnung GMM WS2009/2010

man aufpassen, dass die Einstellung stimmt. Ist der Taschenrechner auf Gradmaß

(DEG) eingestellt, muss man den Winkel in Grad eingeben, steht er auf Bogenmaß

(RAD) so muss er in Bogenmaß angegeben werden.

B Nicht ganz Neues aus der

Differential- und Integralrechnung

WARNUNG: Dieser Anhang ersetzt nicht die Mathematikvorlesungen. Aus

Grunden der Ubersichtlichkeit wird auf Vollstandigkeit, insbesondere auf die

Uberprufung und Angabe von Voraussetzungen verzichtet!

B.1 Kettenregel

Aus der Schule ist bekannt

f(g(x))′ = f ′(g(x)) · g′(x)

Was bedeutet dies nun? f und g sind Funktionen, wobei f(g(x)) die Hinter-

einanderausfuhrung der Funktionen f und g bedeutet. Die Funktion f wird

hierbei auch aussere Funktion ung g innere Funktion genannt.

f ◦ g(x) = f(g(x))

Das heisst, dass die unabhangige Variable x zuerst mittels der Funktion g

abgebildet wird und das Ergbnis dann mittels f abgebildet wird.

Die Kettenregel gibt nun an, wie die Ableitung der verketteten Funktion f ◦ g

zu bilden ist:

df(g(x))

dx=

df(g(x))

dg(x)· dg(x)

dx

Der erste Term auf der rechten Seite heisst außere Ableitung, der zweite Term

auf der rechten Seite innere Ableitung.

106

GMM WS2009/2010 Anhang B – Differential- und Integralrechnung

Setzt man nun y=g(x), so sieht das Ganze schon einfacher aus:

df(g(x))

dx=

df(y)

dy· dg(x)

dxoder kurz

df

dx=

df

dy· dy

dx

In der Modellierung ist wird die unabhangige Variable haufig mit t fur die Zeit

und die innere Funktion mit x bezeichnet:

f ◦ x(t) = f(x(t))

Fur die Ableitung gilt dann

df

dt=

df

dx· dx

dt

B.2 Substitutionsregel


∫ b

a

f(ϕ(x)) · ϕ′(x) dx =

∫ ϕ(b)

ϕ(a)

f(y) dy

Was bedeutet dies nun? Die Substitutionsregel ist ein Hilfsmittel um kompli-

zierte Funktionen in einfachere zu uberfuhren, die sich dann leichter integrieren

lassen. Diese Uberfuhrung heisst auch Substitution. Im obigen Beispiel wird

ϕ(x) durch y substituiert, also y = ϕ(x).

Letzlich ist die Substitutionsregel nichts anderes als die Kettenregel:

Ist F eine Stammfunktion von f , so gilt

F (ϕ(b)) − F (ϕ(a)) =

∫ ϕ(b)

ϕ(a)

f(y) dy

Es gilt aber auch

F (ϕ(b)) − F (ϕ(a)) =

∫ b

a

[F (ϕ(x))]′ dx

Die Ableitung von F (ϕ(x))′ ist nach der Kettenregel mit y = ϕ(x) aber

dF (ϕ(x))

dx= f(ϕ(x)) · ϕ′(x)

Schreibt man nun ϕ′(x) = dϕ(x)dx

so gilt also

107

B.3. Partielle Integration GMM WS2009/2010

∫ b

a

f(ϕ(x)) · dϕ(x)

dxdx =

∫ ϕ(b)

ϕ(a)

f(y) dy

Und nun sieht man auch, was die Physiker meinen, wenn sie das dx einfach

kurzen.

Haufig findet man auch andere Notationen, z.B.

∫ t

0

f(N) N dt =

∫ N(t)

N(0)

f(N) dN

Streng genommen darf das so nicht, da Integrationsgrenze und unabhangige

Variable verschiedene Bezeichner brauchen (hier heissen beide t), aber in der

Praxis kommt es ofter mal vor.

B.3 Partielle Integration

Die partielle Integration ergibt sich direkt aus der Produktregel fur Ableitun-

gen.


[f(x) · g(x)]′ = f ′(x) · g(x) + f(x) · g′(x)

Integration auf beiden Seiten ergibt

∫

[f(x) · g(x)]′ dx =

∫

f ′(x) · g(x) dx +

∫

f(x) · g′(x) dx

damit gilt

[f(x) · g(x)] =

∫

f ′(x) · g(x) dx +

∫

f(x) · g′(x) dx

durch Umstellen ergibt sich die Regel zur partiellen Integration

∫

f(x) · g′(x) dx = [f(x) · g(x)] −∫

f ′(x) · g(x) dx

108

GMM WS2009/2010 Anhang C – Differentialgleichungen

C Differentialgleichungen

Ersetzt nicht die Vorlesung zu Differentialgleichungen !!!

In einer Differentialgleichung (DGL) wird eine (unbekannte) Funktion mit ih-

ren Ableitungen verknupft, z.B.

dy(x)

dx= y(x) (C.1)

Gewohnliche Differentialgleichungen 1. Ordnung

Eine Gleichung der Form

dy(x)

dx= f(x, y(x)) (C.2)

heisst gewohnliche Differentialgleichung 1. Ordnung.

1.Ordnung: es tauchen keine hoheren Ableitungen nach x auf

gewohnlich: es gibt nur die Ableitung nach der unabhangigen Variablen x

Modellierungstechnisch gesprochen:

Die DGL ist eine Vorschrift, die die Anderung einer Zustandsvariablen y bezuglich

x angibt.

Eine Funktion, die die DGL erfullt, heisst Losung der DGL.

Eine DGL zu losen bedeutet, Funktionen y zu finden, die der Vorschriftdy(x)

dx= f(x, y(x)) genugen.

Ist die unabhangige Variable die Zeit, schreibt man anstelle von x gerne t:

dy(t)

dt= f(t, y(t)) oder kurzer

dy

dt= f(t, y) oder auch y = f(t, y)

(C.3)

Eine DGL heisst autonom, wenn f nicht explizit von der unabhangigen Variable

x (bzw. t) abhangt:

dy

dx= f(y) (C.4)

109

Anhang C – Differentialgleichungen GMM WS2009/2010

z.B.

dy

dx= α · y (C.5)

sonst heisst sie nicht-autonom, z.B.

dy

dx= −x · y (C.6)

Was kann man einer DGL ansehen?

Beispiel: y = α · y Die Steigung (Ableitung) der Funktion y an der Stelle t ist

proportional zum Funktionswert an der Stelle t.

Um dies zu visualisieren betrachtet man ein t, y-Koordinatensystem. Zu belie-

big vielen Paaren (t, y)-Paaren zeichnet man nun die Steigung an diesem Punkt

ein, indem man an (t, y) ein Linienelement mit der Steigung y′ = f(t, y) ein-

zeichnet (siehe Abbildung C.1). Ohne die Losungen der DGL zu kennen, weiss

man, dass eine Losung, die durch einen Punkt (t, y) geht, das Linienlement

als Tagente haben muss. Zeichnet man viele Linielemente, so kann man am

Richtungsfeld den Verlauf der Losungen”sehen“.

0

0.5

1

1.5

2

2.5

3

y(t)

0.5 1 1.5 2 2.5 3

t

Abbildung C.1: Richtungsfeld von y = y . Das Richtungsfeld ist t-translationsinvariant, da die DGL autonom ist.

Frage: Wie sieht eine Funktion aus, die in dieses Richtungsfeld passt?

Beobachtung: Es gibt viele Funktionen, die”passen“.

Die Losungen stimmen entweder uberein oder beruhren sich nicht.

Zu jedem Anfangswert y(t0) = y0 gibt es genau eine Losung.

Frage: Ist das bei allen DGLn so ?

Nein, aber unter gewissen Voraussetzungen, ja! Weiter hilft hier der Satz von

110


Picard-Lindelof (siehe Heuser: Differentialgleichungen ): f muss Lipschitzbe-

dingung fur y erfullen

Anfangswertproblem (AWP):

Gesucht ist eine Funktion y, die dy(t)dt

= f(t, y(t)), y(t0) = y0 genugt.

Das AWP zu losen bedeutet, eine Losung der DGL zu finden, die von (t0, y0)

ausgehend in ihrem gesamten Verlauf auf das Richtungsfeld”passt“.

Analytische Losung des AWP y = α · y , y(t0) = y0

Separation der Variablen:

y

y= α

Integration beider Seiten:

∫ t

t0

˙y(s)

yds =

∫ t

t0

α ds

Substitution auf der linken Seite x = y(s):

∫ y(t)

y(t0)

1

xds = α(t − t0)

[ln |x|]y(t)y(t0) = α(t − t0)

ln |y(t)| − ln |y(t0)| = α(t − t0)

ln | y(t)

y(t0)| = α(t − t0)

| y(t)

y(t0)| = eα(t−t0) > 0 ,

also gilt | y(t)y(t0)

e−α(t−t0)| = 1.

Da der Betrag eine stetige Abbildung ist und der Ausdruck in den Betragsstri-

chen nur -1 oder 1 sein kann, muss fur y(t0) > 0 auch y(t) > 0 f.a t gelten.

Damit ist

y(t) = y(t0)eα(t−t0)

die Losung des Anfansgwertproblems.

111

Anhang C – Differentialgleichungen GMM WS2009/2010

Die allgemeine Losung der DGL lautet

y(t) = C eα(t−t0) ,

wobei C eine freie Konstante ist. C ergibt sich fur eine spezielle Losung aus

dem Anfansgwert. Man kann zeigen, dass alle Losungen der DGL diese Form

besitzen.

Eine Losung einer DGL mit einer freien Konstanten heisst allgemeine Losung,

wenn sich jede Losung der DGL in dieser Form schreiben lasst.

Beispiel fur eine DGL, deren Anfangswertproblem keine eindeutige Losung hat

y = 3 3√

y 2 = 3 y2/3 , y(0) = 0

y1(t) = 0 ist Losung

y2(t) = t3 ist Losung, denn y2(t) = 3t2 und 3 ·(

3√

t3)2

= 3t2

Lineare Differentialgleichungen

y + a(t)y = 0 heisst homogene, lineare DGL (in y)

y + a(t)y = b(t) heisst inhomogene, lineare DGL (in y)

Bemerkung: Die inhomogene Gleichung ist streng genommen affin.

Den Term b(t) nennt man auch Storung.

Beispiele:

y = αy ist linear

y = αy2 ist nicht linear

Grundprinzipien bei linearen Differentialgleichungen

1. Uberlagerungsprinzip (Superpositionsprinzip) bei homogenen linearen DGL:

Sind y1, y2 Losungen der homogenen DGL, so ist auch

y3 = αy1 + βy2, α, β ∈ R Losung.

Die Losungsmenge ist ein Vektorraum !!

112


2. Ist y Losung der inhomogenen DGL und z eine beliebige Losung der zu-

gehorigen homogenen DGL, dann ist y+z Losung der inhomogenen DGL.

3. Sind y und w Losungen der inhomogenen DGL, dann ist y-w Losung der

zugehorigen homogenen DGL.

C.1 Numerisches Losen von Differentialgleichungen

Die wenigsten Differentialgleichung, die man in der Modellierung aufstellt, las-

sen sich analytisch losen. Man lost die Gleichungen daher numerisch. Dies

bedeutet, das man die Losung Stuck fur Stuck bestimmt. Gegeben sei das

Anfangswertproblem:

y′ = f(t, y(t)), y(t0) = y0

Euler Verfahren

Zuerst wird an der Stelle (t0, y0) die Steigung bestimmt. Es wird nun ange-

nommen, dass sich die Steigung in einem kleinen Schritt h (Schrittweite) nicht

andert. Die Losung wird also in dem Intervall [t, t + h] durch eine Gerade ap-

proximiert. Man erhalt den approximativen Wert y1 zur Zeit t1. Nun wird die

Steigung an dieser Stelle bestimmt und so weiter:

y(t0 + h) = y(t) + f(t0, y(t0)) · h

allgemein wird also zu einer beliebigen Zeit t der neue Wert zur Zeit t + h

durch

y(t + h) = y(t) + f(t, y(t)) · h

bestimmt. Ausgedruckt in Zeitschritten bedeutet dies:

yk+1 = yk + f(tk, yk) · h k = 0, 1, ..

113

C.1. Numerisches Losen von Differentialgleichungen GMM WS2009/2010

Dieses Verfahren ist recht ungenau.

Beispiel C.1.1 Gegeben ist das Anfangswertproblem y′ = y y(0) = 1, dessen

explizite Losung y(t) = et bekannt ist. Bestimmt man die Losung zur Zeit

T=10 mit dem Euler-Verfahren und einer Schrittweite von h = 1, so erhalt

man als approximative Losung 1024. Der exakte Wert betragt (auf 2 Nach-

kommastellen gerundet) 22026,47.

Fehlerabschatzung

Wie groß ist der Fehler, den man beim numerischen Losen macht? Hierzu

betrachtet man das Taylorpolynom von y:

y(t + h) =n∑

k=0

1

k!y(k)(t) · hk +

1

(n + 1)!· y(n+1)(ξ) · hn+1 ξ ∈ [t, t + h]

= pn(t) + Rn(t)

Wird anstelle der exakten Losung y(t) das Taylorpolynom n-ten Grades be-

trachtet, so ist der Abbruchfehler in jedem Schritt proportional zu hn+1.

Im Falle des Eulerverfahrens wird das Taylorpolynom 1. Grades betrachtet:

y(t + h) = y(t) + y′(t) · h

Der Abbruchfehler ist also in jedem Schritt proportional zu h2.

Runge-Kutta-Verfahren 2. Ordnung (Verbessertes Euler-Verfahren ,

Heun-Verfahren)

Statt der Steigung an der Stelle t als Geradensteigung zu verwenden, werden

die Steigungen in (t, y) und (t+h, y +f(t, y) ·h) gemittelt. Dies bedeutet, dass

man mit dem Funktionswert an der Stelle t+h, der sich nach einem Euler-

Schritt ergibt, die Steigung an t+h noch einmal neu berechnet. Das Mittel

der beiden Steigungen wird dann als endgultige Steigung fur den Zeitschritt

angenommen:

y(t + h) = y(t) +f(t, y(t)) + f(t + h,

Funktionswert nach

einem Euler-Schritt︷︸︸︷

y(t) + f(t, y(t)) · h)

2· h

114

GMM WS2009/2010 Anhang D – Zweidimensionale lineare Abbildungen

Ausgedruckt in Zeitschritten

yk+1 = yk +f(tk, yk) + f(tk+1, yk + f(tk, yk) · h)

2· h k = 0, 1, ...

Runge-Kutta-Verfahren 4. Ordnung

Das bekannteste Verfahren approximiert die Losung bis zur 4. Ordnung. Hier-

bei wird an einer weiteren Stutzstelle bei t + h2

die Steigung berechnet. auch

dies geschieht wieder mit dem verbesserten Verfahren. Die endgultige Steigung

fur den Zeitschritt wird aus 4 Einzelsteigungen gewichtet gemittelt, wobei die

Steigungen in der Mitte bei t+ h2

ein hoheres Gewicht bekommen. Es wird wie

folgt berechnet:

y(t + h) = y(t) +s1 + 2s2 + 2s3 + s4

6· h

s1 = f(t, y(t))

s2 = f(t + h2, y(t) + s1 · h

2)

s3 = f(t + h2, y(t) + s2 · h

2)

s4 = f(t + h, y(t) + s3 · h)

Ausgedruckt in Zeitschritten

yk+1 = yk +s1 + 2s2 + 2s3 + s4

6· h k = 0, 1, ...

s1 = f(tk, yk)

s2 = f(tk + h2, yk + s1 · h

2)

s3 = f(tk + h2, yk + s2 · h

2)

s4 = f(tk + h, yk + s3 · h)

Dieses Verfahren ist in vielen Programmen implementiert.

C.1.1 Verbesserungen

Bessere Verfahren sind implizite Verfahren, Verfahren mit einer automatischen

Zeitschrittanpassung resp. Prediktor-Korrektor-Verfahren.

115

Anhang D – Zweidimensionale lineare Abbildungen GMM WS2009/2010

D Zweidimensionale lineare Abbildungen

Betrachten wir zunachst die Abbildung

(u

v

)

=

(

a b

c d

)

·(

x

y

)

= A ·(

x

y

)

Die quadratische Matrix A bildet den R2 in sich selbst ab.

Eine solche Abbildung setzt sich aus Drehungen, Scherungen, Streckungen,

Projektionen und Spiegelungen zusammen:

Drehung

(u

v

)

=

cos ϕ − sin ϕ

sin ϕ cos ϕ

·(

x

y

)

Scherung

(u

v

)

=

1 b

0 1

·(

x

y

)

Streckung

(u

v

)

=

a 0

0 d

·(

x

y

)

Projektion

(u

v

)

=

1 0

0 0

·(

x

y

)

Spiegelung

(u

v

)

=

−1 0

0 1

·(

x

y

)

Die Determinate det A = ad − bc der Matrix A gibt an, um welchen Faktor die

Flache des Einheitsquadrats bei Abbildung vergroßert bzw. verkleinert wird.

Enthalt die Abbildung eine Spiegelung, wechselt das Vorzeichen. Enthalt die

Matrix eine Projektion, verschwindet die Determinante.

Betrachtet man die Bilder verschiedener Vektoren unter einer Matrix, so stellt

man fest, dass sie i.a. um verschiedene Winkel gedreht und um verschiedene

Faktoren gestreckt werden. Die Frage ist, inwieweit es Vektoren gibt, die nur

gestreckt, aber nicht gedreht werden. Wird ein Vektor nicht gedreht, so wird

jeder Vektor, der in dieselbe Richtung (oder in die entgegengesetzte Richtung)

zeigt, ebenfalls nicht gedreht. Diese Vektoren laufen entlang der sogenannten

Eigenrichtung der Matrix. Jeder Vektor, der dies erfullt (ausser der Nullvek-

tor) heisst Eigenvektor.

116

GMM WS2009/2010 Anhang D – Zweidimensionale lineare Abbildungen

Fur jeden Vektor −→v in Eigenrichtung kann die Matrix A also durch einen

Faktor λ ersetzt werden, der den jeweiligen Streckungsfaktor angibt. Dieser

Faktor heisst Eigenwert zu dieser Eigenrichtung.

Fur einen Eigenvektor −→v zum Eigenwert λ gilt

A · −→v = λ · −→v

oder anders mit der Einheitsmatrix I geschrieben:

(A − λI) · −→v =−→0

Diese Gleichung ist auch der Schlussel, um die Eigenwerte zu finden, sie heisst

Eigenwertgleichung.

Aus der Algebra ist bekannt: Ein lineares, homogenes Gleichungssystem mit

einer quadratischen Matrix B hat genau dann eine nicht-triviale Losung, wenn

det B=0 gilt, also hat die Eigenwertgleichung nur dann nicht-triviale Losungen,

wenn

p(λ) := det(A − λI) = 0

gilt. p(λ) heisst auch charakteristisches Polynom der Matrix A. Man findet

die Eigenwerte, indem man die Nullstellen des charakteristischen Polynoms

bestimmt. Im zweidimensionalen Fall gibt es hochstens zwei verschiedene Ei-

genwerte und man kann sie einfach bestimmen:

p(λ) = det

(

a − λ b

c d − λ

)

= (a−λ)(d−λ)−bc = λ2−(a+d)λ+ad−bc = 0

Mit det A = ad − bc und spur A = a + d erhalt man die Losungen

λ1,2 =spur A

2±√

(spur A)2

4− det A

Die Eigenvektoren erhalt man, indem man fur ein fest vorgegebenen Eigenwert

die Eigenwertgleichung lost:

(a − λ b

c d − λ

)

·(

v1

v2

)

= 0

Da die Eigenvektoren bis auf einen Faktor bestimmt werden konnen, kann

117

Anhang D – Zweidimensionale lineare Abbildungen GMM WS2009/2010

man eine Komponente festlegen, und die andere bestimmen. Fur b 6= 0 kann

man v1 = b setzen und erhalt den Eigenvektor(

bλ−a

). Fur c 6= 0 kann man

v2 = c setzen und erhalt den Eigenvektor(

λ−dc

). Fur b = c = 0 ist A eine

Streckmatrix mit den Eigenwerten a und d. Die zugehorigen Eigenvektoren

sind die Einheitsvektoren in x- bzw. y-Richtung. Gilt a = d, so ist A Vielfaches

der Einheitsmatrix und jeder Vektor (ausser dem Nullvektor) ist Eigenvektor.

118

GMM WS2009/2010 Anhang E – Stabilitatsanalyse autonomer Systeme (2D)

E Stabilitatsanalyse autonomer Systeme (2D)

E.1 System

X = f(X, Y ) (E.1)

Y = g(X, Y )

E.2 Stationare Zustande

Man erhalt die stationaren Zustande des ) Systems E.1, indem

f(x, y) = 0 (E.2)

g(x, y) = 0

gesetzt wird, und die Losungen (X∗, Y ∗) dieses Gleichungssystems bestimmt

werden.

E.3 Geschlossene Bahnen

Im Gegensatz zu eindimensionalen Systemen konnen im zweidimensionalen

Fall geschlossene Losungsbahnen auftreten. Im Normalfall ist die geschlossene

Bahn das Bild einer periodischen Losung. Man nennt eine solche periodische

Bahn stabil, wenn alle Bahnen, die in ihrer Nahe starten, zu dieser hin laufen.

Man nennt solche stabilen Bahnen auch Grenzkreise. Periodische Bahnen

konnen auch instabil sein.

E.4 Satz von Poincare

Bei autonomen Systemen mit zwei Zustandsvariablen lauft jede Losung ent-

weder

• auf einen stationaren Zustand zu, oder

• nahert sich einer geschlossenen Kurve (Grenzkreis) an, oder

• lauft nach unendlich

119

E.5. Teilgleichgewichte GMM WS2009/2010

E.5 Teilgleichgewichte

Betrachten wird die Menge aller Punkte fur die X = 0 bzw. Y = 0 gilt:

Teilgleichgewicht von X = {X = 0} = {(X, Y )|f(X, Y ) = 0}


die Bewegung nach rechts bzw. links erfolgt. Analog gilt

Teilgleichgewicht von Y = {Y = 0} = {(X, Y )| g(X, Y ) = 0}


die Bewegung nach oben bzw. unten erfolgt.

Die Schnittpunkte der Teilgleichgewichte sind gerade die stationaren Zustande.

E.5.1 Stabilitat stationarer Zustande

Die Stabilitat eines stationaren Zustands kann wie im eindimensionalen Fall

untersucht werden. Man betrachtet hierzu das System als Vektor

(X

Y

)

=

(f(X, Y )

g(X, Y )

)

(E.3)

Analog zum eindimensionalen Fall entwickelt man die rechte Seite des Differen-

tialgleichungssystems am stationaren Zustand (X∗, Y ∗) nach Taylor bis zum

ersten Glied und erhalt so das linearisierte System (ist fur jeden stationaren

Zustand unterschiedlich!!):

(x

y

)

= J(X∗, Y ∗) ·(

x

y

)

(E.4)

Die Matrix J ist die Jacobi-Matrix:

J(X∗, Y ∗) =

∂f(X∗, Y ∗)∂X

∂f(X∗, Y ∗)∂Y

∂g(X∗, Y ∗)∂X

∂g(X∗, Y ∗)∂Y

(E.5)

Das linearisierte System ist am stationaren Zustand (X∗, Y ∗) entwickelt und

gibt die Abweichung (Storung) von diesem an. Zu linearen Abbildungen, Ei-

genwerten etc. siehe Anhang D.

120


Hat die Matrix J zwei verschiedene Eigenwerte λ1 und λ2, so kann man eine

allgemeine Losung des linearen Systems E.4 angeben:

~z(t) = C1 · eλ1t · ~v1 + C2 · eλ2t · ~v2

Hierbei ist ~v1 Eigenvektor zum Eigenwert λ1 und ~v2 Eigenvektor zum Eigenwert

λ2, C1 und C2 sind freie Konstanten.

Sind beide Eigenwerte reell so handelt es sich bei den Summanden entweder

um exponentielles Wachstum oder exponentiellen Zerfall. Ist also einer der

Eigenwerte großer Null, wachst die Storung, der Zustand ist instabil. Sind beide

Eigenwerte negativ, so verschwindet eine anfangliche Storung, der Zustand ist

stabil.

Fur einen komplexen Eigenwert λ = a + ib gilt

eλt = e(a+ib)t = eat · eibt = eat · (cos bt + i sin bt)

Das System vollzieht also Schwingungen. Ist der Realteil a kleiner Null, so ver-

schwindet die Storung, ist der Realteil a großer Null, so explodieren Storungen.

Bei a = 0 bleiben Storungen erhalten. Uber die Stabilitat des Ursprungssy-

stems kann dann nicht entschieden werden.

Bezuglich der Stabilitat des Systems E.1 gilt der folgende Satz:

Satz E.5.1 Stabilitat von stationaren Zustanden autonomer Systeme

1. Der stationare Zustand (X∗, Y ∗) von E.1 ist asymptotisch stabil, falls

alle Eigenwerte von negative Realteile besitzen.

2. Der stationare Zustand (X∗, Y ∗) ist instabil, falls mindestens ein Ei-

genwert einen positiven Realteil besitzt.

3. Sind alle Realteile aller Eigenwerte kleiner oder gleich Null und min-

destens ein Realteil gleich Null, so kann nicht uber die Stabilitat des

stationaren Zustands (X∗, Y ∗) entschieden werden.

Der Satz gilt auch fur hoherdimensionale autonome Systeme!

121

E.6. Veranschaulichung des Systemverhaltens GMM WS2009/2010

E.6 Veranschaulichung des Systemverhaltens

Sind beide Eigenwerte reell und verschieden, so gibt es zwei Eigenrichtungen.

Diese Richtungen sind selbst Losungen und konnen von anderen Losungen

nicht uberquert werden. Entlang der Eigenrichtungen bewegen sich die Losun-

gen in gerader Linie auf den stationaren Zustand hin (λ < 0) oder weg (λ > 0).

Die absolute Große des Eigenwerts bestimmt dabei die Geschwindigkeit. Fur

Losungen, die nicht auf einer Eigenrichtung starten gibt es Einflugschneisen.

Fur z.B. λ1 = −0.1 und λ2 = −1 bewegt sich das System zuerst schnell in

Richtung von v2, und dann langsam in Richtung von v1.

Ist ein Eigenwert 0 so findet auf der zugehorigen Eigenrichtung keine Bewegung

statt.

Sind die Eigenwerte komplex (genauer: nicht reell) gibt es keine reellen Eigen-

richtungen. Das System bewegt sich an keiner Stelle radial auf den stationaren

Zustand zu.

Fur die Eigenwerte des Systems gilt

λ1,2 =spur A

2±√

(spur A)2

4− det A =:

µ

2±√

µ2

4− D

Eigenwerte µ D Name Stabilitat Beispiel

λ2 < λ1 < 0 µ < 0 0 < D < µ2

4Knoten stabil

(

−1 1

0.5 −1

)

λ2 < 0 < λ1 D < 0 Sattel instabil

(

0.5 −1

−1 0.5

)

0 < λ2 < λ1 µ > 0 0 < D < µ2

4Knoten instabil

(

0.5 0.5

−1 2

)

λ1,2 ∈ C, Re(λ) < 0 µ < 0 D > µ2

4Focus stabil

(

−1 −1

1 −1

)

λ1,2 ∈ C, Re(λ) = 0 µ = 0 D > µ2

4Zentrum marginal stabil

(

0 −1

1 0

)

λ1,2 ∈ C, Re(λ) > 0 µ > 0 D > µ2

4Focus instabil

(

1 −1

1 1

)

Bei zusammenfallenden Eigenwerten und bei Eigenwerten, die Null werden,

gibt es weitere Formen des Systemverhaltens, auf die hier nicht naher einge-

gangen wird.

122


0spur

det

reelle Eigenwerte reelle Eigenwerte

(spur)2/4=det

komplexe Eigenwerte

stabiler Focus instabiler Focus

Sattel

Zentrumstabiler Knoten instabiler Knoten

Sattel

Abbildung E.1: Stabilitat stationarer Zustande im zweidimensionalen System.Dargestellt sind nur die wesentlichen Falle.

123

Grundlagen der mathematischen Modellierung · Grundlagen der Mathematischen Modellierung Skript zur...

Documents

Transcript of Grundlagen der mathematischen Modellierung · Grundlagen der Mathematischen Modellierung Skript zur...