Quantitative Systemwissenschaften - Karl-Franzens

Quantitative Systemwissenschaften
G.W. Desch und G. Propst
Institut für Mathematik und Wissenschaftliches Rechnen
Karl-Franzens-Universität Graz
Quantitative Systemwissenschaften 1:
1. Grundsätzliches zu mathematischen Modellen
2. Empirische Modelle und Regressionsgerade
3. Parameteranpassung durch kleinste Quadrate
4. Mengenbilanzen
5. Dimensionen und Einheiten
6. Simulation dynamischer Systeme
7. Fragestellungen zu dynamischen Systemen
8. Simulationsdiagramme
Quantitative Systemwissenschaften 2:
9. Qualitative Analyse mathematischer Modelle
10. Verteilte Parameter
11. Diskrete deterministische Modelle
12. Diskrete probabilistische Modelle
1
2
1. Grundsätzliches zu mathematischen Modellen
Viele Fragestellungen der Umweltwissenschaften werden mit Hilfe quantitativer
Methoden bearbeitet. Gemeint ist hier insbesondere die Arbeit mit mathematischen
Modellen von realen Systemen. Im Rahmen dieser Vorlesung verstehen wir unter einem
“System” einen Teilbereich der Wirklichkeit, für den ein mathematisches Modell erstellt
wird. Das ist eine Sammlung gekoppelter Gleichungen und Formeln, deren Lösungen
einige der Eigenschaften des Systems wiedergeben oder sein Verhalten beschreiben. Das
Modell wird mit Hilfe mathematischer und computergestützter Methoden untersucht.
Sofern das Modell die Eigenschaften des Sytems gültig wiedergibt, ist die mathematische
Untersuchung des Modells eine “Analyse” des Systems selbst. Daher gehört zur
Systemanalyse in dieser Bedeutung insbesondere auch die Erstellung eines
mathematischen Modells, die verständnisvolle Arbeit damit, seine Verbesserung,
Anpassung oder - gegebenenfalls - Verwerfung. Ein genauerer Titel der vorliegenden
Lehrveranstaltung wäre
Modellbildung und Simulation.
Populationen, mechanische Systeme, chemische Reaktoren, Biotope, Handelsbetriebe,
Volkswirtschaften, Körperorgane, Klimasysteme: dies sind einige Beispiele für
Teilbereiche der Realität, für deren Modellierung Methoden, die in dieser Vorlesung
besprochen werden, angewendet werden können.
Ein wesentlicher Aspekt mathematischer Modelle ist die vollständige Definiertheit ihrer
Bestandteile. Dies macht sie erstens einer mathematischen Analyse zugänglich und kann
zweitens zur quantitativen Berechnung konkreter Beispiele verwendet werden. Folgende
Ziele können Zweck der mathematischen Modellbildung sein:
• Modellbildungsprozess selbst: Aufdeckung von Wirkungszusammenhängen
aufgrund systematischer Vorgangsweise, Erweiterung des Verständnisses des
modellierten Systems
• quantitative numerische Simulation und Prognose möglicher Entwicklungen
• qualitative mathematische Analyse des Modells
• Untersuchung von Auswirkungen der Änderung von Modell-Parametern
• rechnerunterstützter Entwurf von Systemen (spart Zeit, vermeidet Risiko)
• Identifikation von nicht oder nur sehr schwer messbaren Größen (z. B. Seismik
(Erdöl), Computertomographie (Gewebe))
• Verbesserung oder Optimierung von Systemeigenschaften.
Typische methodische Elemente der Bildung eines mathematischen Modells sind
• Vereinfachung, Idealisierung
• Festlegung des Wesentlichen, Weglassung von Unwesentlichem
• schrittweises Vorgehen von einfach zu komplex
• Versuch und Irrtum.
In den Modellen können alle möglichen mathematischen Objekte verwendet werden,
jedoch beeinflusst der jeweilige Zweck des Modells die Auswahl unter möglichen
1. GRUNDSÄTZLICHES ZU MATHEMATISCHEN MODELLEN
3
Modell-Typen. Grundsätzlich braucht ein Modell nicht komplizierter sein als für seinen
Zweck erforderlich, selbst wenn das betrachtete System sehr komplex ist.
Als Beispiel für ein mathematisches Modell mit dem Zweck einer Prognose betrachten
wir die Verunreinigung eines Sees durch einen Schadstoff. Die Ursache der
Verunreinigung spielt keine Rolle, jedenfalls wird nach dem Zeitpunkt t = 0 kein
Schadstoff mehr zugeführt. Die zu beantwortende Frage lautet: Wie lange wird es dauern,
bis – allein auf Grund der Strömung – die Schadstoffkonzentration im See auf 10% der
Anfangskonzentration reduziert ist? Der See hat das konstante Volumen V , pro Tag
strömen r m3 Wasser zu und ab. Zur Vereinfachung nehmen wir an, dass der Schadstoff
immer gleichmäßig im See verteilt ist. Dann ändert sich die Schadstoffkonzentration S(t)
gemäß der Mengenbilanz (dies wird später mehrmals genauer erläutert)
Änderung der Menge Schadstoff im See = Schadstoff Zufluss - Schadstoff Abfluss,
d
(V S(t)) = 0 − rS(t).
dt
Die Lösung dieser Differentialgleichung (darüber wird später mehr gesagt) mit der
Anfangsbedingung S(0) = S0 ist
S(t) = S0 e−rt/V .
Der gesuchte Zeitpunkt T ist jener, für den S(T ) = S0 /10, also T = ln(10) · V /r. Für den
Lake Michigan ist V = 4.871 × 1012 m3 und r = 4.331 × 108 m3 /Tag, also T=25897 Tage
= 71 Jahre. Natürlich enthält das Modell Vereinfachungen, z.B. die Annahme der totalen
Vermischung, die aber den Abfluss an Schadstoff eher beschleunigen würde; andererseits
sind Schadstoff vermindernde Prozesse wie chemische Reaktionen oder Sedimentation
ausser Acht gelassen.
Der springende Punkt der Brauchbarkeit eines mathematischen Modells für den
angepeilten Zweck ist die Güte des Modells. Das System wird beobachtet, es gibt
qualitative und ev. quantitative Daten. Die Güte des mathematischen Modells wird an
Hand des Grads der Übereinstimmung der “Vorhersagen” des Modells mit den Daten
bewertet. Diese Übereinstimmung wird mitunter quantifiziert (z.B. Summe der FehlerQuadrate) und durch geignete Wahl der Modellparameter optimiert. Wenn das Modell
allerdings strukturell falsch ist, dann gibt es keine geeigente Wahl von Parametern und
die Modellbildung ist (vorläufig) gescheitert.
Im Laufe der Vorlesung werden die dargestellten Methoden der mathematischen
Modellbildung und Analyse an Hand von Beispielen erläutert. Dabei kommt es nicht
darauf an, dass die Modelle möglichst ausgefeilt sind, vielmehr werden allgemeine
Modellierungsprinzipien an Hand einfacher Beispiele besprochen. Ferner wird
demonstriert, wie man mathematische Modelle am Computer verwendet und untersucht.
Es ist ja die Verfügbarkeit von Rechenleistung, bedienungsfreundlicher Software und
graphischer Darstellung ein Hauptfaktor der zunehmenden Anwendung von
Systemanalyse, die auf mathematischer Modellbildung beruht.
4
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
2. Empirische Modelle und Regressionsgerade
2.1. Graphische Darstellung der Daten.
Beispiel 2.1. Die folgende Tabelle zeigt die Bevölkerung der Vereinigten Staaten im
Zeitraum von 1790 bis 1950 (in Einheiten zu tausend Einwohnern).
1790
3929
1800
5308
1810
7240
1820
9638
1830 12866
1840 17069
1850 23192
1860 31443
1870 38558
1880 50156
1890 62948
1900 75995
1910 91972
1920 105711
1930 122775
1940 131669
1950 150697
Gesucht ist ein einfaches Modell, das erlaubt, aus diesen Daten das weitere Wachstum
der Population vorauszusagen.
Für die folgenden Untersuchungen benennen wir mit t die Zeit, in Jahren, und mit P (t)
die Population zum Zeitpunkt t, in tausend Einwohnern.
Abbildung 2.1 zeigt die Daten graphisch. Wir haben die Datenpunkte hier mit
Geradenstücken verbunden, um die Graphik deutlicher zu machen.
Wesentlich deutlicher als aus der Tabelle werden die Wachstumstrends sichtbar. Die
Steigung der Geradenstücke repräsentiert die Bevölkerungszunahme in den einzelnen
Jahrzehnten. Wir sehen, daß die Bevölkerung im gesamten betrachteten Zeitraum
zunimmt, und daß auch die Wachstumsrate zunächst zunimmt. Etwa ab 1900 wird aber
der Zuwachs wieder flacher.
Als eine mögliche Erklärung bietet sich an: Im ungeheuer großen Lebensraum der
Vereinigten Staaten entwickelt sich die Population zunächst ungehemmt, je größer die
Population, desto mehr Nachkommen. Mit immer dichterer Bevölkerung werden aber
Mechanismen wirksam, die das Wachstum einbremsen, geringere Kinderfreudigkeit
und/oder erhöhte Sterblichkeit bewirken. Ein gutes Modell müßte diese Effekte
wiedergeben können.
Wenn wir das Phänomen nicht grob vereinfachend als Beispiel, sondern ernsthaft
analysieren wollten, sollten wir jetzt Bevölkerungsstatistiken heranziehen und
herausfinden, ob der Wachstumsrückgang auf eine höhere Sterblichkeit oder eine
5
6
Abbildung 2.1. Population der USA
x10 4
Population der USA
16
o
14
o
o
12
o
Population
10
o
8
o
o
6
o
4
o
o
o
2
0
1780
o
o
o
1800
o
o
1820
o
1840
1860
1880
1900
1920
1940
1960
Jahr
geringere Geburtenrate zurückzuführen ist. Wir wollen für dieses Beispiel annehmen, daß
uns zur Modellbildung nicht mehr als die obigen Daten und unser Hausverstand zur
Verfügung steht.
Wir werden zunächst nicht auf die Vorgänge und Zusammenhänge innerhalb der
Population, also Geburts- und Sterbeprozesse, eingehen, sondern nur versuchen, welche
Formeln möglichst gut zu den gemessenen Daten passen.
Merksatz 2.2. Ein empirisches mathematisches Modell ist eine Gleichung, die das
Verhalten eines Systems mit möglichst wenig Aufwand möglichst treffend wiedergibt. Zur
Ableitung eines empirischen Modells werden nicht Naturgesetze und die Funktionsweise
des Systems herangezogen, sondern nur die Form des Datenmaterials und eine Funktion,
von der man weiß, daß sie ähnlich verläuft.
Der Vorteil eines empirischen Modells liegt darin, daß es einfach in der Handhabung und
mit wenig Aufwand zu erstellen ist. Der Nachteil liegt darin, daß es zwar die
Systemeigenschaften beschreibt, aber keine Erklärung liefert, welche Ursachen diese
Eigenschaften bewirken.
2.2. Regressionsgerade.
Das einfachste empirische Modell ist eine Geradengleichung. Man spricht von der
Regressionsgeraden an einen Datensatz.
P (t) ≈ kt + d.
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
7
Abbildung 2.2. Regressionsgerade
x10 4
Population USA - lineare Regression
16
o
14
o
o
12
o
10
Population
o
8
o
o
6
o
4
o
o
o
2
0
o
-2
1780
o
1800
o
o
1820
o
o
1840
1860
1880
1900
1920
1940
1960
Jahr
Die Parameter k und d werden so bestimmt, daß die Gerade möglichst genau den
Datensatz wiedergibt. Abbildung 2.2 zeigt den Datensatz mit eingezeichneter
Regressionsgerade:
Wir sind mit diesem Modell nicht zufrieden, sondern beobachten die folgenden Fehler:
• Die Datenpunkte liegen nicht genau auf der Geraden. Aber das sollte uns nicht
stören, denn perfekte Daten, und schon gar perfekte Modelle gibt es nicht. Jedes
Modell ist eine Vereinfachung und kann den Datensatz nicht genau wiedergeben.
• Die Datenpunkte liegen nicht “wie zufällig” um die Gerade gestreut, sondern der
Datensatz zeigt deutliche Trends, die das Modell nicht wiedergibt, und die wir
oben beschrieben haben. Das Modell versagt in der Wiedergabe wesentlicher
Eigenschaften des Systems und wird daher verworfen.
An einem Beispiel mit weniger Rechenaufwand soll gezeigt werden, wie man eine
Regressionsgerade berechnet.
Beispiel 2.3. Bestimme die Regressionsgerade
y ≈ kx + d
zum folgenden Datensatz:
x 0 1 2 3 4
y 1 1 2 6 5
Lösung. Gegeben sind n = 5 Datenpaare (x1 , y1 ), · · · , (x5 , y5 ). Wir erstellen eine
Tabelle aus den Werten von xi , yi , xi yi und x2i und summieren die Spalten:
8
Abbildung 2.3. Regressionsgerade zu Beispiel 2.3
Regressionsgerade
6
o
5
o
y
4
3
2
o
1o
0
0
o
0.5
1
1.5
2
2.5
3
3.5
4
x
x y xy x2
0 1 0 0
1 1 1 1
2 2 4 4
3 6 18 9
4 5 20 16
10 15 43 30
n = 5 ist die Anzahl der Datenpaare. Wir bestimmen zunächst die Mittelwerte nach den
Formeln
n
1X
10
x=
xi =
= 2,
n i=1
5
n
1X
15
y=
yi =
= 3.
n i=1
5
Die Steigung der Geraden errechnet sich nach der Formel
Pn
1
xi yi − x · y
k = n1 Pi=1
n
2
2
i=1 xi − (x)
n
=
43
−3·2
5
30
− 22
5
= 1.3.
Der Parameter d wird dadurch bestimmt, daß die Gerade durch die Mittelwerte geht:
kx + d = y
1.3 · 2 + d = 3
d = 0.4.
Abbildung 2.3 zeigt den Datensatz und seine Regressionsgerade. Auf die graphische
Darstellung als Probe sollte man nicht verzichten.
¤
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
9
Merksatz 2.4. Gegeben sei ein Satz von n Datenpaaren (x1 , y1 ), · · · , (xn , yn ). Die
Regressionsgerade y = kx + d errechnet man nach folgenden Formeln:
n
1X
x=
xi ,
n i=1
n
1X
yi ,
n i=1
Pn
1
xi y i − x · y
k = n1 Pi=1
,
n
2
2
i=1 xi − (x)
n
d = y − kx.
y=
Wissenschaftliche Taschenrechner und Statistikprogramme auf Computern berechnen
Regressionsgeraden automatisch.
2.3. Exponentialfunktion und logarithmische Darstellung.
Wie schon bemerkt, eignet sich die Regressionsgerade nicht zur Beschreibung der
Population aus Beispiel 2.1. Kleine Populationen, deren Wachstum nicht durch die
Beschränkungen ihres Lebensraumes behindert werden, kann man durch
Exponentialfunktionen modellieren:
P (t) ≈ Cekt .
Dabei sind die Parameter C und k so zu wählen, daß der Datensatz möglichst gut
wiedergegeben wird. Ein Blick auf Abbildung 2.1 zeigt, daß die Exponentialgleichung
nicht den gesamten Datensatz gut wiedergeben wird, denn die Exponentialfunktion ist
konvex, d.h. ihre Steigung nimmt ständig zu, während die Steigung der beschriebenen
Population ab etwa 1900 abnimmt. Es wäre aber denkbar, daß sich die Population in den
Jahren davor, solange sie noch klein ist, durch eine Exponentialfunktion hinreichend gut
beschreiben läßt. Leider ist das Beurteilen einer Kurve, ob sie etwa eine Parabel, eine
Exponentialkurve, oder vielleicht eine Hyperbel ist, mit dem freien Auge so gut wie
unmöglich. Wir verwenden einen Trick, der die Exponentialfunktion auf eine Gerade
reduziert:
Wir gehen von der Gleichung einer Exponentialfunktion aus und nehmen (natürliche)
Logarithmen:
P (t) ≈ Cekt
ln(P (t)) ≈ ln(Cekt ) = ln(C) + ln(ekt )
= ln(C) + kt.
Setzen wir
y(t) = ln(P (t)), d = ln(C),
so erhalten wir eine Geradengleichung
y(t) ≈ kt + d.
10
Abbildung 2.4. Logarithmischer Plot
Population USA - log-Plot
14
13
12
log(P)
o
o
o
o
o
o
o
11
o
o
o
o
10
o
o
o
9
o
o
o
8
1780
1800
1820
1840
1860
1880
1900
1920
1940
1960
Jahr
Gegeben sind die Daten von t1 , t2 , · · · und P (t1 ), P (t2 ), · · ·. Wir berechnen die Werte
y1 = ln(P (t1 )), y2 = ln(P (t2 )), usw.
und tragen y gegen t auf. Wenn sich eine Gerade ergibt, genügt P einem
Exponentialgesetz (Abbildung 2.4).
t
P ln(P )
1790
3929 8.28
1800
5308 8.58
1810
7240 8.89
1820
9638 9.17
1830 12866 9.46
1840 17069 9.75
1850 23192 10.05
1860 31443 10.36
1870 38558 10.56
1880 50156 10.82
1890 62948 11.05
1900 75995 11.24
1910 91972 11.43
1920 105711 11.57
1930 122775 11.72
1940 131669 11.79
1950 150697 11.92
Die Graphik zeigt eine sehr gute Übereinstimmung mit einer Geraden für etwa die ersten
8 Datenpaare. Wir haben für diese Daten die Regressionsgerade berechnet und
eingezeichnet. Diese Rechnung führt auf
k = 0.0295, d = −44.55.
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
11
Abbildung 2.5. Population der USA - Exponentialapproximation
x10 5
Population USA - Fit durch Exponentialfunktion
3
2.5
Population
2
1.5
o
o
o
o
1
o
o
o
0.5
o
o
0
1780
o
o
1800
o
o
1820
o
o
1840
o
o
1860
1880
1900
1920
1940
1960
Jahr
Der Parameter C errechnet sich aus d:
d = ln(C),
ed = C,
C = e−44.55 = 4.48 · 10−20
Als empirisches Modell für die Jahre zwischen 1790 und 1860 erhalten wir
P (t) ≈ 4.48 · 10−20 · e0.0295t = e0.0295t−44.55
Zur Überprüfung tragen wir jetzt P (t) und die Approximation gegen t auf. Bei der
Bewertung der Güte der Näherung ist zu beachten, daß nur die ersten 8 Datenpaare
gefittet wurden. Die weitere Entwicklung wird durch die Exponentialfunktion völlig
falsch wiedergegeben (Abbildung 2.5). Man sieht daraus, daß auch ein Modell, das
gegenwärtige Zustände exzellent beschreibt, die Zukunft nicht immer richtig vorhersagt.
Merksatz 2.5. Gegeben sei ein Satz von Datenpaaren (x1 , y1 ), · · · , (xn , yn ). Einen
logarithmischen Plot von y über x erstellt man, indem man die Logarithmen der Daten
zi = ln(yi ) über xi aufträgt. Ergibt der logarithmische Plot eine Gerade z ≈ kx + d, so
genügen die Originaldaten einem Exponentialgesetz: y ≈ Cekx , dabei ist C = ed .
Logarithmische Plots kann man auch verwenden, wenn die Zahlenwerte eines
Datensatzes sich um viele Zehnerpotenzen unterscheiden und im linearen Plot die kleinen
Zahlen nur mehr auf der Nullinie “kleben”.
12
2.4. Andere Transformationen.
Auch andere Modelle als die Exponentialkurve lassen sich durch geschickte
Transformation der Variablen auf Geraden reduzieren.
Beispiel 2.6. In einer Lösung zerfällt ein chemischer Stoff X. Alle Viertelstunden wird
die Konzentration von X gemessen. Mit ti bezeichnen wir die Zeitpunkte (in Minuten),
mit ui die gemessenen Konzentrationen (in Mol pro Liter).
t
u
0 0.0190
15 0.0130
30 0.0091
45 0.0069
60 0.0057
75 0.0052
90 0.0043
105 0.0039
120 0.0035
Die Zerfallsreaktion heißt eine Reaktion m-ter Ordnung, wenn zum Eintreten der
Reaktion m Moleküle von X zusammentreffen müssen. Bei einer Reaktion erster
Ordnung genügt die Konzentration einem Exponentialgesetz
u(t) = Ce−kt .
Bei einer Reaktion höherer Ordnung entspricht die Konzentration dem Gesetz
1
u(t) = (kt + d)− m−1 .
Welche Ordnung hat die beschriebene Reaktion? Wie lauten die Parameter für das
Zerfallsgesetz?
(Die Daten dieses Beispieles sind fiktiv.)
Lösung. Wir beginnen mit einer graphischen Darstellung der Daten
(Abbildung 2.6). Man sieht deutlich das Abnehmen der Konzentration, aber es ist
schwer, zu entscheiden, ob es einer Exponentialkurve oder einer Hyperbel folgt.
Um das Exponentialgesetz zu überprüfen, versuchen wir einen logarithmischen Plot
(Abbildung 2.7):
t
u ln(u)
0 0.0190 -3.96
15 0.0130 -4.34
30 0.0091 -4.70
45 0.0069 -4.98
60 0.0057 -5.17
75 0.0052 -5.26
90 0.0043 -5.45
105 0.0039 -5.55
120 0.0035 -5.65
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
13
Abbildung 2.6. Daten aus Beispiel 2.6
Konzentration
0.02
o
0.018
0.016
0.014
Mol/Liter
o
0.012
0.01
o
0.008
o
0.006
o
o
o
0.004
0.002
0
20
40
60
80
o
100
o
120
Minuten
Abbildung 2.7. Log-Plot zu Beispiel 2.6
Konzentration - log-Plot
-3.5
log(Konzentration)
-4 o
o
-4.5
o
o
-5
o
o
o
-5.5
o
o
-6
0
20
40
60
80
100
120
Minuten
Wir haben zum Vergleich die Regressionsgerade eingezeichnet, aber man sieht mit freiem
Auge, daß die Daten durchaus einer gekrümmten Kurve und keiner Geraden folgen. Die
Reaktion ist bestimmt nicht erster Ordnung.
14
Abbildung 2.8. Plot der Reziprokwerte zu Beispiel 2.6
Konzentration - transformiert
300
o
o
250
1/Konzentration
o
200
o
o
150
o
o
100
o
50 o
0
20
40
60
80
100
120
Minuten
Um eine Reaktion zweiter Ordnung ins Auge zu fassen, verwenden wir eine andere
Transformation:
u(t) ≈
1
kt + d
1
≈ kt + d.
u(t)
1
kommen wir auf eine Geradengleichung. Wir tragen also den
Mit y(t) = u(t)
Reziprokwert von u gegen t auf (Abbildung 2.8):
1
t
u
u
0 0.0190 52.63
15 0.0130 76.92
30 0.0091 109.89
45 0.0069 144.93
60 0.0057 175.44
75 0.0052 192.31
90 0.0043 232.56
105 0.0039 256.41
120 0.0035 285.71
Die Übereinstimmung mit einer Geraden ist sehr gut. Wir nehmen daher an, daß eine
Reaktion zweiter Ordnung stattfindet. Die Formeln der Regressionsgeraden liefern die
Parameter für k und d:
k = 1.96,
d = 52.08,
u(t) ≈
1
.
1.96t + 52.08
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
15
Abbildung 2.9. Daten und Approximation als Reaktion zweiter Ordnung
Konzentration
0.02
o
0.018
0.016
0.014
Mol/Liter
o
0.012
0.01
o
0.008
o
0.006
o
o
o
0.004
0.002
0
20
40
60
80
o
100
o
120
Minuten
Als Probe tragen wir noch einmal die Originaldaten und die Approximation auf
(Abbildung 2.9)
¤
Wir geben noch einige Transformationen an, mit denen sich nichtlineare Modelle auf
Geraden reduzieren lassen:
Merksatz 2.7. Die folgenden Transformationen führen Funktionen y(x) in Geraden
über:
1
(s = x)
z = y1
z = ax + b
1) y(x) = ax+b
ax
1
1
2) y(x) = x+b
s= x
z=y
z = ab s + a1
y
ax
s= x
(z = y)
y = −bs + a
3) y(x) = x+b
4) y(x) = Ceλx (s = x)
z = ln(y)
z = λx + ln(C)
λx
5) y(x) = Ce
(s = x)
z = log10 (y) z = λ log10 (e)x + log10 (C)
6) y(x) = Cxα
s = ln(x)
z = ln(y)
z = αs + ln(C)
α
7) y(x) = Cx
s = log10 (x) z = log10 (y) z = αs + log10 (C)
8) y(x) = Cxeλx (s = x)
z = ln( xy )
z = λx + ln(C)
Manche dieser Transformationen haben Namen:
2) Lineweaver-Burk Transformation
3) Eadie-Hofstee Transformation
4,5) Logarithmische Transformation
6,7) Doppelt logarithmische Transformation
Merksatz 2.8. Verwendung von Merksatz 2.7: Gegeben sind Datenpaare
(x1 , y1 ), · · · , (xn , yn ) und ein empirisches Modell y ≈ f (x). Die Parameter des
empirischen Modells sind an die Daten anzupassen.
16
• Berechne aus den gegebenen Daten xi und yi die Daten si und zi .
• Trage z über s auf. Wenn der Plot annähernd eine Gerade ergibt, paßt das
empirische Modell zu den Daten.
• Bestimme die Parameter der Regressionsgeraden z ≈ ks + d.
• Berechne die Parameter des Modells aus den Parametern der
Regressionsgeraden.
• Graphische Probe mit xi , yi .
Zum Beispiel für die Lineweaver-Burk Transformation sieht der vorletzte Schritt so aus:
seien k und d die Parameter der s, z-Geraden, dh. k = ab , d = a1 . Daraus folgt
a = d1 , b = ak = kd , also
y(x) =
1
x
d
x+
k
d
=
x
.
dx + k
Wir haben hier nur sehr einfache empirische Modelle behandelt, insbesondere solche, die
nur zwei Größen in Beziehung setzen. Die mathematische Statistik bietet unter dem
Namen (lineare und nichtlineare) Regression Methoden zur Anpassung und Validierung
(dh. Beurteilung der Gültigkeit) viel komplexerer Modelle an umfangreichere Datensätze
an.
2.5. Rechnen und Erstellen von Graphiken.
Während der Vorlesung werden die besprochenen Rechnungen und graphischen
Darstellungen am Computer demonstriert. Es gibt viele Nummerik-Pakete, die für diese
Zwecke geeignet sind. Eines davon ist Matlab (professionell aber teuer), das auf vielen
Rechnern in den Benutzerzentren der KFU (bzw. am Terminalserver des ZID) installiert
ist.
Auf der Seite http://www.uni-graz.at/georg.propst/teaching propst.html sind
Verknüpfungen zu Einführungen in Matlab abrufbar. Ausserdem bietet Matlab selbst
ausführliche Hilfe und eine Einführung. In diesem Skriptum sind einige der verwendeten
Matlab Kommandos abgedruckt, um ein Gefühl für die Logik bei der Arbeit am
Computer zu vermittlen. Die Beherrschung der Syntax einer Programmiersprache ist
aber kein Lehrziel dieser Vorlesung.
Matlab ist eine interaktiv verwendbare Software, das heisst der Benützer gibt im
command window (Eingabeaufforderung >>) Befehle ein, die bei Betätigung der
Eingabetaste (←-) sofort ausgeführt werden. Mehrmals zu verwendende Befehle schreibt
man in sogenannte script files (mit der Erweiterung .m), deren Befehlszeilen bei Aufruf
des filenamens ausgeführt werden. Also schreibt man eine Datei, namens chemdata.m,
sagen wir, die die Zeilen
t=[0:15:120];
u=[0.019 .013 .0091 .0069 .0057 .0052 .0043 .0039 .0035];
lu=log(u);
2. EMPIRISCHE MODELLE UND REGRESSIONSGERADE
17
enthält, sodass im command window bei der Anweisung
>>chemdata ←den Vektoren t und p die Werte des obigen Beispiels zugewiesen werden. lu ist ein
beliebig gewählter Variablen-Name für einen Vektor, der die Ergebnisse von log(u)
enthalten soll. Die Matlab function log berechnet komponentenweise den logarithmus
naturalis (zur Basis e) des Vektors u. Ein Strichpunkt am Ende einer Befehlszeile
unterdrückt die manchmal unerwünschte Ausgabe der Ergebnisse des Befehls im
command window. Mit der Matlab function polyfit,
>>c=polyfit(t,lu,1) ←werden die Koeffizienten des Polynoms ersten Grades, das die Punkte
(ti , ln ui ), i = 1, . . . 9, interpoliert, unter c gespeichert. Das heisst der erste Koeffizient
c(1) ist die Steigung k der Regressionsgeraden, der zweite ist ihr Wert d auf der
senkrechten Achse. Mit
>>rg=c(1)*t+c(2); ←werden daher die Werte rg(i) der Regressionsgeraden im Vektor rg gespeichert. Die
Graphik in Abbildung 2.7 erhält man dann mit
>>plot(t,lu,’o:’,t,rg,’-’) ←Die Punkte (t(i),lu(i)) werden als kleine Kreise eingezeichnet und mit punktierten
Geradenstücken verbunden (Linientyp ’o:’), und die Punkte (t(i),rg(i)) werden nur
mit durchgezogenen Geradenstücken verbunden (Linientyp ’-’).
18
3. Parameteranpassung durch kleinste Quadrate
3.1. Logistisches Modell für Populationswachstum.
Wir haben die Populationsdaten der USA aus Beispiel 2.1 bis jetzt nicht zufriedenstellend
modelliert, die Exponentialfunktion eignet sich nur für die Phase des lawinenartigen
Wachstums. Um den gesamten Datensatz zu modellieren, brauchen wir eine Funktion,
die zunächst fast exponentiell anwächst, die aber später einen Wendepunkt durchläuft
und konkav (also mit Rückgang des Wachstums) fortfährt. Ein Beispiel einer solchen
S-förmigen Funktion ist die logistische Funktion (siehe Abbildung 3.1):
et
.
1 + et
Auf der Basis dieser Funktion formulieren wir als Modellgleichung
f (t) =
(3.1)
e
P (t) = K
t−t0
τ
t−t0 .
1+e τ
Die Größen K > 0, t0 ∈ R und τ > 0 sind dabei freie Parameter, wir werden sie so
wählen, daß dadurch der Datensatz möglichst genau wiedergegeben wird.
Wir werden später sehen, daß es auch ein einfaches Strukturmodell gibt, das auf genau
diese Funktion führt. Vorläufig sehen wir die Gleichung einfach als ein empirisches
Modell an, eine Funktion, die das augenfällige Erscheinungsbild des Datensatzes mehr
oder weniger gut wiedergeben kann. Es gibt natürlich viele andere Funktionen mit
ähnlichen Kurven, die wir ebensogut heranziehen könnten. Bestechend an der
logistischen Funktion ist zunächst die verhältnismäßig einfache Formel.
Wir werden in den nächsten beiden Unterabschnitten die Modellgleichung mit Papier
und Bleistift untersuchen. Wir wollen damit herausfinden, ob sie zur Modellierung des
Abbildung 3.1. Logistische Funktion
Logistische Funktion
1
0.9
0.8
0.7
logit(x)
0.6
0.5
0.4
0.3
0.2
0.1
0
-3
-2
-1
0
x
1
2
3
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
19
Abbildung 3.2. Der Parameter K in Modell 3.1
2
1.8
1.6
1.4
P
1.2
1
0.8
0.6
0.4
0.2
0
-4
-3
-2
-1
0
1
2
3
4
t
t0 = 0, τ = 1; durchgezogen: K = 1, strichliert: K = 2, strichpunkt: K = 0.5.
Systems grundsätzlich geeignet scheint, und welche Eigenschaften ein System hat, das
durch diese Gleichung beschrieben wird. In den späteren Abschnitten werden wir eine
Methode zur Bestimmung der Parameter kennenlernen.
3.2. Auswirkung der Parameter.
Merksatz 3.1. Ein empirisches Modell sollte nur mit sovielen freien Parametern
ausgestattet werden, wie unbedingt nötig sind, damit es dem Datensatz vernünftig
angepaßt werden kann. Die Parameter sollten möglichst unabhängige und
verschiedenartige Auswirkungen auf den Verlauf der Funktion nehmen. Besonders
günstig sind Parameter, die eine anschauliche Interpretation zulassen. Zuviele Parameter
(“Überparametrisierung”) führen zu folgenden Problemen:
• Lange Parametertabellen sind unübersichtlich, und die enthaltene Information
ist nur schwer abzulesen.
• Fast jeder Datensatz kann irgendwie dargestellt werden, man findet nicht heraus,
ob das Modell die Situation wirklich treffend karikiert.
• Wenn mehrere Parameter fast dieselbe Auswirkung auf den Verlauf der Funktion
haben, können annähernd gleiche Datensätze durch sehr verschiedene Parameter
gefittet werden: Wir sagen, die Parameter hängen instabil von den Daten ab. An
den Zahlenwerten der Parameter kann man dann nicht mehr sehen, ob man das
Modell zwei sehr verschiedenen oder sehr ähnlichen Datensätzen angepaßt hat.
• Der Rechenaufwand zum Anpassen vieler Parameter ist groß, die Programme
können durch die unvermeidbaren Rechenungenauigkeiten unzuverlässig werden,
vor allem dann, wenn die Parameter instabil von den Daten abhängen.
Wir wollen, bevor wir uns auf die Technik der Parameteranpassung einlassen, beurteilen,
welche Rolle die Parameter im Modell 3.1 spielen. Alle drei Parameter sind hier leicht zu
erklären.
20
Abbildung 3.3. Der Parameter t0 in Modell 3.1
1
0.9
0.8
0.7
P
0.6
0.5
0.4
0.3
0.2
0.1
0
-5
-4
-3
-2
-1
0
1
2
3
4
5
t
K = 1, τ = 1; durchgezogen: t0 = 0, strichliert: t0 = 1, strichpunkt: t0 = −1.
K skaliert die Größe der Population, bei Verdoppelung von K liefert das Modell doppelt
so große Populationen. Die Population liegt stets zwischen 0 und K, denn
e(t−t0 )/τ
< 1,
1 + e(t−t0 )/τ
der Nenner ist immer um 1 größer als der Zähler. Abbildung 3.2 zeigt das Modell bei
verschiedenen Werten von K.
0<
Eine Änderung von t0 bewirkt eine Parallelverschiebung der Kurve in t-Richtung um den
Betrag t0 . Abbildung 3.3 zeigt das Modell bei verschiedenen Werten von t0 . Am
deutlichsten erkennt man die Auswirkung von t0 , wenn man die Verschiebung des
Wendepunktes beobachtet.
Der Parameter τ skaliert die Zeit, er gibt an, wie schnell sich die Population entwickelt,
wie schnell also die Kurve der Funktion durchlaufen wird. Die Zeit t geht mit dem
0
in die Modellgleichung ein. Bei Verdoppelung von τ halbiert sich dieser
Exponenten t−t
τ
Exponent. Die Entwicklung der Population erfolgt nur halb so schnell. Abbildung 3.4
zeigt das Modell bei verschiedenen Werten von τ .
3.3. Eigenschaften des logistischen Modells.
Merksatz 3.2. Ein empirisches Modell steht und fällt damit, ob es das Verhalten des
dargestellten Systems treffend wiedergibt. Das qualitative Verhalten eines einfachen
Modells kann man mit Hilfe von Grenzwert- und Differentialrechnung, zum Beispiel mit
einer Kurvendiskussion beurteilen.
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
21
Abbildung 3.4. Der Parameter τ in Modell 3.1
1
0.9
0.8
0.7
P
0.6
0.5
0.4
0.3
0.2
0.1
0
-5
-4
-3
-2
-1
0
1
2
3
4
5
t
K = 1, t0 = 0; durchgezogen: τ = 1, strichliert: τ = 2, strichpunkt: τ = 0.5.
Wir fassen einige wichtige Eigenschaften des logistischen Modells zusammen:
(1) P (t) ist immer positiv.
(2) P (t) ist immer kleiner als K.
(3) limt→∞ P (t) = K, limt→−∞ P (t) = 0.
0
(4) Ist t−t
tief im negativen Bereich (also t viel kleiner als t0 ), so gilt
τ
t−t0
näherungsweise P (t) ≈ Ke τ .
(5) P (t) ist monoton wachsend.
(6) Für t < t0 ist P konvex, für t > t0 ist P konkav. An t = t0 , P = K2 liegt der
Wendepunkt.
Diese Eigenschaften des logistischen Modells (3.1) lassen sich aus Abbildung 3.1 ablesen
(man beachte, dass K und τ positiv sind), und sind mit Hilfe der ersten und der zweiten
Ableitung von P (t) mathematisch beweisbar.
Welche Bedeutung haben diese Eigenschaften im Zusammenhang mit der Modellierung
einer Population?
(1) Wenn das Modell nicht positive Bevölkerungszahlen voraussagen würde, wäre es
unsinnig.
(2) Das logistische Modell sagt voraus, daß der Lebensraum keine größere
Population als K tragen kann. Wir bezeichnen K auch als die Kapazität des
Lebensraumes. Aus dem Datensatz können wir nicht ablesen, ob das für das
System “Population der USA” wirklich zutrifft, aber es ist plausibel.
(3) Davon ist vor allem der Grenzwert für t → ∞ interessant: Das Modell sagt
voraus, daß sich auf Dauer die Populationsdichte der Kapazität annähert.
22
(4) Das logistische Modell sagt voraus: Solange die Population noch viel kleiner als
K ist, folgt sie ungefähr einem Exponentialgesetz. Wir haben schon gesehen,
daß das sehr gut auf unsere Daten paßt.
(5) Eine Population mit einem signifikanten Bevölkerungsrückgang läßt sich nicht
mit Modell (3.1) darstellen.
(6) Konvexität bedeutet, daß die Wachstumsrate zunimmt, Konkavität, daß die
Wachstumsrate abnimmt. Der Wechsel von konvex zu konkav im Wendepunkt
simuliert unsere Beobachtung, daß das Bevölkerungswachstum der USA in den
letzten beobachteten Jahrzehnten zurückgegangen ist.
Mit Hilfe der Information aus Punkt (6) und (4) kann man graphisch die Parameter K, τ
und t0 schätzen. Betrachtet man Abbildung 2.1 und sucht nach der geeignetsten Stelle
für den Wendepunkt, drängt sich die Stelle t = 1910 auf. Zu diesem Zeitpunkt war die
Größe von P etwa 90000. Wegen Punkt (6) wählen wir daher
t0 = 1910 und K = 180000.
Sind t0 und K einmal fixiert, kann man sie und eines der Datenpaare (t, P (t)) aus
Beispiel 2.1 in das Modell (3.1) einsetzen. Dies ergibt eine Gleichung aus der man eine
Schätzung für τ ausrechnen kann. Für jedes der Datenpaare würde man voneinander
geringfügig abweichende Werte für τ bekommen. Eine kürzere Variante benützt die
Information aus Punkt (4): für die Anfangsjahre wird die logistische Funktion durch die
Exponentialfunktion P (t) = Ke(t−t0 )/τ approximiert. Setzt man das Jahr 1790, P ≈ 4000
und die bereits bestehenden Schätzungen für t0 und K ein, ergibt sich
4000 = 180000e(1790−1910)/τ , also
τ = −120/ ln
4
≈ 31.52
180
In diesem Beispiel kommt man relativ schnell zu einer Schätzung der Modellparameter
aus den Daten. In anderen Situationen wird man andere Ideen und Tricks suchen und
probieren müssen. Insbesondere kann man versuchen, mit mathematischen Methoden
Parameter-abhängige Merkmale des Modells zum Vorschein zu bringen.
In unserem Beispiel könnten wir nun die geschätzten Parameter verwenden und eine
graphische Probe machen. Durch ’Drehen’ an den Parametern könnten wir eventuell
noch eine Verbesserung der Anpassung erreichen. Es gibt aber eine systematische
Vorgangsweise, die optimalen Parameter zu suchen:
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
23
Abbildung 3.5. Daten und Approximationsfehler
x10 5
Approximationsfehler
1.8
*
*
1.6
o
*
1.4
o
*
Population
1.2
o
o
1
*
o
0.8
*o
o
0.6
o
0.4
o
*
o
o
0.2
0
1780
o
*
o
*
1800
o
*
o
*
1820
o
*
o
*
1840
*
*
*
*
1860
1880
1900
1920
1940
1960
Jahr
Kreise: Daten, Sterne: Modellvoraussage
3.4. Methode der kleinsten Quadrate.
Nur sehr einfache Modelle erlauben eine Schätzung der Parameter mit graphischen
Methoden wie oben. Die folgende Methode funktioniert auch für komplizierte Modelle
und läßt sich automatisieren.
Um sie zu erklären, zeigen wir in Abbildung 3.5 die Daten und den Verlauf des
Modells (3.1) mit Parametern, die wir absichtlich so gewählt haben, daß die Daten
schlecht wiedergegeben werden. (Die Parameter sind K = 200000, t0 = 1910, τ = 20.)
Kreise bedeuten die Felddaten, die Sterne sind die Populationswerte, die das Modell
voraussagt. Jedes Felddatum weicht vom Fit, dem vorausgesagten Wert, um einen
gewissen Wert ab, den Approximationsfehler. Durch senkrechte Linien haben wir den
Approximationsfehler deutlich gemacht. Für das Jahr 1950 liefert dieser Parametersatz
einen Fehler von ungefähr -25000, für das Jahr 1880 etwa 15000. Das negative Vorzeichen
soll andeuten, daß die Felddaten unter den Schätzwerten aus dem Modell liegen.
Natürlich ist zu erwarten, daß bei jeder Wahl von Parametern ein Approximationsfehler
auftritt, das Modell kann nicht alle Details der Daten perfekt wiedergeben. Wenn die
Parameter aber gut gewählt sind, werden die Approximationsfehler insgesamt sehr klein
ausfallen.
Drücken wir noch einmal dasselbe mathematisch aus: Gegeben sind ein Satz von n
Datenpaaren (t1 , y1 ), · · · , (tn , yn ). Das sind die beobachteten Daten, die wir als Tabelle in
Beispiel 2.1 angegeben haben. Gegeben ist ferner eine Funktion f (t; α1 , · · · , αm ) als
Modell, die außerdem von m Parametern α1 , · · · , αm abhängt. In unserem Fall ist das die
logistische Funktion (3.1), abhängig von t und den drei Parametern K, t0 und τ :
f (t; K, t0 , τ ) = K
e
t−t0
τ
1+e
t−t0
τ
.
24
Das Modell sagt anstelle der y-Daten die Werte
f (t1 ; α1 , · · · , αm ) · · · f (tn ; α1 , · · · , αm )
voraus. Damit entsteht an jedem Datenpaar ein Approximationsfehler
ei = yi − f (ti , α1 , · · · , αm ).
In unserem Fall ist das
ei = yi − K
e
ti −t0
τ
ti −t0 .
1+e τ
(Der Buchstabe e hat sich für “error”=Fehler eingebürgert. Bitte nicht mit der
Eulerschen Zahl aus der Exponentialfunktion verwechseln.) Ein gutes Modell mit guten
Parametern liefert kleine Werte für die Approximationsfehler ei .
Um den Gesamtfehler zu bewerten, dürfen wir nicht einfach die einzelnen
Approximationsfehler summieren: Dabei könnten sich große positive und negative Fehler
wegheben und einen kleinen Gesamtfehler vortäuschen. Es gibt verschiedene Arten, den
Gesamtfehler zu bewerten. Die bekannteste besteht darin, die einzelnen Fehler erst zu
quadrieren. Damit werden alle Vorzeichen positiv.
n
X
E=
e2i = e21 + e22 + · · · + e2n .
i=1
Wir betrachten diesen Ausdruck als den Gesamtfehler, nach dem wir die Wahl der
Parameter beurteilen können. Die Parameter passen umso besser, je geringer der
Gesamtfehler ist. Damit stehen wir vor der folgenden Aufgabe: Finde jene Parameter, für
die der Gesamtfehler am kleinsten ist.
Merksatz 3.3 (Kleinste Quadrate). Gegeben seien n Datenpaare (t1 , y1 ), · · · , (tn , yn )
und eine Modellfunktion
f (t; α1 , · · · , αm ),
abhängig von Parametern α1 , · · · , αm .
Die Parameteranpassung nach der Methode der kleinsten Quadrate ist die Suche nach
jenen Parametern α1 , · · · , αm , für die die Quadratsumme der Approximationsfehler
n
X
E(α1 , · · · , αm ) =
(yi − f (ti ; α1 , · · · , αm ))2
i=1
am kleinsten wird.
Für die Bewertung der Güte der Anpassung eignet sich noch besser der relative
quadratische Approximationsfehler
Pn
(yi − f (ti ; α1 , · · · , αm ))2
Pn 2
.
Erel (α1 , · · · , αm ) = i=1
i=1 yi
Er unterscheidet sich von E nur durch Division durch einen konstanten Faktor. Wenn
man den Datensatz einfach durch f = 0 approximiert, erhält man Erel = 1. Erel = 0
bedeutet völlig fehlerfreie Darstellung der Daten durch das Modell.
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
25
Auch die Regressionsgerade an einen Datensatz ist in Wirklichkeit eine Gerade
f (t) = kt + d, deren Parameter k und d nach der Methode der kleinsten Quadrate an den
Datensatz angepaßt sind.
3.5. Optimumsuche mit Computer.
Wir stehen jetzt vor einer Extremwertaufgabe, einem Optimierungsproblem:
Merksatz 3.4. Ein Optimierungsproblem hat die folgende Gestalt: Gegeben ist eine
Menge Ω von Vektoren in Rm (können auch, im einfachsten Fall, Skalare in R1 sein),
und eine Funktion J, die jedem Vektor x aus Ω eine reelle Zahl J(x) zuordnet. Gesucht
ist jener Vektor x ∈ Ω, für den die Funktion J(x) am kleinsten (oder am größten) ist.
Die Menge Ω heißt der zulässige Bereich, J heißt die Zielfunktion, und die Lösung x des
Optimumproblems heißt die optimale Lösung.
 
K
In unserem Fall der Parameteridentifizierung besteht Ω aus allen Tripeln x =  t0  mit
τ
K ≥ 0, t0 ∈ R und τ > 0. Die Zielfunktion J ist die Quadratsumme E der
Approximationsfehler. Gesucht ist das Minimum von E.
Um Mißverständnisse zu vermeiden, betonen wir, daß sich Optimierungsaufgaben
durchaus nicht nur im Zusammenhang mit Parameteranpassung ergeben müssen.
Überall, wo Mathematik eingesetzt wird, um unter verschiedenen Möglichkeiten die beste
zu bestimmen, kommt man letztlich auf ein Optimierungsproblem.
Wenn J nur von einer reellen Zahl, nicht einem Vektor, abhängt, läßt sich oft durch
Nullsetzen der ersten Ableitung der optimale Parameter bestimmen. Auch bei
Funktionen von mehreren Parametern kann das Optimum gelegentlich durch Nullsetzen
aller partiellen Ableitungen bestimmt werden, aber meistens ergeben sich dabei
Gleichungen, die mit Bleistift und Papier nicht mehr zu lösen sind. Man ist dann auf
sogenannte numerische Verfahren angewiesen.
Es gibt viele verschiedene Algorithmen (Rechenverfahren) zur numerischen
Minimumsuche. Grundsätzlich arbeiten sie nach dem folgenden Prinzip: Der Benützer
gibt einen Startwert, eine erste Schätzung für die optimale Lösung, an. Das Programm
tastet ab, in welche Richtung die Parameter zu verändern sind, damit sich die
Zielfunktion verkleinert. In diese Richtung geht es ein Stück und verbessert dadurch den
Schätzwert der optimalen Lösung. Vom neuen Schätzwert tastet es sich weiter. Wenn die
Schritte, die von Schätzung zu Schätzung gemacht werden, sehr klein geworden sind, und
sich der Wert der Zielfunktion nur (mehr) wenig ändert, wird die letzte Schätzung als
endgültige Näherung für die optimale Lösung angenommen.
26
Abbildung 3.6. Ein lokales Minimum fängt einen Optimumsucher
Wege eines Optimumsuchers von zwei Startwerten aus:
Von Startpunkt 1 wird das Minimum, von 2 aus nur ein lokales Minimum gefunden.
Merksatz 3.5. Ein numerisches Verfahren zur Optimumsuche sucht eine Stelle, an der
eine gegebene Zielfunktion ein Optimum hat. Ausgehend von einem Startwert geht es
systematisch und schrittweise zu immer besseren Stellen weiter. Wenn ein
Abbruchkriterium erfüllt ist, indem sowohl die Schrittweite als auch die Änderung der
Funktionswerte kleiner als eine vorgegebene Toleranz wird, bricht das Verfahren die
Suche ab und gibt die letzte erreichte Stelle als Lösung aus.
Der Benützer eines numerischen Optimumsuchers muß bereitstellen
(1) ein Unterprogramm, das die Zielfunktion (und, bei Verwendung eines
Gradientenverfahrens, auch ihre Ableitungen) aus den Argumenten (das sind
hier die Parameter) berechnet,
(2) einen Startwert (nach Möglichkeit in der Nähe der Lösung),
(3) eine Toleranz in Bezug auf die Ungenauigkeit der Lokalisation der
Optimal-Stelle.
Aus der algorithmischen Vorgangsweise ergeben sich auch die Probleme: Das Programm
tastet zwar die Umgebung einer gefundenen Schätzung ab, aber es kann nicht beurteilen,
ob nicht weit entfernt eine noch bessere Lösung existiert. Deshalb sind Zielfunktionen
mit mehreren lokalen Minima schwierig zu minimieren: Der Optimumsucher läuft das
lokale Minimum an, auf das er ausgehend von seinem Startpunkt als erstes trifft, und
erkennt nicht, ob er damit das globale Minimum erreicht hat. Im Zweifelsfall hilft Starten
von verschiedenen Startpunkten aus, und ein Vergleich der verschiedenen erreichten
Lösungen. Natürlich gibt es Zielfunktionen, die gar kein Minimum besitzen. Das
Programm läuft dann theoretisch bis ins Unendliche weiter, wenn es nicht eigens dagegen
abgesichert ist. Auch bei einer Zielfunktion mit Minimum kann sich ein Programm ins
Unendliche verlaufen, wenn das Minimum vom Startwert weit entfernt ist.
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
27
Nach der Art, wie das Verfahren seine Suchrichtung wählt, unterscheiden sich die
verschiedenen Optimierungsmethoden. Wenn die partiellen Ableitungen der Zielfunktion
greifbar sind, weiß man auch die Richtung, in die der schnellste Abfall zu erwarten ist
(nämlich den Gradienten). Gute Verfahren folgen allerdings nicht stur dem Gradienten,
sondern verwerten auch Information aus den bisher zurückgelegten Schritten. Sind die
Ableitungen der Zielfunktion nicht greifbar, müssen immer genügend viele Punkte in der
Nähe der aktuellen Schätzung abgetastet werden, um zu beurteilen, welche Richtung den
meisten Erfolg verspricht. Ferner gibt es sogenannte ’genetische Algorithmen’, die eine
Gruppe von Schätzungen nach den Prinzipien Kreuzung, Mutation und Selektion
verändern und dadurch - wie bei der biologischen Evolution - an die Daten anpassen. Es
gibt kein allgemein bestes Verfahren. Wenn bei einem Problem ein Verfahren kein
Minimum finden will, kann man mit einem anderen Verfahren noch immer Erfolg haben.
Es kann aber auch vorkommen, dass das Parameteranpassungsproblem gar keine
befriedigende Lösung hat, weil das Modell nicht zu den Daten passt. Bei vergeblicher
Suche nach einem geeigneten Minimum ist das Modell in Frage zu stellen, und eine
Abänderung des Modells ins Auge zu fassen.
3.6. Lösung des Parameteranpassungsproblems in MATLAB.
Zunächst muss man eine Matlab function schreiben, nennen wir sie zielfkt, die zu
gegebenen Parametern K, t0 , τ den Approximationsfehler, das ist die Summe der
Fehlerquadrate, berechnet. Eine solche function wird durch einen file namens zielfkt.m
definiert, der folgende erste und letzte Zeilen enthält:
function fehler = zielfkt(x)
.
.
.
fehler = · · ·
Beim Aufruf dieser function enthält der Argumentvektor x die aktuellen Parameter und
der output-Variablen fehler muss die zugehörige Summe der Quadrate zugewiesen
werden. Im command window gibt man dann die Befehle
>>options=optimset(’TolX’,1e-4,’MaxFunEvals’,1000); ←>>xstart=[1 2 3]; ←>>xend=fminsearch(’zielfkt’,xstart,options) ←um (zum Beispiel) die Toleranz auf 10−4 , die maximal zugelassene Anzahl von
Auswertungen von zielfkt auf 1000 und den Startvektor auf [1 2 3] zu setzen. Die
Matlab Routine fminsearch beginnt bei xstart mit der Suche nach einem Minimum
und ruft dazu die vom Benützer programmierte function zielfkt immer wieder auf. Der
Such-Algorithmus bestimmt systematisch die Stellen an denen die Auswertung von
zielfkt benötigt wird, diese sollte entsprechend flexibel programmiert sein. Findet der
Algorithmus Parameter, die zur geforderten Toleranz passen, dann werden diese (hier)
dem Vektor xend zugewiesen.
28
Abbildung 3.7. Schritte des Nelder-Mead-Verfahrens
Dick ausgezogen: Neues Tetraeder
3.7. Das Verfahren von Nelder-Mead.
Nur für Interessierte erklären wir, wie das in fminsearch implementierte Suchverfahren von
Nelder-Mead funktioniert.
Nehmen wir an, daß die Zielfunktion name als Argument (so wie im obigen Beispiel) einen
dreidimensionalen Vektor hat. Wir geben einen dreidimensionalen Startvektor x0 vor, und
rufen damit fminsearch(’name’,x0) auf (das Weglassen des dritten Arugments beim Aufruf
von fminsearch ist gleichbedeutend mit der Wahl der default Optionen).
Was läuft im Computer ab?
Die Idee des Programms ist, die optimale Lösung durch ein Tetraeder von 4 Schätzpunkten
einzukreisen. (Wäre name eine Funktion eines zweidimensionalen Vektors, würde das Optimum
mit Dreiecken eingekreist. In höheren Dimensionen m nennt man das Analogon zu Dreieck und
Tetraeder ein Simplex. Es hat immer m + 1 Ecken.) Der erste Schritt besteht darin, in der
Umgebung des Startvektors 4 Punkte als Ecken des ersten Tetraeders zu fixieren. An jedem
Eckpunkt wird der Wert der Zielfunktion ausgerechnet. Alle weiteren Schritte bestehen darin,
das einmal erreichte Tetraeder zu verbessern.
Einer der Eckpunkte des gerade aktuellen Tetraeders hat den größten, also schlechtesten Wert
der Zielfunktion. Wir konstruieren ein neues Tetraeder, indem gerade dieser Eckpunkt durch
einen anderen ersetzt wird, an dem die Zielfunktion einen möglichst guten (=kleinen) Wert
annimmt. Die verschiedenen Strategien, die wir jetzt erklären werden, sind in Abbildung 3.7
zusammengefaßt.
Der erste Versuch besteht darin, daß man den schlechtesten Eckpunkt am Schwerpunkt der
gegenüberliegenden Seitenfläche spiegelt. Dieser Schritt erscheint sinnvoll, weil er von der
schlechtesten Ecke weg in die Richtung besserer Werte der Zielfunktion führt. Am neuen
Eckpunkt wird der Wert der Zielfunktion berechnet. Ob die Spiegelung ein Erfolg war,
entscheidet sich an Hand der Zielfunktion.
3. PARAMETERANPASSUNG DURCH KLEINSTE QUADRATE
29
Ist die Zielfunktion am neuen Eckpunkt sehr gut, nämlich kleiner als die Werte an allen
anderen Eckpunkten, dann ist die Suchrichtung erfolgversprechend. Wir strecken den Eckpunkt
doppelt so weit von der Fläche weg, an der er gespiegelt wurde, dadurch wird die Suche in die
erfolgversprechende Richtung vorangetrieben. Wenn auch dieser Eckpunkt besser als alle
anderen ist, wird er beibehalten. Ist aber der Eckpunkt des ausgedehnten Simplex schlechter als
eine der anderen Ecken, behalten wir den gespiegelten Eckpunkt bei.
Wir behalten auch den gespiegelten Eckpunkt bei, wenn die Zielfunktion dort zwar nicht besser
als an allen anderen Ecken, aber nicht die schlechteste des ganzen verbleibenden Tetraeders ist.
Wenn aber der gespiegelte Eckpunkt schlechter als alle Ecken außer höchstens der
weggelassenen ist, ist die Spiegelung offenbar die falsche Strategie zur Optimumsuche. Man
geht davon aus, daß man bereits so nahe am Optimum ist, daß jeder große Schritt nur mehr
Verschlechterungen bringt. Daher reduziert man die Schrittweite. Man wählt vom gespiegelten
Eckpunkt und seinem Urbild den Punkt mit der besseren Zielfunktion, und halbiert die Strecke
von diesem Punkt zum Schwerpunkt der gegenüberliegenden Seitenfläche (an dem vorhin
gespiegelt wurde). Wenn der neue Eckpunkt nicht schlechter als alle anderen verbleibenden
Ecken ist, behält man ihn bei, das Tetraeder wurde kontrahiert.
Wenn auch die Kontraktion nichts hilft, wählt man aus dem gesamten ursprünglichen Tetraeder
die beste Ecke, und ersetzt alle anderen Ecken durch die Mittelpunkte zwischen ihnen und der
besten Ecke. Damit werden alle Seitenlängen des Tetraeders halbiert, es schrumpft.
Nach der Bestimmung der neuen Ecke und damit des neuen Tetraeders prüft man, wie weit die
Ecken voneinander entfernt sind, also die Seitenlängen des Tetraeders. Sind alle Seiten kürzer
als die geforderte Toleranz, betrachtet man die Ecke mit dem besten Wert der Zielfunktion als
Lösung und beendet das Verfahren.
3.8. Die optimalen Parameter des logistischen Modells.
Wir berichten kurz von den Erfahrungen und Ergebnissen bei der Optimumsuche für das
Populations-Modell. Wie wir schon gesehen haben, sind geeignete Parameter von sehr
unterschiedlicher Größenordnung, nämlich K ∼ 2 × 105 , t0 ∼ 1.8 × 103 , τ ∼ 3 × 10. Dies
ist eine für die Optimumsuche sehr nützliche Vor-Information, denn wenn man keine
Ahnung hat, wo man suchen soll, dann ist es schwierig etwas Passendes zu finden.
Ausserdem ist es für das automatische Such-Verfahren günstig, wenn die Auswirkung der
Variation für alle Parameter ähnlich groß sind (die Tetraeder im Nelder-Mead Verfahren
werden sonst sehr spitz und flach).
Daher verwenden wir statt K, t0 , τ die skalierten Parameter
K̃ = K/105 , t˜0 = t0 /103 , τ̃ = τ /10. Wir schreiben also ein Unterprogramm, das die
Summe der Fehlerquadrate aus K̃, t˜0 , τ̃ berechnet. Wenn man dieses Unterprogramm
dem Suchverfahren zur Verfügung stellt, dann variiert es die Argumente K̃, t˜0 , τ̃ und
zwar alle drei etwa zwischen 1 und 4.
Startet man das Nelder-Mead Verfahren bei (1,1,1), dann konvergiert es ohne
Fehlermeldungen oder Warnungen nach (0.5536, 1.1090, 1.1022). Graphisch ist allerdings
sofort klar, dass dies keine guten Parameter sind: die zugehörige Funktion ist beinahe
konstant und zwar der Mittelwert der Daten. Dort befindet sich scheinbar ein lokales
Minimum der Summe der Fehlerquadrate.
30
Abbildung 3.8. Logistisches Modell für Beispiel 2.1
x10 4
16
o
14
o
o
12
o
Population
10
o
8
o
o
6
o
4
o
o
o
2
0
1780
o
o
o
1800
o
o
1820
o
1840
1860
1880
1900
1920
1940
1960
Jahre
Kreise: Felddaten. Durchgezogen: Modell
Startet man bei (1,2,3), dann konvergiert das Verfahren nach (2.007, 1.9158, 3.2496). Die
zugehörige Graphik sieht gut aus. Beim weiteren Experimentieren mit den Startwerten
bekommt man entweder sinnlos schlechte Ergebnisse (lokale Minima) oder sehr gute
Ergebnisse ganz in der Nähe von K = 20070, t0 = 1915.8, τ = 32.496. Dort scheint sich
also das globale Minimum der Summe der Quadrate zu befinden. Den Vergleich der
Daten und dem Modell mit diesen Parametern zeigt Abbildung 3.8. Die
Übereinstimmung zwischen Modell und Daten ist hervorragend. Vergleichen Sie auch die
Ergebnisse mit den graphisch geschätzten Werten. Nur den Parameter K haben wir
deutlich unterschätzt.
Trotzdem wäre es falsch zu glauben, das empirische Modell für das
Bevölkerungswachstum in den USA gefunden zu haben. Setzt man zum Beispiel in das
Modell mit den soeben bestimmten Parametern das Jahr 2000 ein, erhält man
P (2000) ∼ 186 Millionen, während laut Census 2000 tatsächlich 281 Millionen Menschen
in den USA lebten. Bei der Bestimmung der Parameter unter der Berücksichtigung von
Daten bis 2000, würden die optimalen Werte sicher anders ausfallen.
4. MENGENBILANZEN
31
4. Mengenbilanzen
4.1. Strukturmodelle.
Wir haben bisher unsere Modellgleichungen an Hand existierender Felddaten ausgesucht,
die Modellfunktion sollte eine Kurve ergeben, die die Felddaten möglichst genau
wiedergeben kann. Ein völlig anderer Weg, ein Modell zu einem System zu gewinnen,
besteht darin, die Einzelteile des Systems und ihre Wechselwirkungen untereinander zu
betrachten, und diese Wechselwirkungen durch mathematische Gleichungen, sozusagen
durch Naturgesetze, wiederzugeben. Wenn man alle Teile zusammenfügt, erhält man ein
System von Gleichungen, das die Einzelteile des modellierten Systems und ihre
Zusammenhänge wiederspiegelt. Ein solches Modell heißt ein Strukturmodell.
Merksatz 4.1. Ein Strukturmodell bildet die Einzelteile eines Systems und ihre
Wechselwirkungen ab. Es verwendet und liefert daher nicht nur Information über die
Reaktion des Systems als Ganzes, sondern auch über die inneren Wirkungsweisen, die zu
dieser Reaktion führen.
Vorteile des Strukturmodells:
• Information über die Wirkungsweise des Systems wird gewonnen.
• Die Parameter haben (oft) eine natürliche Bedeutung.
• Der Modellbildungsvorgang trägt zum Verständnis des Systems bei.
Vorteile des empirischen Modells:
• Die Gleichungen sind einfach, übersichtlich, und rechnerisch leicht zu behandeln.
• Wenige Parameter reichen meist.
• Selbst Systeme, von deren Wirkungsweise man gar keine Ahnung hat, lassen sich
modellieren.
Es hängt vom modellierten System und dem Zweck der Modellbildung ab, welche
Hilfsmittel man zur Erstellung eines Strukturmodells heranzieht. Da es sehr
unterschiedliche Systeme gibt, kommen alle Zweige der Mathematik zum Tragen. Eine
der wichtigsten Modellierungsstrategien ist aber das Erstellen von Mengenbilanzen, das
wir an den folgenden Beispielen erörtern werden.
4.2. Modell eines Kessels.
Beispiel 4.2. Ein zylindrischer Kessel (Abbildung 4.1) besitzt einen Zufluß, durch den
ständig mit gleichbleibender Zuflußrate Wasser zugeführt wird. Am Boden des Kessels
befindet sich ein Abfluß, der durch ein enges Rohr führt. Auf Grund des
Strömungswiderstandes ist die abfließende Wassermenge proportional zum
hydrostatischen Druck am Kesselboden. Dieser ist wieder proportional zur Höhe der
Wassersäule im Kessel.
(Die Zuflußrate, die Größe des Kessels, und den Strömungswiderstand des Abflusses
setzen wir als gegeben voraus.)
(1) Welcher Wasserstand wird sich im Kessel einstellen?
32
Abbildung 4.1. Kessel aus Beispiel 4.2
(2) Wie entwickelt sich der Wasserstand im Verlauf der Zeit, wenn wir von einem
leeren Kessel ausgehen und dann den Zufluß voll aufdrehen?
Wir sehen einen wesentlichen Unterschied zwischen den beiden Fragestellungen:
Frage (1) bezieht sich auf einen Dauerzustand, den Wasserstand, der sich auf Dauer im
Kessel einstellen wird. Wir erwarten, daß sich der Wasserstand nicht mehr (wesentlich)
ändern wird, wenn dieser Zustand einmal erreicht ist — die Zeit spielt dann keine Rolle
mehr. Später werden wir einen solchen Zustand einen Gleichgewichtszustand nennen.
Der Frage unterliegt eine statische Betrachtungsweise. Ihre Antwort wird eine einzige
Zahl sein: Der Wasserstand ist so und so hoch. Dagegen geht Frage (2) von einer
dynamischen Betrachtungsweise aus: Wir fragen nach dem Verlauf eines Vorgangs im
Lauf der Zeit. Die Antwort könnte zum Beispiel durch eine Kurve gegeben werden, die
den Wasserstand in Abhängigkeit von der Zeit darstellt, die seit dem Einschalten des
Zuflusses verstrichen ist.
Merksatz 4.3. Eine statische Betrachtungsweise bezieht sich auf einen Dauerzustand
des Systems, der sich im Laufe der Zeit nicht ändert. Eine dynamische Betrachtungsweise
untersucht, wie sich der Zustand des Systems im Lauf der Zeit verändert.
Statische Modellierung. Wir wenden uns zunächst der Beantwortung der
statischen Fragestellung zu. Wir benennen alle Größen, die im System eine Rolle spielen,
und legen ihre Maßeinheiten fest. Die Form des Kessels ist durch seinen Radius r und
seine Gesamthöhe H gegeben. Als Einheit kommt eine Längeneinheit in Frage. Im
4. MENGENBILANZEN
33
SI-Einheitensystem ist die Längeneinheit Meter (m). Mit z bezeichnen wir die Zuflußrate,
also die Wassermenge, die pro Zeiteinheit durch den Zufluß strömt. Damit wird also eine
Volumsmenge auf Zeit bezogen, die Einheit ist (im SI-System) Kubikmeter pro Sekunde
(m3 /s). Mit a bezeichnen wir die Abflußrate, also die Wassermenge, die pro Zeiteinheit
durch den Abfluß fließt. Die Einheit ist wieder m3 /s. Während aber z gegeben ist,
müssen wir a erst aus anderen Gesetzmäßigkeiten herleiten. Letztlich suchen wir nach
dem Wasserstand im Kessel, den wir mit h bezeichnen, die Einheit ist wieder m.
Die Abflußrate ergibt sich aus dem hydrostatischen Druck d am Kesselboden, die Einheit
des Drucks im SI-System ist ein Pascal (1Pa = 1kg/ms2 ). Der Zusammenhang zur
Abflußrate ergibt sich durch ein Proportionalitätsgesetz:
(4.1)
a = d/w
Dabei ist w der Strömungwiderstand des Abflußrohres in kg/m4 s. (Wie man zu dieser
Einheit kommt, erklären wir in einem späteren Kapitel.)
Der hydrostatische Druck errechnet sich wieder aus der Höhe der Wassersäule über dem
Boden, multipliziert mit dem spezifischen Gewicht des Wassers:
(4.2)
d = γh
Dabei ist γ das spezifische Gewicht von Wasser, in kg/m2 s2 . Der Zahlenwert ist bekannt:
kg
γ ≈ 104 2 2 .
ms
Tabelle 4.1 faßt noch einmal alle Größen und Gesetze des statischen Modells zusammen.
(Die Mengenbilanz werden wir gleich erklären.)
Unsere bisherige Arbeit hat darin bestanden, alles Wissen aus den Angaben zu
formalisieren, indem wir den wesentlichen Größen geeignete Namen gegeben haben, und
die verbal dargestellten Beziehungen durch Gleichungen ausgedrückt haben. Der
springende Punkt zur Beantwortung der Frage (1) ist aber eine einfache Mengenbilanz:
Wenn der Wasserstand gleichbleiben soll, muß genauso viel zufließen wie abfließen. Im
Gleichgewicht gilt also
z = a.
Wir versuchen, alle Hilfsgrößen durch die gesuchte Gröse h und die gegebenen Parameter
auszudrücken:
d = γh nach (4.2),
d
γ
a = = h nach (4.1).
w
w
Damit ergibt sich aus der Mengenbilanz
γ
z = h.
w
Wir können damit den Wasserstand aus den bekannten Größen berechnen.
wz
(4.3)
h=
γ
34
Tabelle 4.1. Statisches Modell zu Beispiel 4.2
Größe
r
H
z
w
γ
h
a
d
Einheit
m
m
m3
s
kg
m4 s
kg
104 m2 s2
m
m3
s
kg
ms2
Modellgrößen
Benennung
Radius
Höhe
Zuflußrate
Kommentar
bekannt, unnötig
bekannt, unnötig
bekannt
Str.widerstand
bekannt
spez. Gew. H2 O bekannt
Wasserstand
gesucht
Abflußrate
hydr. Druck
Modellgleichungen
d
,
w
d = γh.
a=
statische Mengenbilanz
z=a
γ
z= h
w
Überzeugen wir uns noch, daß die Formel ein plausibles Ergebnis liefert. Wir sehen aus
der Formel: Je größer z, also die Zuflußrate, desto größer der Wasserstand, der sich
einstellen wird. Je größer w, der Strömungswiderstand, also je enger der Abfluß, desto
höher der Wasserstand. Beide Zusammenhänge erscheinen intuitiv durchaus
einleuchtend. Interessant ist auch, daß die Abmessungen des Kessels überhaupt keinen
Einfluß auf den Wasserstand nehmen, natürlich nur, solange er nicht überläuft.
¤
Merksatz 4.4. Eine statische Mengenbilanz beruht auf dem folgenden Grundsatz: Die
vorhandene Menge ändert sich nicht, befindet sich also im Gleichgewicht, wenn die
Summe aller Zuflüsse Z gleich der Summe aller Abflüsse A ist:
0 = Z − A.
Dynamische Modellierung. Wir wollen jetzt nicht das statische Gleichgewicht
des Systems betrachten, sondern die Art, wie sich der Wasserstand, also der Zustand des
Systems, im Lauf der Zeit ändern kann. Benennen wir die Zeit mit t (in Sekunden s).
Der Zeitpunkt t = 0 sei der Zeitpunkt des Einschaltens des Zuflusses. Einige unserer
Modellgrößen bleiben immer gleich, dazu gehören die Abmessungen des Kessels r, H, das
spezifische Gewicht γ des Wasser, der Abflußwiderstand w und, sobald der Zufluß
aufgedreht ist, die Zuflußrate z. Das sind also konstante (= von der Zeit unabhängige)
4. MENGENBILANZEN
35
Tabelle 4.2. Dynamisches Modell zu Beispiel 4.2
Größe
t
V (t)
r
H
z
w
γ
π
h(t)
a(t)
d(t)
λ
Einheit
s
m3
m
m
m3
s
kg
m4 s
kg
104 m2 s2
3.14
m
m3
s
kg
2
ms
γ
1
2
r πw s
Modellgrößen
Benennung
Zeit
Volumen
Radius
Höhe
Zuflußrate
Kommentar
gesucht
bekannt
bekannt, unnötig
bekannt
Str.widerstand bekannt
spez. Gew.
Kreiszahl
Wasserstand
Abflußrate
bekannt
bekannt
gesucht
hydr. Druck
Hilfsparameter
Modellgleichungen
d(t)
,
w
d(t) = γh(t),
a(t) =
V (t) = r2 πh(t).
dynamische Mengenbilanz
d
V (t) = z − a(t) = z − λV (t).
dt
Parameter des Systems. Dagegen wird sich der Wasserstand im Lauf der Zeit ändern,
und mit ihm alle Größen, die von ihm abhängen. Dazu gehören der hydrostatische Druck
und die Abflußrate, außerdem eine Größe, auf die wir bisher gar nicht geachtet haben,
und die jetzt eine zentrale Rolle einnehmen wird, das Volumen V , also die Gesamtmenge
des Wassers im Kessel (in m3 ). Diese Größen betrachten wir jetzt als Funktionen, die von
t abhängen: h(t), d(t), a(t) und V (t).
An den Zusammenhängen zwischen Wasserstand, Druck und Abflußrate ändert sich
nichts im Vergleich zum statischen Modell, außer daß diese Größen jetzt zu
verschiedenen Zeitpunkten verschiedene Werte annehmen. Das Volumen errechnet sich
aus Grundfläche und Höhe der Wassersäule im Kessel:
(4.5)
d(t)
,
w
d(t) = γh(t),
(4.6)
V (t) = r2 πh(t).
(4.4)
a(t) =
Tabelle 4.2 zeigt das gesamte dynamische Modell, die Mengenbilanz erklären wir gleich.
36
Abbildung 4.2. Erste Ableitung als Wachstumsrate
Weil sich jetzt die Wassermenge im Kessel ändern kann, müssen zufließende und
abfließende Menge nicht mehr gleich sein. Beide tragen zur Gesamtänderung des
Wasservolumens im Kessel bei. Das Volumen V (t) ist eine Funktion von t. Abbildung 4.2
zeigt schematisch, wie das Volumen in Abhängigkeit von der Zeit verlaufen könnte. In
einem kurzen Zeitintervall τ wächst das Volumen um den Betrag V (t + τ ) − V (t) an.
Wenn τ sehr klein ist, können wir anstatt der Volumskurve selbst ihre Tangente zur
Abschätzung des Wachstums heranziehen. Die Steigung k der Tangente an die
Volumskurve beschreibt die Wachstumsrate des Volumens: In einem sehr kurzen
Zeitintervall der Länge τ Sekunden wächst das Volumen ungefähr um kτ Kubikmeter an.
Rechnerisch erhält man die Steigung der Tangente durch Differenzieren der
Volumsfunktion nach t:
d
k = V (t).
dt
Zwei Faktoren tragen zu dieser Volumsänderung im Zeitraum von τ Sekunden bei: Der
Zufluß liefert zτ Kubikmeter Wasser, der Abfluß saugt a(t)τ Kubikmeter Wasser ab.
Daher ist die gesamte Volumsänderung
kτ = zτ − a(t)τ.
Durchkürzen und Umschreiben von k als erste Ableitung von V liefert:
d
V (t) = z − a(t).
dt
(4.7)
Versuchen wir, alle unbekannten Größen durch V (t) und die bekannten Parameter
auszudrücken:
V (t)
nach (4.6),
r2 π
γ
d(t) = γh(t) = 2 V (t) nach (4.5),
r π
d(t)
γ
a(t) =
= 2 V (t) nach (4.4).
w
r πw
h(t) =
4. MENGENBILANZEN
37
Wir setzen in die Mengenbilanz ein und erhalten
d
γ
(4.8)
V (t) = z − 2 V (t).
dt
r πw
Es ist gut, eine Gleichung mit möglichst wenigen Parametern aufzuschreiben. Fassen wir
die vorkommenden Parameter zusammen
γ
λ= 2 ,
r πw
so erhalten wir die einfachere Form
d
(4.9)
V (t) = z − λV (t).
dt
Diese Gleichung setzt die unbekannte Funktion V (t) in Beziehung zu ihrer Ableitung
d
V (t), sie ist also eine Differentialgleichung. Damit ist die Modellbildung beendet. Was
dt
man mit einer Differentialgleichung anfangen kann, werden wir später sehen.
¤
Merksatz 4.5. Eine dynamische Mengenbilanz beruht auf dem Grundsatz: Die
Änderung der Gesamtmenge ist die Summe aller Zuflüsse abzüglich der Summe aller
Abflüsse. Die Modellbildung besteht darin, eine Funktion f zu erstellen, die die Differenz
der Zu- und Abflussraten Z(t) und A(t) als Funktion der Zeit t und der Gesamtmenge
M (t) angibt:
Z(t) − A(t) = f (t, M (t)).
Änderungsrate bedeutet Änderung pro Zeiteinheit. Ist M eine differenzierbare Funktion
von t, so ist die Änderungsrate von M zum Zeitpunkt t gegeben durch die erste Ableitung
von M (t) nach t, dtd M (t). Dynamische Mengenbilanzen führen daher auf
Differentialgleichungen der Form
d
M (t) = f (t, M (t)).
dt
Merksatz 4.6. Eine Differentialgleichung beschreibt eine Funktion, indem sie sie in
Beziehung zu einer oder mehrerer ihrer Ableitungen setzt. Die Suche nach Lösungen
einer Differentialgleichung ist also die Suche nach Funktionen, die die
Differentialgleichung erfüllen.
Beispiel 4.7. Wir beziehen uns auf das dynamische Modell zu Beispiel 4.2 in
Tabelle 4.2 und die dazugehörige Differentialgleichung (4.9). Sei V0 eine beliebige
Konstante. Zeigen Sie, daß die Funktion
z
z
(4.10)
V (t) = + e−λt (V0 − )
λ
λ
die Differentialgleichung (4.9) löst, und für t = 0 die Anfangsbedingung
V (0) = V0
erfüllt.
Wie lautet die Lösung, wenn zu Beginn (wie in Frage (2) angegeben) der Kessel leer ist?
38
Lösung. Obwohl es meistens schwierig ist, eine Lösung einer Differentialgleichung zu
finden, ist es leicht, nachzuprüfen, ob eine Funktion die Gleichung löst, wenn man sie
einmal “erraten” hat. Man muß dazu nur differenzieren:
d
d hz
z i
V (t) =
+ e−λt (V0 − )
dt
dt λ
λ
z
−λt
= 0 − λe (V0 − )
λ
z
−λt
= −λe (V0 − ).
λ
Zum Vergleich berechnen wir die rechte Seite der Differentialgleichung:
z
z
z − λV (t) = = z − λ − λe−λt (V0 − )
λ
λ
z
−λt
= −λe (V0 − ).
λ
Wir sehen, daß die Ableitung und die rechte Seite der Differentialgleichung
übereinstimmen. Daher erfüllt die Funktion V (t) aus Formel (4.10) die
Differentialgleichung (4.9).
Wir überprüfen noch die Anfangsbedingung:
z
z
V (0) = + e−λ0 (V0 − )
λ
λ
z
z
= + V0 −
λ
λ
= V0 .
Nun soll zum Zeitpunkt t = 0 der Kessel leer, also das enthaltene Wasservolumen gleich
0 sein. Wir setzen in Formel (4.10) also V0 = 0 und erhalten
z
z
z
V (t) = − e−λt = (1 − e−λt ) .
λ
λ
λ
Sie sehen, daß Formel (4.10) viele verschiedene Lösungen beschreibt, für jeden
Anfangswert V0 eine andere. Das versteht man leicht: Das Volumen im Kessel läßt sich
nur dann voraussagen, wenn man nicht nur weiß, was zuströmt und was abfließt (das
wird durch die Differentialgleichung ausgedrückt), sondern auch, was zu Beginn im
Kessel vorhanden ist. Die Bedingung V (0) = 0 ist eine Anfangsbedingung für die
Differentialgleichung: Sie beschreibt den Zustand des Systems “Kessel” zu Anfang des
Beobachtungszeitraumes.
¤
Merksatz 4.8. Ob eine vorgelegte Funktion eine gegebene Differentialgleichung erfüllt
oder nicht erfüllt, kann man durch Einsetzen überprüfen. Eine Differentialgleichung hat
im Allgemeinen unendlich viele Lösungen. Unter diesen wird durch Zusatzbedingungen
(etwa Anfangsbedingungen) eine passende eindeutig ausgewählt.
4. MENGENBILANZEN
39
4.3. Das logistische Populationsmodell als Strukturmodell.
Beispiel 4.9. In einem beschränkten Lebensraum entwickelt sich eine große Population.
Konstruieren Sie ein einfaches mathematisches Modell, das folgenden Gesetzmäßigkeiten
Rechnung trägt:
(1) Die Geburtenrate (Anzahl der Geburten pro Jahr) ist ungefähr proportional zur
Anzahl der Individuen.
(2) Je größer die Population, desto größer ist auch der Anteil der Individuen, die
jährlich sterben: Durch die Erschöpfung der Ressourcen des Lebensraumes
werden die Überlebensmöglichkeiten jedes Individuums eingeschränkt.
(3) Der Lebensraum ist geschlossen, Einwanderung und Auswanderung finden nicht
(oder nur in vernachlässigbar kleinem Ausmaß) statt.
Das Modell soll die Entwicklung der Bevölkerungszahl im Lauf der Zeit beschreiben.
Modellbildung. Wir bezeichnen mit t die Zeit (in Jahren), und mit P (t) die
Bevölkerung zum Zeitpunkt t. P (t) ist die gesuchte Funktion, die durch das Modell
beschrieben werden soll.
Wir formalisieren zunächst das Geburtengesetz. Die Anzahl der Neugeburten pro
Zeiteinheit (Jahr) bezeichnen wir mit B(t). (B ist natürlich auch eine Funktion von t.
Wenn sich im Lauf der Zeit die Gesamtbevölkerungsdichte ändert, wird sich damit auch
die laufende Anzahl der Geburten ändern.) Laut Annahme (1) gibt es einen konstanten
Parameter β, die Fertilität, sodaß
(4.11)
B(t) = βP (t).
Der Faktor β ist die Anzahl der Nachkommen, die ein Individuum pro Jahr im Mittel
hervorbringt.
Das Mortalitätsgesetz ist in Annahme (2) nur sehr vage formuliert. Wir haben
Handlungsfreiheit, eine mathematische Beschreibung zu wählen, die diese Eigenschaften
hat und dabei auf keine allzu komplizierte Formel führt. Für jedes Individuum besteht
eine gewisse Wahrscheinlichkeit, im laufenden Jahr zu sterben, die Mortalität. Wir
bezeichnen sie mit µ(t). (Weil nach Annahme (2) die Überlebenswahrscheinlichkeit von
der Populationsdichte abhängt, ist auch µ keine Konstante, sondern kann sich mit der
Zeit ändern.) Mit M (t) bezeichnen wir die Anzahl der Todesfälle pro Zeiteinheit. Weil
jedes Individuum dem Todesrisiko unterliegt, erhalten wir
(4.12)
M (t) = µ(t)P (t).
Nun soll laut Annahme µ von der Populationsdichte abhängen: Je größer die Population,
desto größer µ. Eine einfache Funktion mit dieser Eigenschaft ist zum Beispiel
(4.13)
µ(t) = γ + δP (t)
mit konstanten Parametern γ und δ. Fassen wir (4.12) und (4.13) zusammen, erhalten wir
(4.14)
M (t) = (γ + δP (t))P (t).
Nach Annahme (3) sind Geburt und Tod die einzigen Vorgänge, die den
Bevölkerungsstand ändern können. Daher sind wir jetzt in der Lage, die dynamische
Mengenbilanz für den Gesamtzuwachs der Population aufzustellen: Das
40
Gesamtwachstum pro Jahr ist die jährliche Anzahl der Geburten abzüglich der jährlichen
Anzahl der Todesfälle:
d
P (t) = B(t) − M (t)
dt
= βP (t) − (γ + δP (t))P (t)
δ
= (β − γ)P (t)(1 −
P (t)).
β−γ
Es ist sinnvoll, ein Modell mit möglichst wenigen Parametern aufzustellen. Wir können
β, γ und δ auf verschiedene Weise anders zusammenfassen, sodaß das Modell nur zwei
Parameter enthält. Zum Beispiel bewährt sich die Einführung der folgenden Parameter:
1
,
β−γ
β−γ
.
K=
δ
τ=
Damit erhalten wir die folgende Differentialgleichung, die sogenannte logistische
Gleichung:
(4.15)
d
1
P (t)
P (t) = P (t)(1 −
).
dt
τ
K
Damit haben wir ein dynamisches Modell für die Population aufgestellt.
¤
Merksatz 4.10. Anders als bei der Modellierung von mechanischen oder anderen
physikalischen oder chemischen Systemen besteht bei der Beschreibung physiologischer,
biologischer und ökonomischer Vorgänge oft große Freiheit in der Wahl der Gleichungen,
weil die zugrundeliegenden Gesetzmäßigkeiten oft nur qualitativ erfaßt sind. Man wählt
dann natürlich möglichst einfache Modelle, die die entsprechenden qualitativen Merkmale
wiedergeben. Dafür muß man auch in vielen Fällen damit rechnen, daß eher die
qualitativen Merkmale der Lösung als die gewonnenen Zahlenwerte für die Beschreibung
des Systems signifikant sind.
Beispiel 4.11. Zeigen Sie: Wenn t0 eine beliebige Konstante ist, ist die Funktion
P (t) = K
e(t−t0 )/τ
1 + e(t−t0 )/τ
eine Lösung der logistischen Gleichung (4.15).
Wie muß t0 bestimmt werden, wenn die Population zu einem Zeitpunkt t1 bekannt ist ?
4. MENGENBILANZEN
41
Tabelle 4.3. Logistisches Populationsmodell
Modellgrößen
Größe
Einheit
Benennung
t
Jahr
Zeit
P (t)
Indiv.
Gesamtpopulation
B(t)
Indiv./Jahr Geburten
M (t) Indiv./Jahr Todesfälle
β
1/Jahr
Fertilität
µ(t)
1/Jahr
Mortalität
γ
1/Jahr
Parameter für µ
δ
1/(Ind. Jahr) Parameter für µ
1
τ
Jahr
Zeitkonstante
β−γ
β−γ
K
Indiv. Kapazität
δ
Kommentar
gesucht
bekannt
bekannt
bekannt
Modellgleichungen
B(t) = βP (t),
M (t) = µ(t)P (t),
µ(t) = γ + δP (t).
dynamische Mengenbilanz
d
1
P (t)
P (t) = B(t) − M (t) = P (t)(1 −
).
dt
τ
K
Lösung. Wir überprüfen die Differentialgleichung: Zuerst differenzieren wir P (t) mit
der Bruchregel:
·
¸
d
d
e(t−t0 )/τ
P (t) =
K
dt
dt
1 + e(t−t0 )/τ
1 (t−t0 )/τ
=Kτ
=
e
(1 + e(t−t0 )/τ ) − e(t−t0 )/τ ( τ1 e(t−t0 )/τ )
(1 + e(t−t0 )/τ )2
e(t−t0 )/τ
K
τ (1 + e(t−t0 )/τ )2
Zum Vergleich berechnen wir die rechte Seite der Differentialgleichung (4.15):
1
P (t)
1
e(t−t0 )/τ
P (t)(1 −
)= K
τ
K
τ
1 + e(t−t0 )/τ
µ
K e(t−t0 )/τ
1−
K 1 + e(t−t0 )/τ
K e(t−t0 )/τ (1 + e(t−t0 )/τ − e(t−t0 )/τ )
τ
(1 + e(t−t0 )/τ )2
K
e(t−t0 )/τ
=
.
τ (1 + e(t−t0 )/τ )2
=
¶
42
Ein Vergleich mit der Ableitung von P (t), die wir oben berechnet haben, zeigt
d
1
P (t)
P (t) = P (t)(1 −
).
dt
τ
K
Also erfüllt die angegebene Funktion P (t) tatsächlich die logistische Gleichung. t0 ist
dabei irgendein konstanter Parameter, je nach Wahl von t0 erhalten wir viele
verschiedene Lösungen.
Wenn jetzt t1 und P1 vorgegeben sind, bestimmen wir den Parameter t0 so, dass
P (t1 ) = P1 ist, und zwar durch Einsetzen von t1 in die Lösungsfunktion P , was P1
ergeben soll:
e(t1 −t0 )/τ
K
= P1
1 + e(t1 −t0 )/τ
Ke(t1 −t0 )/τ = P1 (1 + e(t1 −t0 )/τ )
(K − P1 )e(t1 −t0 )/τ = P1
P1
e(t1 −t0 )/τ =
K − P1
µ
¶
P1
t1 − t0
= ln
,
τ
K − P1
¶
µ
P1
.
t0 = t1 − τ ln
K − P1
Wieder läßt sich der freie Parameter mit Hilfe einer Zusatzbedingung bestimmen.
¤
Dieselbe Funktion P (t) haben wir im vorigen Kapitel sehr erfolgreich als empirisches
Modell an Felddaten für die Population der USA angepaßt. Wir sehen jetzt, daß es ein
Strukturmodell gibt, das gerade auf diese Funktion führt. Gemessen an den groben
Verallgemeinerungen, die zum Strukturmodell geführt haben, ist die hervorragende
Übereinstimmung zwischen Daten und Modell trotzdem überraschend.
Bemerkung. Genaugenommen ist die Anzahl der Individuen einer Population eine
ganze Zahl. Diese ändert sich unregelmäßig und sprunghaft zu den Zeitpunkten von
Geburt oder Tod einzelner Individuen. Also beruht das Modell mit nicht-ganzzahliger
sich stetig ändernder Bevölkerungszahl P (t) auf einer Idealisierung, die aber die
Modellbildung wesentlich vereinfacht. Für viele Zwecke ist das stetige Modell ausserdem
völlig ausreichend; eine ganzzahlige Modellierung ist zwar grundsätzlich möglich (siehe
Teil 2), für große Populationen aber nicht angebracht.
4. MENGENBILANZEN
43
4.4. Das Punktreaktormodell.
Beispiel 4.12. Die Vorgänge in einem Kernspaltungsreaktor lassen sich am besten im
Hinblick auf die freien Neutronen im Reaktor beschreiben. Das Schicksal der Neutronen
wird von den folgenden Vorgängen geprägt:
(1) Freie Neutronen können von Atomkernen eingefangen werden, ohne daß eine
weitere Reaktion eintritt.
(2) Freie Neutronen können von spaltbaren Atomkernen eingefangen werden, wobei
eine Kernspaltung ausgelöst wird. Bei der Kernspaltung werden wiederum
Neutronen frei (die sogenannten prompten Neutronen), außerdem entstehen
instabile radioaktive Isotope.
(3) Wenn gewisse instabile Isotope zerfallen, werden Neutronen frei. (Diese
Neutronen heißen verzögerte Neutronen, weil sie indirekt durch die
Kernspaltung, aber nicht sofort, sondern über den Umweg der entstandenen
Isotope etwas später freigesetzt werden.)
(4) Abgesehen von diesen Reaktionen gibt es eine Neutronenquelle, die —
unabhängig von den laufenden Kernspaltungen — ständig eine gewisse kleine
Menge an Neutronen freisetzt.
Die Aktivität eines Reaktors läßt sich durch die Anzahl der freien Neutronen im Reaktor
beschreiben (wegen der oben beschriebenen Beziehungen ist diese Zahl stark zur Anzahl
der ablaufenden Kernspaltungen korreliert). Entwickeln Sie ein einfaches dynamisches
Modell zur Vorhersage der Neutronenpopulation im Reaktor.
Modellbildung. Wir bezeichnen mit t die Zeit (in Sekunden), und mit N (t) die
Anzahl der freien Neutronen im Reaktor zum Zeitpunkt t. Außerdem müssen wir den
Vorrat an instabilen Isotopen im Auge behalten. Der Einfachheit halber nehmen wir
(nicht realistisch) an, daß in der Reaktionskette nur eine Sorte instabile Isotope
vorkommt. Sei I(t) die Anzahl der instabilen Atomkerne.
Betrachten wir jetzt den Reaktor vom Standpunkt eines freien Neutrons aus: Die
Wahrscheinlichkeit, daß das Neutron in der nächsten Sekunde eingefangen wird, ohne
eine weitere Reaktion auszulösen, sei β. Die Wahrscheinlichkeit, daß das Neutron in der
nächsten Sekunde eingefangen wird und eine Kernspaltung auslöst, sei α. Da jedes
Neutron diesen Risiken ausgesetzt ist, haben wir die folgende Verlustbilanz für die
Neutronen:
(4.16)
(4.17)
αN (t) + βN (t) Neutronen werden pro Sekunde eingefangen.
αN (t) Kernspaltungen finden pro Sekunde statt.
Sei ν die Anzahl der Neutronen und γ die Anzahl der instabilen Kerne, die bei einer
Kernspaltung freigelassen werden. Aus (4.17) folgt dann:
(4.18)
ναN (t) prompte Neutronen werden pro Sekunde frei.
(4.19)
γαN (t) instabile Kerne werden pro Sekunde frei.
Mit λ bezeichnen wir die Wahrscheinlichkeit, die für einen instabilen Kern besteht,
innerhalb der nächsten Sekunde zu zerfallen. Bei jedem Zerfall entstehen µ verzögerte
44
Tabelle 4.4. Punktreaktormodell mit einer instabilen Spezies
Modellgrößen
Größe Einheit Benennung
t
s
Zeit
N (t)
Ind.
freie Neutronen
I(t)
Ind.
instabile Kerne
β
1/s
Wahrsch.: Einfangen ohne Reaktion
α
1/s
Wahrsch.: Einfangen mit Reaktion
λ
1/s
Wahrsch.: radioaktiver Zerfall
γ
instabile Kerne pro Spaltung
ν
prompte Neutronen pro Spaltung
µ
verzögerte Neutronen pro Zerfall
S
Ind./s Neutronen aus der Quelle
Kommentar
gesucht
gesucht
bekannt
bekannt
bekannt
bekannt
bekannt
bekannt
bekannt
dynamische Mengenbilanzen
d
N (t) = (να − α − β)N (t) + µλI(t) + S,
dt
d
I(t) = γαN (t) − λI(t).
dt
Neutronen. Weil jeder instabile Kern diesem Risiko unterliegt, erhalten wir:
(4.20)
(4.21)
λI(t) instabile Kerne zerfallen pro Sekunde.
µλI(t) verzögerte Neutronen werden pro Sekunde frei.
Letztich benennen wir den Beitrag der Neutronenquelle
(4.22)
S Neutronen pro Sekunde werden von der Quelle freigesetzt.
Wir stellen jetzt zwei Mengenbilanzen auf, eine für die freien Neutronen, eine für die
instabilen Kerne. Der Gesamtzuwachs der freien Neutronen pro Sekunde besteht aus den
freiwerdenden prompten Neutronen und den verzögerten Neutronen, der Verlust aus den
eingefangenen Neutronen. Der Gewinn an instabilen Atomkernen kommt von den
Kernspaltungen, der Verlust vom radioaktiven Zerfall. Wir erhalten daher die folgenden
Bilanzen:
d
(4.23)
N (t) = (να − α − β)N (t) + µλI(t) + S,
dt
d
(4.24)
I(t) = γαN (t) − λI(t).
dt
Das ist ein System von zwei Differentialgleichungen für zwei unbekannte Funktionen N
und I.
Es ist nicht schwierig, das Modell so auszuweiten, daß mehrere verschiedene Isotope mit
verschiedenen Halbwertszeiten im Reaktor beschrieben werden können, und auch bei
Zerfall instabiler Isotope wieder andere instabile Kerne freiwerden, sodaß in der Folge
einer Kernspaltung eine ganze Kette radioaktiver Zerfallsreaktionen abläuft.
¤
4. MENGENBILANZEN
45
Weil wir in diesem Modell nicht die räumlichen Verhältnisse im Reaktor berücksichtigt haben
(etwa, wie weit ein Neutron fliegen muß, bis es am nächsten Brennelement eingefangen werden
kann, um eine Kernspaltung auszulösen), sondern gerechnet haben, als würden alle Vorgänge
an einem Punkt ablaufen, heißt diese Art von Neutronenbilanz ein Punktreaktormodell.
Für das technische Design eines speziellen Reaktors, das ja darin besteht, Brennelemente,
Bremsstäbe, und viele andere Einrichtungen möglichst wirkungsvoll und betriebssicher im
Raum anzuordnen, ist das Punktreaktormodell natürlich viel zu grob. Man kann aber an Hand
des Modells Einiges über die Eigenheiten von Reaktoren lernen, etwa, daß diese Einrichtungen
überhaupt nur technisch möglich sind, weil ein Teil der Neutronenbilanz relativ langsam über
die verzögerten Neutronen abläuft, die prompten Neutronen reagieren viel zu schnell, als daß
man den Reaktor rechtzeitig steuern könnte. Andere Eigenheiten sind zum Beispiel, daß Start
und Abschaltung, und vor allem ein Neustart knapp nach einer Abschaltung, besonders
kritische Betriebsphasen sind.
Beispiel 4.13. Wir betrachten das Punktreaktormodell aus Beispiel 4.12. Nehmen wir
an, daß der Reaktor stabil ist und sich ein Gleichgewicht einspielt. Welche
Neutronenpopulation wird sich im Reaktor als Gleichgewicht einstellen?
Lösung. Das ist eine statische Mengenbilanz. Wir können also im
Punktreaktormodell die Zuflüsse und Abgänge von Neutronen und instabilen Isotopen
jeweils gleichsetzen. Anders ausgedrückt, müssen die Gesamtzuwachsraten von Neutronen
und Isotopen gleich Null sein: In den Differentialgleichungen (4.23) und (4.24) setzen wir
d
d
N = 0,
I = 0.
dt
dt
Damit erhalten wir ein lineares Gleichungssystem
0 = (να − α − β)N (t) + µλI(t) + S,
0 =
γαN (t)
− λI(t).
Das System besitzt die Lösung
S
,
α + β − γµα − να
γαS
I=
.
λ(α + β − γµα − να)
N=
N und I hängen nun natürlich nicht mehr von t ab: Das Gleichgewicht besteht ja gerade
darin, daß sich der Zustand des Systems mit der Zeit nicht ändert.
¤
Auf analoge Art und Weise können Sie die Gleichgewichtslagen des dynamischen
Kessel-Modells oder des Populationsmodells bestimmen.
Merksatz 4.14. Die Gleichgewichtslagen eines Systems von Differentialgleichungen
erhält man, indem man die Ableitungen mit Null ansetzt und das entstehende
Gleichungssystem löst. Das Null-Setzen der Ableitungen macht aus dem dynamischen
Modell ein statisches. Falls dieses Lösungen besitzt, sind diese die Gleichgewichtslagen
des Differentialgleichungsmodells.
46
Viel wichtiger ist allerdings die Frage, ob der Reaktor überhaupt stabil ist. Die Antwort
darauf liefert die Theorie der linearen Differentialgleichungen mit Hilfe von
Matrizenrechnung. Stabilitätsuntersuchungen gehören zur qualitativen Analyse
mathematischer Modelle und werden in einem späteren Kapitel besprochen.
4.5. Kinetik chemischer Reaktionen.
Beispiel 4.15. Stickoxid (NO) spielt nach neuesten Erkenntnissen eine wesentliche
Rolle in der Blutdruckregulation. Pharmakologen haben Versuche zu folgender Reaktion
durchgeführt: Einer Lösung wird eine gewisse Menge eines Stoffes (“Donor”) beigesetzt,
aus dem durch Zerfall NO entsteht. NO selbst wird, zusammen mit Sauerstoff, weiter
abgebaut. Um konstante Reaktionsgeschwindigkeiten zu sichern, wurde die Temperatur
der Lösung konstant gehalten. Die NO-Konzentration in der Lösung wurde etwa 40
Minuten lang durch ständige Messungen verfolgt und die Daten mittels Computer
aufgenommen. Finden Sie ein Modell, das diese Daten wiedergeben kann, unter
Berücksichtigung der folgenden chemischen Sachverhalte:
(1) Der Zerfall des Donors ist eine Reaktion erster Ordnung, das heißt, in den
Zerfall ist jeweils nur ein Molekül des Donors involviert. Die Anzahl der pro
Sekunde zerfallenden Moleküle ist proportional zur Konzentration des Donors in
der Lösung. Die Reaktion ist irreversibel.
(2) Der Zerfall von NO ist, von NO aus gesehen, eine Reaktion zweiter Ordnung. Je
zwei Moleküle von NO reagieren, zusammen mit einem Molekül O2 , in der
Abbaureaktion. Die Anzahl der pro Sekunde zerfallenden Moleküle ist
proportional zum Quadrat der NO-Konzentration und proportional zur
O2 -Konzentration. Die Reaktion ist irreversibel.
(3) Sauerstoff ist in der Lösung so reichlich vorhanden, daß die Konzentration von
O2 trotz der Zerfallsreaktion als annähernd konstant angenommen werden darf.
(4) Außer den beiden angeführten Reaktionen und der Beigabe von Donor zur
Lösung am Anfang des betrachteten Zeitraumes finden keine Reaktionen statt,
die die Konzentration von NO beeinflussen.
(5) Vor der Beigabe des Donors befindet sich kein NO in der Lösung.
Modellbildung. Die Problemstellung verlangt ein dynamisches Modell. t
bezeichne die Zeit, die seit Versuchsbeginn verstrichen ist (in Sekunden). Wir bezeichnen
mit cN O (t) die Konzentration von NO und mit cD (t) die Konzentration des Donors zum
Zeitpunkt t. Die Einheiten sind Mol pro Liter. (Ein Mol ist eine Gewichtseinheit, die sich
auf das Molekulargewicht der betreffenden Substanz bezieht: Die Masse eines Stoffes in
Mol ist die Gesamtmasse in Gramm, gebrochen durch das Molekulargewicht. Chemische
Berechnungen vereinfachen sich durch diese Konvention: Zum Beispiel liefert der Zerfall
eines Mols vom Donor genau ein Mol von NO.) c0 bezeichne die Donorkonzentration in
der Lösung unmittelbar nach Zugabe des Donors, also zur Zeit t = 0 (in Mol/l). Anfangs
befindet sich in der Lösung noch gar kein Stickoxid.
(4.25)
(4.26)
cD (0) = c0 ,
cN O (0) = 0.
4. MENGENBILANZEN
47
Mit o2 bezeichnen wir die Konzentration an Sauerstoff. Nach Annahme (3) können wir o2
als konstant voraussetzen. Die beiden Reaktionen laufen dann nach den folgenden
Gesetzmäßigkeiten ab:
Die Masse NO, die in einem Liter Lösung wird pro Sekunde durch Donorzerfall
freigesetzt wird, ist proportional zur Donorkonzentration. In Mol gerechnet, wird
gleichzeitig der Donorbestand um den gleichen Betrag verringert. Der
Proportionalitätsfaktor k1 (in 1/s) beschreibt die Reaktionsgeschwindigkeit des Zerfalls:
(4.27)
k1 cD (t) = Menge zerfallender D = Menge gebildetes NO pro Sek.
Die Masse NO, die in einem Liter Lösung abgebaut wird, ist proportional zum Quadrat
der NO-Konzentration und zur Sauerstoffkonzentration:
(4.28)
k2 o2 c2N O (t) = Menge abgebautes NO pro Sek.
Die Mengenbilanz für den Donor enthält nur Verluste durch den Zerfall. Die
Mengenbilanz des Stickoxids enthält Zuwachs vom Donorzerfall und Verlust vom
NO-Abbau.
d
(4.29)
cD (t) = −k1 cD (t),
dt
d
(4.30)
cN O (t) = k1 cD (t) − k2 o2 cN O (t)2 .
dt
Wir erhalten ein System von 2 Differentialgleichungen in den zwei unbekannten
Funktionen cD und cN O . Obwohl sich dieses System mit Bleistift und Papier nicht durch
eine geschlossene Formel lösen läßt, kann die Lösung mit Hilfe numerischer Verfahren auf
dem Computer erfolgen, wenn man die Parameter k1 , k2 , c0 und o2 kennt.
¤
Die einzige Gleichgewichtslage ist die Null-Lage, denn aus ċD = 0 folgt mit k1 6= 0 dass
cD (t) ≡ 0 und dann aus ċN O = 0 dass cN O (t) ≡ 0 (ein Punkt über einer Funktion ist eine
Kurzschreibweise für die Ableitung dieser Funktion). Die Gleichgewichtslage des Modells
wird ausgehend von jedem positiven Anfangszustand asymptotisch im Lauf der Zeit
angenähert, aber nie ganz erreicht.
Bemerkung: B. M. Mayer und K. Schmidt am Institut für Pharmakologie und Toxikologie,
K. F. Universität Graz, haben diese Versuche durchgeführt, und das Modell dazu aufgestellt.
Das Modell war bei der Validierung der Versuche nützlich: Die anfängliche Diskrepanz zwischen
Modell und Daten hat die Aufmerksamkeit darauf gelenkt, daß bei diesem Prozeß der Einfluß
der Temperatur auf die Reaktionsgeschwindigkeit sehr groß ist. Nach Modifikation der Versuche
mit genauer regulierter Temperatur ist die Übereinstimmung zwischen Modell und Meßdaten
bei geeigneter Wahl der Parameter (Reaktionsgeschwindigkeiten) sehr überzeugend.
Diese Vorgehensweise ist eine Standardmethode zur Modellierung der Kinetik chemischer
Reaktionen. Geht man vom dynamischen Modell dann zur statischen Mengenbilanz über,
erhält man für reversible Reaktionen das bekannte Massenwirkungsgesetz.
48
Tabelle 4.5. Stickoxidbildung durch Zerfall einer Donorsubstanz
Größe
Einheit
t
s
cD (t)
Mol/l
cN O (t)
Mol/l
c0
5.6 · 10−6 Mol/l
o2
2 · 10−4 Mol/l
k1
1.4 · 10−3 1/s
k2
9 · 106 l2 /(Mol2 s)
Modellgrößen
Benennung
Zeit
Konzentration Donor
Konzentration NO
Anfangskonzentration Donor
Konzentration O2
Reaktionskonstante Donorzerfall
Reaktionskonstante NO-Abbau
Kommentar
gesucht
gesucht
bekannt
bekannt
bekannt
bekannt
dynamische Mengenbilanzen
d
cD (t) = −k1 cD (t),
dt
d
cN O (t) = k1 cD (t) − k2 o2 cN O (t)2 .
dt
Anfangsbedingungen
cD (0) = c0 ,
cN O (0) = 0.
4.6. Ein kleines Weltmodell.
1972 wurden vom Club of Rome Simulationsergebnisse veröffentlicht, die etwa für das
Jahr 2000 einen dramatischen Bevölkerungszusammenbruch voraussagten. Diese
Prognosen beruhten auf mathematischen Modellen der globalen Entwicklung, wie sie
damals von J. W. Forrester, D. Meadows u. a. entwickelt wurden. Sie werden oft kurz als
Weltmodelle bezeichnet. Wir betrachten hier eine einfache Variante, die sich aus der
Bilanz des Zusammenspiels der drei Kenngrößen Bevölkerung, Wirtschaft und Umweltund Ressourcenbelastung ergibt. Ausgangspunkt ist folgende grobe Beschreibung der
”Weltdynamik”(nach H. Bossel, Modellbildung und Simulation, Vieweg 1994):
Beispiel 4.16. Wir beobachten heute weltweit eine zunehmende Belastung der
natürlichen Ressourcen (Rohstoffreserven) und der natürlichen Umwelt. Die Gründe
hierfür liegen in einer ständigen Zunahme der Bevölkerung, damit auch der Verbräuche
der verschiedensten Rohstoffe und der Abgabe von Abfallstoffen jeder Art an die Umwelt.
Eine wichtige Bestimmungsgröße dieser Ressourcen- und Umweltbelastung ist der
spezifische Verbrauch an Rohstoffen und Energie pro Kopf. Dieser spezifische Verbrauch
steigt noch tendenziell mit der wachsenden Umweltbelastung (durch wachsende
Aufwendungen für den Umweltschutz und schwieriger werdende Abbaubedingungen).
Andererseits ist festzustellen, dass sich mit wachsendem Einsatz von Rohstoffen und
Energien pro Kopf, also mit wachsendem spezifischen Konsum, zunächst die
Versorgungsmöglichkeiten verbessern, so dass mehr Menschen versorgt werden können.
Festzustellen ist aber auch, dass auf Grund der wachsenden Umweltbelastungen mit
Schadstoffen wie auch der schwindenden natürlichen Ressourcenbasis Rückwirkungen auf
4. MENGENBILANZEN
49
die Gesundheit und die Lebenserwartung der Bevölkerung bestehen. Die
Umweltbelastungen lassen ein zunehmendes politisches Handeln erwarten, um
schädlichen Entwicklungen zu begegnen.
Dieser Beschreibung folgend (wie immer man dazu stehen mag), ergeben sich drei für
den Zustand und die zeitliche Entwicklung maßgebliche Größen, nämlich die
Bevölkerung, der spezifische Konsum (Kurzbegriff stellvertretend für den gesamten
Verbrauch an natürlichen Ressourcen pro Kopf) und die Belastung der Umwelt.
Dementsprechend definieren wir die drei Zustandsgrößen
V Größe der Bevölkerung
L Umweltbelastung
K spezifischer Konsum
und entnehmen der obigen Beschreibung dynamische Mengenbilanzen für diese Größen.
Der Kürze halber werden hier Überlegungen bzgl. adäquater Maßeinheiten beiseite
gelassen.
Die Zunahme der Bevölkerung ist proportional zur Größe der Bevölkerung, die
Zuwachsrate wird mit wachsendem Konsum größer und mit wachsender
Umweltbelastung kleiner. Eine Formel mit dieser Eigenschaft ist
K
Zuwachsrate von V = β V.
L
Die Abnahme der Bevölkerung pro Zeiteinheit ist proportional zu V und soll sich bei
steigender Umweltbelastung verstärken; wir setzen
Abnahmerate von V = µLV.
Darin sollen β und µ positive Konstanten sein. Die Umweltbelastung soll mit wachsender
Bevölkerung steigen, ebenso mit wachsendem Konsum, was z. B. in der Formel
Zuwachsrate von L = λV K
mit λ > 0 der Fall ist. Um zum Ausdruck zu bringen, dass die Umweltbelastung durch
natürliche Prozesse – wie den Zerfall von Schadstoffen oder die CO2 -Aufnahme der
Pflanzen – vermindert wird, setzen wir die Formel einer chemischen Reaktion erster
Ordnung an, nämlich
Abnahmerate von L = αL,
mit einer positiven Konstanten α. Ökologische Abbauprozesse haben Kapazitätsgrenzen,
die z. B. durch Nährstoff-, Licht- und Wasserbeschränkungen gegeben sind. Bei
Überlastung kann der Abbau bestenfalls an dieser Grenze operieren; es gibt daher einen
kritischen Wert L0 so, dass oberhalb dieses Wertes der Abbau bei einem konstanten
Wert hängt, also
Abnahmerate von L = αL0 , für L > L0 .
Schließlich suchen wir eine Formel für die Änderung des Konsums. Ungehemmtes
Wirtschaftswachstum wird durch exponentielles Wachstum K̇ = κK modelliert.
Realistischer ist es, eine Kapazitätsgrenze K0 zu modellieren; wie im Populationsbeispiel
ist dies z. B. mit der logistischen Bilanz K̇ = κK(K0 − K) der Fall. Um zu
berücksichtigen, dass laut obiger Beschreibung zunehmende Umweltbelastung den
50
Verbrauch zunächst erhöht, andererseits bei zunehmender Belastung der Verbrauch
durch politisches Handeln gedrosselt wird, ersetzt man auf der rechten Seite der
logistischen Bilanz K durch das Produkt KL, und erhält
Änderungsrate von K = κKL(σ − KL)
mit positiven Konstanten κ, σ. Zusammengefasst ergibt sich das folgende mathematische
Modell für die zeitliche Änderung der Größen V (t), L(t), K(t) zur Zeit t:
(4.31)
(4.32)
(4.33)
d
K(t)
V (t) = β
V (t) − µL(t)V (t),
dt
L(t)
½
d
αL(t), für L(t) ≤ L0
L(t) = λV (t)K(t) −
αL0 , für L(t) > L0
dt
¡
¢
d
K(t) = κK(t)L(t) σ − K(t)L(t) .
dt
Dies ist ein System von 3 nichtlinearen Differentialgleichungen, das man für numerische
Simulationen verwenden kann. Die Dynamik der Lösungen hängt stark von den
(relativen) Werten der Parameter β, µ, λ, α, L0 , κ und σ ab. Es muss versucht werden,
wenigstens einige der Parameter in vom Modell unabhängigen
wirtschaftswissenschaftlichen Untersuchungen zu bestimmen. Einige der Parameter sind
durch gesellschaftliches und politisches Handeln beeinflussbar, etwa β
(Geburtenkontrolle), λ (Umweltschutz) oder σ (Drosselung des Konsums). Daraus erklärt
sich der Appell-Charakter des Einsatzes von Weltmodellen: diese oder jene Entwicklung
wird vorausgesagt, wenn nicht rechtzeitig diese oder jene Maßnahmen greifen.
Das kleine Modell ist zwar zu simpel, um ernstzunehmende Prognosen liefern zu können,
es wird hier nur zur Demonstration von Mengenbilanzen herangezogen. Trotz seiner
Kleinheit zeigt das Modell jedoch die auch für ausgebaute Weltmodelle typischen
Verhaltensmuster: gedämpfte Schwingungen – oder aber Aufschaukelung und
Zusammenbruch – je nach Parametersatz.
An den Beispielen dieses Kapitels sieht man die Allgemeinheit und breite Anwendbarkeit
der Methode der dynamischen Mengenbilanz; sie führt zu Modellen mit
Differentialgleichungen, mit denen die Änderung der Mengen beschrieben wird. Das Wort
’Mengen’ wird hier bewusst in seiner ungenauen Bedeutung verwendet, z. B. kann
Umweltbelastung sowohl Schadstoffe in Gewässern als auch Defekte der
Atmosphärenzusammensetzung umfassen.
Später wird gezeigt, dass man am Computer die Lösungen (eines Systems) von
Differentialgleichungen recht einfach näherungsweise berechnen und graphisch darstellen
kann. Zusätzlich zum Verständnisgewinn bei der Modellbildung selbst, hat man also
Möglichkeiten durch (systematisches) Probieren und Simulieren qualitative und
quantitative Eigenschaften des Modells zu untersuchen.
5. DIMENSIONEN UND EINHEITEN
51
5. Dimensionen und Einheiten
5.1. Das SI-System.
Jeder Meßgröße kommt eine Maßeinheit zu, etwa Meter (für Längen), Kubikmeter
(Volumen), Newton (Kraft), Watt (Leistung) u.a.. Bekanntlich gibt es verschiedene
Einheiten für dieselbe Größe, man kann etwa Leistung in Watt oder PS angeben. Das
SI-System (système international) erzielt eine Vereinheitlichung, indem alle Maßeinheiten
auf einige wenige Grundeinheiten zurückgeführt werden, und zwar
Zeit
Länge
Masse
el. Stromstärke
Temperatur
Stoffmenge
Lichtstärke
s
m
kg
A
K
mol
cd
Sekunde
Meter
Kilogramm
Ampere
Kelvin
Mol
Candela
Wir geben einige abgeleitete Maßeinheiten der Mechanik an:
Frequenz
Geschwindigkeit
Beschleunigung
Kraft
Druck
Arbeit
Leistung
Hz s−1
m s−1
m s−2
N kg m s−2
Pa N m−2 = kg m−1 s−2
J
N m = kg m2 s−2
W N m s−1 = kg m2 s−3
Hertz
Newton
Pascal
Joule
Watt
Merksatz 5.1. Achten Sie darauf, daß die Zahlenangaben in Ihren Modellen auf ein
einheitliches Maßsystem bezogen sind. Ständige Umrechnungen zwischen
Einheitensystemen machen die Arbeit undurchsichtig und sind Fehlerquellen.
5.2. Dimensionen.
Jede Maßeinheit gehört zu einer Dimension, das ist die Art der Größe, die durch die
Einheit beschrieben wird. So entspricht das Meter (ebenso wie die Meile, Seemeile oder
ein Angström) der Dimension Länge, ein Kilogramm, Pfund, oder eine Tonne gehören
zur Dimension Masse. In den beiden SI Tabellen stehen in der ersten Spalte Dimensionen
und in der letzten Spalte die zugehörigen SI-Einheiten. In korrekt formulierten Modellen
unterliegen dimensionierte Zahlengrößen bestimmten Regeln:
Merksatz 5.2. Die folgenden Vorschriften regeln den Gebrauch von dimensionierten
Größen in Formeln:
(1) Man darf nur Größen derselben Dimension gleichsetzen.
(2) Man darf nur Größen derselben Dimension addieren oder subtrahieren. Die
Summe (Differenz) hat dann wieder dieselbe Dimension.
52
(3) Multipliziert (dividiert) man eine Größe der Dimension A mit einer Größe der
Dimension B, so hat das Produkt die Dimension AB (der Quotient die
Dimension A/B).
(4) Differenziert man eine Größe der Dimension A nach einer Größe der
Dimension B, so hat die Ableitung die Dimension A/B. Integriert man eine
Größe der Dimension A nach einer Größe der Dimension B, so hat das Integral
die Dimension AB.
(5) In alle anderen Funktionen (Winkelfunktionen, Exponentialfunktion,
Logarithmus usw.) dürfen nur dimensionslose Größen eingesetzt werden. Die
Resultate sind wieder dimensionslose Größen.
Ein System von Formeln, das alle diese Vorschriften erfüllt, nennen wir homogen in der
Dimension.
Beispiel 5.3. Welche Dimensionen und Einheiten haben die Größen im Modell zu
Beispiel 4.15?
Lösung. Die Größe t bezeichnet Zeit. Als Einheit kommt nach dem SI-System die
Sekunde in Frage.
cD und cN O sind Konzentrationen, also Teilchenzahlen pro Volumen. Als Einheit käme
nach dem SI-System mol/m3 in Frage. Wir verwenden im Modell mol/l, das sind 1000
mol/m3 , denn 1 m3 umfaßt 1000 l. (Zwar ist Liter keine Grundeinheit des SI-Systems.
Weil aber in diesem Modell weder Längen noch Flächen vorkommen, bringt es nichts,
das Volumen in die dritte Potenz einer Länge zu zerlegen.)
Die Einheiten der Modellparameter ergeben sich daraus, daß alle obigen Regeln erfüllt
sein müssen:
Wir untersuchen die Gleichung
d
cD (t) = −k1 cD (t)
dt
auf Dimension. Dabei lassen wir die in Klammern stehenden t ausser Betracht, diese
bedeuten ja nur “zum Zeitpunkt t”. dtd cD ist die Ableitung einer Konzentration nach der
Zeit und hat nach Regel (4) als Dimension Konzentration/Zeit, d.i. Masse/(Volumen ·
Zeit). Daher muß auch k1 cD dieselbe Dimension haben. Weil cD eine Konzentration ist,
bleibt nach Regel (3) für den Faktor k1 die Dimension 1/Zeit, die Einheit 1/s.
Ebenso untersuchen wir die Gleichung
d
cN O (t) = k1 cD (t) − k2 o2 cN O (t)2 .
dt
Wieder hat dtd cN O die Dimension Konzentration/Zeit. Die beiden Summanden auf der
rechten Seite müssen dieselbe Dimension haben. Da k1 als Dimension 1/Zeit hat und cN O
eine Konzentration ist, paßt der erste Summand bereits. cN O und o2 sind
Konzentrationen. Daher hat o2 c2N O als Dimension Konzentration3 . Da k2 o2 cN O als
Dimension Konzentration/Zeit hat, bleibt für den Faktor k2 die Dimension
1/(Zeit.Konzentration2 ). Als Einheit ergibt sich l2 /(mol2 s).
¤
5. DIMENSIONEN UND EINHEITEN
53
Zur Bestimmung von Dimensionen ist es zweckmäßig mit fest gewählten Einheiten zu
arbeiten. Wir verwenden eckige Klammern, um die Einheit einer Größe zu bezeichnen,
[Größe] = algebraische Kombination der Grundeinheiten.
Wenn g im Modell z. B. eine Beschleunigung bezeichnet, schreiben wir
m
[g] = 2 .
s
Man “rechnet” dann mit den Abkürzungen für die Einheiten, als ob sie algebraische
Variable wären. In den zusammenfassenden Tabellen zu den in früheren Kapiteln
gebildeten Modellen finden Sie die Einheiten der Modellgrößen. Mit deren Hifle können
Sie die Modelle hinsichtlich der Dimensionsregeln überprüfen.
Merksatz 5.4. Die Überprüfung von Modellgleichungen auf Dimensionen ist eine gute
Fehlerprobe, ob die Gleichungen wenigstens formal richtig sind.
Wenn man sich zum Beispiel nicht sicher ist, ob das Gesetz von Torricelli für die
Ausflussgeschwindigkeit v aus einem
√ Behälter mit Spiegelhöhe h aufgrund der
Erdbeschleunigung g wirklich v = 2gh lautet, kann man zur Kontrolle eine
Dimensionsbetrachtung durchführen
[v 2 ] =
m2
s2
und [2gh] = [2][g][h] =
m
m,
s2
also ist die Formel v 2 = 2gh tatsächlich homogen in der Dimension. Die Zahl 2 ist
dimensionslos und wird daher einfach weggelassen. Übrigens ergibt sich der
Zusammenhang mit dem Kessel-Modell a = d/w aus Tabelle 4.1 mit a = Av, γ = ρg,
w = ρv/2A, mit der Dichte der Flüssigkeit ρ und dem Querschnitt des Abflussrohrs A.
Prüfen Sie dies nach, insbesondere auch die Dimensionen (die Dichte ist die Masse einer
Substanz pro Einheitsvolumen, [ρ] = kg/m3 ).
Die Überprüfung (vorläufiger) Modelle auf Dimensionsregeln mag auf den ersten Blick
wie nutzlose Pedanterie erscheinen, in der Praxis der Modellbildung ist sie aber ein
äusserst wirksames und einfach einsetzbares Werkzeug. Dabei darf allerdings nicht der
Versuchung nachgegeben werden, zusätzliche (und überflüssige) Modell-Parameter
einzuführen, nur zum Zweck des Zurechtbiegens der Dimensionen.
Beispiel 5.5. Der Gebrauchsanweisung eines Fallschirms entnehmen Sie, dass sein
Widerstandsbeiwert ζ in Quadratmetern angegeben ist; insbesondere schließen Sie
daraus, dass [ζ]=m2 . Schätzen Sie ab, wie schwere Lasten der Masse m Sie damit
abwerfen können, wenn die Fallgeschwindigkeit einen vorgegebenen Wert v nicht
überschreiten darf.
54
Lösung. Wir stellen zunächst eine statische Kräftebilanz Gewicht= Reibung auf,
also mg = R mit der Erdbeschleunigung g. Ein plausibles Modell für die Reibung
aufgrund des Luftwiderstands des Fallschirms sollte sowohl mit ζ und v, als auch mit der
Atmosphärendichte ρ zunehmen. Der einfachste Ansatz mit dieser Eigenschaft ist das
Produkt R = ζρv. Die entsprechende Kräftebilanz lautet mg = ζρv. Wir bestimmen die
Dimensionen dieses Modells:
m
kg m
kg
.
[mg] = kg 2 und [ζρv] = [ζ][ρ][v] = m2 3 =
s
m s
s
Links und rechts vom Gleichheitszeichen stehen also verschiedene Dimensionen! Daher
müssen wir dieses Modell verwerfen.
Wir sehen, dass in der Dimension von R im Zähler noch eine Länge und im Nenner eine
Zeit fehlt, insgesamt also ein Faktor Geschwindigkeit. Dies bringt uns auf die Idee den
Ansatz R = ζρv 2 zu versuchen. Er ist ebenso plausibel und hält sogar der
Dimensionsbetrachtung stand:
kg m2
m kg
[R] = m 3 2 = 2 = [mg].
m s
s
Die Formel R = ζρv 2 ist somit ein brauchbarer Kanditat für R. Dieses Modell liefert für
die maximal abwerfbare Masse den Wert m = ζρv 2 /g. Es kann aber ohne genauere
Angaben nicht eindeutig geklärt werden, ob in der Gebrauchsanweisung wirklich diese
Formel gemeint ist, oder ob in ζ vielleicht noch dimensionslose Faktoren versteckt sind.
Zum Beispiel werden in Gleichungen für den Druck in Rohrströmungen Reibungsterme
in der Form 12 ζρv 2 angesetzt mit einem dimensionslosen Widerstandsbeiwert ζ: der
Faktor 12 ist eine altehrwürdige Konvention und erinnert an die kinetische Energie der
Strömung.
¤
2
5.3. Dimensionslose Modelle.
Das Stickoxid-Modell enthält 4 Modellparameter (k1 , k2 , c0 und o2 ). Wir sehen aber, daß
manche der Parameter nur gemeinsam vorkommen. Zum Beispiel treten k2 und o2 nur
gemeinsam als Produkt k2 o2 auf. Wir könnten sehr gut die beiden Parameter in einen
einzigen zusammenfassen. Wenn man sich eine Übersicht über das Verhalten des Modells
für verschiedene Parameterwerte verschaffen will (und das kann eine sehr wichtige
Fragestellung sein), ist es ein Vorteil, wenn in den Modellgleichungen möglichst wenige
Parameter vorkommen. Eine systematische Methode, ein Modell auf möglichst wenige
Parameter zu reduzieren, besteht darin, es in ein Modell mit dimensionslosen Größen
umzuformen.
Es gibt in sinnvollen Modellen immer dimensionslose Produkte, das sind Kombinationen
von Modellgrößen derart, dass sich die Dimensionen durch die Kombination wegkürzen.
Man spricht dann von einer dimensionslosen Größe und schreibt für die zugehörige
Einheit 1.
Beispiel 5.6. Bilden Sie einige dimensionslose Produkte aus den Modellgrößen im
Modell zu Beispiel 4.15.
5. DIMENSIONEN UND EINHEITEN
55
Lösung. Basteln wir uns zum Beispiel ein dimensionsloses Produkt, das k2
enthält,wobei [k2]=l2 /(mol2 s). Um die Sekunde wegzukürzen, dividieren wir durch k1 .
Die Konzentrationen kürzen wir, indem wir mit c20 multiplizieren. Wir erhalten das
Produkt
k2 c20
.
k1
Rechnen wir noch einmal die Einheit des Produktes nach:
mol mol
l2
·
· s = 1.
2 ·
l
l
mol s
Alle Einheiten haben sich weggekürzt, wir haben ein dimensionsloses Produkt gefunden.
Überzeugen Sie sich selbst, daß auch die folgenden Produkte dimensionslos sind:
cD c3N O cD
,
, k1 t.
c0 c20 o22
Wir können leicht unendlich viele weitere dimensionslose Produkte erzeugen, indem wir
die, die wir bereits haben, in verschiedener Weise miteinander multiplizieren.
¤
Die dimensionslosen Produkte können als neue Modellgrößen herangezogen werden. Man
erhält dann dimensionslose Gleichungen. Wenn man genügend viele genügend
unabhängige Produkte einführt, erhält man ein dimensionsloses Modell. Dieses ist
gleichwertig mit dem Ausgangsmodell in dem Sinn, dass die Größen der beiden Modelle
eindeutig ineinander umgerechnet werden können.
Merksatz 5.7 (Pi-Theorem). Jedes Formelsystem mit dimensionsbehafteten Größen,
das homogen in der Dimension ist, läßt sich in ein Formelsystem mit dimensionslosen
Produkten dieser Größen umschreiben.
Unter dem Stichwort Dimensionsanalyse findet man in Lehrbüchern über mathematische
Modellierung eine systematische Methode, möglichst wenige dimensionslose Produkte
aufzufinden, aus denen sich alle anderen dimensionslosen Produkte durch
Multiplikationen und Divisionen zusammensetzen lassen. Diese wenigen Produkte
können dann als Größen eines dimensionslosen Modells herangezogen werden. Obwohl
die Methode ziemlich einfach ist, sie beruht auf der Lösung eines linearen
Gleichungssystems, behandeln wir sie hier nicht weiter. In vielen Fällen genügt Probieren
und Geschick ohne viel Theorie.
Beispiel 5.8. Machen Sie das logistische Modell der Populationsdynamik (Beispiel 4.9)
dimensionsfrei. Das Modell lautet:
µ
¶
1
P (t)
d
P (t) = P (t) 1 −
.
dt
τ
K
56
Dabei ist t die Zeit (in Jahren), P (t) die Bevölkerungszahl zur Zeit t (in Individuen). K
ist die Kapazität des Lebensraumes (in Individuen), und τ ist eine Zeitkonstante (in
Jahren), die die Geschwindigkeit der Populationsentwicklung beschreibt.
Lösung. t und τ haben beide die Dimension Zeit. Daher führen wir als neue
dimensionslose Variable
t
s=
τ
ein und betrachten alle Funktionen in Abhängigkeit von s statt t. Wir merken vor:
ds
1
= .
dt
τ
P (t) und K haben dieselbe Dimension (nämlich Bevölkerungszahl). Als neue
dimensionslose Zustandsgröße führen wir
P (t)
u(s) =
K
ein. Also bezeichnet u(s) den Bruchteil der Kapazität des Lebensraumes, der zum
Zeitpunkt t = τ s von der Bevölkerung ausgeschöpft wird. Bei u(s) = 1 ist der
Lebensraum voll gesättigt.
Wir formen jetzt die Differentialgleichung um:
µ
¶
d
1
P (t)
P (t) = P (t) 1 −
dt
τ
K
d
1
Ku(s) = Ku(s)(1 − u(s))
dt
τ
d
ds
K
Ku(s) ·
= u(s)(1 − u(s))
ds
dt
τ
K d
K
u(s) = u(s)(1 − u(s))
τ ds
τ
d
u(s) = u(s)(1 − u(s))
ds
¤
Im dimensionslosen Modell
d
u(s) = u(s)(1 − u(s))
ds
gibt es gar keine Parameter mehr. In dieser Version des Modells ist eine Anpassung an
Daten (die entsprechend transformiert werden müssten) nicht möglich. Parameter, deren
Werte noch nicht festgelegt sind oder zur Anpassung des Modells an verschiedene
Datensätze verwendet werden sollen, wird man also nicht wegtransformieren, es sei denn
vorübergehend, um die Struktur des Modells kennenzulernen.
Die Größen eines dimensionslosen Modells haben meist eine weniger anschauliche
Bedeutung als diejenigen des Ausgangsmodells. Aber der Übergang zu einer möglichst
einfachen Form eines Modells bringt dessen essenzielle Struktur am klarsten zum
Vorschein und eliminiert redundante Größen, die bei der Arbeit mit dem Modell
hinderlich sind. Daher ist es für jede Modelliererin von Interesse, gleichwertige aber
vereinfachte Formen ihres Modells zu erarbeiten.
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
57
6. Simulation deterministischer dynamischer Systeme
In diesem Abschnitt betrachten wir einige Aspekte der Simulation dynamischer Systeme
an Hand von Differentialgleichungs-Modellen.
Merksatz 6.1. Wenn der Anfangszustand und alle externen Größen und Parameter
eines deterministischen dynamischen Modells gegeben sind, lassen sich - zumindest
näherungsweise - zukünftige Zustände am Computer berechnen. Diese Rechnung für eine
Folge von zukünftigen Zeitpunkten heißt Simulation oder Vorwärtsrechnung. Das
Ergebnis bezeichnet man als numerische Lösung.
6.1. Simulationswerkzeuge.
Für die Simulation von dynamischen Systemen steht eine Vielfalt von Hilfsmitteln zur
Verfügung. Viele Kriterien spielen bei der Auswahl mit (die Reihenfolge der Aufzählung
sagt nichts über die Wichtigkeit aus):
(1) Flexibilität: Wird nur ein bestimmtes System in verschiedenen
Parameterkonfigurationen simuliert? Wird ein Programm zur Modellierung
verschiedener dynamischer Systeme gebraucht? Wird ein Universalwerkzeug für
verschiedenste mathematische Operationen gebraucht?
(2) Kapazität und Rechengeschwindigkeit: Wird die Methode auf ein sehr großes
System angewendet? Soll die Simulation im Online-Betrieb eingesetzt werden
und daher fast augenblicklich auf die Eingaben reagieren können?
(3) Bequemlichkeit: Welche Programmierkenntnisse hat der Benützer? Wird ein
Simulationswerkzeug nur einmal eingerichtet und oft benützt, sodaß auch ein
erheblicher Arbeitsaufwand bei der Programmierung gerechtfertigt ist? Soll das
Programm interaktiv betrieben werden?
(4) Vorhandene Hardware.
(5) Kosten.
Low-Level Programmiersprachen, zum Beispiel BASIC, FORTRAN, C oder PASCAL,
sind darauf ausgerichtet, universell alle Funktionen eines Computers programmieren zu
können. Compilersprachen setzen das Programm erst in Maschinensprache um und
haben anschließend einen sehr schnellen Code zur Verfügung. Interpreter (viele Versionen
von BASIC) übersetzen jeden eingegebenen Schritt und führen ihn sofort aus. Das
braucht länger, dafür ist ein interaktiver Betrieb möglich: Man kann die Rechnung sofort
modifizieren, wenn die Ergebnisse der ersten Nebenrechnungen es nahelegen. Low-Level
Programmierung kann die Eigenheiten der verwendeten Hardware optimal ausschöpfen
und wird daher für sehr rechenintensive Probleme eingesetzt. Der Vorbereitungsaufwand
für Low-Level-Programmierung ist sehr hoch, weil sich der Programmierer um viele
Details wie die Dateneingabe und Ausgabe oder die Sicherung des Programmes bei
unvorhergesehenen Vorfällen (logisch inkonsistente Dateneingabe) selbst kümmern muß.
Für die Lösung von Differentialgleichungen, so wie für viele andere Routineaufgaben,
stehen allerdings fertige, ausgetestete Unterprogramme in professionellen
Programmbibliotheken zur Verfügung. Sie müssen nur mehr ins Hauptprogramm
eingelinkt werden.
58
Tabellenkalkulationsprogramme (etwa EXCEL) sind zwar keine Simulationswerkzeuge,
erlauben aber die Programmierung einfacher mathematischer Prozeduren mit viel
weniger Arbeitsaufwand als Programmiersprachen. Die Bearbeitung erfolgt interaktiv.
Für jede Größe wird am Bildschirm eine Zelle eingerichtet, in die die Formel eingetragen
wird, nach der sich die Größe aus den Einträgen der anderen Zellen errechnet. Die
Grenzen sind erreicht, wenn komplexe Rechenprozesse oder große Datenmengen
bewältigt werden sollen.
Interaktive Mathematikpakete sind Mehrzweckprogramme für alle Arten von
mathematischen Problemen. Es handelt sich um Pakete, die sowohl numerisch rechnen
als auch symbolisch Formeln umformen können (zweiteres wird in dieser Vorlesung nicht
verwendet); Beispiele sind MATLAB, MATHEMATICA, MAPLE. Diese
Mathematikpakete stellen unter vielem anderen auch Lösungsverfahren für
Differentialgleichungen zur Verfügung. Der Vorteil dieser Werkzeuge besteht in ihrer
Vielseitigkeit und Bequemlichkeit. Erst bei sehr großen Problemen stößt man auf die
Grenzen der Kapazität.
Simulationssprachen (wie SIMULINK, VENSIM, STELLA) sind auf die Simulation
dynamischer Systeme spezialisiert. Für diese Bedürfnisse bieten sie optimale
Bequemlichkeit, nämlich eine besonders einfache graphische Programmierung zur
Eingabe von Simulationsdiagrammen, vorbereitete Ein- und Ausgaberoutinen mit
Graphik, und vorbereitete Lösungsverfahren, die man nur mehr aus demMenu auswählen
muß. Sie sind flexibel genug, um verschiedenste Systeme zu simulieren. Die hohe
Bequemlichkeit schlägt sich gelegentlich in deutlichen Einschränkungen der Kapazität
und Rechengeschwindigkeit nieder.
In dieser Vorlesung werden bei Computer-Demonstrationen Matlab und dessen
Simulationspaket Simulink verwendet.
Spezialprogramme sind für bestimmte Anwendungen geschrieben, die immer
wiederkehren (Analyse von Grundwasserströmungen, Schadstoffausbreitung, Spannungsund Dehnungsverhalten in elastischen Strukturen, chemische Reaktionen). Sie sind auf
die Anwendung optimal zugeschnitten und verlangen vom Benützer am wenigsten
Vorbereitungsarbeit.
Wir erwähnen noch ein Simulationswerkzeug, das auf Grund der hohen Leistung der
digitalen Computer jetzt mehr und mehr in den Hintergrund tritt, den Analogrechner.
Im Analogrechner werden die simulierten Systeme als elektrische Schaltkreise
nachgebaut. Man entwirft also aus Transistoren, Widerständen und anderen
Bauelementen Schaltungen, in denen Spannungen auftreten und Ströme fließen, welche
genau die simulierten Differentialgleichungen erfüllen. Die Rechenergebnisse lassen sich
dann als Spannungen an bestimmten Stellen der Schaltung messen, zum Beispiel mit
Oszillographen. Ein Vorteil der Analogrechner ist ihre Geschwindigkeit: Während der
Digitalrechner eine Sequenz von Rechenschritten abarbeitet, “geschieht” das simulierte
System gleichzeitig in allen Stromkreisen des Analogrechners. Die Formulierung
dynamischer Syteme durch Simulationsdiagramme, also Blockschaltbilder, geht
traditionell auf den Einsatz der Analogrechner zurück.
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
59
6.2. Das Eulersche Polygonzugverfahren.
Bevor eine Simulation mit einem Simulationspaket gestartet werden kann, muß der
Benützer Verfahrens-Parameter wie die Rechengenauigkeit oder die Schrittweite
festlegen, und manchmal auch eines unter mehreren Rechenverfahren auswählen. Um die
Rolle dieser Entscheidungen zu beurteilen, müssen wir die Funktionsweise einer
Simulationsrechnung zumindest in groben Zügen verstehen.
Ein kontinuierliches Modell wird auf einem ganzen Zeitintervall t0 < t < tend betrachtet,
sein Zustand numerisch aber nur zu diskreten Zeitpunkten ti zwischen t0 und tend
berechnet. Den Zustand zur Zeit t bezeichnen wir mit x(t), wieder ist x(t) ein Vektor,
der alle Zustandsgrößen zusammenfaßt. Die Entwicklung des Systems wird durch eine
Differentialgleichung beschrieben:
d
x(t) = f (t, x(t)).
dt
Dabei ist dtd x(t) der Vektor, der aus den Wachstumsraten aller Zustandsgrößen besteht,
und f ist eine Funktion, deren Formel bekannt ist. Wir nehmen an, daß der
Anfangszeitpunkt t0 und der Anfangszustand x0 = x(t0 ) bekannt sind, und wollen
zukünftige Zustände zu gewissen Zeitpunkten ti+1 > ti > t0 berechnen.
Wir nehmen der Einfachheit halber zunächst an, daß das System nur eine Zustandsgröße
beinhaltet. Wenn die Zukunft des Systems bekannt wäre, könnten wir also eine Kurve
zeichnen, indem wir waagrecht die Zeit und senkrecht den Zustand auftragen. Genau
diese Kurve kennen wir aber noch nicht, wir sollen sie berechnen.
An Hand der Differentialgleichung
d
(6.1)
x(t) = −x2 (t), x(0) = 2,
dt
erklären wir ein Verfahren, das diese Aufgabe (zunächst schlecht und recht) bewältigt. In
diesem Fall hängt die rechte Seite f nur von x und nicht von t ab, und es ist f (x) = −x2 .
Abbildung 6.1 zeigt das Schema des Verfahrens, und Abbildung 6.2 zeigt maßstabgerecht
die exakte Lösung und zwei Näherungslösungen von (6.1). Die Zahlenwerte sind in
Tabelle 6.1 zusammengestellt.
Gleichung (6.1) ist besonders einfach und läßt sich auch auf dem Papier lösen. Die exakte
Lösung ist
1
.
x(t) =
0.5 + t
(Machen Sie die Probe.) Deshalb kann man an diesem Beispiel gut sehen, wie weit die
Näherungslösungen von der exakten Lösung abweichen.
Wir suchen die Lösung x(t) für 0 ≤ t ≤ 1. Wir unterteilen das Intervall [0, 1] in n
Teilschritte. Der Einfachheit halber, und weil kein Grund besteht, eine andere
Unterteilung zu wählen, wählen wir die Schritte in gleichem Abstand (das wäre aber
nicht unbedingt notwendig):
0 = t0 < h = t1 < 2h = t2 < · · · < nh = 1.
60
Abbildung 6.1. Schema des Eulerschen Polygonzugverfahrens
Tabelle 6.1. Eulerverfahren für (6.1)
ti
0.0
0.25
0.5
0.75
1.0
xi
2.0
1.0
0.75
0.6094
0.5165
x2i
4.0
1.0
0.5625
0.3713
0.2668
hx2i
1.0
0.25
0.1406
0.0928
0.0667
x(ti )
2.0
1.3333
1.0
0.8
0.6667
Näherungslösungen für Schrittweite h = 0.25 und exakte Lösung
Den Abstand bezeichnen wir auch als Schrittweite, und die Stellen t0 , · · · , tn als
Stützstellen. Am liebsten wüßten wir den exakten Wert x(ti ) an jeder der Stützstellen.
Wir werden dafür nur Näherungswerte berechnen, die wir mit x0 , · · · , xn bezeichnen.
Damit unser Beispiel nicht zu lang ausfällt, wählen wir nur n = 4, also nur 5 Stützstellen
und eine Schrittweite von 0.25.
Der Wert x(t0 ) = x(0) = 2 ist bekannt. Wenn wir diesen Wert in f einsetzen, erhalten
wir f (x0 ) = −22 = −4. Das ist aber noch nicht der nächste Zustand, sondern nur die
Wachstumsrate, also die Steigung der Tangente an die Lösungskurve am Punkt t = 0.
Weil wir im Augenblick nichts genaueres über die Lösungskurve wissen, folgen wir der
Tangente bis zur nächsten Stützstelle und hoffen, daß wir uns dabei von der
Lösungskurve nicht allzuweit entfernen. Die Gleichung der Tangente ist
x(t) = x0 + (t − t0 )f (x0 ) = 2 − (t − t0 ) · 4.
An der nächsten Stützstelle t1 = 0.25 nimmt die Tangente den Wert
x1 = 2 − 0.25 · 4 = 1
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
61
an. Wir akzeptieren diesen Wert als Näherungswert für die Lösung. Weil wir aber nicht
der Lösungskurve selbst, sondern nur der Tangente gefolgt sind, müssen wir mit einem
Fehler rechnen, wir liegen nicht mehr genau auf der exakten Lösung: Dieser Fehler heißt
Abbruchfehler.
Wir kennen also jetzt, zumindest näherungsweise, einen Wert der Lösung zur Zeit
t1 = 0.25, nämlich x1 = 1. Wir nehmen an, daß dieser Wert die Lösung zumindest
einigermaßen genau wiedergibt. Wir wiederholen jetzt genau dieselbe Prozedur: Nach der
Differentialgleichung ist die Steigung der Tangente der Lösungskurve im Punkt t1
d
x(t1 ) = f (x(t1 )) ≈ f (x1 ).
dt
In unserem Beispiel ist das
d
x(t1 ) ≈ −x21 = −1.
dt
Wir folgen also der Geraden
x = x1 + (t − t1 )f (x1 ) = 1 − (t − 0.25) · 1
bs zur nächsten Stützstelle t2 = 0.5 und erhalten dort
x2 = x1 + (t2 − t1 )f (x1 ) = 1 − 0.25 · 1 = 0.75.
Selbst wenn x1 der exakte Wert der Lösung an der Stelle t1 gewesen wäre, müßten wir
einen Abbruchfehler in Kauf nehmen, denn wieder sind wir einer Tangente statt der
Lösungskurve selbst gefolgt.
Wir wiederholen diesen Schritt solange, bis wir alle Stützstellen erreicht haben. Die
Formel zur Berechnung des Näherungswertes xi+1 an der Stützstelle ti+1 aus dem
Näherungswert xi an der Stützstelle ti ist, wie wir an der bisherigen Rechnung gesehen
haben,
xi+1 = xi + (ti+1 − ti )f (ti , xi ).
Tabelle 6.1 zeigt, welche Zahlenwerte dabei herauskommen.
Merksatz 6.2. Numerische Lösungsverfahren für Differentialgleichungen bestimmen
vom Anfangszustand aus mit Hilfe der Differentialgleichung die Richtung, in die sich der
Zustand weiterentwickeln wird. Sie gehen ein kurzes Geradenstück in diese Richtung und
bestimmen dann eine neue Richtung durch Einsetzen des ersten Näherungszustands in
die rechte Seite der Differentialgleichung - und so weiter bis zum gewählten Endzeitpunkt.
Je länger wir den Tangenten folgen, ohne unseren Kurs zu korrigieren, desto weiter
entfernen wir uns (im Allgemeinen) von der Lösungskurve. Daher wird das Verfahren
eine bessere Näherung ergeben, wenn wir kürzere Schritte machen. Abbildung 6.2 zeigt
die exakte Lösung von (6.1) und die Näherungslösungen für die Schrittweiten 0.25 und
0.1. Wir sehen, daß bei der kürzeren Schrittweite die Lösung besser wiedergegeben wird.
Um den Abbruchfehler gering zu halten, muß man also eine kleine Schrittweite wählen.
Leider entsteht außer dem Abbruchfehler ein zweiter Rechenfehler: Weil auch der beste
Computer nur endlich viele Ziffern speichern kann, wird bei jedem Rechenschritt der
Zahlenwert des Ergebnisses auf eine gewisse Anzahl von Dezimalstellen gerundet. Dabei
62
Abbildung 6.2. Eulersches Polygonzugverfahren für (6.1)
2 *o
1.8
1.6
*
1.4
x
*
1.2
*
1
*
o
*
*
0.8
*
o
*
*
o
0.6
*
o
0.4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t
Exakte Lösung und zwei Näherungslösungen:
Kreise: Schrittweite 0.25; Sterne: Schrittweite 0.1
entsteht ein Rundungsfehler. Dieser Rundungsfehler wirkt besonders gravierend, wenn
Zahlen von sehr verschiedener Größe zueinander addiert werden. Zum Beispiel
verschwindet bei Rundung auf 5 signifikante Stellen der zweite Summand in der Addition
1 + 0.00001 völlig, er wirkt sich erst auf der sechsten Stelle aus: 1.00001. (Noch
schlimmer wirkt der Rundungsfehler bei Subtraktion von fast gleichgroßen Zahlen.)
Wenn die Schrittweite sehr klein gewählt wird, wird in jedem Schritt zur bisherigen
Näherungslösung xi nur eine winzige Zahl addiert, und der Rundungsfehler kann sich
deutlich auswirken. Dazu kommt, daß bei einer kleinen Schrittweite sehr viele Schritte
durchgerechnet werden müssen. Das macht die Rechnung nicht nur langsam, vor allem
können sich die Rundungsfehler aufschaukeln.
Merksatz 6.3. Zwei Fehler treten bei Anwendung von numerischen Lösungsverfahren
für Differentialgleichungen auf. Der Abbruchfehler stammt daher, daß man mit jedem
Schritt einer Geraden statt der tatsächlichen Lösungskurve folgt (die man ja nicht kennt).
Der Rundungsfehler entsteht, weil das Ergebnis jedes Rechenschrittes gerundet wird.
Die Wahl der Schrittweite muß daher einen Kompromiß schließen: Durch die Wahl einer
kleinen Schrittweite wird der Abbruchfehler verringert, aber der Rundungsfehler
vergrößert.
Gute Verfahren sind so gebaut, daß sie schon bei relativ großen Schrittweiten nur geringe
Abbruchfehler machen.
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
63
6.3. Bessere Verfahren zur numerischen Lösung von
Differentialgleichungen.
Wir haben das Eulersche Polygonzugverfahren eingeführt, weil es besonders einfach ist,
und sich daran das Prinzip der numerischen Lösung von Differentialgleichungen schön
erklären läßt. Im Hinblick auf die Relation des Abbruchfehlers zur Schrittweite ist es
leider ein schlechtes Verfahren. Wir werden jetzt überlegen, wie man das Verfahren
effizienter gestalten könnte.
Im Zeitintervall von t0 zu t1 ändert die Lösungskurve kontinuierlich ihre Richtung,
während das Näherungsverfahren mangels besserer Information seine Richtung beibehält.
Erst zum Zeitpunkt t1 wird eine neue Richtung bestimmt. Aus der neuen Richtung läßt
sich aber der Trend ermessen, nach dem die Lösungskurve ihre Richtung geändert hat,
und diese Information kann auf zwei Arten genützt werden:
Aus der Richtungsänderung der Tangenten zwischen ti und ti+1 kann man auf die Größe
des Abbruchfehlers schließen. Wenn man den geschätzten Abbruchfehler mit der Anzahl
der Schritte multipliziert, hat man eine Schätzung, wie weit die Näherungslösung von der
exakten Lösung entfernt sein wird. Wenn die Rechnung zu ungenau ausfällt, muß man
den Zeitschritt in kürzere Intervalle teilen, um den Abbruchfehler zu verringern. Anders
ausgedrückt: Der Schritt von ti bis ti+1 wird nur akzeptiert, wenn |f (xi+1 ) − f (xi )| eine
gewisse Schranke nicht überschreitet. Wenn der Unterschied zu groß ist, wird der Schritt
von ti aus mit einer kleineren Schrittweite wiederholt. Wenn dagegen der Unterschied
besonders gering ausfällt, kann man für den nächsten Schritt eine größere Schrittweite
ins Auge fassen. Auf Grund solcher Überlegungen kann man Verfahren programmieren,
die automatisch eine geeignete Schrittweite wählen.
Merksatz 6.4. Programmen mit automatischer Schrittweitensteuerung gibt der
Benützer eine größte und eine kleinste erlaubte Schrittweite vor, und eine
Größenordnung der Genauigkeit, mit der die Näherungslösung die exakte Lösung
wiedergeben soll. Das Verfahren versucht dann, sich innerhalb dieses Rahmens von
Schritt zu Schritt die geeignete Schrittweite auszuwählen.
Die Gerade, entlang der man den Punkt xi+1 zur Zeit ti+1 erreicht hat, hat im
Eulerverfahren die Steigung f (xi ). Wir haben damit eine erste Schätzung des nächsten
Punktes der Lösung (Predictor – Voraussage). Die Steigung der Lösungskurve hat sich
aber kontinierlich verändert und liegt näherungsweise bei f (xi+1 ), wenn der Zeitpunkt
ti+1 erreicht ist. Wir können diese Trendinformation ausnützen und den Schritt von ti
auf ti+1 mit einer neuen Steigung wiederholen, zum Beispiel mit 12 (f (xi+1 ) + f (xi )).
Dieser Wert wird insgesamt der Lösungskurve besser entsprechen, weil er die Steigung zu
Beginn und Ende des Schrittes berücksichtigt. Der Predictor dient also nur als erster
Versuch zur Erkundung des Trends, der endgültige Schritt wird im Sinne des Trends
korrigiert (Corrector). Professionelle Lösungsverfahren, wie etwa das
Runge-Kutta-Verfahren, tasten jedes Schrittintervall erst vier- bis fünfmal ab, bevor die
endgültige Richtung des Schrittes festgelegt wird. Indem der einzelne Schritt sehr
sorgfältig abgewogen wird, erfolgt er mit einem viel kleineren Abbruchfehler als beim
Eulerverfahren. Der Preis dafür ist, daß bis zur Durchführung des Schrittes die Funktion
64
Abbildung 6.3. Lösung einer steifen Differentialgleichung
x10 -6
1
0.9
0.9
0.8
0.8
0.7
0.7
Konzentration von NO
Konzentration des Donors
x10 -6
1
0.6
0.5
0.4
0.3
0.6
0.5
0.4
0.3
0.2
0.2
0.1
0.1
0
0
5000
Zeit t
0
0
5000
Zeit t
k1 = 5 · 10−3 , k2 = 107 , o2 = 2.5 · 10−4 , c0 = 10−6 .
f mehrmals für verschiedene vorläufige Werte berechnet werden muß. Das kostet
Rechenzeit. Um die genauen Formeln und eine Aufzählung von verschiedenen Methoden
wollen wir uns hier nicht kümmern.
Merksatz 6.5. Es gibt sehr viele verschiedene Lösungsverfahren für
Differentialgleichungen, die bei verschiedenen Problemen verschieden gut arbeiten.
Simulationspakete stellen meist mehrere Verfahren zur Verfügung. In den Handbüchern
bzw. in der online Hilfe findet man Hinweise, welche Bedingungen für die Anwendung
der jeweiligen Methoden sprechen. Wenn ein Lösungsverfahren mit einer Fehlermeldung
abstürzt oder sehr langsam rechnet, sollte man andere Verfahren versuchen.
Es gibt eine bestimmte Sorte von Differentialgleichungen, die sich der numerischen
Lösung besonders hartnäckig widersetzen, sogenannte steife Differentialgleichungen.
Steife Differentialgleichungen beschreiben Systeme, in denen Prozesse mit sehr
verschiedenen Geschwindigkeiten ablaufen. Stellen Sie sich zum Beispiel einen Kessel vor,
dem langsam Chemikalien zugeführt und entzogen werden. Zwischen den Chemikalien,
die dort aufeinandertreffen, können sehr schnelle Reaktionen ablaufen. Die
Konzentrationen im Kessel werden vom langsamen Materialtransport und den schnellen
Reaktionen beeinflußt.
Merksatz 6.6. Steife Differentialgleichungen beschreiben Systeme, in denen Prozesse
mit sehr unterschiedlichen Geschwindigkeiten ablaufen. Solche Modelle sind numerisch
besonders fehleranfällig. Simulationspakete bieten oft spezielle Lösungsverfahren für steife
Differentialgleichungen an.
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
65
Abbildung 6.4. Ein Differentialgleichungslöser stürzt ab
x10 -6
1
0.8
0.8
0.6
0.6
0.4
0.4
Konzentration von NO
Konzentration des Donors
x10 -6
1
0.2
0
-0.2
-0.4
0.2
0
-0.2
-0.4
-0.6
-0.6
-0.8
-0.8
-1
0
5000
Zeit t
-1
0
5000
Zeit t
k1 = 5 · 10−3 , k2 = 107 , o2 = 2.5 · 10−4 , c0 = 10−6 .
Ein Runge-Kutta-Verfahren wird mit zu grober Genauigkeitsvorgabe angesetzt.
Das Verfahren versucht viel zu lange Schrittweiten und verliert völlig die Orientierung.
Als Beispiel für ein steifes System zeigen wir in Abbildung 6.3 den Verlauf der
Zustandsgrößen des NO-Modells aus Beispiel 4.15 unter Bedingungen, bei denen der
Donor im Verhältnis sehr viel schneller zerfällt als das NO. (Weil die Zerfallsrate von NO
proportional zum Quadrat der Konzentration ist, während die Zerfallsrate vom Donor
proportional zur Konzentration ist, treten solche Bedingungen bei sehr kleinen
Konzentrationen ein.) Der rapide Anstieg zu Beginn kommt davon, daß fast
augenblicklich der gesamte Donor zu NO zerfällt. Die langsam abfallende Kurve zeigt
den allmählichen Abbau von NO. Wir sehen, warum steife Differentialgleichungen so
schwer zu lösen sind: Der schnelle Prozeß bestimmt die Anfangsrichtung der Tangente.
Wir haben strichliert die Tangenten zum Anfangszeitpunkt eingezeichnet. Der schnelle
Zerfall des Donors bewirkt, daß diese Tangenten sehr steil sind. Wenn das Verfahren
diesen Tangenten zu lange folgt, schießt es gleich mit dem ersten Schritt weit über das
Ziel hinaus und berechnet von dort aus völlig unsinnige Werte weiter. Die Schrittweite
muß also sehr kurz sein. Um den langsamen Prozeß zu beobachten, muß aber ein
ausreichend langes Zeitintervall mit diesen kurzen Schritten überdeckt werden. Dabei
können die Rundungsfehler dramatisch kumulieren.
Die Gleichung in unserem Beispiel ist noch nicht sehr steif. Wenn man die
Rechentoleranz genügend klein vorgibt, reichen auch Runge-Kutta-Verfahren aus.
Als Beispiel für das Versagen eines Differentialgleichungslösers lassen wir das Programm
ode23 aus Matlab das NO-Modell lösen, ohne die Genauigkeiten vorzugeben. Das
Verfahren zielt dann von sich aus auf relative und absolute Genauigkeit von 10−3 und
10−6 ab, das ist zu grob im Vergleich mit dem Anfangswert 10−6 . Bilder wie
Abbildung 6.4 erhält man mitunter beim ersten Versuch, eine Differentialgleichung
66
numerisch zu lösen. Die Abhilfe kommt durch kleinere Fehlertoleranzen und
Schrittweiten.
Merksatz 6.7. Wenn ein Verfahren bei der Lösung einer Differentialgleichung auf
offensichtlich unsinnige Werte kommt oder mit einer Fehlermeldung abbricht, zwingt
man es durch Vorgabe kleinerer Toleranzen und kleinerer Schrittweiten dazu, sich
langsamer vorzutasten, oder man wählt ein anderes Verfahren.
VORSICHT! Es gibt aber auch Differentialgleichungen, die gar nicht bis in alle Zukunft
Lösungen besitzen, weil ihre Lösungen schon in endlicher Zeit gegen unendlich streben
(“explodieren”). Auch ein gutes Verfahren wird bei solchen Gleichungen aufgeben
müssen, wenn die Werte zu groß werden.
Bemerkung 6.8. Zwar werden bei der Beschreibung der Minimumsuche (Seite 25 f)
und der obigen Beschreibung der numerischen Lösung von Differentialgleichungen gleiche
Worte verwendet (etwa ’Startwert’, ’Schrittweite’, ’Toleranz’), dennoch handelt es sich
um völlig unabhängige Verfahren zur Lösung völlig unterschiedlicher Probleme. Verfallen
Sie nicht in den Fehler, diese Dinge gedanklich miteinander zu vermischen, sondern
machen Sie sich die Unterschiedlichkeit der beiden Problemstellungen und die
Verschiedenheit der beiden Verfahren klar bewusst.
6.4. Lösung von Differentialgleichungen in MATLAB.
Wir beschreiben hier, wie man die Differentialgleichung dtd x(t) = f (t, x(t)) in MATLAB
lösen kann.
Es stehen mehrere Verfahren mit automatischer Schrittweitensteuerung zur Verfügung:
Die gebräuchlichsten sind das Runge-Kutta-Verfahren ode45 und das Verfahren variabler
Ordnung ode15s für steife Differentialgleichungen (das Kürzel ode kommt von ’ordinary
differential equation’). Sie unterscheiden sich in den verwendeten Näherungen für die
Integration und damit auch in der Anzahl der Auswertungen von f in jedem einzelnen
Zeitschritt. x(t) ist dabei ein Vektor, es können also Systeme von Differentialgleichungen
mit mehreren Zustandsgrößen gelöst werden.
Man schreibt zunächst einen MATLAB-Funktionsfile, der die Funktion f (t, x) erklärt,
und speichert ihn. Wir zeigen als Beispiel den Funktionsfile, der die rechte Seite des
Donorproblems beschreibt:
function y=funk(t,x);
% rechte Seite des Donorproblems
% 1. Komponente von x: Donor
% 2. Komponente von x: NO
global k1 k2 o2
dzr=k1*x(1); % Donorzerfallsrate
6. SIMULATION DETERMINISTISCHER DYNAMISCHER SYSTEME
67
nor=k2*o2*x(2)*x(2); % NO-Oxidationsrate
y(1)=-dzr; % Zuwachsrate Donor
y(2)=dzr-nor; % Zuwachsrate NO
y=y’; % das Resultat muss ein Spaltenvektor sein
Diese function berechnet aus den im Argument x vom rufenden Programm übergebenen
Konzentrationen die rechte Seite des Gleichungssystems. Das Argument t wird in diesem
Beispiel nicht verwendet, muss aber an erster Stelle der Argumenteliste stehen. Wie
schon erwähnt ist % das Kommentarzeichen, was rechts davon steht ist bei der
Ausführung der function irrelevant. Wir sichern diesen File unter dem Namen funk.m.
Zur Steuerung der Rechengenauigkeit und anderer Verfahrensparameter verwendet man
die MATLAB-Funktion odeset. In odeset wird die gewünschte absolute Genauigkeit
mit der Bezeichnung AbsTol und die gewünschte relative Genauigkeit mit der
Bezeichnung RelTol gesetzt. Die Namen und Werte sind Argumente von odeset, ihr
Ergebnis wird als Variable gespeichert und an ode45 als Argument übergeben. Dabei
meint AbsTol den Abstand der berechneten Lösung x(t) zur (unbekannten) wahren
Lösung x̂(t) und muss Problem-spezifisch gewählt werden. Dagegen muss RelTol nicht
der absoluten Größe von x̂(t) entsprechend gewählt werden, denn gefordert wird
|x(t) − x̂(t)|
< RelTol.
|x̂(t)|
Es besteht allerdings keine Garantie, dass die berechnete Lösung überall näher bei der
exakten Lösung liegt als die gewählten Toleranzen.
Wir legen das Zeitintervall (hier tspan), die Anfangsdaten (hier in x0) und die
Verfahrensparameter (hier unter options) fest und rufen dann ode45 auf. Anschließend
sehen wir uns das Ergebnis mit plot am Bildschirm an:
>>global k1 k2 o2 ←>>k1=5e-3; ←>>k2=1e+7; ←>>o2=2.5e-4; ←>>tspan=[0 1000]; ←>>x0=[1e-6;0]; ←>>options=odeset(’AbsTol’,1e-10,’RelTol’,1e-4); ←>>[t,x]=ode45(’funk’,tspan,x0,options); ←>>plot(t,x(:,1),’-’,t,x(:,2),’--’) ←10−10 klingt sehr genau, aber bedenken Sie, daß der Anfangswert 10−6 ist. Der String
’funk’ ist der Name des Funktionsfiles, der die rechte Seite der Differentialgleichung aus
den jeweils vom rufenden Programm (hier ode45) übergebenen Argumenten t und
x=[x(1);x(2)] berechnet. Als Ausgabe von ode45 erhält man den Spaltenvektor t, das
ist eine Liste aller Stützstellen, und die Matrix x. Jede Zeile dieser Matrix entspricht
einer Stützstelle ti und enthält die Werte des Lösungsvektors x(ti ),
 


t1
x1 (t1 ) · · · xn (t1 )
t = t2  , x = x1 (t2 ) · · · xn (t2 ) ,
..
..
..
.
.
.
(im Donorproblem ist n = 2).
68
Logisch wesentlich an dieser Vorgangsweise ist also: Der Benützer des im Numerik-Paket
bereitgestellten Dgl-Lösers (hier ode45) muss eine Funktion programmieren (hier funk),
die aus einem Zeitpunkt t und einem Zustandsvektor x den zugehörigen Vektor f (t, x)
der rechten Seite des Dgl-Systems berechnet. Diese Funktion sollte möglichst flexibel
programmiert sein, bleibt es doch dem Dgl-Löser überlassen, an welchen Stellen und wie
oft er das Benutzerprogramm aufruft. Da dies typischerweise mehrere tausend mal
geschieht, ist bei großen Problemen der Benützer auch an der Rechenzeit-Effizienz seines
Programms (hier funk) interessiert.
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
69
7. Fragestellungen zu dynamischen Systemen
7.1. Dynamische Systeme.
Wir verzichten, wie immer, auf eine exakte mathematische Definition und begnügen uns
mit einem Appell an die Anschauung:
Merksatz 7.1. In einem dynamischen System ändert sich ein Zustand in Abhängigkeit
von einer Zeitgröße. Dabei bestehen Gesetzmäßigkeiten, die Voraussagen auf zukünftige
Zustände an Hand des gegenwärtigen Zustands gestatten.
Sehr oft bezeichnet die Zeitgröße tatsächlich die Zeit, wie in allen Modellen, die wir
bisher formuliert haben. Allerdings ist es vom rein mathematischen Standpunkt natürlich
auch möglich, daß die Zeitgröße eine völlig andere physikalische (biologische, . . . )
Interpretation hat. Um Sie nicht unnötig zu verwirren, nehmen wir von solchen
Beispielen Abstand.
Beispiel 7.2. An einer Feder hängt ein Gewicht. Wenn das Gewicht aus der Ruhelage
gebracht wird, zum Beispiel durch Anziehen, wird es auf- und abschwingen. Wir
beobachten die Bewegung des Gewichtes. Im Lauf der Zeit ändert sich der Zustand des
Systems: Der Zustand besteht aus zwei unabhängigen Größen, nämlich der Höhe, in der
sich das Gewicht im Augenblick befindet, und der Geschwindigkeit, mit der es sich bewegt.
Beispiel 7.3. Ein Wald ist von parasitären Insekten befallen. Alle Jahre wird durch
Stichproben der Bestand an Larven ermittelt. Der Zustand des Systems, der sich von
Jahr zu Jahr ändert, ist die Anzahl der Larven.
Beide Beispiele, 7.2 und 7.3, passen in das Konzept eines dynamischen Systems. Je nach
Beschaffenheit und vor allem nach Betrachtungsweise können wir aber dynamische
Systeme in sehr unterschiedliche Klassen einteilen.
In Beispiel 7.2 betrachten wir am besten die Zeit als ein fließendes Kontinuum. Die
Beobachtungszeitpunkte bilden ein Intervall, ein ganzes Geradenstück auf der
Zahlengeraden. Die Bewegung kann man durch eine Kurve darstellen. Das
mathematische Modell besteht aus einer Differentialgleichung.
Dagegen bietet sich für Beispiel 7.3 eine diskrete Betrachtungsweise an: Weil es zu
bestimmten Jahreszeiten gar keine Larven gibt, erfaßt man sinnvollerweise die Anzahl
der Larven jedes Jahr zu bestimmten Stichtagen. Die Beobachtungszeitpunkte bilden
eine Folge. Um in diesem System eine Gesetzmäßigkeit zu formulieren, sucht man eine
Formel, die den Larvenbestand des nächsten Jahres aus dem heurigen Bestand
berechnet, eine sogenannte Differenzengleichung.
Es ist nicht zwingend vorgegeben, ob ein natürliches System diskret oder kontinuierlich
modelliert wird. Wenn man die Bewegung des Gewichtes filmt, empfindet man die
70
gefilmte Bewegung zwar als kontinuierlich, aber genau genommen wird nur zu diskreten
Zeitpunkten fotografiert. Man könnte auch die Entwicklung der Larven innerhalb eines
Monats im Rahmen eines Räuber-Beute-Modells durch eine Differentialgleichung
beschreiben wollen.
Merksatz 7.4. Ein diskretes dynamisches System wird nur zu einzelnen, fest
vorgegebenen Zeitpunkten betrachtet, die eine Folge bilden. Im kontinuierlichen System
läuft die Betrachtung über ein ganzes Zeitintervall.
Die Simulation von deterministischen diskreten dynamischen Systemen ist
vergleichsweise einfach. Das System wird zu den Zeitpunkten
t0 < t1 < t2 < · · ·
betrachtet. Zu diesen Zeitpunkten durchläuft das System die Zustände
x0 , x 1 , x 2 , · · · .
Im Vektor xi werden alle Zustandsgrößen zum Zeitpunkt ti zusammengefaßt. Der
Zusammenhang zwischen Gegenwart und Zukunft wird durch eine Differenzengleichung
bestimmt:
xi+1 = F (ti , xi ).
Dabei ist F das Modell, eine Funktion, deren Formel man kennt. Um vom bekannten
Zustand x0 aus den nächsten Zustand x1 zu erhalten, muß nur in die Formel eingesetzt
und F (t0 , x0 ) berechnet werden. Wenn man diese Prozedur wiederholt, erreicht man nach
und nach alle zukünftigen Zustände.
In Beispiel 7.2 kann man auf Grund einfacher physikalischer Gesetze die Bewegung des
Gewichtes sehr genau vorhersagen, wenn die Masse des Gewichtes und die Stärke der
Feder bekannt sind. Ein System, das erlaubt, (zumindest theoretisch) die Entwicklung in
der Zukunft exakt aus dem gegenwärtigen Zustand zu bestimmen, heißt es ein
deterministisches System.
In Beispiel 7.3 sind die Voraussagen weniger sicher. Der Larvenbestand wird von vielen
Zufälligkeiten beeinflußt, die sich der Voraussage entziehen (z.B. Wetter). Die
Zusammenhänge zwischen Gegenwart und Zukunft sind nur Gesetzmäßigkeiten im Sinne
von Wahrscheinlichkeitsaussagen. Man kann zum Beispiel voraussagen, mit welcher
Wahrscheinlichkeit der Larvenbestand im nächsten Jahr größer sein wird als heuer. Ein
System, in dem die Gegenwart die Zukunft nicht eindeutig bestimmt, sondern nur
Wahrscheinlichkeiten für das Eintreffen zukünftiger Zustände vorgibt, heißt ein
stochastisches System. Statt stochastisches System sagt man gleichbedeutend auch
probabilistisches System.
Auch die Wahl zwischen deterministischer und stochastischer Modellierung liegt eher an
der Betrachtungsweise als am beobachteten System selbst. Die Quantenmechanik geht
davon aus, daß alle Bewegungen nur stochastisch beschrieben werden können. Es wäre
wahrscheinlich überzogen, das Feder-Masse-System aus Beispiel 7.2 mit
Quantenmechanik anzusetzen. Aber Einflüsse durch Luftzug, oder Vibrationen der
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
71
Aufhängung durch vorbeifahrende Sechsachser könnten durchaus stochastisch modelliert
werden. Andererseits kann man nach dem Gesetz der großen Zahlen Systeme, in denen
Zufallsgesetze gelten, gut durch deterministische Modelle beschreiben, wenn nur die
betrachteten Populationen groß genug sind, sodaß sich die Zufälligkeiten gegeneinander
aufheben.
Merksatz 7.5. Die Zusammenhänge zwischen dem gegenwärtigen Zustand und den
zukünftigen Zuständen eines stochastischen Systems sind durch
Wahrscheinlichkeitsaussagen beschrieben. Im deterministischen System wird jeder
zukünftige Zustand eindeutig durch den gegenwärtigen Zustand bestimmt.
Lassen Sie sich bitte von der Wahl der Beispiele nicht irreleiten: Es gibt auch
stochastische kontinuierliche, und deterministische diskrete Systeme. Alle Modelle, die
wir in diesem Semester behandeln, sind deterministisch und kontinuierlich. Im zweiten
Semester werden wir auch Beispiele für diskrete und für probabilistische Modelle
kennenlernen.
7.2. Anatomie eines dynamischen Systems.
Betrachten wir zur Illustration eine Modifikation von Beispiel 4.15.
Beispiel 7.6. In einer Lösung zerfällt, wie schon in Beispiel 4.15 beschrieben, eine
Donorsubstanz zu NO. Dieses wird oxidiert und dadurch aus der Lösung entfernt. Bei
beiden Prozessen wird Wärme frei. Ein Sensor mißt laufend die NO-Konzentration der
Lösung und gibt auf einem Plotter die Kurve aus. Der Analytiker kann kontinuierlich
Donor zugeben, wobei er die Zuflußgeschwindigkeit über eine regelbare Pumpe steuern
kann. Außerdem kann er die Temperatur der Lösung durch ein Thermometer
überwachen, und der Lösung durch eine regelbare Kühleinrichtung Wärme entziehen,
indem durch eine Kühlschleife Kühlflüssigkeit geleitet wird.
Abbildung 7.1 zeigt ein Wirkungsdiagramm dieses Systems. Versuchen Sie selbst, die
dargestellten Zusammenhänge zu rechtfertigen.
Wir wollen an diesem Beispiel erklären, welche Arten von Größen in einem dynamischen
System vorkommen können.
Im Zentrum des Systems stehen die Zustandsgrößen. Sie ändern sich im Verlauf der Zeit,
und geben zu jedem Zeitpunkt den augenblicklichen Zustand des Systems wieder. In
Beispiel 7.6 kann man als Zustandsgrößen die Konzentrationen von Donor, NO und O2
sowie die Temperatur der Lösung heranziehen. Man kann aber auch einen anderen Satz
von Zustandsgrößen aufbauen, etwa die Masse des gesamten gelösten Donors, das
gesamte gelöste NO, das gesamte gelöste O2 , und die in der Lösung verfügbare freie
Wärme. Offensichtlich beschreiben beide Datensätze den Zustand des Systems gleich gut.
Dem Analytiker bleibt die Wahl der bequemsten Formulierung. Man wählt aber
jedenfalls so wenig Zustandsgrößen wie möglich. Die Zeit wird durch die Zeitgröße
beschrieben, in unserem System eine kontinuierliche Größe, in Sekunden.
Die Gesetzmäßigkeiten der chemischen Reaktionen und der Wärmezufuhr werden mit
Hilfe der Zuwachsraten der Chemikalien und der Temperatur, also generell über die
72
Abbildung 7.1. Wirkungsdiagramm zu Beispiel 7.6
Kursiv: Steuerungs- und Beobachtungsgrößen
Zuwachsraten der Zustandsgrößen beschrieben. Das ist typisch für kontinuierliche
Systeme, die auf Differentialgleichungen führen.
Die weitere Entwicklung des Systems wird durch den gegenwärtigen Zustand, aber auch
durch exogene Größen bestimmt. Exogene Größen beschreiben Effekte, die von
außerhalb auf das System einwirken. Wenn diese Größen gezielt beeinflußt werden
können, sprechen wir auch von Steuergrößen. In stochastischen Systemen könnte man
auch die Zufallseinflüsse als externe Größen auffassen. In unserem System gibt es zwei
Steuergrößen: Den Zufluß an Donor und den Durchsatz an Kühlflüssigkeit.
Nicht alle Systemgrößen sind dem Beobachter direkt zugänglich. Sehr oft sind gerade die
Zustandsgrößen nicht direkt meßbar, sondern nur durch Rückschlüsse aus ihren
Auswirkungen errechenbar. Zugängliche Größen bezeichnen wir als Beobachtungsgrößen.
In unserem System können zwei Zustandsgrößen direkt gemessen werden, sind also auch
Beobachtungsgrößen: die NO-Konzentration und die Temperatur.
Die Wechselwirkungen zwischen allen diesen Größen werden durch verschiedene
Gesetzmäßigkeiten geregelt. In unserem Modell hängen viele Größen voneinander ab: die
Zerfallsgeschwindigkeiten der Chemikalien von der Temperatur, der Wärmeentzug durch
die Kühlanlage vom Kühlmitteldurchsatz und der Lösungstemperatur, die Temperatur
von der freien Wärme, die freiwerdende Wärme von der Anzahl der oxidierten
NO-Moleküle, und so weiter. In den Formeln dieser Gesetze kommen Konstanten vor
(etwa die spezifische Wärme des Lösungsmittels, die Zerfallsgeschwindigkeit des Donors
bei 25◦ C), die Eigenschaften des Systems beschreiben, die vom augenblicklichen Zustand
oder von der Zeit nicht abhängen. Diese Konstanten heißen Parameter.
Auf dem Weg der Berechnung der Zuwachsraten der Zustandsgrößen aus den
Zustandsgrößen, den exogenen Größen und den Parametern treten viele
Nebenrechnungen auf. Die Resultate dieser Nebenrechnungen heißen Zwischengrößen. Zu
den Zwischengrößen in unserem Beispiel gehören der augenblickliche Wärmeentzug durch
die Kühlung, oder wieviel NO pro Sekunde oxidiert wird.
Merksatz 7.7. Größen folgender Art können in einem dynamischen System auftreten:
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
73
• Die Zeitgröße.
• Zustandsgrößen. Sie hängen von der Zeitgröße ab und beschreiben den
augenblicklichen Zustand.
• Exogene Größen. Sie quantifizieren Einflüsse, die von außen an das System
herangetragen werden. Wenn man eine exogene Größe gezielt beeinflussen kann,
spricht man auch von einer Steuergröße.
• Beobachtungsgrößen. Sie beschreiben Auswirkungen des Zustands, die einer
direkten Beobachtung oder Messung zugänglich sind.
• Parameter. Sie beschreiben Eigenschaften des Systems oder seiner Teile, die vom
augenblicklichen Zustand nicht abhängig sind und sich nicht mit der Zeit ändern.
• Zwischengrößen. Das sind Größen, die sich aus den bisher beschriebenen
Größen berechnen lassen. Im Allgemeinen treten sie als Produkte von
Nebenrechnungen bei der Berechnung der Zuwachsraten auf.
• Wachstumsraten der Zustandsgrößen (nur in kontinuierlichen Modellen). Sie
lassen sich aus den oben beschriebenen Größen berechnen und bestimmen die
Weiterentwicklung des Zustandes.
Beispiel 7.8. Betrachten Sie den Kessel aus Beispiel 4.2. Wir nehmen an, daß der
Kessel ein geschlossener Stahltank ist. Der Zufluß wird über ein Ventil geregelt. Im
Abfluß befindet sich ein Strömungsmeßgerät, das die pro Sekunde abfließende
Wassermenge anzeigt. Klassifizieren Sie die Größen in diesem System.
Wie so oft, bedeutet die Zeitgröße wirklich Zeit. Sinnvollerweise werden wir hier die Zeit
kontinuierlich modellieren. Der Zustand ist beschrieben, wenn wir wissen, wieviel Wasser
sich im Augenblick im Kessel befindet. Wahlweise können wir als Zustandsgröße
entweder die Höhe des Wasserstandes oder das enthaltene Wasservolumen heranziehen.
Wir wählen aber nur soviele Zustandsgrößen, wie unbedingt gebraucht werden, also eine
von beiden. Die andere ist dann bei Bedarf eine Hilfsgröße.
Der Wasserzufluß ist eine exogene Größe, eine Steuergröße, die wir durch Drehen am
Ventil beeinflussen können. Der Wasserabfluß wird gemessen, und kann daher als einzige
Größe direkt beobachtet werden, denn ins Innere des Stahltanks können wir nicht sehen.
Glücklicherweise läßt sich der Wasserstand aus der Abflußgeschwindigkeit berechnen.
(Wie?)
Das Modell enthält viele Parameter, nämlich die Abmessungen des Kessels, den
Strömungswiderstand des Abflußrohres, und das spezifische Gewicht von Wasser.
Wir bemerken noch, daß für dieses, durch strenge physikalische Gesetzmäßigkeiten
bestimmte System ein deterministisches Modell durchaus adäquat ist. Das ändert sich
sofort, wenn wir eine Risikoanalyse durchführen: “Mit Wahrscheinlichkeit von 30% ist
einmal im Jahr der Abfluß verstopft . . . ”.
74
7.3. Verschiedene Fragestellungen.
Man untersucht dynamische Systeme aus zwei Gründen: Entweder braucht man eine
neue Veröffentlichung auf der Publikationsliste, oder man will an Hand der Studie einen
Komplex von Fragen beantworten. Welche Fragen beantwortet werden sollen, sollte man
sich schon vor der Modellbildung gründlich überlegen. Wir geben eine kurze
Klassifizierung von Fragen, die im Zusammenhang mit einem dynamischen System
auftreten können.
Simulation oder Vorwärtsrechnung: Die Systemeigenschaften, also Struktur und
Parameter, sind bekannt. Man kennt auch den gegenwärtigen Zustand und die externen
Größen, die zu erwarten sind. Zu berechnen ist der zukünftige Verlauf der
Zustandsgrößen. Die Vorwärtsrechnung liefert also die direkten Voraussagen, die sich
ableiten lassen, wenn ein Modell einmal etabliert ist.
Beobachtungsproblem: Die Systemeigenschaften sind bekannt, aber der Zustand des
Systems nicht. Dagegen kennen wir den Verlauf der Beobachtungsgrößen über ein
gewisses Zeitintervall. Daraus soll der gegenwärtige Zustand (oder vielleicht der Zustand
zu Beginn der Beobachtung) errechnet werden. Wenn der Zustand aus verrauschten,
unvollkommenen Beobachtungen geschätzt werden soll, reden wir auch von einem
Filterungsproblem.
Parameteridentifikation oder -anpassung: Die Struktur des Modells ist festgelegt, aber
die Parameter müssen erst bestimmt werden. Das geschieht natürlich so, daß reale Daten
mit den Daten verglichen werden, die das Modell liefern würde, und solange an den
Parametern gedreht wird, bis sich Modell und Meßdaten möglichst gut decken.
Parameteroptimierung: Die Struktur des Systems ist bekannt, und die Parameter können
innerhalb eines gewissen Bereichs eingestellt werden (durch geeignete Konstruktion des
Systems, geeignete Materialen, oder Einstellen von Regeleinrichtungen). Die Parameter
sollen so eingestellt werden, daß das System ein möglichst günstiges Verhalten zeigt (zum
Beispiel möglichst stabil ist, oder möglichst kostengünstig operiert). Obwohl die
Zielsetzung in der Anwendung völlig anders ist, ist dieses Problem mathematisch mit der
Parameteridentifizierung sehr eng verwandt.
Optimale Steuerung: Die Systemeigenschaften sind bekannt. Über eine oder mehrere
Steuergrößen soll das System so beeinflußt werden, daß es ein möglichst günstiges
Verhalten zeigt (zum Beispiel einen gewünschten Zustand in möglichst kurzer Zeit oder
möglichst kostengünstig erreicht und daran festhält).
Qualitative Untersuchung: Nicht die Zahlenwerte der Lösung, sondern charakteristische
Eigenschaften stehen im Mittelpunkt. Stellt sich im Lauf der Zeit eine Gleichgewichtslage
ein? Erholt sich das Gleichgewicht, wenn man es stört? Kann das System periodische
Schwingungen ausführen?
7.4. Parameteridentifikation.
Als Beispiel einer komplexeren Fragestellung behandeln wir die Parameteridentifikation.
Wir werden sehen, daß die Vorwärtsrechnung hier als ein wichtiges Hilfsmittel auftritt.
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
75
Abbildung 7.2. Ablauf einer Parameteridentifikation
Abbildung 7.2 zeigt die typische Vorgangsweise bei einer Parameteridentifikation. Man
erfaßt Felddaten, also reale Daten über das modellierte System. Andererseits läßt man
mit verschiedenen Sätzen von Parametern Simulationen laufen. Wenn die Simulationen
die Felddaten ausreichend gut wiedergeben, nimmt man die Parameter an. Sonst ändert
man die Parameter systematisch und versucht weitere Simulationen. Natürlich kann es
vorkommen, daß sich gar keine Parameter finden, die die Felddaten genügend gut
nachvollziehen. In diesem Fall müssen die Modellgleichungen selbst verändert werden. Im
Ablauf des Identifikationsverfahrens wird die Vorwärtsrechnung oft und oft wiederholt.
Weil die Vorwärtsrechnung selbst im allgemeinen langwierig ist, ist eine
Parameteridentifikation mit sehr viel Rechenaufwand verbunden.
Wir nehmen nun wieder das Stickoxidmodell aus Beispiel 4.15, und werden anhand
dieses Beispiels eine Parameteridentifikation durchführen.
Beispiel 7.9. Wir modellieren die chemischen Reaktionen aus Beispiel 4.15 mit den
Differentialgleichungen
d
cD (t) = −k1 cD (t),
dt
d
cN O (t) = k1 cD (t) − k2 o2 cN O (t)2 ,
dt
cD (0) = c0 ,
cN O (0) = 0.
Die Größen k1 und k2 , die die Reaktionsgeschwindigkeiten bestimmen, sollen durch einen
Versuch in Erfahrung gebracht werden. Dazu setzen wir der Lösung zunächst
c0 = 10−5 mol/l Donor zu. Wir messen eine (im Versuchsverlauf annähernd konstante)
Sauerstoffkonzentration der Lösung von o2 = 2.5 · 10−4 mol/l. Alle 100 Sekunden wird die
Konzentration von Donor und NO gemessen, die Meßwerte finden sich in Tabelle 7.1. Zu
bestimmen sind die Parameter k1 und k2 . Erfahrungsgemäß dürften die
Größenordnungen dieser Parameter etwa k1 ≈ 10−3 und k2 ≈ 107 betragen.
76
Tabelle 7.1. Daten zu Beispiel 7.9
t
000
100
200
300
400
500
600
700
800
900
1000
1100
1200
cD
1.00 · 10−5
7.40 · 10−6
5.57 · 10−6
3.90 · 10−6
2.87 · 10−6
1.95 · 10−6
1.46 · 10−6
1.03 · 10−6
7.79 · 10−7
5.66 · 10−7
4.21 · 10−7
2.95 · 10−7
2.24 · 10−7
cN O
0.00 · 10−6
2.34 · 10−6
2.79 · 10−6
2.74 · 10−6
2.50 · 10−6
2.30 · 10−6
1.82 · 10−6
1.75 · 10−6
1.64 · 10−6
1.29 · 10−6
1.10 · 10−6
9.64 · 10−7
9.30 · 10−7
7.5. Parameteridentifikation mit freiem Auge.
Die Methode scheint im Computerzeitalter antiquiert zu sein, aber sie hat einen Vorteil.
Indem man selbst viele Parameter ausprobiert, sammelt man Erfahrung über die Rolle,
die die Parameter im System spielen. Bevor wir die Identifikation durchführen, müssen
wir die Simulation vorbereiten und die Daten einlesen.
In Matlab geschieht dies am besten in einem script, sagen wir dondata.m, das die Daten
unter tdat, dondat, nodat speichert und die Variablennamen c0, o2, k1, k2 global
setzt, sodass sie auch in funk und in parziel verwendet werden können. Der folgende
M-File führt eine Vorwärtsrechnung durch und plottet die Ergebniskurven und die
Meßdaten zum Vergleich. Im Zentrum der Rechnung steht wieder die MATLAB-Routine
ode45. Sie löst die Differentialgleichung x0 (t) = f (t, x(t)) mit Hilfe eines
Runge-Kutta-Verfahrens. Dabei ist die Funktion f durch einen MATLAB-Funktionsfile
gegeben (hier funk.m).
% Vorwaertsrechnung Donorproblem
% Name der Funktion: funk
% Anfangswerte: 0,c0
% Zeitintervall: 0-1200
% Genauigkeit: 1e-10
options=odeset(’AbsTol’,1e-10);
[t,x]=ode45(’funk’,[0,1200],[c0;0],options);
%Plotten und Beschriften:
plot(t,x(:,1),’--’,t,x(:,2),’-’,...
tdat,dondat,’*’,tdat,nodat,’o’)
xlabel(’t’),ylabel(’Donor: --,*; NO: -,o’)
title(sprintf(’k1=%8.4e, k2=%8.4e’, k1,k2))
%(Titelzeile druckt die verwendeten Parameter)
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
77
Abbildung 7.3. Simulation des Donorproblems: erster Versuch
x10 -5
1*
k1=1.0000e-003, k2=1.0000e+007
0.9
0.8
*
Donor: --,*; NO: -,o
0.7
0.6
*
0.5
0.4
*
0.3
o
o
o
*
o
0.2
o
*
o
o
*
0.1
0o
0
200
400
600
o
o
*
*
800
*
o
o
*
*
1000
o
*
1200
t
Simulation: Donor strichliert, NO durchgezogen.
Meßdaten: Donor Sternchen, NO Kreise.
Wir sichern diesen File als simu.m.
Mit diesen Werkzeugen können wir die Parameteridentifikation beginnen. Wir stehen
jetzt im Workspace und laden zuerst die Daten. Anschließend machen wir den ersten
Versuch, mit k1 = 10−3 und k2 = 107 . Wir geben diese Parameter im Workspace ein und
starten die Vorwärtsrechnung. Abbildung 7.3 zeigt das Ergebnis. Die Übereinstimmung
ist noch schlecht, obwohl zumindest die Form der Kurven ungefähr paßt:
>>dondata ←>>k1=1e-3; ←>>k2=1e7; ←>>simu ←Wir können jetzt die Auswirkungen der Parameter ausprobieren. Es ist gut, immer nur
einen Parameter zu verändern, damit man weiß, welchen Grund die Änderung der
Kurven hat, die man dann beobachtet. Anscheinend fällt die Donorkonzentration in der
Simulation zu langsam ab. Wir verdoppeln die Zerfallsgeschwindigkeit k1 :
>>k1=2e-3; ←>>simu ←und versuchen weiter, die Parameter zu verändern. Wie wirkt sich k1 auf die
Konzentrationskurven aus? Auf welche Kurve wirkt k2 ? Wenn man einige Erfahrung mit
den Parametern gesammelt hat, könnte man etwa auf folgende Simulation kommen:
>>k1=3e-3; ←>>k2=8e6; ←>>simu ←-
78
Abbildung 7.4. Simulation des Donorproblems: nach einiger Erfahrung
x10 -5
1*
k1=3.0000e-003, k2=8.0000e+006
0.9
0.8
*
Donor: --,*; NO: -,o
0.7
0.6
*
0.5
0.4
*
0.3
o
o
o
*
o
0.2
o
*
o
o
*
0.1
0o
0
200
400
600
o
o
*
*
800
*
o
o
*
*
1000
o
*
1200
t
Simulation: Donor strichliert, NO durchgezogen.
Meßdaten: Donor Sternchen, NO Kreise.
Abbildung 7.4 zeigt eine recht gute Anpassung an die Meßdaten.
Beobachten Sie, daß die Donorkurve nur durch k1 beeinflußt werden kann. Das hat einen
Grund: Die Differentialgleichung für cD enthält nur cD und k1 , es besteht kein Einfluß
von cN O auf die Donorkonzentration. Je größer k1 , desto schneller fällt die
Donorkonzentration ab, und desto schneller steigt zu Beginn die NO-Konzentration. Je
größer k2 , desto schneller fällt später wieder die NO-Konzentration ab. Wenn man diese
Zusammenhänge verstanden hat, ist es leicht, sich zu geeigneten Modellparametern
vorzutasten.
7.6. Kleinste Quadrate.
Wie für das empirische logistische Modell können wir auch die Parameter eines
dynamischen Systems an die Felddaten mit Hilfe der Methode der kleinsten Quadrate
anpassen. Die einzige Neuerung ist jetzt, daß jede Berechnung der Zielfunktion eine
Vorwärtsrechnung enthält.
Wir wiederholen noch einmal das Prinzip der kleinsten Quadrate:
(1) Gegeben ist ein Satz von Zeitpunkten t1 , t2 , · · · , tn und ein Satz von Meßdaten
x1 , x2 , · · · , xn . (Dabei können, wie in unserem Beispiel, die Meßdaten Vektoren
sein, also aus mehreren Komponenten bestehen.
(2) Wenn wir einen Satz von Parametern vorgeben, können wir durch
Vorwärtsrechnung bestimmen, welche Zustandswerte das Modell für die
Zeitpunkte t1 , · · · , tn voraussagt. Nennen wir diese simulierten Werte
x̃1 , x̃2 , · · · , x̃n .
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
79
(3) Der quadratische Fehler für diesen Satz von Parametern ist dann
E = |x1 − x̃1 |2 + |x2 − x̃2 |2 + · · · + |xn − x̃n |2 .
(Sind die Daten Vektoren, dann bezeichnet |a| die Norm des Vektors a.) Weil die
Simulationsdaten von den vorgegebenen Parametern abhängen, hängt auch der
quadratische Fehler von den Parametern ab. Der quadratische Fehler ist also
eine Funktion der Parameter.
(4) Durch ein Minimumsuchprogramm bestimmen wir jene Parameter, für die der
quadratische Fehler am kleinsten wird. Die Zielfunktion der Minimumaufgabe
ist also der quadratische Fehler.
Wir kennen bereits den Minimumsucher fminsearch in MATLAB. Die Zielfunktion muß
als MATLAB-Funktionsfile programmiert werden. Vorab ein technisches Detail: Der
Minimumsucher sucht einen Vektor von Parametern. Er arbeitet am besten, wenn alle
Parameter dieselbe Größenordnung haben. Wir packen daher die beiden Parameter
skaliert in einen Vektor:
µ
¶
k1 · 103
kk =
.
k2 · 10−7
Wir beginnen mit dem Editieren der Zielfunktion:
function fehl=parziel(kk);
%parziel: Zielfunktion Parameteridentifikation
global k1 k2 o2 c0 tdat dondat nodat
%Der Vektor kk enthaelt die beiden Parameter
%Auspacken der Argumente
k1=kk(1)*1e-3;
k2=kk(2)*1e7;
options=odeset(’AbsTol’,1e-10);
[t,x]=ode45(’funk’,tdat,[c0;0],options);
Jetzt bilden wir die Differenz der Meßdaten von den Simulationsdaten und errechnen
daraus den quadratischen Fehler:
%Berechnung des Fehlers
donfehl=x(:,1)-dondat;
nofehl=x(:,2)-nodat;
fehl=donfehl’*donfehl+nofehl’*nofehl;
Damit wäre die Zielfunktion fertig programmiert, denn ihre Aufgabe ist die Berechnung
des quadratischen Fehlers aus den Parametern. Nur damit wir dem Programm beim
Arbeiten zusehen können, programmieren wir auch eine graphische Ausgabe der
neuesten Daten. Das verlangsamt natürlich den Ablauf der Rechnung, und dient hier nur
zu didaktischen Zwecken.
%Beobachtung des Verlaufs
plot(t,x(:,1),’--’,t,x(:,2),’-’,...
tdat,dondat,’*’,tdat,nodat,’o’)
xlabel(’t’),ylabel(’Donor: --,*; NO: -,o’)
80
Abbildung 7.5. Simulation des Donorproblems: kleinste Quadrate
−5
1.2
k1=3.1379e−003, k2=7.7613e+006, fehler=1.6686e−013
x 10
1
Donor −−,*
NO −,o
0.8
0.6
0.4
0.2
0
0
200
400
600
t
800
1000
1200
Simulation: Donor strichliert, NO durchgezogen.
Meßdaten: Donor Sternchen, NO Kreise.
title(sprintf(...
’k1=%8.4e, k2=%8.4e, fehl=%8.4e’,k1,k2,fehl))
pause(0.1)
Wir sichern die Zielfunktion unter dem Namen parziel.m. Beachten Sie dass das zweite
Argument im Aufruf von ode45 nicht wie früher ein Vektor aus zwei Komponenten für
Start- und Endzeit ist, sondern mehr als 2 Komponenten besitzt, nämlich hier 13. Dies
bewirkt, dass ode45 in [t,x] 13 Zeilen erzeugt, genau eine für jeden Meßzeitpunkt
t = 0, . . . , 1200. Dies geschieht durch automatische Interpolation, während ode45 mit viel
kleinerer variabler Schrittweite arbeitet.
Wir führen jetzt im Workspace die Parameteridentifikation durch. Der Minimumsucher
benötigt den Namen der Zielfunktion, einen Startwert, und eventuell die Toleranz. Wir
haben schon durch die händische Identifikation ziemlich genaue Vorstellungen, wo die
Parameter liegen, nämlich k1 ≈ 3 · 10−3 , k2 ≈ 8 · 106 , das wäre nach der Skalierung:
kk ≈ (3, 0.8)T . Absichtlich geben wir einen ungünstigeren Startwert vor, um dem
Programm beim Arbeiten zuzusehen.
>>dondata ←>>kkstart=[1;1]; ←>>kk=fminsearch(’parziel’,kkstart) ←kk =
3.1379
0.7761
Wir beobachten auf dem Bildschirm, wie der Minimumsucher zunächst scheinbar ziellos
Parameter ausprobiert und bessere und schlechtere Simulationen erzeugt, aber nach
7. FRAGESTELLUNGEN ZU DYNAMISCHEN SYSTEMEN
81
Abbildung 7.6. Hierarchie der Programme
einiger Zeit geht er ziemlich gezielt auf immer bessere Kurven los. Wenn die Rechnung
abbricht, erhalten wir die optimalen Parameter. Vergessen wir nicht die Skalierung, die
wir vorgenommen haben. Wir erhalten k1 = 3.1379 · 10−3 und k2 = 0.7761 · 107 .
Abbildung 7.5 zeigt die Kurven der Simulation mit diesen Parametern im Vergleich zu
den Meßdaten. Abbildung 7.6 fasst die Struktur des Rechenvorgangs zusammen. Im
workspace werden mit dondata die Daten und die globalen Größen zugewiesen und
fminsearch aufgerufen, das die Minimumsuche durchführt und dabei parziel für
verschiedene Parameter k1 , k2 aufruft. parziel berechnet nach einer Simulation mit
ode45 die Summe der Fehler-Quadrate. Die von den rufenden Routinen übergebenen
Argumente stehen in der linken Spalte von Abbildung 7.6, die von den gerufenen
Routinen gelieferten Ergebnisse in der rechten. In der Rechnung zu Abbildung 7.5 wurde
parziel 84 mal aufgerufen, soviele Vorwärtsrechnungen wurden also durchgeführt.
Dabei wurde funk insgesamt 7644 mal ausgewertet; dies kann man durch Einbau
einfacher Zähler (etwa z=z+1) feststellen. Die gesamte Rechenzeit beträgt ca. 1 Sekunde
auf einem Standard PC.
An dieser Stelle sei an Bemerkung 6.8 erinnert, denn hier wurde eine Problemstellung
beschrieben, in der beide Verfahren zum Einsatz kommen. Trotzdem sind sie nicht
miteinander vermischt: die Minimumsuche steht in der obigen Hierarchie über dem
Differentialgleichungslöser und die beiden Verfahren arbeiten intern unabhängig.
Das Modell wurde somit erfolgreich an Daten angepasst, wobei allerdings alle hier
vorhandenen Daten verwendet wurden, sodass die Parameter k1 , k2 genau für diesen
Datensatz optimiert wurden. Es besteht die Hoffnung, dass diese k1 , k2 vom speziellen
Datensatz unabhängige Reaktionskonstanten sind. Dies müsste allerdings erst mit Hilfe
anderer, zusätzlicher Messungen untersucht werden, z.B. mit anderen Donorzugaben
(anderes c0 im Modell) oder anderem Sauerstoffgehalt (anderes o2 im Modell).
Grundsätzlich sollte man nicht alle verfügbaren Daten in die Parameteridentifikation
stecken, sondern einen unabhängigen Anteil davon für Kontroll-Zwecke aussparen.
82
8. Simulationsdiagramme
Wir geben in diesem Abschnitt eine Einführung in blockorientierte
Simulationsdiagramme. Simulation ist nichts anderes als numerische Lösung der
Modellgleichungen. Simulationssprachen sind Computerprogramme, die darauf ausgelegt
sind, daß man die Modellgleichungen auf möglichst anschauliche Art eingeben kann, und
mit den Techniken der numerischen Gleichungslösung so wenig wie möglich belastet
wird. Blockorientierte Simulation hat den Vorteil, daß das Modell auf der Basis einer
graphischen Darstellung entwickelt und eingegeben wird.
8.1. Wirkungsdiagramm.
Ein Wirkungsdiagramm ist eine Darstellung eines Systems und der Wechselwirkungen
seiner Komponenten. Dabei werden die zu modellierenden Zustandsgrößen und Prozesse
als Worte so angeordnet, dass die Wechselwirkungen mit Pfeilen symbolisiert werden
können. Dies ist ein nützlicher erster Schritt bei der Bildung eines mathematischen
Modells, unterstützt das Nachdenken über die Wirkungen im System und stellt eine
Vorstufe bei der Entwicklung eines Simulationsdiagramms dar. Insbesondere fällt schon
hier eine wichtige - vielleicht vorläufige - Entscheidung über die Auswahl der
Zustandsgrößen und die berücksichtigten oder weggelassenen Wechselwirkungen.
Beispiel 8.1. In Beispiel 4.15 wird einer wässerigen Lösung eine sogenannte
Donorsubstanz zugesetzt, die in einer Reaktion erster Ordnung zu NO zerfällt. In einer
Reaktion zweiter Ordnung wird NO oxidiert und dadurch aus der Lösung entfernt. Stellen
Sie die wechselseitige Beeinflussung der Systemgrößen Donorkonzentration,
Stickoxidkonzentration, Sauerstoffkonzentration, Donorzerfallsrate und
NO-Oxidationsrate graphisch dar.
Diese Vorgänge werden in Abbildung 8.1 graphisch zusammengefaßt: Drei Substanzen in
der Lösung stehen in Wechselwirkung, nämlich der Donor, NO und O2 . Die zwei
chemischen Vorgänge werden durch die Zerfallsrate des Donors und die Abbaurate von
NO beschrieben. Durch Pfeile bezeichnen wir die Art der Wechselwirkung: Je mehr
Donor vorhanden ist, desto mehr Donor zerfällt auch, daher wirkt sich der Donorbestand
positiv auf die Donorzerfallsrate aus. Wir zeichnen einen mit + markierten Pfeil vom
Donorbestand zur Donorzerfallsrate. Je mehr Donor zerfällt, desto mehr NO wird
gebildet, aber desto weniger Donor bleibt übrig. Die Donorzerfallsrate wirkt sich also
positiv auf den NO-Bestand und negativ auf den Donorbestand aus. Wir zeichnen einen
positiv markierten Pfeil vom Donorzerfall zum NO-Bestand, und einen negativ
markierten Pfeil vom Donorzerfall zum Donorbestand. Je mehr O2 und NO vorhanden
sind, desto mehr NO wird oxidiert. Der NO-Abbau wirkt sich negativ auf den Bestand
von NO und O2 aus. Weil wir annehmen, daß soviel Sauerstoff gelöst ist, daß der
geringfügige Verlust durch die NO-Oxidation nicht ins Gewicht fällt, markieren wir die
negative Auswirkung der NO-Oxidation auf den Sauerstoffbestand nur mit einem
strichlierten Pfeil, wir werden sie im Modell nicht berücksichtigen.
8. SIMULATIONSDIAGRAMME
83
Abbildung 8.1. Wirkungsdiagramm zu Beispiel 4.15
Die eben erstellte Abbildung ist ein Wirkungsdiagramm.
Merksatz 8.2. Ein Wirkungsdiagramm stellt durch Pfeile dar, welche Komponenten
eines Systems einander beeinflussen. Durch Vorzeichen wird dargestellt, ob der Einfluß
positiv oder negativ ist. Das Wirkungsdiagramm ist ein qualitativer Überblick über die
Struktur des Systems, es enthält keine quantitative Information.
8.2. Erstellen eines Simulationsdiagramms.
Wir entwickeln in diesem Abschnitt ein Simulationsdiagramm, das ist eine graphische
Beschreibung des Systems, die auch quantitative Information enthält, und zwar so
detailliert, daß sie sich direkt in Gleichungen umsetzen läßt. Wieder nehmen wir als
Beispiel das Stickoxidmodell.
Die Zustandsgrößen unseres Systems sind die Bestände an Donor und NO, sie
beschreiben zu jedem Zeitpunkt t den augenblicklichen Zustand der Lösung, alle anderen
Größen (etwa die Donorzerfallsrate und die Oxidationsrate) lassen sich daraus
berechnen. Den Sauerstoffbestand betrachten wir in unserem Modell als eine bekannte
Konstante. Wir zeichnen zwei Rechtecke, sogenannte Blöcke, für die Zustandsgrößen
(Abbildung 8.2). Unter die Blöcke schreiben wir die Namen der dargestellten Größen.
Wir stellen uns vor, daß am Ausgang der Blöcke die Zahlenwerte cD und cN O der Größen
abgelesen werden können. (Wie, das sehen wir später.) Am Eingang werden Einflüsse
stehen, die auf diese Größen einwirken.
Merksatz 8.3. Den Aufbau eines Simulationsdiagrammes beginnt man mit den Blöcken
für die Zustandsgrößen.
Wir modellieren jetzt die Wirkungsweise des Donorzerfalls (Abbildung 8.3): Die
Zerfallsrate DZR des Donors wird ebenfalls durch einen Block repräsentiert,
ausnahmsweise zeichnen wir diesen Block dreieckig — warum, das wird gleich erklärt.
Am Ausgang des Blocks steht also, wieviel Mol Donor in der Sekunde zerfallen. Weil der
Donorbestand auf die Zerfallsrate Einfluß nimmt, zeichnen wir einen Pfeil vom Block
Donor zum Eingang des Blocks DZR. Soweit sieht die Modellierung nur qualitativ aus.
Wir stellen aber auch den quantitativen Zusammenhang dar:
DZR(t) = k1 · cD (t).
Die Form des Blocks symbolisiert die Art der Gleichung: Ein dreieckiger Block
bezeichnet einen sogenannten Verstärker. Die Spitze des Dreiecks zeigt zum Ausgang.
84
Abbildung 8.2. Blöcke beschreiben Systemgrößen
Abbildung 8.3. Wirkungsweise des Donorzerfalls
Der Ausgang eines Verstärkers entsteht aus dem Eingang durch Multiplikation mit einer
bekannten Konstanten. Die Konstante, hier k1 , tragen wir in das Dreieck ein.
Der Donorzerfall selbst wirkt auf den NO-Bestand, das modellieren wir später, und auf
den Donorbestand. Weil er darauf negativ wirkt, drehen wir das Vorzeichen um.
Mathematisch gesehen ist das eine Multiplikation mit −1, kann also durch einen
Verstärker mit der Konstanten −1 dargestellt werden. Am Ende des zweiten Verstärkers
steht also die negative Donorzerfallsrate, und diese wirkt direkt auf den Donorbestand.
Wir führen also einen Pfeil vom zweiten Verstärker an den Eingang des Blocks
Donorbestand.
Beginnend mit einem Anfangswert für die Donorkonzentration integriert die Schleife in
Abb. 8.3 die negative Donorzerfallsrate. Das können wir wahlweise mit Hilfe der
Ableitung oder eines Integrals schreiben:
(8.1)
(8.2)
d
c (t)
dt D
= −DZR(t), cD (t0 ) = c0
Rt
cD (t) = c0 + t0 (−DZR(s)) ds.
t0 bezeichnet dabei den Anfangszeitpunkt der Simulation, den Zeitpunkt, zu dem Donor
in die Lösung eingebracht wurde. Diese beiden Zeilen sind gleichwertig in dem Sinn, dass
cD (t) dann und nur dann eine Lösung von (8.1) ist, wenn cD (t) auch die
Integralgleichung (8.2) erfüllt (Hauptsatz der Differential- und Integralrechnung).
Der Ausgang des Blocks Donor ist also das Integral des Eingangs plus Anfangswert
cD (t0 ) = c0 . Der Block Donor ist ein sogenannter Integrationsblock. Wir symbolisieren
das, indem wir in den Block das Symbol 1/s schreiben. (Die Schreibweise 1/s für einen
Integrationsblock können wir hier nicht erklären. Sie stammt aus der Regeltechnik.) Der
Integrationsblock beinhaltet einen Anfangswert, nämlich den Donorbestand c0 zu Beginn.
8. SIMULATIONSDIAGRAMME
85
Abbildung 8.4. Einige Typen von Blöcken
Beachten Sie, daß die Graphik sich direkt in Gleichungen umsetzen läßt. In
Abbildung 8.3 haben wir die folgenden Beziehungen symbolisiert:
DZR = k1 cD ,
Z t
cD (t) = c0 +
(−DZR(s)) ds.
t0
Merksatz 8.4. In einer blockorientierten Simulationssprache werden die
mathematischen Zusammenhänge innerhalb eines Systems anstatt durch Gleichungen
graphisch dargestellt. Dabei wird jede Systemgröße durch einen Block symbolisiert. Am
Ausgang des Blocks steht die dargestellte Größe. Pfeile zum Eingang des Blocks zeigen
die Größen an, die auf den Block Einfluß nehmen. Die Art des Blocks zeigt an, wie die
Ausgangsgröße aus den Eingangsgrößen errechnet wird.
Bevor wir das Simulationsdiagramm weiterentwickeln, stellen wir in Abbildung 8.4 die
Blocktypen zusammen, die wir noch brauchen werden. Integrationsblock und Verstärker
wurden bereits beschrieben. Ein Additionsblock hat zwei oder mehrere Eingänge, die
addiert oder subtrahiert werden. Jeder Eingang hat ein Vorzeichen, das festlegt, ob er
positiv oder negativ zur Summe beiträgt. Auch ein Multiplikationsblock hat mehrere
Eingänge. Am Ausgang steht das Produkt der Eingänge. Eine feste Konstante wird
durch ein kleines Rechteck symbolisiert. Weil die Konstante eben von keinen Einflüssen
verändert wird, hat dieser Block keine Eingänge. Der Oszillographenblock stellt keine
Gleichung dar: Er ist die Anweisung, während der Simulation die Kurve derjenigen
Größe zu zeichnen, die am Eingang hängt.
Wir stellen jetzt das Simulationsdiagramm fertig. Zwei Größen beeinflussen den
NO-Bestand, nämlich die Donorzerfallsrate (positiv), und die NO-Oxidationsrate
(negativ). An den Eingang des Blocks NO hängen wir daher einen Additionsblock, der
86
Abbildung 8.5. Simulationsdiagramm des Stickoxidmodells
die Oxidationsrate von der Donorzerfallsrate subtrahiert. Die Donorzerfallsrate ist
bereits modelliert und wird an den positiven Eingang des Addierers angeschlossen. Der
negative Eingang wird für die Oxidationsrate reserviert. Die Differenz der beiden
Reaktionsraten ist die Wachstumsrate des NO-Bestandes. Wieder steht am Eingang des
Blocks NO die Wachstumsrate, also die Ableitung: Auch NO ist ein Integrationsblock,
was wir durch das Symbol 1/s darstellen. Der Anfangswert ist 0, zu Beginn befindet sich
kein Stickoxid in der Lösung.
Beachten Sie ein wichtiges Detail: Die Größe DZR wurde mehrfach weiterverarbeitet, die
Ausgangsleitung von DZR wurde sozusagen angezapft. Wo sich die Leitungen
verzweigen, haben wir einen deutlichen Punkt eingezeichnet. In komplizierten
Simulationsdiagrammen ist man gezwungen, Leitungen zu überkreuzen, die miteinander
nichts zu tun haben. Der Punkt zeigt an, daß nicht eine Kreuzung von zwei
Abbildung 8.6. Kreuzungen und verbundene Leitungen
8. SIMULATIONSDIAGRAMME
87
unabhängigen Leitungen aus zeichentechnischen Gründen erfolgt ist, sondern tatsächlich
eine Leitung in mehrere Richtungen verzweigt wurde (Abbildung 8.6).
Die NO-Oxidation ist eine Reaktion zweiter Ordnung nach NO. Die Oxidationsrate ist
das Produkt von vier Faktoren
N OR = k2 o2 c2N O = k2 o2 cN O cN O .
Dabei sind k2 und o2 zwei bekannte Konstanten, die wir durch Konstantenblöcke
symbolisieren. NOR selbst ist dann ein Multiplikationsblock mit vier Eingängen.
Damit sind die mathematischen Beziehungen zwischen den Systemgrößen fertig
beschrieben:
Z t
cD (t) = c0 +
(−DZR(s)) ds,
0
DZR(t) = k1 cD (t),
Z t
cN O (t) = 0 +
(DZR(s) − N OR(s)) ds,
0
N OR(t) = k2 o2 c2N O (t).
und das Simulationsdiagramm in Abb. 8.5 ist die blockorientierte Darstellung des
Differentialgleichungsmodells (4.29),(4.30). Wir schließen noch zwei Oszillographenblöcke
an: In der Simulation wollen wir die Kurven der Donorkonzentration und der
NO-Konzentration beobachten.
Schreiben Sie, zur Übung, ein kleines System von Differentialgleichungen Ihrer Wahl an
und zeichnen Sie das zugehörige Simulationsdiagramm auf. An die Eingänge der
Integratorblöcke werden die Wachstumsraten angeschlossen. Diese werden gemäß den
rechten Seiten der Differentialgleichung gezeichnet, wobei für die Zustandsgrößen die
Ausgänge der Integratorblöcke verwendet werden. Damit haben Sie Ihr System - im Sinn
von Simulationen mit geeigneter software - auch schon gelöst.
8.3. Implementieren eines Simulationsdiagramms in SIMULINK.
Bisher kennen wir das Simulationsdiagramm nur als eine graphische Alternative zur
Darstellung mathematischer Beziehungen. Nun wird demonstriert, wie man in Simulink,
das ist ein Teilpaket von Matlab, ein Simulationsdiagramm eingeben und dann damit
Simulationen durchführen kann. Anders als im Workspace gibt man in erster Linie nicht
Formeln, sondern Graphik ein. Als Ergebnis der Simulation erhält man die Kurven der
Lösungen des Modells. Beachten Sie, dass die Simulationen dieses Abschnitts in sich
selbstständig sind und nicht auf Matlab-Programmen früherer Kapitel aufbauen.
Die Eingabe eines Simulationsdiagramms in Simulink umfaßt die folgenden 8 Schritte:
(1) Öffnen eines Fensters für das Simulationsdiagramm.
(2) Beschaffung der benötigten Blöcke aus der SIMULINK-Blockbibliothek.
(3) Anordnung der Blöcke.
88
(4)
(5)
(6)
(7)
(8)
Eintragen der Verbindungen.
Zurechtrücken und Feinkorrekturen.
Eintragen der Parameter und Anfangswerte.
Einstellen der Oszilloskope.
Speichern des Simulationsdiagramms.
Schritt (1) geschieht durch das Matlab-Kommando simulink oder durch Klicken des
Simulink Symbols. Die Blockbibliothek ist ein Fenster, das die Blöcke (ohne
Verbindungsleitungen) enthält. Mit der Maus zieht man daraus die benötigten Blöcke ins
Diagramm-Fenster. Die Schritte (3) und (5) werden auch mit der Maus und/oder mit
Tasten-Kombinationen durchgeführt. Die Verbindungsleitungen zeichnet man ebenfalls
mit der Maus ein, dies erfordert etwas handwerkliches Geschick und Übung. Zur
Festlegung der Block-Parameter und der Anfangswerte bedient man sich des pop-up
Formulars, das man durch Doppelklicken des Blocks erhält. Die Oszilloskope kann man
(nach der Simulation) auch automatisch skalieren lassen. Abbildung 8.5 zeigt das fertige
Simulink-Diagramm zum Stickoxidmodell.
Wir wollen jetzt Beispiele rechnen mit dem Modell, das durch das Simulationsdiagramm
dargestellt wird. Dazu müssen wir dem Programm noch mitteilen, welchen Zeitraum die
Simulation durchschreiten soll, wie genau und mit welchem Verfahren es rechnen soll.
Dazu öffnen wir in der Menuzeile des Fensters, in dem unser Simulationsdiagramm steht,
unter dem Menupunkt Simulation die Option Parameter. Es öffnet sich ein Formular
mit mehreren Zeilen für die Verfahrensparameter. (Die Verfahrensparameter steuern das
numerische Verfahren und haben mit den Modellparametern nichts zu tun.)
Zunächst wählen wir durch Anklicken die Lösungsmethode. In Frage kommt etwa das als
default angebotene Variable-step Runge-Kutta-Verfahren ode45. Weil der Computer die
Gleichungen nur näherungsweise lösen kann, braucht er ein geeignetes
Näherungsverfahren, und verschiedene Verfahren bewähren sich bei verschiedenen
Gleichungen.
Start Time und Stop Time umspannen den Simulationszeitraum. Wir wählen Start 0
und Stop 2000, das sind 2000 Sekunden.
Differentialgleichungslöser tasten sich in kleinen Schritten die Lösungskurve entlang (vgl.
Seite 59 f). Je kürzer die Schritte, desto genauer können sie der Kurve folgen, aber um
den Preis von viel Rechenaufwand. Durch die Aufhäufung von vielen Rundungsfehlern
geht wiederum Genauigkeit verloren, wenn die Schrittweite allzu kurz ist. Außerdem
wird bei kurzen Schritten die Rechnung sehr langsam. Weil das gesamte Zeitintervall
2000 Sekunden lang ist, spezifizieren wir als größte erlaubte Schrittweite etwa 10, als
kleinste 1. Zum Vergleich können wir auch die step sizes auf auto setzen.
Als Rechengenauigkeit verlangen wir nicht allzu viel. Im Hinblick darauf, daß die
Konzentrationen in der Größenordnung von 10−5 liegen, verlangen wir 1e-8 absolute
8. SIMULATIONSDIAGRAMME
89
Abbildung 8.7. Simulationsergebnisse des Stickoxidmodells
x10 -5
1
0.9
0.9
0.8
0.8
0.7
0.7
Konzentration von NO
Konzentration des Donors
x10 -5
1
0.6
0.5
0.4
0.3
0.6
0.5
0.4
0.3
0.2
0.2
0.1
0.1
0
0
500
1000
Zeit t
1500
2000
0
0
500
1000
1500
2000
Zeit t
tolerance. Das Programm ist so ausgelegt, daß die berechnete Lösung von der
tatsächlichen Lösung um ungefähr diese Größenordnung abweichen darf.
Wir haben jetzt alle Verfahrensparameter eingegeben. Wir schließen das zugehörige
Formular und starten die Simulation durch Klicken der Start-Taste. Eine laufende
Simulation kann man durch Anklicken von Stop abbrechen, die Rechenzeit für dieses
Beispiel beträgt aber nur Bruchteile einer Sekunde.
Das Lösungsverfahren arbeitet also in diskreten Zeitschritten, die wie eine tickende Uhr
wirken. Obwohl es nicht tatsächlich so ist, kann man sich vorstellen, dass diskrete Signale
als Pulse im Takt der Uhr durch die Leitungen laufen. Startwerte sind die Anfangswerte
in den Integratorblöcken. Die Amplitude eines Pulses entspricht dem Wert der
zugehörigen (Zwischen-)Größe. In den Blöcken werden die Pulse transformiert bzw.
miteinander verknüpft. In den Oszilloskopen werden die diskreten Simulations-Zeitpunkte
waagrecht, die zugehörigen Amplituden der Eingangspulse senkrecht aufgetragen (weil
die Schritte klein sind, sieht die durchgezogene Verbindungslinie glatt aus). Gerade die
Anschaulichkeit dieser Vorstellung erleichtert die Interpretation eines Diagramms.
Ausserdem erleichtern die Block-Bibliothek sowie die Anschaulichkeit der graphischen
Programmierung das Erlernen und Verwenden eines derartigen Simulationswerkzeugs.
Abbildung 8.7 zeigt die Kurven der Donor- und NO-Konzentration, wie sie für c0=1e-5,
o2=2.5e-4, k1=5e-3, k2=1e7 auf den Oszilloskopen entstehen. Allerdings müssen dazu
die Achsen der Oszilloskope geeignet skaliert sein, was man am einfachsten durch
Autoscale erreicht. Der Verlauf der Kurven ist sehr plausibel: Der Donor unterliegt nur
einem exponentiellen Zerfallsgesetz, hat am Anfang den vorgegebenen Wert und
verschwindet langsam aus der Lösung. Zunächst befindet sich kein NO in der Lösung,
durch den anfänglich raschen Zerfall von Donor steigt die Konzentration von NO an.
90
Damit steigt aber auch die Oxidationsrate, und der Nachschub durch den Donor klingt
allmählich ab. Daher fällt nach einem Maximum auch die NO-Konzentration gegen Null
ab.
8.4. Ein Infektionsmodell.
Beispiel 8.5. Eine Population entwickelt sich auf der Basis von Zuwanderung und Tod.
In dieser Population bricht eine Infektionskrankheit aus. Gesucht ist ein einfaches
Simulationsmodell, an dem demonstriert werden kann, wie sich die Infektion innerhalb
der Bevölkerung entwickeln könnte, ob sie ausstirbt, oder ob sie sich endemisch einnistet,
sodaß auf lange Sicht immer ein gewisser Anteil der Bevölkerung erkrankt sein wird. Wir
treffen dabei die folgenden qualitativen Annahmen:
(1) Ansteckung erfolgt durch Kontakt infizierter Personen mit noch nicht infizierten.
Sobald eine Person infiziert ist, ist sie in der Lage, die Krankheit zu übertragen.
(2) Die Krankheit ist soweit ungefährlich, als die Mortalität von infizierten und
gesunden Personen gleich ist.
(3) Die Anzahl der Infizierten nimmt keinen Einfluß auf die Anzahl der
Zuwanderungen und auf die Heilungschancen der einzelnen Infizierten.
(4) Infizierte Personen können geheilt werden, und sind nach der Heilung völlig und
auf immer immun gegen die Krankheit.
(5) Das Infektionsrisiko ist für alle Personen der Population gleich.
(6) Neu zugewanderte Personen sind weder infiziert noch immun.
Stellen Sie diese Vorgänge durch ein Wirkungsdiagramm graphisch dar.
Wir bemerken, daß diese Annahmen natürlich eine sehr grobe Vereinfachung jedes realen
Infektionsgeschehens sind. Überlegen Sie selber, welche Phänomene dazu beitragen, daß
die Annahmen 1–6 unrealistisch sein könnten. Die Aufgabe des Modells, das wir
entwickeln werden, soll nicht sein, den Verlauf einer einzelnen Infektion präzise
vorauszusagen, sondern wir wollen uns nur einen Überblick verschaffen, welchen Verlauf
eine Infektion generell nehmen könnte.
Da die Population nach Punkt (5) homogen ist, müssen wir nur drei Klassen von
Personen unterscheiden:
S Gesunde noch nicht immune, also noch anfällige Personen (in englischer
Terminologie Susceptible).
I Infizierte Personen, sie sind zugleich Krankheitsüberträger (englisch Infective).
R Immune Personen, die dadurch aus dem Wirkungskreis der Infektion entfernt
sind (englisch Removed).
Ein Modell, das auf der Basis dieser drei Populationsklassen aufbaut, heißt ein
SIR-Modell.
Abbildung 8.8 zeigt das Wirkungsdiagramm. Im Zentrum stehen die Bestände der drei
Klassen S, I und R. Die Populationsdynamik der drei Klassen wird durch folgende
Vorgänge geregelt:
8. SIMULATIONSDIAGRAMME
91
Abbildung 8.8. Wirkungsdiagramm zu Beispiel 8.5
Jedes Individuum jeder Klasse unterliegt einer gewissen Sterbewahrscheinlichkeit. Die
Anzahl der täglichen Verluste jeder Klasse durch Tod steigt, je mehr Individuen in der
Klasse sind. Diese Verluste sind natürlich dem Bestand der Klasse abträglich. Wir haben
daher eine positive Wirkung vom Bestand der Klasse zur Häufigkeit der Todesfälle aus
der Klasse, und eine negative Wirkung von der Todesrate auf den Bestand.
Die Zuwanderung bewirkt einen Zuwachs der Susceptible-Klasse. Sie wird selbst durch
keinen Vorgang des Infektionsgeschehens beeinflußt.
Je mehr Anfällige und je mehr Infizierte in der Population vorhanden sind, desto eher
kommt es zu Kontakt zwischen einem Anfälligen und einem Infizierten, und desto
häufiger tritt auch Ansteckung auf. Die Häufigkeit der Ansteckungen wird daher sowohl
vom Bestand von S als auch von I positiv beeinflußt. Jede Ansteckung befällt eine
Person aus der Klasse S und führt sie in die Klasse I über: Die Ansteckungshäufigkeit
wirkt sich auf den Bestand von S negativ, auf I positiv aus.
Je mehr Personen infiziert sind, desto mehr können auch geheilt werden. Der Bestand
von I wirkt sich positiv auf die Häufigkeit der Heilungen aus. Jedes geheilte Individuum
wird der Klasse I entzogen und als immun der Klasse R zugezählt.
Beispiel 8.6. Entwickeln Sie ein Simulationsmodell, das mit möglichst einfachen
mathematischen Beziehungen die oben beschriebenen Vorgänge im SIR-Modell
beschreiben kann. Führen Sie verschiedene Simulationen durch, um eine Vorstellung von
den Möglichkeiten zu bekommen, wie sich eine Infektion in einer Bevölkerung ausbreiten
kann.
Wir beginnen unser Simulationsdiagramm mit den Beständen S(t), I(t), und R(t). Diese
Funktionen geben an, wieviel Individuen die drei Klassen zur Zeit t enthalten.
Zuwanderung und Tod, sowie Ansteckung und Heilung tragen zu den Wachstumsraten
der einzelnen Klassen bei. Wir modellieren jede Klasse durch einen Integrationsblock,
dem ein Additionsblock vorgelagert ist, in dem die einzelnen Beiträge zur
Wachstumsrate aufsummiert werden.
In Abbildung 8.9 modellieren wir zunächst Zuwanderung und Tod. Da die
Zuwanderungsrate (Zuwanderungen pro Tag) unbeeinflußt von allen anderen Vorgängen
bleibt, ist sie eine Konstante β > 0. Sie wird zur Wachstumsrate der S-Klasse zugezählt.
92
Abbildung 8.9. Die Wirkung von Zuwanderung und Tod in Beispiel 8.5
Jede Klasse hat ihre eigene Häufigkeit der Todesfälle. Wir nehmen an, daß jedes
Individuum derselben konstanten Sterbewahrscheinlichkeit pro Tag unterliegt, diese
Wahrscheinlichkeit bezeichnen wir mit µ > 0. Die Anzahl der täglichen Todesfälle aus
der Klasse S ist daher µS und wird aus S durch einen Verstärker errechnet. Sie trägt
negativ zur Wachstumsrate von S bei. Dieselbe Konstruktion wird auch für die Klassen I
und R durchgeführt.
Wir modellieren jetzt die Ansteckung. Die einfachste mathematische Beziehung, die die
qualitativen Eigenschaften der Infektion nachahmt, beschreibt die Anzahl der täglichen
Ansteckungen durch λIS mit einem konstanten Faktor λ > 0. Ein solches Modell ist
gerechtfertigt, wenn man davon ausgeht, daß die Anzahl der Kontakte zwischen
Infizierten und Anfälligen proportional sowohl der Anzahl der Infizierten als auch der
Anfälligen ist. Im Simulationsdiagramm errechnen wir das Produkt IS durch einen
Multiplikationsblock. Der anschließende Verstärker trägt den Faktor λ bei. Die
Häufigkeit der Ansteckungen trägt negativ zur Wachstumsrate von S und positiv zur
Wachstumsrate von I bei.
Wir modellieren die Heilung unter der Annahme, daß für jede infizierte Person die
Wahrscheinlichkeit, innerhalb eines Tages geheilt zu werden, γ > 0 beträgt. Die Anzahl
der täglichen Heilungen ist dann γI, was durch einen Verstärker modelliert wird. Die
Heilungen tragen negativ zum Wachstum von I und positiv zum Wachstum von R bei.
Damit sind die mathematischen Beziehungen fertig modelliert. Abbildung 8.10 enthält
das zugehörige Simulationsdiagramm (in Simulink).
Wir müssen noch Ausgabeblöcke festlegen, die die Simulationsergebnisse anzeigen, und
die Parameter in die Blöcke eintragen. Wir wollen an Hand dieses Modells
demonstrieren, wie eng Simulink mit dem Matlab-Workspace zusammenarbeitet. Anstatt
die Kurven auf Oszilloskopen anzuzeigen, schreiben wir die Resultate in den Workspace.
In der Bibliothek befindet sich der Block “to Workspace” (ein längliches Rechteck);
dieser Block schreibt einfach alle Daten, die an seinen Eingang gelangen, in einen langen
Spaltenvektor im Workspace. Den Namen des Vektors können wir festlegen, indem wir
den Block öffnen. Wir schreiben die Bestände aller drei Klassen in den Workspace, als
Vektoren susc, infec und remov. Diese Daten können wir nach der Simulation im
Workspace weiterverarbeiten. Damit wir wissen, zu welchen Zeitpunkten die Daten
8. SIMULATIONSDIAGRAMME
93
Abbildung 8.10. Simulationsdiagramm des SIR-Modells
gelten, schreiben wir auch die Simulationszeit als time in den Workspace. Die diskreten
Simulationszeitpunkte erhält man aus dem Block Clock, der wie eine Uhr aussieht.
Neuere Versionen von Simulink speichern die Simulationszeitpunkte auch ohne Clock in
einem Vektor (default tout).
Auch bei der Festlegung der Parameter in den Blöcken verwenden wir das enge
Wechselspiel zwischen dem Workspace und dem Simulationsdiagramm.
Merksatz 8.7. Ein Parameter in einem SIMULINK-Block kann durch einen
Zahlenwert festgelegt werden, aber auch durch jede Formel, die der Workspace zu Beginn
der Simulation ausrechnen kann.
Wir schreiben daher in die Verstärker- und Konstantenblöcke einfach die Ausdrücke
lambda, mu, gamma, beta und müssen nur bei Beginn der Simulation sicherstellen, daß die
Zahlenwerte dafür im Workspace aufliegen. Das ist so gut und so schlecht wie die
Eingabe der Zahlenwerte in die Blöcke selbst, nur für µ ergibt sich ein echter Vorteil:
Wenn wir im Lauf der Simulationen verschiedene Werte für µ probieren wollen, müssen
wir nur den Wert im Workspace ändern und nicht drei Blöcke aufmachen.
Als Anfangswert für S geben wir 100000 im Integrationsblock S ein. (Das ist der Wert,
der sich mit den Parametern, die wir verwenden werden, ohne Infektion in der
Bevölkerung als Gleichgewicht einstellen würde.) Als Anfangswert für I geben wir 100
ein, dh. 100 Personen werden unvermittelt krank. Wir beginnen mit 0 Individuen in der
Klasse R.
94
Wir führen jetzt eine Simulation durch. In den Workspace schreiben wir die
Wertzuweisungen
>>beta=100; ←>>mu=0.001; ←>>gamma=0.4; ←>>lambda=5e-6; ←Als Verfahrensparameter (einzusetzen im Menupunkt Simulation, Option Parameters)
verwenden wir:
Runge-Kutta Verfahren,
Startzeit 0,
Stopzeit 4000,
kleinste Schrittweite 1,
größte Schrittweite 10,
Genauigkeit 1.
Wir sollten später auch ausgiebig mit Ergebnissen bei anderen Verfahrensparametern
(z.B. default) vergleichen, um Vertrauen (oder Misstrauen) in die Simulationsergebnisse
zu bekommen. Wir starten die Simulation und hören dann einen kurzen Piepton.
Jetzt liegen die Resultate der Simulation im Workspace, und wir können sie ansehen und
weiterverarbeiten. Von Interesse ist zum Beispiel, wie sich die Anzahl der Infizierten im
Lauf der Zeit entwickelt:
>>plot(time,infec) ←Die Kurve zeigt zunächst einen rapiden Anstieg der Infizierten bis zu einem Maximum,
und im folgenden ein langsames Auf- und Abschwingen der Anzahl der Infizierten mit
langsam abnehmender Amplitude.
Wir könnten auch die Anzahl der Infizierten gegen die Anzahl der Anfälligen auftragen:
>>plot(susc,infec) ←Wir erhalten die linke Kurve in Abbildung 8.11.
Wir ändern jetzt die Ansteckungsrate und betrachten eine weniger ansteckende
Krankheit. Dazu verändern wir λ im Workspace:
>>lambda=3e-6; ←Die anderen Parameter lassen wir unverändert und starten eine neue Simulation. Nach
Beenden des Laufes können wir die Resultate wieder plotten:
>>plot(susc,infec) ←Wir erhalten die rechte Kurve in Abbildung 8.11. Auf diese Weise könnten wir die
verschiedensten Kombinationen von Parametern durchprobieren.
Beispiel 8.8. Interpretieren Sie die Simulationsresultate für die Parameter
β = 100, µ = 0.001, γ = 0.4, λ = 5 · 10−6 bzw. λ = 3 · 10−6 ,
die in Abbildung 8.11 dargestellt sind.
8. SIMULATIONSDIAGRAMME
95
Abbildung 8.11. Simulationsergebnisse des SIR-Modells
lambda = 5e-6, 4000 Tage
lambda = 3e-6; 4000 Tage
2500
10
9
2000
8
7
Infective
6
Infective
1500
1000
5
4
3
500
2
1
0
6
o
7
8
Susceptible
9
*
10
x10 4
0
9.997
9.998
9.999
Susceptible
*
10
x10 4
Parameter:β = 100, µ = 0.001, γ = 0.4, λ = 5 · 10−6 bzw. λ = 3 · 10−6 , Simulationszeitraum: 4000.
Wir bemerken, daß die Zahlenwerte des Modells nicht realistisch gewählt wurden,
sondern so, daß das Verhalten, das ein SIR-Modell zeigen kann, besonders deutlich
sichtbar wird.
Die Kurven liefern einen Vergleich der Ergebnisse in Abhängigkeit vom Parameter λ, der
mißt, wie ansteckend die Krankheit innerhalb der Population wirkt. Alle anderen
Parameter sind in beiden Simulationen gleich. In den Parameter λ gehen sowohl
Eigenschaften der Krankheit als auch die hygienischen Verhältnisse ein. (Masern sind
zum Beispiel besonders ansteckend, sodaß λ für diese Krankheit relativ groß ist. Durch
sorgfältige Quarantäne von Infizierten kann man λ klein halten.)
Ist λ genügend klein (rechte Kurve), so stirbt die Infektion aus: Die Neuansteckungen
werden von den Heilungen überwogen. Zunächst nimmt die Anzahl der Anfälligen wegen
der Ansteckung ab. Aber auch die Anzahl der Infizierten fällt rasch, sodaß die Infektion
ausstirbt, und sich letztlich die Anzahl der Anfälligen jenem Wert nähert, der das
Gleichgewicht in der völlig gesunden Population darstellt (mit einem kleinen * rechts
unten eingezeichnet.)
Ist λ groß (linke Kurve), überwiegt zunächst die Ansteckung die Heilung. Die Anzahl der
Infizierten nimmt schnell zu, und durch die Ansteckung nimmt die Anzahl der Anfälligen
ab. Letztlich fehlt der Nachschub an Gesunden, die angesteckt werden könnten. In dieser
96
Tabelle 8.1. SIR-Modell zur Ausbreitung einer Infektion
Größe
Einheit
t
Tag
S(t)
Ind.
I(t)
Ind.
R(t)
Ind.
β
Ind/Tag
µ
1/Tag
γ
1/Tag
λ
1/(Tag.Ind)
Modellgrößen
Benennung
Zeit
Anzahl Anfällige
Anzahl Infizierte
Anzahl Immune
Zuwanderungsrate
Mortalität
Heilungswahrscheinlichkeit
Maß für Ansteckung
Kommentar
gesucht
gesucht
gesucht
bekannt
bekannt
bekannt
bekannt
dynamische Mengenbilanzen
d
S(t) = β − µS(t) − λS(t)I(t),
dt
d
I(t) = λS(t)I(t) − µI(t) − γI(t),
dt
d
R(t) = γI(t) − µR(t).
dt
Phase der Epidemie fällt sowohl die Anzahl der Anfälligen als auch die der Infizierten.
Die Infektion scheint fast auszusterben. Anschließend steigt langsam die Anzahl der
Anfälligen wieder an. Wenn aber genug Nachschub an Gesunden besteht, die angesteckt
werden können, bricht die nächste Welle der Epidemie aus, nicht ganz so heftig wie die
erste. Dieser Zyklus wiederholt sich in immer weniger ausgeprägten Wellen. Wir können
uns vorstellen, daß sich auf die Länge ein Gleichgewicht einstellt, indem stets ein
bestimmter Prozentsatz der Population infiziert ist: Die Infektion ist endemisch. (Das
endemische Gleichgewicht ist durch einen Kreis angedeutet.)
Natürlich könnte man jetzt ähnliche Untersuchungen über die Auswirkung der anderen
Parameter anstellen. Man könnte etwa durch Simulationen herausfinden, wie groß λ
mindestens sein muß, daß bei den gegebenen anderen Parametern die Infektion nicht
ausstirbt (nämlich λ = 4.01 · 10−6 , wie wir auch durch Untersuchung der Gleichgewichte
finden werden).
Anmerkung. Wir haben uns die Klasse R als Geheilte gedacht und nicht, was ja auch
mit removed gemeint sein könnte, als auf Grund der Infektionskrankheit Verstorbene.
Aber auch die erste Variante hat etwas Unmenschliches an sich: wenn die Geheilten
immun sind, dann sollten doch vor allem die Geheilten bei der Bekämpfung der Epidemie
mitwirken! Wie könnten wir die Heilungsrate γ in diesem Sinn von R abhängigig
machen? Und wie würde das entsprechende Simulationsdiagramm ausschauen?
8. SIMULATIONSDIAGRAMME
97
Es wird Ihnen sicher aufgefallen sein, dass wir Simulationen durchführen konnten, ohne
das mathematische Modell anzuschreiben. Dies ist deshalb möglich, weil das
Simulationsdiagramm inklusive festgelegter Anfangswerte und Modell-Parameter
eindeutig und detailliert genug ist, um von Simulink als ein System von
Integralgleichungen interpretiert zu werden. Insbesondere sind an die Eingänge der
Integratorblöcke die Änderungsraten der Zustandsgrößen angeschlossen (vgl.
(8.2), (8.1)). Und diese wurden modelliert, ohne alle Formeln explizit anzuschreiben.
Wir können das Diagramm als System von Differentialgleichungen interpretieren und
erhalten die Gleichungen in Tabelle 8.1. Auf der linken Seite eines Gleichheitszeichens
steht die Ableitung einer Zustandsgröße, rechts davon der Eingang des zugehörigen
Integratorblocks, nämlich die Zusammenfassung der Additions- und
Multiplikationsblöcke, in die S, I, R und β eingehen.
Es ist immer vorteilhaft, mit beiden Formen des Modells zu arbeiten, schon zur
gegenseitigen Kontrolle. Die Differentialgleichungs-Variante ist zwar nicht direkt für
Simulationen zu gebrauchen, sie ist aber kompakter, übersichtlicher, schneller
anzuschreiben und kann - wie wir sehen werden - einer qualitativen Analyse unterzogen
werden.
8.5. Trajektorien im Zustandsraum.
In diesem Abschnitt besprechen wir Konzepte, die von der speziellen Simulationsmethode
völlig unabhängig sind. Diese Konzepte haben also nichts mit Simulationsdiagrammen zu
tun, es ist hier aber eine geeignete Stelle dieses wichtige Thema zu besprechen.
Wir können die zeitliche Entwicklung eines Systems graphisch darstellen, indem wir jede
Zustandsgröße als Funktion der Zeit zeichnen. Im Donor-Problem hat man z. B. die
beiden Konzentrationen cD (t) und cN O (t), t ≥ 0, also zwei Zeichnungen, wie in
Abbildung 8.7. Aus diesen beiden Kurven können wir eine einzige machen, wenn wir auf
der waagrechten Achse nicht die Zeit sondern die erste Zustandskomponente cD (t)
auftragen und auf der senkrechten die zweite, cN O (t). Wir zeichnen also die Kurve
{(cD (t), cN O (t)), t ≥ 0}. Das Ergebnis ist in Abbildung 8.12 zu sehen. (In Matlab
würde man dazu in Abschnitt 6.4 einfach das Kommando plot(x(:,1),x(:,2))
eingeben.) Diese Abbildung ist wie folgt zu interpretieren: jeder Punkt der Kurve ist ein
Zustand des Modells, der im Laufe der Simulation angenommen wurde. Und zwar ist der
Anfangszustand der Punkt (10−5 , 0) rechts unten. Von dort bewegt sich der Punkt
(cD (t), cN O (t)) zunächst nach links oben, in Richtung höherer NO-Konzentration und
abnehmender Donor-Konzentration. Wenn der Maximalwert der NO-Konzentration
überschritten ist, bewegt sich das System nach links unten, beide Konzentrationen
nehmen ab. Die Punkte auf der Kurve könnten wir mit der Hand zeichnen, in dem wir zu
jedem Simulations-Zeitpunkt t die Größe von cD (t) aus Abbildung 8.7 entnehmen und
als waagrechte Koordinate auftragen und als senkrechte Koordinate die Größe von
cN O (t) für dieselbe Zeit t. Die Kurve in Abbildung 8.12 heisst eine Trajektorie des
Donor-Modells. Sie besteht aus allen Punkten (cD (t), cN O (t)), die das Modell, ausgehend
von einem speziellen Anfangswert, im Lauf der Zeit durchläuft. Obwohl cD auf der
waagrechten und cN O auf der senkrechten Achse aufgetragen wird, ist die Trajektorie
nicht so zu interpretieren, dass cN O eine Funktion von cD wäre. Die Trajektorie ist eine
98
Abbildung 8.12. Simulationsergebnisse des Stickoxidmodells
−6
3.5
x 10
3
Konzentration von NO
2.5
2
1.5
1
0.5
0
0
0.2
0.4
0.6
0.8
Konzentration des Donors
1
1.2
−5
x 10
Teilmenge des Zustandsraums, das ist der Vektorraum aller möglichen Zustände; für das
Donorproblem ist der “Raum” zweidimensional, eine Ebene.
Definition 8.9. Die Werte der n Zustandsgrößen eines Modells fassen wir zu einem
Zustandsvektor mit n Komponenten zusammen, kurz Zustand genannt. Der Raum aller
Zustandsvektoren heisst der Zustandsraum (oder Phasenraum), das ist ein Vektorraum
mit der Dimension n. Die zeitliche Entwicklung eines kontinuierlichen dynamischen
Systems erzeugt eine Trajektorie
{(x1 (t), . . . , xn (t)) : t ≥ t0 } ⊂ Rn
im Zustandsraum, das ist eine Kurve, die aus den Punkten besteht, die der Zustand ausgehend von einem Anfangszustand zum Zeitpunkt t0 - durchläuft. (In der
englischsprachigen Literatur wird statt trajectory auch das Wort orbit verwendet.)
Beachten Sie, dass die zwei Kurven in Abbildung 8.7 mehr Information beinhalten als die
Kurve in Abbidung 8.12 allein. Man kann nämlich, wie geschildert, von der Abbidung 8.7
eindeutig zur Trajektorie gelangen aber nicht umgekehrt: die Zeit ist in Abbidung 8.12
nicht explizit enthalten. Daher sieht man einer Trajektorie ohne Zusatzinformation nicht
an, wo sie beginnt und endet. Auch über die Geschwindigkeit mit der sich der Zustand
ändert ist in Abbidung 8.12 nichts enthalten. Aus Abbidung 8.7 wissen wir jedoch, dass
der Zustand sich mit abnehmender Geschwindigkeit dem Nullpunkt nähert. Wir könnten
in Abbidung 8.12 Zusatzinformation einzeichnen indem wir den Anfangszustand
markieren, oder die Geschwindigkeit, mit der die Trajektorie durchlaufen wird, dadurch
andeuten, dass wir in regelmäßigen zeitlichen Abständen ein Symbol auf die Trajektorie
setzen. Dann würde der Abstand der Symbole in Richtung links unten immer kleiner
8. SIMULATIONSDIAGRAMME
99
Abbildung 8.13. Die drei Zeit-Plots der endemischen Epidemie
4
10
x 10
S(t)
9
8
7
6
0
500
1000
1500
2000
2500
3000
3500
4000
0
4
x 10
500
1000
1500
2000
2500
3000
3500
4000
0
500
1000
1500
2000
Zeit t
2500
3000
3500
4000
2000
I(t)
1500
1000
500
4
R(t)
3
2
1
0
werden. Oft deutet man wenigstens die Richtung, in der die Trajektorie durchlaufen
wird, mit Pfeilen an.
Ein Beispiel mit einem dreidimensionalen Zustandsraum ist das Epidemie-Modell: die
Trajektorien bestehen aus Punkten mit den Koordinaten (S(t), I(t), R(t)).
Abbildung 8.13 zeigt die wellenförmige zeitliche Entwicklung der Epidemie: immer
wieder verschwindet die Krankheit fast vollkommen aus der Population, aber nie ganz.
Durch das wiederholte Ausbrechen der Infektionskrankheit (schneller Anstieg von I(t))
geht die Zahl der Ansteckbaren S(t) zurück, gleichzeitig steigt die Zahl der Geheilten
Abbildung 8.14. Die Trajektorie der endemischen Epidemie
4
x 10
4
R(t)
3
2
1
0
2500
2000
1500
10
9.5
1000
9
8.5
8
500
7.5
7
I(t)
0
6.5
6
S(t)
4
x 10
100
R(t). Die Parameter und der Anfangszustand sind hier gleich gewählt wie in Beispiel 8.8.
Die dadurch erzeugte Trajektorie sehen wir in Abbildung 8.14. Tatsächlich ist die linke
Kurve in Abbildung 8.11 nichts anderes als die Projektion der Trajektorie auf die
S − I-Ebene. Wieder wissen wir nur auf Grund der Zeitplots, dass die immer enger
werdende Spirale in Abbildung 8.14 rechts unten beginnt und ein im Zentrum liegendes
endemisches Gleichgewicht im Gegenuhrzeigersinn umkreist.
Es gibt also zwei fundamental unterschiedliche Möglichkeiten der graphischen
Darstellung der zeitlichen Entwicklung eines dynamischen Systems: die Zeitplots der
einzelnen Komponenten des Zustands und die Trajektorien des Zustands im
Zustandsraum. Beide Darstellungen haben ihre Vor- und Nachteile; verwendet werden
alle beide je nach Zweck und Anwendung. Als eine reduzierte Variante können wir auch
nur einige der Komponenten der Trajektorie zeichnen. Dies ist bei mehr als drei
Zustandsgrößen unumgänglich, kann aber auch einfach die Interpretation erleichtern
(vgl. Abbildung 8.14 und Abbildung 8.11).
Die geometrischen Merkmale der Trajektorien bringen besonders anschaulich qualitative
Eigenschaften des Systems zum Ausdruck. Wenn das System sich etwa in einer
Gleichgewichtslage befindet, sich also der Zustand nicht mit der Zeit ändert, dann ist die
entsprechende Trajektorie ein einziger Punkt im Zustandsraum. Viele
Differentialgleichungs-Modelle haben die Eigenschaft, dass die Lösung zu gegebenem
Anfangswert eindeutig bestimmt ist. Dies bedeutet aber, dass sich zwei Trajektorien
eines solchen Modells nicht überkreuzen können, sonst wäre der Schnittpunkt ein
Anfangswert mit zwei möglichen Lösungen. Viele Differentialgleichungs-Modelle haben
sogar die Eigenschaft, dass je zwei verschiedene Trajektorien gar keinen gemeinsamen
Punkt haben können. Dann kann eine Trajektorie, die sich einer Gleichgewichtslage
nähert, dies nur mit nach Null gehender Geschwindigkeit tun. Sie erreicht die
Gleichgewichtslage nie ganz, sonst hätte sie mit dieser einen gemeinsamen Punkt. Zum
Beispiel in unserem Stickoxid-Modell gibt es für beliebig lange Zeitintervalle noch
Restkonzentrationen größer als Null. Im Epidemie-Modell wird das endemische oder auch
das gesunde Gleichgewicht nie erreicht (vgl. Abbildung 8.11), obwohl die Trajektorie im
Lauf der Zeit beliebig nahe herankommt.
Dieses - streng genommen - unrealistische Verhalten der Modelle kommt durch die
kontinuierliche und stetige Entwicklung der Modell-Trajektorien zustande. Für sehr
kleine Konzentrationen müssten wir die Reaktionen einzelner Moleküle modellieren, denn
irgendwann zerfällt das letzte Molekül Donor und oxidiert das letzte Molekül NO. Und
in der Nähe des endemischen Gleichgewichts schwanken die Anzahlen von Individuen in
den Klassen S, I, R ganzzahlig und unregelmäßig, bleiben aber in der Nähe des
Gleichgewichts. Diese unstetigen Änderungen in den realen Systemen kann ein
kontinuierliches Modell nicht nachahmen (ganzzahlige probabilistische Modelle werden
wir später kennenlernen). Wie wir gesehen haben gibt es aber sehr wohl wichtige
Aspekte der Entwicklungen in realen Systemen, die von Modellen mit
Differentialgleichungen gut und effizient wiedergegeben werden.
8. SIMULATIONSDIAGRAMME
101
Zwischenbilanz und Ausblick:
Hier endet der erste Teil des Zyklus Quantitative Systemwissenschaften. Dem
Studienplan der Umweltsystemwissenschaften an der Universität Graz entsprechend
wurde versucht, im ersten Teil schon die wichtigsten Elemente der mathematischen
Modellbildung und Simulation unterzubringen. Didaktisch wurde ein flacher Einstieg
gewählt, mit Themen, die schon in anderen Lehrveranstaltungen vorbereitet wurden. Im
Laufe der Vorlesung wird die Dichte des Dargebotenen größer, am Ende werden wichtige
Konzepte schon recht komprimiert dargestellt. Dies geschieht unter der Annahme, dass
die HörerInnen sich gewöhnen an die Sprechweise in Modellen, Funktionen, Lösungen
und Numerik. Es ist empfehlenswert die mathematisch ausformulierten Teile und die
konkreten Programmierbeispiele nicht einfach zu überspringen. Lässt man sich einmal
ernstlich darauf ein, wird man feststellen, dass es sich um nachvollziehbare Teile handelt,
deren Verständnis viel zur Sicherheit gegenüber der zu erlernenden Materie beiträgt.
Ein Charakteristikum dieser Vorlesung ist die Gleichzeitigkeit von Beispielen, allgemeinen
Konzepten und konkreter Vorgangsweise. Dies birgt die Gefahr eines verwirrenden
Durcheinanders, kann aber die mitunter spröde anmutenden mathematischen Techniken
durch anwendungsorientierte Interessen motivieren. Versuchen Sie, die allgemeinen
Konzepte nicht - vor lauter Beispielen und Details - aus den Augen zu verlieren.
Der zweite Teil des Zyklus ist nach den gleichen Prinzipien entworfen. Unter anderem
baut er auf den Modellen des ersten Teils auf, insbesondere wird qualitative Analyse an
Hand des Epidemie-Modells demonstriert. Nach einem kleinen Ausflug in die nichtlineare
Dynamik (Bifurkationen, Chaos) werden räumlich verteilte Systeme modelliert.
Durch Grenzübergang gelangt man von einem Vielzellen-Modell zu partiellen
Differentialgleichungen. Anhand der Schadstoffausbreitung in einer Flussströmung
werden die wichtigsten Typen solcher Gleichungen plausibel gemacht. Die dazu
notwendigen Vorkenntnisse entsprechen - selbstverständlich - dem Studienplan USW
Bakkalaureat. An Hand der Matlab software Femlab wird demonstriert, wie partielle
Differentialgleichungen zwar recht bequem numerisch simuliert werden können, aber
auch, welche Fehler dabei auftreten.
Dann werden zeit-diskrete Modelle behandelt. Ein klassisches Epidemie-Modell wird
erstellt und qualitativ untersucht (Schwellenwert). Der Bezug zum SIR Modell erfolgt
durch Grenzübergang. Als nächstes werden Zellularautomaten als räumlich verteilte
diskrete Modelle erläutert. An Hand des legendären Game of Life wird sichtbar, wie
komplex das Verhalten eines äusserst einfach definierten Modells sein kann. Hier beginnt
die Beschäftigung mit Modellen - ich nenne sie Simulationsmodelle - die nicht versuchen
ein real existierendes System mathematisch abzubilden. Vielmehr werden solche Modelle
nur numerisch implementiert und qualitative Eigenschaften des Modell-Verhaltens in
Simulationen zum Vorschein gebracht.
Schließlich wenden wir uns stochastischen diskreten Modellen zu, und zwar zunächst dem
Langzeitverhalten linearer Markov-Ketten. Dann betrachten wir ein ganzzahliges
probabilistisches Epidemie-Modell: an Hand einer hoch nichtlinearen Markov-Kette wird
102
erläutert und demonstriert, wie man mit Monte-Carlo-Simulationen zu statistischen
Prognosen kommt. Am Ende der Vorlesung mit immanenten Übungen steht ein
probabilistischer Zellularautomat: Forest Fire. Dies ist ein Prototyp von
Simulationsmodellen, die Systemeigenschaften wie Phasenübergang, Selbstorganisation
oder Musterbildung veranschaulichen und verstehen helfen.
Bücher zu den Vorlesungen Quantitative Systemwissenschaften 1 und 2:
• P. Doucet and P.B. Sloep, Mathematical Modeling in the Life Sciences, Ellis
Horwood 1992.
Eine gut lesbare Einführung in die Methoden der mathematischen
Modellbildung, der einige Ideen, die wesentlichen Konzepte zu erläutern,
entnommen wurden.
• MT. Hütt, Datenanalyse in der Biologie, Springer 2001.
Eine kompakte und verständlich geschriebene Einführung in die nichtlineare
Dynamik, fraktale Geometrie und Informationstheorie, der die
Bifurkations-Beispiele und der Waldbrand-Automat entnommen wurden.
• D. Kaplan and L. Glass, Understanding Nonlinear Dynamics, Springer 1995.
Eine populäre, leicht verständliche, für Nicht-Mathematiker geschriebene
Einführung in Chaos, Fraktale, Zellular-Automaten und nichtlineare Dynamik
mit vielen realen Beispielen.
• J.B. Snape, I.J. Dunn, J. Ingham and J.E. Prenosil, Dynamics of Environmental
Bioprocesses, VCH 1995.
Eine umfassende Sammlung von Grundlagen, Methoden und Beispielen der
Modellbildung zu Umweltprozessen.
• D.G. Luenberger, Introduction to Dynamic Systems, Wiley 1979.
Klassische Einführung in die Mathematik linearer und nichtlinearer dynamischer
Systeme, der die Darstellung der Markov-Ketten entnommen wurde.
• F.C. Hoppensteadt and C.S. Peskin, Modeling and Simulation in Medicine and
the Life Sciences, Springer 1992, 2002.
Naturwissenschaftlich anspruchsvoll geschriebene Darstellung der Anwendung
mathematischer Modelle in den Biowissenschaften, der die diskreten
Infektionsmodelle entnommen wurden.
• H. Bossel, Modellbildung und Simulation, Vieweg 1994.
Eine lehrmeisterlich geschriebene Anleitung zur Bildung mathematischer
Modelle und Simulation (inklusive Software), der das kleine Weltmodell
entnommen wurde.

Download Report