3) Kontinuierliche Modellbildung und Simulation
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 1 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.1.
Populationsdynamik: Modelle und ihre numerische Behandlung 3.1.1.
Modelle der Populationsdynamik
• die numerische Simulation stützt sich typischerweise auf kontinuierliche Modelle (reellwertige Größen, Apparat der Analysis) • zwei große Klassen:
Populationsdynamik: . . .
– Probleme mit Berücksichtigung des Ortes bzw. Raumes führen auf partielle Differentialgleichungen oder kurz PDE
Regelungstechnik: . . .
– Probleme ohne Berücksichtigung des Ortes bzw. Raumes führen auf gewöhnliche Differentialgleichungen oder kurz ODE
Wärmeleitung: Modell . . .
Verkehrsfluss: . . .
• Standardbeispiel für letztere: Populationsdynamik – Entwicklung bzw. Wachstum von Populationen, * entweder als isolierte Spezies (ohne externe Einflüsse) * oder in (friedlicher / feindlicher) Koexistenz verschiedener Spezies – Modellierung hat hier große Tradition (Biologie und Mathematik) – ein Beispiel mit zwei „Spezies“ schon in Kapitel 1: Wettrüsten! – klassischer, einfacher Vertreter eines 1-Spezies-Modells: Modell nach Maltus (1798)
Page 2 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Modell nach Maltus • betrachtet wird nur eine einzige Spezies P : – konstante Geburtenrate γ pro Zeiteinheit und Individuum – konstante Sterberate δ pro Zeiteinheit und Individuum – dadurch auch konstante Wachstumsrate λ = γ − δ (positiv/negativ) • Studium der Entwicklung von p(t), der Anzahl von Individuen von P : p(t + ∆t) = p(t) + λ · p(t) · ∆t
Populationsdynamik: . . .
(Wachstum ist proportional zu Populationsgröße und Zeit)
Regelungstechnik: . . . Verkehrsfluss: . . .
• führt auf die gewöhnliche Differentialgleichung
Wärmeleitung: Modell . . .
p(t) ˙ = λ p(t) mit der Lösung p(t) = p0 eλt , wobei p(0) = p0 • Anmerkung: – exponentielles Verhalten (Wachstum oder Abnahme extrem schnell) – zunächst diskrete Realität (Anzahl von Individuen ist ganzzahlig!), aber dennoch kontinuierliches Modell!
Page 3 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskrete Herleitung? • obige Herleitung war kontinuierlich (Analysis, Grenzwert ∆t → 0 ) • es geht aber auch rein diskret : – Modellhypothese einer konstanten Wachstumsrate entspricht Annahme der Verdopplung des Bestands alle τ Jahre – Start zur Zeit t = 0 mit der Anfangsstärke p0 – nach k · τ Jahren folgende Populationsstärke: ln 2·k· ττ
p(k · τ ) = 2 · p0 = p0 · e k
– mit λ = ln 2/τ also
= p0 · e
Populationsdynamik: . . . ln 2 τ ·k·τ
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
λ·(k·τ )
p(k · τ ) = p0 · e
– „wagt“ man jetzt noch den Schritt vom ganzzahligen k · τ zum reellwertigen t von vorhin, so ist unser bisheriges kontinuierliches Modell (die ODE) auch mittels diskreter Argumente hergeleitet • beruhigend – wir bleiben deshalb bei kontinuierlichen Ansätzen
Page 4 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Einschub: Bezüge zur Verzinsung • Startkapital K, Jahresnominalzins R% – nach n Jahren (jährliche Verzinsung): K · (1 + r)n , r = R/100 12·n r – nach n Jahren (monatliche Verzinsung): K · 1 + 12 365·n r – nach n Jahren (tägliche Verzinsung): K · 1 + 365 – im Grenzwert (unendlich fein): N ·n N ·n limN →∞ K · 1 + Nr = limN →∞ K · 1 + Nr·n = K · er·n ·n – exponentielles Wachstum ( r = λ, n = t , Lob den Banken!)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Verdopplung des Kapitals nach τ Jahren (vgl. Folie 4): 2K = K · (1 + r)τ ,
also
τ=
ln 2 0.7 70 ≈ = ln(1 + r) r R
2K = K · er·τ ,
also
τ=
ln 2 0.7 70 ≈ = r r R
bzw. Page 5 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Erste Verfeinerung: Sättigung • Ist exponentielles Wachstum realistisch? – Bevölkerung der Erde zwischen 1700 und 1961: ja! * Wachstumsrate von ungefähr 0.02 * das heißt: Verdopplung etwa alle 34.67 Jahre – im Allgemeinen jedoch: nein! * begrenzte Aufnahmekapazität der Erde, begrenzte Ressourcen * zunehmender Wettbewerb/Kampf um Nahrung, Wasser oder Luft wirkt bremsend auf das Bevölkerungswachstum
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• Verfeinerung gemäß Verhulst und anderen (19. Jahrhundert):
Wärmeleitung: Modell . . .
– Bevölkerungszahl strebt gegen Sättigungsgrenzwert – lineare Geburten- und Sterberaten (jetzt nur pro Zeiteinheit): γ(t) = γ0 − γ1 p(t),
δ(t) = δ0 + δ1 p(t),
γ0 > δ0 > 0, γ1 , δ1 > 0
(ein hoher Bestand dämpft Geburten- und erhöht Sterberate) – Für t → ∞ existiert ein Grenzwert 0 < p∞ < ∞ Aussterben noch Bevölkerungsexplosion!
– also weder Page 6 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Sättigung nach Verhulst (2) • zugrunde gelegte Differentialgleichung: p(t) ˙ = γ(t) − δ(t) = γ0 − δ0 − (γ1 + δ1 ) · p(t) = −m · (p(t) − p∞ ) mit m = γ1 + δ1 > 0,
p∞ =
γ0 − δ0 >0 γ1 + δ1
• Startwert zur Zeit Null:
Populationsdynamik: . . .
p(0) = p0
Regelungstechnik: . . . Verkehrsfluss: . . .
• als Lösung errechnet man:
Wärmeleitung: Modell . . . −m·t
p(t) = p∞ + (p0 − p∞ ) · e • Eigenschaften der Lösung: – Grenzwert gegen unendlich: lim p(t) = p∞
t→∞
– beschränktes Wachstum, in vielen Szenarien realistischer
Page 7 of 165
– allerdings: gebremste Zunahme / Abnahme von Anfang an
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation
– oft tritt Sättigungseffekt jedoch erst nach einiger Zeit ein
Hans-Joachim Bungartz
Zweite Verfeinerung: Logistisches Wachstum • Ist Verhulsts erster (Sättigungs-) Ansatz realistisch? – zweite Ableitung (Krümmung) von p(t) ändert ihr Vorzeichen nicht! – S-Form (logistisches Wachstum) kommt jedoch oft vor • neues Modell Verhulsts, um Wendepunkt zu erlauben: p(t) ˙ = a − b · p(t) · p(t) = a · p(t) − b · p2 (t), a >> b > 0,
p(0) = p0 Populationsdynamik: . . .
mit der Lösung a · p0 p(t) = b · p0 + (a − b · p0 ) · e−at
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Diskussion: – Grenzwert für t gegen unendlich existiert: lim p(t) =
t→∞
a b
– p0 < a/b : es gibt einen Wendepunkt, p(t) < a/b ∀t – p0 > a/b : kein Wendepunkt, p(t) monoton fallend, p(t) > a/b ∀t – Beispiel USA 1790-1950: a = 0.03134, b = 1.5587 · 10−10 – Modell mächtiger als zuvor, aber immer noch ohne äußere (und insb. agierende bzw. wechselwirkende) Einflüsse
Page 8 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Logistisches Wachstum (2) • Typischerweise ist a wesentlich größer als b : – quadratischer Term schlägt dadurch erst bei großem p(t) zu – Warum überhaupt quadratisch und nicht kubisch? Dies ist eben eine Modellthese! – Rechtfertigung des quadratischen Modellansatzes: * Störung eines Individuums ist proportional zu p(t) * somit Störung der Allgemeinheit proportional zum Quadrat von p(t) • logistisches Wachstum auch anderweitig zu erzielen; beim Tumorwachstum etwa folgende Ansätze:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– p(t) ˙ = λ(t) · p(t) * dabei λ(t) positiv, monoton fallend und stetig * Ausgangspunkt: empirisch begründete Annahme zeitlich abklingender Wachstumsraten – p(t) ˙ = f (p(t)) · p(t) * dabei f (p) positiv, monoton fallend und für t → ∞ gegen 0 strebend * dies ist Oberklasse zu beiden bisher genannten logistischen Modellen * Empirie: reicht aus für geschlossene Populationen
Page 9 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Weitere Möglichkeit: Oszillationen • Idee dieses qualitativ anderen Modellansatzes: – Grenzwert limt→∞ p(t) = p∞ soll wieder existieren – Überschwingen und allmähliches Einpendeln auf Grenzwert – Modell: Schwingungen (konkret: linearer gedämpfter Oszillator) • Differentialgleichung: p¨(t) + µ · p(t) ˙ + ω 2 · p(t) − p∞ = 0 µ > 0,
2
Populationsdynamik: . . . Regelungstechnik: . . .
2
4·ω −µ ≥0
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• neu: – zweite Ableitung von p(t) tritt auf: somit also Wechselspiel aus Populationsstärke, Wachstum und Wachstumsänderung • Lösung bei Anfangsbedingung: p(0) = p0 − µt 2
p(t) = (p0 − p∞ ) · e
· cos
µ2 ω2 − ·t 4
+ p∞ Page 10 of 165
• Interpretation (insbesondere der zweiten Ableitung)?
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Mehr als eine Spezies • nächster Schritt in Richtung mehr Realismus: betrachte zwei Spezies P und Q, z.B. p(t) ˙ = f (p(t), q(t)) · p(t), q(t) ˙ = g(p(t), q(t)) · q(t), • f und g: Wachstumsraten, sinnvoll für positive Werte von p und q • es gebe p¯, q¯ mit f (¯ p, q¯) = g(¯ p, q¯) = 0 (nicht immer der Fall): p(t) p¯ – = wird Gleichgewichtslösung genannt, falls p¯, q¯ > 0 q(t) q¯ f (¯ p, q¯) = g(¯ p, q¯) = 0 ⇒ p(t) ˙ = q(t) ˙ =0
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Gleichgewicht beschreibt stabilen bzw. stationären Zustand des Systems (kein Wachstum); die Lösung (p)t = p¯, q(t) = q¯ heißt stationäre Lösung • entscheidende Fragen: – Gibt es eine Gleichgewichtslösung? – Falls ja, ist sie attraktiv (d.h. wird sie irgendwann oder im Limes angenommen, und verharren p und q dann auch in ihr)? – vgl. Situation beim Beispiel des Wettrüstens in Kapitel 1
Page 11 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Einfaches Beispiel: Wettrüsten • schon in Kapitel 1: Wettrüst-Modell x(t) ˙ = f (x, y) = −mx(t) + ay(t) + c y(t) ˙ = g(x, y) = bx(t) − ny(t) + d, a, b, c, d, m, n > 0
x(t) ˙ −m = y(t) ˙ b
a x(t) c · + −n y(t) d
• beachte: einfacheres Modell (linear!) als auf Folie 11! • Interpretation:
Populationsdynamik: . . .
– Rüstungsausgaben x(t), y(t) von zwei Großmächten X und Y – Abrüstraten m und n , Aufrüstraten a und b
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– konstante Aufrüstbeiträge c und d (Abrüstbeiträge, falls negativ) • man kann hier zeigen: – wenn m · n = a · b , dann existiert (¯ x, y¯) mit f (¯ x, y¯) = g(¯ x, y¯) = 0 – somit hat man die stationäre Lösung x(t) x ¯ = y(t) y¯
∀t
(Gleichgewicht, wenn beide Komponenten positiv) – wenn m · n > a · b , dann liegt zudem Attraktivität vor (Deutung?)
Page 12 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Wettrüsten: Gleichgewicht
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
x(t) ˙ −1 1/2 x(t) 3/2 = · + y(t) ˙ 1/2 −1 y(t) 0
• Eigenwerte der Matrix: −1/2, −3/2 • es gilt m · n > a · b 1/2 x ¯ 2 x0 = , Gleichgewicht (attraktiv) bei = • Start bei 3/2 y¯ 1 y0
Page 13 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Wettrüsten: Explosion
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
x(t) ˙ −3/4 1 x(t) 1/2 = · + y(t) ˙ 1 −3/4 y(t) −5/4
• Eigenwerte der Matrix: 1/4, −7/4 • es gilt m · n < a · b x0 5/4 x ¯ 2 • Start bei = , Gleichgewicht (nicht attraktiv) bei = 9/4 y¯ 1 y0
Page 14 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Wettrüsten: noch ’was Nettes
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
x(t) ˙ −3/4 = y(t) ˙ −1
1 x(t) 0 · + −3/4 y(t) 5/2
• Eigenwerte der Matrix: −3/4 ± i 1/2 x ¯ 8/5 x0 = , Gleichgewicht (attraktiv) bei = • Start bei y0 11/4 y¯ 6/5
Page 15 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Wettrüsten (die letzte) • d.h.: Y rüstet auch bei großem X(t) ab und fährt mit dieser „pazifistischen“ Marschroute sicherheitspolitisch nicht schlechter als mit „konventioneller“ Strategie! • entscheidende Frage: – Ist das Gleichgewicht wirklich stabil? Bedeutet nicht z.B. x(t) x ¯ 100 = = y(t) y¯ 1 in der Logik des Kalten Krieges, dass X irgendwann doch angreift – und sich einen Dreck um die theoretische Stabilität schert?
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Modell also noch unvollständig: * große positive Differenz x(t) − y(t) führt wohl zu Aufrüsten bei Y , große negative Differenz zu Abrüsten (analog für X) * vgl. „Two-Power-Standard“ (Royal Navy stärker als zwei beliebige andere europäische Marinen zusammen, British Naval Defence Act 1889) • damit sind Wege zur weiteren Modellverfeinerung aufgezeichnet
Page 16 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Existenz eines Gleichgewichts • Folie 11 sowie Theorie der ODE: hinreichend für Attraktivität eines Gleichgewichts (¯ p, q¯) sind negative Realteile der Eigenwerte der Jacobimatrix zu f (p, q) · p F (p, q) = g(p, q) · q in (¯ p, q¯) , also der Eigenwerte von p, q¯) · p¯ + f (¯ p, q¯) fq (¯ p, q¯) · p¯ p, q¯) · p¯ fq (¯ p, q¯) · p¯ fp (¯ fp (¯ = p, q¯) = JF (¯ gp (¯ gp (¯ p, q¯) · q¯ gq (¯ p, q¯) · q¯ + g(¯ p, q¯) p, q¯) · q¯ gq (¯ p, q¯) · q¯
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• Vielzahl denkbarer Konstellationen (verschiedene Eigenschaften von f und g , wichtig für Gleichgewicht):
Wärmeleitung: Modell . . .
– Wettbewerb / Konkurrenz : beide Spezies von derselben Charakteristik, konkurrieren um denselben Lebensraum – Räuber-Beute : eine Spezies ist natürliche Beute der anderen – Symbiose : beide profitieren voneinander in einer Gemeinschaft
Page 17 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Konkurrenz-Szenario • Die Spezies P und Q sind keine Feinde im Sinne von Nahrung, aber sie konkurrieren um dieselben Ressourcen (etwa Löwen und Leoparden): fp (p, q), fq (p, q), gp (p, q), gq (p, q) < 0 für p, q > 0 • Als hinreichende Eigenwertbedingung ergibt sich fp (¯ p, q¯) · gq (¯ p, q¯) − fq (¯ p, q¯) · gp (¯ p, q¯) > 0 (Deutung: Eigenbehinderung größer als Fremdbehinderung) • einfachstes (konkretes) Modell: lineares f und g f (p, q) = a1 + a2 · p + a3 · q,
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
g(p, q) = a4 + a5 · p + a6 · q
mit (aufgrund der obigen Modellannahmen an f und g ) a1 , a4 > 0,
a2 , a3 , a5 , a6 < 0,
a2 · a6 > a3 · a5
• System f (¯ p, q¯) = g(¯ p, q¯) = 0 ist eindeutig lösbar (lineare Algebra), die zugehörige stationäre Lösung ist attraktiv (Eigenwerte!) • einzig p¯ > 0, q¯ > 0 ist noch zu klären
Page 18 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Konkurrenz-Szenario (2) • Ermittlung des Gleichgewichtszustands: – zu lösen sind zwei lineare Gleichungen in zwei Unbekannten: f (p, q) = a1 + a2 · p + a3 · q = 0 g(p, q) = a4 + a5 · p + a6 · q = 0 – eindeutige Lösung: a3 a4 − a1 a6 , p¯ = a2 a6 − a3 a5
a1 a5 − a4 a2 q¯ = a2 a6 − a3 a5
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
(Nenner größer Null!) – Lösung definiert Gleichgewicht, falls p¯ und q¯ größer Null; dies sichergestellt, falls a2 a1 a3 > > a5 a4 a6 – Attraktivität steckt schon in Modellannahme drin! • im linearen Fall ist auch der (Entartungs-) Fall a2 a6 = a3 a5 behandelbar: – Eigenbehinderung = Fremdbehinderung – die Folge kann beispielsweise das Aussterben einer Spezies sein
Page 19 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Konkurrenz: Gleichgewicht
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
√ √ 3 5 3 5 f (p, q) = + − ·p− ·q 2 24 8 24 √ √ 7 3 3 3 3 7 g(p, q) = + − ·p− ·q 8 2 8 8
• Eigenwerte der Matrix: −1/2, −1 • es gilt a2 a6 − a3 a5 = 0.5 > 0 p0 1/4 p¯ 4 • Start bei = , Gleichgewicht (attraktiv) bei = 3 q¯ 1 q0
Page 20 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Konkurrenz: Aussterben von P
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
71 23 25 f (p, q) = − ·p− ·q 8 12 12 73 25 23 g(p, q) = − ·p− ·q 8 12 12 • Eigenwerte der Matrix: −4, 1/6 • es gilt a2 a6 − a3 a5 = − 23 < 0 1/4 p0 = • Anfangswerte 1/2 q0
Page 21 of 165
• Gleichgewicht (nicht attraktiv, sondern Sattelpunkt) bei
p¯ 3 = q¯ 3/2
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Was macht Q so ganz allein? • P stirbt (und scheidet somit) aus! • Q gehorcht somit nun der ODE 73 23 73 23 2 q˙ = g(0, q) · q = − ·q ·q = ·q− ·q 8 12 8 12 • dies ist ODE für logistisches Wachstum (vgl. Folie 8)! • mit a =
73 8 ,b
=
23 12
Populationsdynamik: . . .
ergibt sich daher folgende Lösung: q(t) =
Regelungstechnik: . . .
a · q0 b · q0 + (a − b · q0 ) · e−at
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
(dabei gibt q0 den Bestand von Q zum Zeitpunkt des Aussterbens von P an – jetzt wieder auf t = 0 gesetzt) • Ablesen zeigt q0 ≈ 4.8 , und damit erhält man q(t) ≈
43.8 , 9.2 − 0.075 · e−9.125·t
lim q(t) =
t→∞
a ≈ 4.761 b Page 22 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel Konkurrenz: Aussterben (2)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Modell:
√ √ 3 3 7 3 3 f (p, q) = 7 − − ·p− ·q 32 4 16 √ √ 13 13 g(p, q) = − + 12 3 − 3 3 · p − ·q 8 4 • es gilt a2 a6 − a3 a5 = 4 > 0 1/4 p0 = • Start bei 1/2 q0 p¯ 4 • stabiler Zustand (attraktiv), aber kein Gleichgewicht = q¯ −1/2
Page 23 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel 1 zur Entartung • jetzt gelte a2 a6 = a3 a5 (Attraktivitätsbedingung verletzt!) • konkret:
f (p, q) = 4 − p − 2q g(p, q) = 6 − 2p − 4q
• es gilt (−1) · (−4) − (−2) · (−2) = 0 • keine Lösung für f = g = 0 ; zwei parallele Geraden; Gleichgewicht: einer stirbt aus
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 24 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel 2 zur Entartung • wieder sei a2 a6 = a3 a5 (Attraktivitätsbedingung verletzt!) • konkret:
f (p, q) = 8 − p − 2q g(p, q) = 16 − 2p − 4q
• es gilt (−1) · (−4) − (−2) · (−2) = 0 • unendliche viele Lösungen für f = g = 0 ; alle sind Gleichgewichtszustände! Populationsdynamik: . . .
• Zielpunkt hängt von Startpunkt ab
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 25 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Räuber-Beute-Szenario • Q ist „Nahrung“ für P , was zu einem unterschiedlichen Wachstumsverhalten führt (Unsymmetrie): fp (p, q), gp (p, q), gq (p, q) < 0,
fq (p, q) > 0
für p, q > 0
• in unserem linearen Modell heißt das a2 , a5 , a6 < 0,
a3 > 0
(Räuber P profitiert natürlich jetzt von starker Population der Beute Q) • Man kann für dieses Szenario zeigen:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Jede Gleichgewichtslösung (sofern existent) erfüllt die Eigenwertbedingung, ist also attraktiv! • Konstruktion von Gleichgewichtslösungen: – löse das 2 × 2-Gleichungssystem (wie zuvor; Nenner aufgrund obiger Bedingung nie Null) – prüfe, ob die berechneten Werte p¯, q¯ positiv sind Page 26 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Räuber-Beute-Szenario (2) • klassischer Vertreter: Modell nach Volterra und Lotka: – keinerlei Eigeneinfluss (P auf P bzw. Q auf Q): a2 = a6 = 0 – unsymmetrischer Fremdeinfluss wie zuvor: a3 > 0,
a5 < 0
Populationsdynamik: . . .
– konstante Terme in f und g (auch hier Unsymmetrie): a1 < 0,
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
a4 > 0
– Gleichgewichtslösung kann berechnet werden: 0 < p¯ = −
a4 , a5
0 < q¯ = −
a1 a3
– Beide Eigenwerte der Jacobimatrix sind rein imaginär (die Realteile der Eigenwerte sind somit Null und nicht negativ)!! – Statt der Attraktivität kann Stabilität gezeigt werden: Entwicklung bleibt gefangen in Umgebung der Gleichgewichtslösung!
Page 27 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel zu Volterra-Lotka • Volterra-Lotka-Szenario: p˙ = f (p, q) · p, q˙ = g(p, q) · q,
1 1 ·q f (p, q) = − + 2 200 1 1 g(p, q) = − ·p 5 50
• Gleichgewichtslösung: p¯ = 10,
q¯ = 100
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 28 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.1.2.
Numerische Behandlung von ODE
• Modelle der Populationsdynamik: eine/mehrere ODE, zusammen mit zusätzlichen Bedingungen zur eindeutigen Lösbarkeit – dort als Anfangswertproblem: Lösung vorgegeben am Start des betrachteten Zeitintervalls (Ausgangsstärke p0 ) – manchmal auch als Randwertproblem: Werte der Lösung an beiden Intervallendpunkten vorgegeben (z.B. optimale Flugbahn eines Space Shuttle) • Prototypen eines Anfangswertproblems (AWP): y(t) ˙ = f (t, y(t)),
y(a) = ya ,
Regelungstechnik: . . . Verkehrsfluss: . . .
t≥a
oder – im Falle eines Systemes – y˙i (t) = fi (t, y1 (t), . . . , yn (t)), yi (a) = yi,a ,
Populationsdynamik: . . .
Wärmeleitung: Modell . . .
t ≥ a,
i = 1, . . . , n
• Falls f nur vom ersten Argument t abhängt: einfaches Integrationsproblem (Trivialsituation)! • i.A.: analytisch nicht lösbar, näherungsweise (also numerische) Berechnung erforderlich! Page 29 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Kurzer Exkurs zur Diskretisierung • Start jedes numerischen Ansatzes: Diskretisierung • Übergang vom Kontinuum zu Diskretem/Endlichem: – statt (überabzählbar vieler) reeller Zahlen nun alternativ * ganzzahlige Arithmetik: feste Auflösung, fester Bereich, Integer * Fixpunkt-Arithmetik: feste Auflösung, fester Bereich, Dezimalpunkt * Gleitpunkt-Arithmetik: variable Auflösung, variabler Bereich Gleitpunktzahlen: IFB,t = {M · B E : M = 0 ∨ B t−1 ≤| M |< B t ; M, E ∈ ZZ} (Basis B, Mantisse M, Exponent E, Stellenzahl t)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Maschinenzahlen: IF = IFB,t,a,b = {f ∈ IFB,t : a ≤ E ≤ b} Auflösung (maximaler Relativabstand): ρ = B 1−t kleinste positive / größte Maschinenzahl: B 1−t · B α , (B t − 1) · B β – statt unendlicher Reihen (Sinus) endliche Polynome – diskrete Punktmenge (Gitter ) statt eines Intervalls etc. – statt Ableitungen nun Differenzenquotienten: y(t) ˙ ≈
y(t + h) − y(t) h
Page 30 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Rundung, Rundungsfehler verschiedene Fehlerquellen in numerischen Algorithmen • Rundungsfehler: durch Verwendung von Gleitpunktzahlen – Nachbarn zu x ∈ IR : fl (x) = max{f ∈ IF : f ≤ x}, fr (x) = min{f ∈ IF : f ≥ x} – Rundung: rd : IR → IF, surjektiv – idempotent – monoton * Aufrunden: rd(x) = fr (x) * Abrunden: rd(x) = fl (x) * korrektes Runden: zum Näheren der beiden (Sonderbehandlung des Mittelpunkts!) * Abhacken: zum näher an Null Gelegenen * Schranke für relativen Rundungsfehler: ρ bzw. ρ/2 – ideale Arithmetik:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
˙ = rd(a ∗ b) ∀a, b ∈ IF ∀∗ ∈ {+, −, ·, /} a∗b (technisch möglich (IEEE!), aber nicht in allen Rechnern realisiert) – Abschwächung: ˙ = (a ∗ b) · (1 + ε(a, b)), a∗b
| ε(a, b) |< ε¯ = O(ρ)
– Rundungsfehleranalyse: akkumulierten Einfluss bei längeren Berechnungen abschätzen (Beispiel: Horner-Schema zur Polynomauswertung)
Page 31 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Weitere Fehlerquellen • Diskretisierungsfehler: durch Verwendung diskreter Punkte anstelle des Kontinuums • Abbruchfehler: durch Abbruch einer Iteration nach N Iterationsschritten, z.B. – Abbruch einer Reihenberechnung nach N Gliedern – Abbruch einer Nullstellensuche nach Newton nach N Schritten – Abbruch einer Nullstellensuche nach Newton, wenn keine signifikante Änderung mehr eintritt • Datenfehler: Eingabedaten sind oft (ungenaue) Messwerte o. Ä. • Alle diese Fehlerquellen sind im Auge zu behalten!
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 32 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Das Euler-Verfahren • Standard-Diskretisierungsansatz für AWP: – Finite-Differenzen-Approximation: nähere die Ableitung durch geeignete Differenzenquotienten an, z.B. im einfachsten Fall y(a + δ t) = ˙ y(a) + δ t · f (t, y(a)) d.h. dann als Rechenvorschrift: yk+1 = yk + δ t · f (tk , yk ), tk = a + kδ t,
k = 0, 1, . . .
– Aus einer kontinuierlichen Gleichung (gültig in allen Punkten) wird somit eine Folge diskret(isiert)er Gleichungen (gültig in jeweils einem Punkt); Abstand zweier Punkte: Schrittweite
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Bereich kann auch endlich sein: t ∈ [a, b] – Vorgehensweise wird Euler-Verfahren genannt – Gretchenfragen: * Was passiert für δt → 0 ? Konvergenz gegen Lösung der ODE? * Wie schnell ist Konvergenz (Fehler bei Halbierung der Schrittweite)? * Gibt es scharfe Maximalgröße für die Schrittweite?
Page 33 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Euler-Verfahren, Diskretisierungsfehler • alternative Herleitung des Euler-Verfahrens: abgehackte Taylor-Entwicklung der exakten Lösung y(t) : y(tk+1 ) = y(tk ) + (tk+1 − tk )y(t ˙ k) + R = ˙ y(tk ) + (tk+1 − tk )f (tk , yk ) • lokaler Diskretisierungsfehler: maximaler lokaler Einfluss durch Verwendung von Differenzen statt Ableitungen; hier (Achtung: y(t) bezeichnet die exakte Lösung des Problems, keine Näherung): l(δ t) = max
(y(t + δ t) − y(t)) δt
[a,b]
− f (t, y(t))
• globaler Diskretisierungsfehler: maximaler Fehler aller berechneten diskreten Näherungswerte:
e(δ t) = max | yk − y(tk ) |
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
[a,b]
• Konsistenz:
l(δ t) → 0
• Konvergenz (stärker): e(δ t) → 0
für δ t → 0
Page 34 of 165
für δ t → 0
Modellbildung und Simulation
• Konsistenz- / Konvergenzordnung k: l(δ t), e(δ t) = O(δ t ) k
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Lokaler und globaler Diskretisierungsfehler im Vergleich
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 35 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Konsistenz- und Konvergenzordnung • Euler (bei Beschränktheit von y¨ und fy ) : – konsistent von erster Ordnung : l(δ t) = O(δ t) – konvergent von erster Ordnung : e(δ t) = O(δ t) – Es gibt aber durchaus konsistente Verfahren, die nicht konvergieren! • schnellere Konvergenz bei höherer Ordnung (k größer) – möglicher Ausgangspunkt Taylor-Entwicklung: führt zu komplizierten Formeln (höhere Ableitungen von f erforderlich)
Populationsdynamik: . . .
– alternativ: benutze zusätzliche Auswertungen von f (neben den Gitterpunkten): Verfahren vom Runge-Kutta-Typ
Verkehrsfluss: . . .
Regelungstechnik: . . .
Wärmeleitung: Modell . . .
• einfachster Vertreter: Methode von Heun δ t yk+1 = yk + f (tk , yk )) + f (tk+1 , yk + δ t · f (tk , yk )) 2 • sowohl konsistent als auch konvergent (2. Ordnung) l(δ t) = O(δ t2 ),
e(δ t) = O(δ t2 ) Page 36 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Verfahren nach Runge und Kutta • berühmtester Vertreter: Runge-Kutta-Verfahren δt (T1 + 2T2 + 2T3 + T4 ) , 6 = f (tk , yk ), δt δt = f tk + , y k + T 1 , 2 2 δt δt = f tk + , y k + T 2 , 2 2 = f (tk+1 , yk + δ t T3 )
yk+1 = yk + T1 T2 T3 T4
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• konsistent und konvergent von vierter Ordnung • Die Verfahren nach Euler bzw. Heun bzw. Runge-Kutta entsprechen übrigens alle einer Quadratur-Regel (Rechtecksregel bzw. Trapezregel bzw. Simpson-Regel ): – Integriere die ODE formal und interpretiere die einzelnen Verfahren als Näherungen für das dann zu bestimmende Integral über f : y˙ = f (t, y(t))
⇒
tk+1
y(tk+1 ) = y(tk ) +
f (t, y(t)) dt
Page 37 of 165
tk Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Alternative: Mehrschrittverfahren • Runge-Kutta-Verfahren sind teuer, da sie viele Auswertungen von f erfordern (f oft nicht in geschlossener Form verfügbar) • Alternative zur Erzielung höherer Ordnung: nutze Historie, also bereits erfolgte Auswertungen (Adams-Bashforth- oder Mehrschrittverfahren) – prominenter Vertreter: Methode 2. Ordnung δt δ t 3f (tk , yk ) − f (tk−1 , yk−1 ) = yk + (3fk − fk−1 ) 2 2
yk+1 = yk +
– allg.: betrachte Polynom P (t) vom Grad p − 1 , welches f in den p diskreten Stützstellen (ti , fi = f (ti , yi )), i = k − p + 1, . . . , k, interpoliert
tk+1
yk+1 = yk +
y(t) ˙ dt = yk + tk
tk+1
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
tk+1
f (t, y(t)) dt = ˙ yk + tk
Populationsdynamik: . . .
P (t) dt tk
– p = 1: Euler – p = 2: obiges Verfahren – allgemeiner Ansatz: Ordnung p – keine zusätzlichen Auswertungen von f (eine pro diskretem Zeitpunkt)
Page 38 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Mehrschrittverfahren (2) • am Anfang: keine oder zu wenige Vorgänger verfügbar • Abhilfe: modifizierte Vorgehensweise – Setze ein Einschrittverfahren oder ein Mehrschrittverfahren mit hinreichend kleiner Schrittzahl ein, um einen Wert für yk+1 und somit für fk+1 zu erhalten. – Ordnung bleibt erhalten! • Nachteil aller bisher vorgestellten Verfahren: – Die Schrittweite δt muss oft sehr klein sein, wenn Instabilitäten (Oszillationen in der berechneten Lösung, die dadurch von grundlegend anderem Charakter als die „echte“ Lösung ist) vermieden werden sollen!
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Eine kleine Schrittweite bedeutet natürlich einen entsprechend hohen Berechnungsaufwand! • Auch hierfür gibt’s eine Abhilfe: implizite Verfahren!
Page 39 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Implizite Verfahren • Alle bisher vorgestellten Verfahren sind explizit: Die jeweilige Vorschrift ermöglicht die direkte Berechnung eines weiteren Zeitschritts. • jetzt: benutze das zu bestimmende yk+1 auch rechts • führt zu Adams-Moulton Mehrschrittverfahren: – benutze Interpolation und Vorgänger wie bei Adams-Bashforth, aber mit dem neuen Zeitpunkt als Stützstelle Populationsdynamik: . . .
– impliziter Euler (1. Ordnung):
Regelungstechnik: . . .
yk+1 = yk + δ t · fk+1
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Variante 2. Ordnung: yk+1 = yk + δ t ·
(fk + fk+1 ) 2
– Variante 4. Ordnung: yk+1 = yk +
δt (fk−2 − 5fk−1 + 19fk + 9fk+1 ) 24
• Problem: Wie ermittelt man yk+1 im impliziten Fall? – Holzhammer: i.A. nichtlineare Gleichung auflösen (Newton o. Ä.) – einfacher und verbreitet: Prädiktor-Korrektor Ansatz
Page 40 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Prädiktor-Korrektor-Ansatz • zweistufiges Verfahren: – Zunächst wird mit einem geeigneten expliziten Verfahren ein vorläufiger Wert für yk+1 ermittelt (Prädiktor). – Dieser wird dann in der impliziten Vorschrift auf der rechten Seite verwendet, um den endgültigen Wert von yk+1 zu bestimmen (Korrektor). • wesentliche Eigenschaft (wenn man alles richtig macht): – vollständig explizites Verfahren (Prädiktor & Korrektor explizit!) – dennoch Charakteristik eines impliziten Verfahrens
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Bewertung impliziter Verfahren: – einzelner Zeitschritt ist teurer (Auflösung einer Gleichung anstelle des einfachen „Fortschaltens“ der Zeit) – Zahl der Zeitschritte i.A. kleiner (größere Schrittweite möglich), dadurch wieder weniger Rechenaufwand
Page 41 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Höhere Ableitungen, Systeme • Modelle der Populationsdynamik: nur erste Ableitung tritt auf • in der Mechanik etwa: zweite Ableitung (Beschleunigung, Kraft) • Vorgehensweise bei höheren Ableitungen: – Reduktion auf System von ODE mit nur erster Ableitung – dazu: Hilfsvariable y1 = y, y2 = y, ˙ y3 = y¨, . . . , yn = y (n−1) ;
Regelungstechnik: . . .
aus einer Gleichung n-ter Ordnung y wird somit
(n)
= f (t; y, y, ˙ y¨, . . . , y
Populationsdynamik: . . .
Verkehrsfluss: . . . (n−1)
)
Wärmeleitung: Modell . . .
y˙ 1 = y2 y˙ 2 = y3 ... y˙ n−1 = yn y˙ n = f (t; y1 , y2 , . . . , yn )
• Systeme: analoge Verfahren wie für eine ODE
Page 42 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schlecht konditionierte Probleme • Kondition eines Problems: Einfluss von Schwankungen in der Eingabe auf das Resultat – Eigenschaft des Problems, nicht eines numerischen Verfahrens! • schlechte Kondition: – Kleinste Trübungen in den Eingabedaten können zu völlig verschiedenen Ergebnissen führen! – Numerische Behandlung wird dann sehr schwer (weil exakte Eingabedaten fast nie verfügbar sind). • ein Beispiel:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– betrachte folgende ODE: y¨(t) − N · y(t) ˙ − (N + 1) · y(t) = 0 – Anfangsbedingungen (zwei, da ODE zweiter Ordnung): y(0) = 1, – exakte Lösung:
y(0) ˙ = −1 Page 43 of 165
y(t) = e−t
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schlecht konditionierte Probleme (2) • jetzt: beliebig kleine Änderung in einer der beiden Anfangsbedingungen: – neuer Wert von y in t = 0 :
yε (0) = 1 + ε
– resultierende neue Lösung: N + 1 −t ε yε (t) = 1 + ε ·e + · e(N +1)t N +2 N +2 – völlig verschiedenes Verhalten für t → ∞ : lim y(t) = 0,
t→∞
lim yε (t) = ∞
t→∞
• Risiko bei numerischen Verfahren:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Eingabedaten liegen nicht exakt vor (z.B. aus Messungen) – Zwischenresultate (Eingabe für nächsten Schritt) nicht exakt wegen Rundungsfehlern • schlechte Kondition: prima Ausrede („da geht eh nix“) oder GAU für die Numerik? Page 44 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Stabilität • Jetzt gehen wir von gut konditionierten Problemen aus (also keine Ausreden mehr!). • Ziel: Lücke zwischen Konsistenz und Konvergenz schließen • betrachte neues AWP: • exakte Lösung:
y(t) ˙ = −2y(t) + 1, e−2t +1 y(t) = 2
• gut konditioniert:
yε (0) = 1 + ε
⇒
y(0) = 1
| yε (t) − y(t) |= εe−2t
• verwende Mittelpunktsregel: yk+1 = yk−1 + 2δ t · fk also
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
yk+1 = yk−1 + 2δ t(−2yk + 1) = yk−1 − 4δt · yk + 2δt, y0 = 1 • 2-Schritt-Verfahren, starte mit y0 und dem exakten y(δt) – Zeitschrittweite
δ t = 1.0 ⇒ y9 = −4945.5,
y10 = 20953.9
– Zeitschrittweite
δ t = 0.1 ⇒ y79 = −1725.3,
y80 = 2105.7
– Zeitschrittweite
δ t = 0.01 ⇒ y999 = −154.6,
y1000 = 158.7 Page 45 of 165
• Fazit: irgendwann geht’s immer schief!
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel zu Stabilität und Mittelpunktsregel • jetzt das Ganze im Schaubild: AWP war y(t) ˙ = −2y(t) + 1, y(0) = 1,
Lösung y(t) =
e−2t+1 2
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 46 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Stabilität (2) • Phänomen Instabilität: – starke Oszillationen – Gestalt von berechneter und exakter Lösung völlig verschieden – berechnete Lösung nicht akzeptabel (nicht interpretierbar als exaktes Ergebnis zu leicht modifizierten Eingabedaten) • Mittelpunktsregel ist
Populationsdynamik: . . .
– konsistent (Differenzenquotient approximiert erste Ableitung) – offensichtlich aber nicht konvergent (sogar unabhängig von δt , wenn das Intervall nach rechts unbeschränkt ist)
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• allgemein gilt (unter Erhalt der Ordnung): Konsistenz + Stabilit¨ at = Konvergenz • Es gibt Stabilitätsbedingungen an ein Verfahren (stellen i.d.R. starke Einschränkungen an die Zeitschrittweite: sehr klein!) • stabile Verfahren: – alle expliziten Einschrittverfahren sind’s, Mittelpunktsregel nicht! – Adams-Bashforth und Adams-Moulton s-Schritt stabil für s > 1
Page 47 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Steifheit • ein weiteres Phänomen: betrachte das AWP y(t) ˙ = −1000y(t) + 1000, • exakte Lösung:
y(0) = y0 = 2
y(t) = e−1000t + 1
• Problem ist gut konditioniert Populationsdynamik: . . .
• explizites Euler-Verfahren liefert
Regelungstechnik: . . .
yk+1 = yk + δ t(−1000yk + 1000) = (1 − 1000δ t)yk + 1000δ t = (1 − 1000δ t)k+1 + 1
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• offensichtlich: Oszillationen und Divergenz für δ t > 0.002 • dies, obwohl – expliziter Euler konsistent – expliziter Euler stabil (wie alle expliziten s-Schritt-Verfahren) – somit auch konvergent!
Page 48 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Steifheit (2) • Phänomen wird Steifheit genannt – eine unwichtige Komponente der Lösung (im Beispiel der negative Exponentialterm e−1000t ) erzwingt Mini-Schrittweite bei der Diskretisierung (und zwar auf gesamtem Bereich!) – dadurch inakzeptabel hoher Berechnungsaufwand – dies, obwohl Lösung fast überall trivial darstellbar Populationsdynamik: . . .
• Erklärung? – Konsistenz und Stabilität sind beides asymptotische Begriffe, treffen also Aussagen für den Grenzfall hinreichend kleiner Schrittweite δt
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Konsistenz/Stabilität/Konvergenz also nicht verletzt, das Betreten der asymptotischen Phase erfolgt allerdings erst für sehr kleine Schrittweite und ist dadurch extrem aufwändig! – gilt für alle expliziten Verfahren – diese taugen nicht für steife Differentialgleichungen • Abhilfe: implizite Verfahren! Page 49 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Steifheit (3) • jetzt impliziter Euler: yk+1 = yk + δ t · f (tk+1 , yk+1 ) = yk + δ t · (−1000 · yk+1 + 1000) yk + 1000 · δ t = 1 + 1000 · δ t 1 +1 = (1 + 1000 · δ t)k+1 • offensichtlich: keine Oszillationen, stets Konvergenz • Erklärung:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Explizite Verfahren approximieren Lösung mit Polynomen- und das funktioniert bei negativ exponentiellem Abklingen eben nicht! – Implizite Verfahren setzen dagegen rationale Funktionen ein! • Ergo: bei steifen ODE stets implizite Verfahren verwenden!
Page 50 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zusammenfassung • Wir haben drei Phänomene betrachtet, die uns bei der numerischen Lösung von ODE das Leben schwer machen: – schlechte Kondition: * bedrohliche Eigenschaft des zugrunde liegenden Problems (hat nichts mit numerischen Näherungsverfahren zu tun) * im Extremfall wenig Spielraum für die Numerik
– Instabilität:
Populationsdynamik: . . .
* bedrohliche Eigenschaft eines numerischen Verfahrens, die i.d.R. zu sehr kleinen Zeitschritten und im Extremfall zur Untauglichkeit des Verfahrens führt * implizite Verfahren sind hier expliziten oft überlegen – Steifheit:
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* bedrohliche Eigenschaft eines Problems, der aber mit impliziten numerischen Verfahren beizukommen ist * implizite Verfahren hier ein Muss
Page 51 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ausblick auf Randwertprobleme • Beispiel: y¨ = f (t, y, y), ˙ t a ≤ t ≤ tb ,
y(ta ) = ya , y(tb ) = yb
• Spezialfall: y¨ = a(t)y(t) ˙ + b(t)y(t) + c(t) gleiche Randbedingungen • a = 0 und, b > 0: RWP hat eindeutige Lösung
Populationsdynamik: . . .
• diskretes Gitter:
Regelungstechnik: . . .
δ t = (tb − ta )/n, t0 = ta , tn = tb , ti = ta + iδ t
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Finite-Differenzen-Approximation für zweite Ableitung: y¨(t) = ˙
y(t + δ t) − 2y(t) + y(t − δ t) δ t2
• diskretes Analogon zur ODE in jedem Gitterpunkt: δ t−2 · (yi+1 − 2yi + yi−1 ) − bi yi = ci ,
i = 1, . . . , n − 1
• Tridiagonalsystem linearer Differenzengleichungen
Page 52 of 165 Modellbildung und Simulation
• Matrixeigenschaften entscheidend für Konvergenz!
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Tridiagonalmatrix im Beispiel • Gestalt des tridiagonalen Gleichungssystems: ⎞ ⎛ ⎛ 2 + b 1 h2 −1 ⎟ ⎜ ⎜ −1 2 + b2 h2 −1 ⎟ ⎜ ⎜ ⎟·⎜ ⎜ .. ⎠ ⎝ ⎝ . −1 −1 −1
2 + bn−1 h2
y1 y2 .. .
⎞
⎛
−h2 c1 + ya −h2 c2 .. .
⎞
⎟ ⎟ ⎜ ⎟ ⎟ ⎜ = ⎟ ⎟ ⎜ ⎠ ⎠ ⎝ 2 yn−1 −h cn−1 + yb
– Schreibweise jetzt: h für die Schrittweite – n − 1 diskrete Gleichungen an den n − 1 Innenpunkten
Populationsdynamik: . . . Regelungstechnik: . . .
– vorgegebene Randwerte in ta , tb wandern auf die rechte Seite
Verkehrsfluss: . . .
– Bedingung b(t) > 0 garantiert strenge Diagonaldominanz der Matrix: | aii |> | aij | ∀ i
Wärmeleitung: Modell . . .
j=i
• Eigenschaften der Matrix: – streng diagonaldominant (folglich nichtsingulär) – symmetrisch und positiv definit – Eigenwerte von A entscheidend für Konvergenz (hier: 2. Ordnung)
Page 53 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Randwertprobleme: Erste Ableitungen • jetzt a(t) nicht notwendig Null, auch erste Ableitungen kommen vor • Diskretisierung der ersten Ableitung: – zentrale Differenz (zweiter Ordnung): y(t) ˙ = ˙
y(t + h) − y(t − h) 2h
– Differenzengleichung ergibt sich in inneren Punkten zu
ai h ai h −1− yi−1 + (2 + bi h2 )yi + − 1 + yi+1 = −h2 ci 2 2
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Nachteil: zur eindeutigen Lösbarkeit sind kleine Zeitschritte nötig! (Bedingung | ai h |≤ 2 etwa garantiert Diagonaldominanz) – deshalb oft Upwind-Diskretisierung: 1 yi+1 − yi y(t) ˙ = ˙ · yi − yi−1 h
für ai < 0 für ai ≥ 0
– jetzt unbedingt diagonaldominant und nichtsingulär – Nachteil: Approximationsgüte nur noch von erster Ordnung!
Page 54 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Randwertprobleme: Randbedingungen • Randbedingungen bisher: – Festschreibung der Werte von y(t) an den Randpunkten des Intervalls: Dirichlet-Randbedingungen • weitere mögliche Randbedingungen: – Neumann-Randbedingungen: * Festschreibung der (ersten) Ableitung im Randpunkt: y(t ˙ a ) = y˙ a * diskretisierte Gleichung am Rand (virtueller Punkt t−1 ) : y˙ a =
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
y1 − y−1 y−1 − 2y0 + y1 2y1 − 2y0 − 2y˙ a h ˙ = ⇒ y¨(ta ) = 2 2h h h2
* Daraus ergibt sich diskrete Gleichung in t0 . * Achtung: y0 , yn sind jetzt auch Unbekannte! – periodische Randbedingungen: * Randwerte stimmen an beiden Enden überein: y(ta ) = y(tb ) = y0
Page 55 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schießverfahren • Ein ganz anderer Ansatz sind Schießverfahren: – Rückführung eines RWP auf eine Folge von AWP! • Ausgangspunkt wieder: y¨ = f (t, y, y), ˙ ta ≤ t ≤ tb , y(ta ) = ya , y(tb ) = yb • gelöst wird nun
Populationsdynamik: . . .
y¨ = f (t, y, y), ˙ ta ≤ t ≤ tb , y(ta ) = ya , y(t ˙ a) = s (mit unbekanntem Parameter s : „Schusswinkel“)
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Aufgabe lautet somit: finde s¯ und AWP-Lösung y(t; s¯) so, dass gilt y(tb ; s¯) = yb • Lösung: – erinnert an inverses Problem – Iteration (Newton-ähnlich, z.B. mit einem AWP pro Iterationsschritt)
Page 56 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zusammenfassung • ein erstes kontinuierliches Beispiel studiert: – Populationsdynamik! • Modellierung: – Herleitung verschiedener und unterschiedlich genauer Modelle für ganz unterschiedliche Szenarien – Rüstzeug: Anfangswertprobleme gewöhnlicher Differentialgleichungen – Diskussion der erhaltenen Modelle (Existenz von Lösungen, stationäre Zustände etc.)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Simulation: – direkte (analytische) Lösung scheitert i.A. – deshalb: numerische Lösung der Modellgleichungen – somit eingehender diskutiert: numerische Lösung von ODE * Anfangswertprobleme * Randwertprobleme
Page 57 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.2.
Regelungstechnik: Deterministische und Fuzzy Logic Ansätze 3.2.1.
Basics zur Regelungstechnik
• Führung technischer Prozesse erfordert und kombiniert – Planung: Konzeption, Definition des Soll-Zustands – Steuerung: Anstreben eines Soll-Zustands durch Einstellung von Steuerparametern (Stellgrößen) * Beispiel: immer mittags zusätzlichen Kühlkreislauf zuschalten – Regelung: Überwachung und Minimierung von Soll-Ist-Abweichungen (Echtzeit-Aufgabe)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* Beispiel: Temperatur im Reaktor konstant halten • äußere und innere Regelschleifen: – äußere: langsame Reaktion, größere Abweichungen erlaubt * Realisierung in Software möglich * Beispiel: Kontostand auf Girokonto – innere: zeitkritisch, nur kleine Abweichungen erlaubt * Realisierung in Hardware erforderlich * Beispiel: Temperatur in chemischem Reaktor konstant halten
Page 58 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Notation • Regelung:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Ziel: Stellgröße u so, dass sich möglichst geringe Ist-Soll-Diskrepanz ergibt • mögliche Optimierungsstrategien / Optimalitätskriterien: – Abweichung schnellstmöglich gegen Null – Abweichung bleibt in vorgegebener Bandbreite – Abweichung im statistischen Mittel Null – geringstmögliche Regelkosten (bspw. Energiebedarf für das Einstellen der Stellgrößen)
Page 59 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Wichtige Typen von Reglern • P-Regler: – enthält nur sogenannte Proportional-Komponente (Stellgröße proportional zur Soll-Ist-Abweichung): u ∼ ∆y = ySoll − yIst – bleibende Regelabweichung ∆y = 0 möglich • PID-Regler:
Populationsdynamik: . . .
– mächtiger, besteht aus Proportional-Komponente (wie zuvor), Integralund Differential-Komponente
Regelungstechnik: . . .
– Integral-Komponente garantiert verschwindende Regelabweichung:
Wärmeleitung: Modell . . .
Verkehrsfluss: . . .
∆y → 0 – Differentialglied berücksichtigt Änderung der Regelabweichung – insgesamt nach Diskretisierung (diskrete Punkte ti , ∆ti = ti − ti−1 ) : u(ti ) = KP · ∆yi + KI · ∆ti ·
∆yj + KD ·
j
∆yi = ySoll (ti ) − yIst (ti ) – entscheidend: Auslegung der Parameter KP , KI , KD (Modell)
∆yi − ∆yi−1 , ∆ti Page 60 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
PID-Regler • Proportional-Komponente: – Regelmaßnahme orientiert sich an bzw. ist proportional zur aktuellen Soll-Ist-Abweichung – soll lokaler Abweichung entgegenwirken • Integral-Komponente: – Regelmaßnahme orientiert sich an bzw. ist proportional zur Summe der bisher aufakkumulierten Soll-Ist-Abweichungen
Populationsdynamik: . . .
– soll Abweichung in der Summe (auf längere Sicht) bekämpfen
Verkehrsfluss: . . .
Regelungstechnik: . . .
Wärmeleitung: Modell . . .
• Differential-Komponente: – Regelmaßnahme orientiert sich an bzw. ist proportional zur Änderung der Soll-Ist-Abweichung von einem Zeitschritt zum nächsten – soll Aufschaukeln und Oszillationen verhindern • Frage: Braucht man alle drei Komponenten? Dazu schauen wir ein Beispiel an!
Page 61 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.2.2.
Beispiel: Lineare Rückführregelung
• Modell: zu regelnder Prozess als homogenes System linearer ODE – – – –
Vektor der Zustandsvariablen x(t) ∈ IRn Koeffizientenmatrix (zeitinvariant) A ∈ IRn,n Vektor der Stellgrößen u(t) ∈ IRm Matrix der Stellkoeffizienten (zeitinvariant) B ∈ IRn,m
• Variante I (ohne Regelkreis): x(t) ˙ = A · x(t) (vgl. Wettrüstmodell auf Folien 12 ff.; dort war x(t) ˙ = A · x(t) + c) • Variante II (mit Regelkreis):
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
x(t) ˙ = A · x(t) + B · u(t) • lineare Rückführregelung: Stellgrößen u linear aus Zustand x – Rückführmatrix K ebenfalls zeitinvariant: u = −K · x, K ∈ IRm,n – Bestimmung von K : Erfüllung von Regelziel und Nebenbedingungen • somit zu lösen: – wie zuvor System linearer ODE (aber mit anderer Systemmatrix) – Ziel: stationärer Grenzzustand (Ableitung 0, o.B.d.A. auch Soll 0) x(t) ˙ = A · x(t) + B · (−K · x(t)) = (A − BK) · x(t)
Page 62 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Lösung des ODE-Systems • zu lösen ist das System x(t) ˙ = M · x, M = A − B · K, x(0) = x0 • Standard-Vorgehensweise bei homogenen linearen ODE (einfach!): – Lösungsansatz x(t) = ν · eλt , ν ∈ IRn konstant – Einsetzen in ODE ergibt λ · ν · eλt = M · ν · eλt ⇔ M ν = λν – λ, ν sind also als Eigenwert und Eigenvektor von M zu wählen – n (komplexe) Eigenwerte λi gibt’s immer; falls auch n linear unabhängige Eigenvektoren νi vorliegen (M symmetrisch oder paarweise verschiedene Eigenwerte, z.B.), ergibt sich als Lösung x(t) =
n
νi · e
λi t
mit x0 =
i=1
n
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
νi
i=1
– sonst: komplizierter, es geht aber immer: x(t) = et·M · x0 – Charakteristik der Eigenwerte hat ablesbare Konsequenzen: * Realteile aller Eigenwerte negativ: Lösung stabil/stationär/gegen 0 * mindestens ein positiver Realteil: exponentiell wachsende Komponente * alle Realteile Null: periodische Schwingungen der Lösung * alle Realteile negativ, alle Imaginärteile Null: Idealfall ohne Oszillationen
Page 63 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Auslegung des Reglers (Matrix K ) regeln heißt jetzt: Rückführmatrix K so bestimmen, dass M = A − BK „besser“ als A, d.h.: • möglichst negative Realteile aller Eigenwerte von M – Wirkung: Abklingen • möglichst große Absolutbeträge der Realteile aller Eigenwerte von M – Wirkung: rasches Abklingen • möglichst kleine Absolutbeträge der Imaginärteile aller Eigenwerte von M – Wirkung: keine hochfrequenten Oszillationen
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Stabilität sollte auch bei Einführung kleiner Störgrößen (ist dann entsprechend komplizierteres Modell) noch gegeben sein! – gibt das Modell der linearen Rückführregelung per se natürlich nicht her (enthält ja keine Störgrößen) • Beispiel: PID-Regler für lineares ODE-System (verfügbar als Maple-Worksheet) Page 64 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.2.3.
Regelung mittels Fuzzy-Logik
• vorneweg ein paar Definitionen: – scharfe Menge (crisp set) X : Menge von Elementen x ∈ X im üblichen Sinn (endlich oder unendlich) – unscharfe Menge (fuzzy set) A˜ über X : charakterisiert über eine sogenannte Zugehörigkeitsfunktion ˜ : X → [0, 1] µ(., X, A) diese gibt für jedes Element x ∈ X die Wahrscheinlichkeit an, mit der es ˜ + µ(x, X, B) ˜ > der unscharfen Menge A˜ angehört (Achtung: µ(x, X, A) 1 ist möglich!) – man schreibt die unscharfe Menge auch
˜ ,x ∈ X , A˜ = x, µ(x, X, A)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
also als Menge von Paaren (Element, Wahrscheinlichkeit) – als Spezialfall gilt dabei: µ(x, X, X) = 1 ∀x ∈ X,
µ(x, X, A) = 1 ∀x ∈ A
˜ = x ∈ X : µ(x, X, A) ˜ >0 – Trägermenge (support) von A˜ : supp(A) – α -Schnitt (-level-set, -cut):
˜ : x ∈ X und µ(x, X, A) ˜ ≥α A˜α = x, µ(x, X, A)
Page 65 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Fuzzy-Logik • noch ein Begriff: – gekappte unscharfe Menge A˜ ↑ α : ˜ ˜ ≤α µ(x, X, A), falls µ(x, X, A) µ(x, X, A˜ ↑ α) = α sonst die Wahrscheinlichkeiten werden also nach oben beschränkt • Modellieraufgabe v.a.: Festlegung der Gestalt (Form und Werte) der Zugehörigkeitsfunktion – 0-1-Sprungfunktion: herkömmliche Mengen – viele andere Formen sind denkbar (stückweise linear, stückweise Polynome höheren Grades, ...)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• verfügbare Fuzzy-Logik-Systeme bieten verschiedene (meist einfache) Varianten an • Wir schauen uns exemplarisch ein paar Möglichkeiten an.
Page 66 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Formen der Zugehörigkeitsfunktion • scharfe Menge: Rechteck! • Dreieck (oft gleichschenklig, aber nicht immer):
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Trapez (meist symmetrisch, aber nicht immer):
Page 67 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Formen der Zugehörigkeitsfunktion • S- bzw. Z-Form (typisch z.B. für einseitig abgegrenzte unscharfe Begriffe):
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Glockenkurve (verwendet bei beidseitig abgegrenzten unscharfen Begriffen):
Page 68 of 165 Modellbildung und Simulation
• verbreitet sind ferner Polygonzüge, Sinus (halbe Periode), etc.
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Operationen auf unscharfen Mengen • üblicherweise benötigt und definiert: Pendants zu den Standard-Mengenbzw. logischen Operationen – Durchschnitt – Vereinigung – Komplement • großer Spielraum bei der konkreten Realisierung (an Anwendung anpassbar) • gewisse Einschränkungen:
Populationsdynamik: . . . Regelungstechnik: . . .
– z.B.: übliche Regeln bzw. Operationen der Logik sollten in den „scharfen“ Grenzfällen noch gelten µ ≡ 1,
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
µ≡0
Page 69 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Definition der Operationen =A ∩B über • Durchschnitt C
= min µ(x, X, A), µ(x, X, B) µ(x, X, C)
=A ∪B über • Vereinigung C
= max µ(x, X, A), µ(x, X, B) µ(x, X, C) Populationsdynamik: . . .
über =A • Komplement C
= 1 − µ(x, X, A) µ(x, X, C)
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 70 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Linguistische Variable • Eine linguistische Variable ist charakterisiert durch ihren Namen v und durch ihre möglichen Werte (Ausprägungen, linguistische Terme). – Jeder linguistischen Variablen ist eine scharfe Menge X zugeordnet. – Die Menge aller Werte von v heißt Termmenge T (v) . – Jeder linguistische Term ist eine unscharfe Menge, definiert über der scharfen Menge X . • Beispiel 1: linguistische Variable „Temperatur“ – zugeordnete scharfe Menge: reelle Achse (Temperaturwerte) – linguistische Terme: eiskalt – kalt – normal – lauwarm – warm – heiß
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Beispiel 2: linguistische Variable „Farbe“ – zugeordnete scharfe Menge: Frequenzbereich 400nm ... 800nm – linguistische Terme: green – yellowish-green – yellow-green – greenishyellow – yellow (vgl. Farbmodell CNS in der Computergraphik)
Page 71 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Bsp.: Linguistische Variable Temperatur • Name: Temperatur • zugehörige scharfe Menge: Temperaturskala (also etwa die reellen Zahlen oder ein Teilintervall) • zugehörige linguistische Terme: siehe vorige Folie • jeweilige Zugehörigkeitsfunktionen (jeder Wert/Term hat als unscharfe Menge eine):
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 72 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Fuzzy-Regelung • heute weit verbreitet, insbesondere bei Konsumgütern wie Waschmaschinen oder Fotoapparaten • Vorteile einer Fuzzy-Regelung : – mathematisches Modell ist oftmals für die zu beschreibenden Prozesse nicht vorhanden oder zu komplex (wir haben ja zuvor nur den einfachen Fall eines Systems homogener linearer ODE betrachtet – die Welt ist jedoch wesentlich komplizierter), erfordert also erheblichen mathematischen Aufwand zu seiner Beherrschung – Fuzzy-Regelung ist dagegen mit Schulmathematik möglich – Fuzzy-Regelung ist intuitiv: linguistische / umgangssprachliche Begriffe zur Nachbildung des WENN-DANN- Regelungsvorgangs
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• beachte: alles, also insbesondere alle Prämissen und alle Folgerungen, wird unscharf ausgewertet! • im Folgenden: weg vom scharfen Messwert über den Umweg des unscharfen Modellapparats schließlich zur scharfen Stellgröße
Page 73 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Entwurfsschritte • 1. Schritt: Fuzzifizierung (!) der gemessenen Größe(n) (also z.B. Temperatur), d.h. Definition der linguistischen Variablen sowie ihrer jeweiligen Terme samt den entsprechenden Zugehörigkeitsfunktionen • 2. Schritt: Erstellung der Regelbasis (Beziehungen auf WENN-DANN-Basis etc.) auf der Grundlage von Expertenwissen • 3. Schritt: Auswahl geeigneter Inferenzoperatoren (Übertragung der Unschärfe von den Messgrößen auf die Stellgrößen)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• 4. Schritt: Defuzzifizierung, d.h. Berechnung der scharfen Stellgrößen • Beispiel nach Bothe, Fuzzy Logic, Springer 1993 – Aufgabe: Regelung einer Kühlventilstellung – Messgröße: Temperatur T – Messwert: gemessene Temperatur t – Stellgröße bzw. Stellwert: Stellung K bzw. k des Kühlventils
Page 74 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 1: Fuzzifizierung • allgemeines Prinzip: – gemessene Werte stammen aus scharfen Mengen (linguistische Variable) – zu jeder linguistischen Variablen gibts mehrere unscharfe Mengen (linguistische Terme, gemäß menschlicher Wahrnehmung) – für jeden Messwert werden zu allen linguistischen Termen die Zugehörigkeitswahrscheinlichkeiten ermittelt • am konkreten Beispiel: – Temperaturwert t wurde gemessen
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– zwei linguistische Variable, zwei scharfe Mengen: * Temperatur T * Stellung K des Kühlventils
Page 75 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 1: Fuzzifizierung (2) • (einige) unscharfe Mengen zu den Werten der linguistischen Variablen „Temperatur“:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• (einige) unscharfe Mengen zu den Werten der linguistischen Variablen „Kühlventilstellung“:
Page 76 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 2: Regelbasis • allgemeines Prinzip: mit linguistischer Variablen X – Elementarbedingungen der Art X = A und Wert (unscharfer Menge) A – logische Verknüpfung von Elementarbedingungen mittels der üblichen Operatoren AND, OR und NOT mit linguistischer Variable Y und Wert – ausgelöste Aktion der Art Y = B (unscharfer Menge) B
Populationsdynamik: . . .
– somit Regeln der Art
Regelungstechnik: . . .
1,j AND X2 = A 2,j THEN Y = B j IF X1 = A
(Regel j)
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
mit entsprechenden linguistischen Variablen und Termen • am konkreten Beispiel: – Regel 1: IF T = niedrig THEN K = halb offen – Regel 2: IF T = mittel THEN K = fast offen
Page 77 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 3: Inferenz • allgemeines Prinzip: – berechne zum gemessenen Wert x ∈ X und zu jeder Elementarbedin die Wahrscheinlichkeit der Zugehörigkeit µ(x, X, A) gung X = A – Verknüpfungen der µ -Werte gemäß: * AND als Durchschnitt unscharfer Mengen (Minimum der Zugehörigkeitswahrscheinlichkeiten) * OR als Vereinigung unscharfer Mengen (Maximum der Zugehörigkeitswahrscheinlichkeiten) * NOT als Komplement einer unscharfen Menge (1 minus Zugehörigkeitswahrscheinlichkeit)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Damit ist für die linke Seite jeder Regel (der IF-Teil) ein µj berechenbar. – Somit ist nun der IF-Teil aller Regeln mit einer unscharfen Zugehörigkeit (Zugehörigkeitswahrscheinlichkeit) belegt.
Page 78 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 3: Inferenz (2) • allgemeines Prinzip (Fortsetzung): – nahe liegende Inferenz-Strategie für jede Regel: der Wert µj des IF-Teils der Regel beschränkt die Zugehörigkeitswahrj des THEN-Teils scheinlichkeit (den µ -Wert) der unscharfen Menge B der Regel – Begründung: Eine nur unwahrscheinliche Prämisse zieht die nur unwahrscheinliche Ausführung der Konklusion nach sich. – Somit wird für jede Regel j die bezüglich µj gekappte unscharfe Menge j ↑ µj berechnet. B • Daraus ist dann im vierten und letzten Schritt mit Hilfe der Defuzzifizierung wieder ein scharfer Ventilstellungswert zu ermitteln.
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 79 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 3: Inferenz (3) • jetzt wieder am konkreten Beispiel: – Regel 1 liefert µ1 = µ(t, T, niedrig) – daraus ergibt sich: halb offen ↑ µ1
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Regel 2 liefert µ2 = µ(t, T, mittel) – daraus ergibt sich: fast offen ↑ µ2
Page 80 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 4: Defuzzifizierung • allgemeines Prinzip: – Zu jeder im THEN-Teil einer Regel auftretenden linguistischen Variablen werden alle entsprechenden Regeln betrachtet. – Aus den hierzu in Schritt 3 ermittelten gekappten unscharfen Mengen wird jeweils die (unscharfe) Vereinigung gebildet. – verschiedene Möglichkeiten zur anschließenden Berechnung scharfer Stellwerte: * * * * *
Mittelwert aus Maximal- und Minimalwert des Trägers der unscharfen Vereinigung Mittelwert aus den Minimal- und Maximalwerten der Träger der einzelnen gekappten unscharfen Mengen Schwerpunkt der Vereinigungsfläche gewichteter Schwerpunkt (bei gewichteten Regeln) ···
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 81 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 4: Defuzzifizierung (2) • Formeln zur Schwerpunktbildung: – Bezeichnungen: * * * * *
x : ein Wert des Ausgabebereichs (Stellwert) µ(x) : Zugehörigkeitsfunktion zur unscharfen Vereinigung gj : Gewicht der Regel j µj : Ergebnis der linken Seite der Regel j : Fläche der Zugehörigkeitsfunktion der unscharfen Menge F (A)
– x-Koordinate xS des Schwerpunkts: x · µ(x) dx xS = µ(x) dx
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– x-Koordinate des gewichteten Schwerpunktmittels: j xS,j · ωj mit xS = j ωj j ↑ µj ) · µj · gj ωj = F (B
Page 82 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schritt 4: Defuzzifizierung (3) • und wieder am konkreten Beispiel: – berechne Schwerpunkt der Vereinigung der beiden gekappten unscharfen Mengen halb offen ↑ µ1 und fast offen ↑ µ2 – Schwerpunkt sei S = (kS , µS ) – Ergebnis: Ventil ist auf den (scharfen) Wert kS einzustellen! Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Anmerkung: Es gibt auch adaptive Fuzzy-Regelungen: – Anpassung der Gewichte bzw. der Zugehörigkeitsfunktionen
Page 83 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.2.4.
Beispiel: Balancierter Stab
• klassisches Beispiel zur Demonstration der Regelung: balancierter Stab bzw. inverses Pendel
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• klassische Modellierung: – vier Gleichungen (Horizontal-, Vertikal- und Drehbewegung des Stabschwerpunkts, Horizontalbewegung des Wagens) – führt auf System von vier linearen ODE in Unbekannten ϕ, ϕ, ˙ s, s˙ Page 84 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Balancierter Stab (2) • auch bei linearer Rückführregelung relativ kompliziert • Alternative: Fuzzy-Logik-Regelung • linguistische Variable: – Winkel ϕ mit linguistischen Termen (z.B) negativ, null, positiv – Winkelbeschleunigung ϕ˙ mit linguistischen Termen negativ , null, positiv – Stellkraft u mit linguistischen Termen groß-negativ, negativ, null, positiv, groß-positiv • neun Regeln, einfache AND-Verknüpfung von Winkel und Winkelbeschleunigung; z.B.: IF ϕ = negativ AND
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
ϕ˙ = positiv THEN u = positiv
Page 85 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Balancierter Stab (3) • alle neun Regeln:
ϕϕ ˙
negativ
null
positiv
negativ
groß-positiv
positiv
negativ/null
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
null
positiv
null
negativ
positiv
positiv/null
negativ
groß-negativ
Wärmeleitung: Modell . . .
Page 86 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Hardware-Realisierungen • Fuzzy Regler als Chips – oft nur Regler mit zwei Eingängen (Variable x und ihre Zeitableitung) sowie einem Ausgang (Stellgröße u ) nötig – hierfür Spezialchips (müssen sehr billig sein für Massenmarkt, z.B. bei Konsumgütern) – zeitkritische Regelungen auf Fuzzy-Basis erfordern Hochleistungs-CoProzessoren
Populationsdynamik: . . .
– typische Eigenschaften (bzw. was heute geht):
Regelungstechnik: . . .
* interne Parallelauswertung von Regeln * bis zu 8 (>256) Eingangsgrößen, bis zu 8 (>64) Ausgangsgrößen * bis zu 7 Terme mit Zugehörigkeitsfunktion pro linguistischer Variable * bis zu 16 Interpolationspunkte pro Zugehörigkeitsfunktion * bis zu 256 (>16384) Regeln * Auswertezeit zwischen 0.1 und maximal 50 Millisekunden (1.25 Mikrosekunden)
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 87 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zusammenfassung • Folgende Prozesskette wurde realisiert: – scharfen Messwert für jede Messgröße ablesen – Fuzzifizierung: Definition der unscharfen Welt – Regelsystem: Folgerungen in der unscharfen Welt aufstellen – Inferenz: * zu jeder Regel mit Hilfe des Messwerts einen Zugehörigkeitswert zum IF-Teil der Regel ermitteln * damit unscharfe Menge des THEN-Teils der Regel kappen
– Defuzzifizierung:
* unscharfe Vereinigung jeweils aller gekappten unscharfen Mengen zu einer linguistischen Variablen * Schwerpunktbildung der unscharfen Vereinigung liefert schließlich einen scharfen Stellwert für jede linguistische „THEN-Variable“
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– scharfen Stellwert für jede Stellgröße einstellen
Page 88 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.3.
Verkehrsfluss: Modellierung über kontinuierliche Größen 3.3.1.
Einführung
• strategisches Verhalten bei hoher Verkehrsdichte: – einzelner Autofahrer: * individuelles Optimum aus Reduzierung der eigenen Fahrzeit und Vermeidung der eigenen Verwicklung in Unfälle * intuitiver positiver Beitrag zu störungsfreiem Verkehrsfluss insgesamt (also auch für die anderen Verkehrsteilnehmer) schwierig – Polizei bzw. Verkehrsleitstelle:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* Vermeidung von Staus und Unfällen (Gesamtsicht) * Mittel: Straßenbau, Verkehrsleitsysteme, flexible Richtungsnutzung von Fahrspuren, Ampelanlagen etc. • erforderlich: Modellierung des Verkehrsflusses – um auch komplexe Situationen verstehen zu können – um sinnvoll steuernd sowie regelnd tätig werden zu können – Wechselwirkungen von Verkehrsdichte und Geschwindigkeit, im Großen (Durchsatz) und im Detail (Fortpflanzung kurzer Stockungen) – wir werden studieren: einfaches Modell nach Lighthill & Whitham
Page 89 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Mögliche Beschreibungsarten • Vielzahl verschiedener möglicher Ansätze: – Straßennetz als Graph bzw. Netzwerk, Dynamik (Autos im System) mit Petri-Netzen – Straßennetz als zellulärer Automat, Dynamik (Verkehr) über geeignete Regelbasis – Straßennetz als Wartenetz im Sinne von Abschnitt 2.3: Ampeln bzw. allgemeine Kreuzungen als elementare Wartesysteme, Zu- und Abfluss von Autos als stochastische Prozesse – Autos als Elementarteilchen (Moleküle wie in Molekulardynamik, oder stochastisch mittels Aufenthaltswahrscheinlichkeiten wie bei Boltzmann, Schrödinger & Co.)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Straßennetz als Kanalsystem, Gesamtverkehr als zähes Fluid, das sich durch die Kanäle quält – Staus, Ampeln: Wellenausbreitung, Akustik – ···
Page 90 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Verkehrsfluss: Analogien • breite Anwendbarkeit solcher Verkehrsmodelle: – Gleichungen zur Verkehrsmodellierung können z.B. auch die Ausbreitung von Verschmutzungen in Strömungsfeldern (Gewässer, Luft) beschreiben – Ausbreitung von Staus analog * zur Ausbreitung von Schockwellen (Düsenflugzeuge, Schüsse etc., Druckwelle nach Explosionen) * zur Fortpflanzung von Feuerfronten bei Waldbränden * zur Versickerung von Wasser im Boden nach starken Regenfällen * zur Bewegung von Eisenbahnwaggons beim Rangieren
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• wir streifen kurz – stationäre Verkehrsströme (im Gleichgewichtszustand) – dynamische Änderungen und singuläre Störungen • dabei durchwegs kontinuierliche Beschreibung („Strömungsgrößen“)
Page 91 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Rolle des Individuums • menschliches (Fahr-) Verhalten ist – variabel und anpassungsfähig – auch von psychologischen Faktoren bestimmt – kaum exakt (im Sinne eines mathematisch-physikalischen Modells) beschreibbar • Einfluss des Individuums unterschiedlich: – einspurige Fahrbahn, wenig Verkehr, Überholverbot: * der Langsamste bestimmt das Tempo * Durchschnittsgrößen wenig sinnvoll – Autobahn, starker Verkehr:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* Einfluss einzelner Autos gering * Verwendung von Durchschnittsgrößen sinnvoll • Schwerpunkt von Modellierung und Simulation: – starker Verkehr – mehrspurige Straßen
Page 92 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.3.2.
Gleichmäßiger Verkehrsfluss
• Studium der Verkehrsströme, nicht individueller Fahrzeugbewegungen – außenstehender Beobachter beschreibt (Euler-Beschreibung) – Alternative wäre Lagrange-Beschreibung (teilchenartig, Bewegungen der Autos werden erfasst) • Szenario (zunächst ziemlich unrealistisch): – Autos alle gleich lang ( l ) und gleich schnell, gleicher Abstand – gleichförmiger, stationärer, einspuriger Verkehrsfluss
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• drei Größen: – Abstand zweier Autos (Straßenlänge normiert auf 1):
Wärmeleitung: Modell . . .
1−N ·l N −1 – Verkehrsgeschwindigkeit V (km/h): Tempo der Autos – Verkehrsdichte N (Autos/km): Anzahl der Autos pro Strecke – Verkehrsfluss F (Autos/h): Durchsatz an einem Kontrollpunkt Page 93 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Eine erste Beziehung: F, N und V • fundamentale Beziehung zwischen Fluss, Geschwindigkeit, Dichte (absolut oder gemittelt): F = N · V – vgl. Formel von Little aus Abschnitt 2.3: Füllung = Verweilzeit mal Durchsatz * Durchsatz entspricht Fluss F * Füllung entspricht Dichte N −1 −1 * Verweilzeit entspricht Geschwindigkeit V – neu ist lediglich die Raumdimension („Strecke“) – Veranschaulichung:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 94 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Eine zweite Beziehung: V und N • charakterisiere zunächst V als Funktion von N : – V streng monoton fallend in N * höhere Geschwindigkeit erfordert höheren Sicherheitsabstand und somit geringere Dichte * im Grenzfall (Maximaldichte, Abstand 0): V →0
für N → Nmax =
1 l
Regelungstechnik: . . .
– Höchstgeschwindigkeit auf leerer Straße: V → Vmax
Populationsdynamik: . . .
Verkehrsfluss: . . .
für N → 0
Wärmeleitung: Modell . . .
– Linearisierung (Vereinfachung): V = Vmax · 1 −
N Nmax
• Verkehrspolitik hätte lieber Abstand und somit N als Funktion von V : – Bremsweg wird maßgeblich von V bestimmt! – vgl. die plakatierte Regel „Abstand halber Tacho“ – Ansatz liefert leider keine vernünftigen Modelle für das Fahrverhalten
Page 95 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Eine dritte Beziehung: F und N • aus dem Gesagten folgt F = F (N ) = V (N ) · N • im Fall der linearen V -N -Abhängigkeit: F (N ) = Vmax · N · 1 −
N Nmax Populationsdynamik: . . .
parabolisches Modell:
Regelungstechnik: . . .
F →0
für N → Nmax , N → 0
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 96 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Gleichgewichtszustände • zwei Möglichkeiten, einen Fluss F < Fmax zu erreichen: – Möglichkeit 1: kleine Dichte, große Geschwindigkeit – Möglichkeit 2: große Dichte, kleine Geschwindigkeit • Möglichkeit 1 ist attraktiver: – für einzelnen Autofahrer, weil schneller am Ziel – für Verkehrsplaner, weil Autofahrer glücklicher
Populationsdynamik: . . . Regelungstechnik: . . .
• tritt Möglichkeit 1 in Praxis ein?
Verkehrsfluss: . . .
– bei wenig Verkehr ja (alle geben Gas)
Wärmeleitung: Modell . . .
– bei viel Verkehr eher nein • jetzt Steuerung/Regelung gefragt: – erzwinge die „bessere“ Alternative – Mittel: Ampeln, Tempohinweise ...
Page 97 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Alternative: kubisches Modell • parabolisches Modell war:
N F (N ) = Vmax · N · 1 − Nmax
• unbefriedigend u.a., weil (oft im Gegensatz zur Beobachtung) – maximaler Fluss starr bei halber Maximaldichte – maximaler Fluss starr bei halber Maximalgeschwindigkeit • daher kubisches Modell:
Populationsdynamik: . . . Regelungstechnik: . . .
F (N ) = Vmax · N · 1 − αN − βN 2
Verkehrsfluss: . . .
Wärmeleitung: Modell . . .
• Bedingungen: – per Ansatz wieder erfüllt: F (N ) → 0 für N → 0
– ebenfalls per Ansatz wieder erfüllt: F (N ) → Vmax
für N → 0
– außerdem soll wieder gelten: F (N ) → 0 für N → Nmax
– schließlich: F (Nopt ) = 0, Nopt · Vopt = Fmax – hierbei Vopt beobachtbar (Geschwindigkeit mit maximalem Fluss) – somit drei zusätzliche Gleichungen (nicht alle linear) für die drei Unbekannten α, β, Nopt
Page 98 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.3.3.
Instationäre Situationen
• jetzt: Dichte variiert in Ort und Zeit (realistischer) • Frage: Wie entwickelt sich der Fluss? – als kontinuierliche Größe: F = F (x, t) – nicht Bewegung einzelner Autos, sondern Verkehrsstrom als Ganzes • Annahmen (zur Erleichterung): – vernünftige Zeit- und Längenskalen seien a priori wählbar
Populationsdynamik: . . .
– weiterhin gelte F = F (N ) = V (N ) · N (quasi-stationär )
Regelungstechnik: . . .
* verzögerungsfreie Anpassung der Autos an Verkehrssituation * damit außen vor: Situation beim Lichtwechsel an Ampeln – kontinuierliche Beschreibung sei möglich (kein Ärger mit Glattheit (Stetigkeit, Differenzierbarkeit) etc.)
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• jetzt: kontinuierliche Beziehungen herleiten
Page 99 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein Erhaltungssatz • geschlossene Teststrecke von x = a nach x = b – keine Auf- oder Abfahrten auf der Teststrecke (Erhaltung) – unidirektionaler Verkehr: Zugänge bei a , Abgänge bei b – kontinuierlich betrachtete kontinuierliche Größen: * Dichte N (x, t) am Ort x zur Zeit t * Fluss F (x, t) am Ort x zur Zeit t * Anzahl der Autos auf der Strecke M (t) zur Zeit t • dann gilt – zu jedem Zeitpunkt:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
x=b
M (t) =
N (x, t) dx x=a
– mit fortschreitender Zeit: ∂ M (t) = Mt (t) = ∂t
x=b
Nt (x, t) dx x=a
– verantwortlich für Änderungen von M (t): Fluss in a und b x=b ∂ Fx (x, t) dx M (t) = Mt (t) = F (a, t) − F (b, t) = − ∂t x=a
Page 100 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein Erhaltungssatz (2) • Visualisierung der Berechnung von M (t) :
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• zwei Wege zur Berechnung der Zeitableitung von M (t) : – via Integral der Zeitableitung der Dichte über ganze Strecke – oder via Zu- und Abgangsbilanz an den Endpunkten a und b
Page 101 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein Erhaltungssatz (3) • insgesamt gilt somit
x=b
Nt (x, t) dx = −
x=a
bzw.
x=b
Fx (x, t) dx x=a
x=b
(Nt (x, t) + Fx (x, t)) dx = 0 x=a
und zwar
Populationsdynamik: . . .
– für jeden Zeitpunkt t – für jede Teststrecke [a, b]
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• deshalb gilt (bei hinreichender Differenzierbarkeit) Nt (x, t) + Fx (x, t) = 0 ∀x, t bzw. mit F = F (N ): Nt (x, t) + FN (N (x, t)) · Nx (x, t) = 0 ∀x, t • Verkehrsgleichung (nichtlineare partielle DGL!) – berücksichtigt Erhaltung von Autos: Dichteänderungen schlagen sich in Flussänderungen nieder und umgekehrt – Unbekannte: Dichte; daraus dann mittels F = F (N ) der Fluss
Page 102 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Verkehrsgleichung • man beachte: – Die Abhängigkeit F = F (N ) ist eine Modellannahme und somit bereitzustellen (z.B. parabolisches oder kubisches Modell). – Die Verkehrsgleichung ist ein einfaches Beispiel einer nichtlinearen hyperbolischen Wellengleichung (beschreiben allgemeine Wellenausbreitungsphänomene). – Die Ableitung des Flusses nach der Dichte, also ∂ F = FN , ∂N wird Signalgeschwindigkeit genannt:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* aus Dimensionsgründen (Einheit!) tatsächlich Geschwindigkeit * Auswirkung einer Dichteänderung (Störung) auf den Fluss, d.h. Ausbreitung der Information einer Änderung oder einer Störung * Signalgeschwindigkeit ist ungleich Autogeschwindigkeit V * Signalgeschwindigkeit kann bei V ≥ 0 positiv, negativ oder Null sein. • Verkehrsgleichung kann in einfachen Fällen analytisch gelöst werden (Charakteristiken-Methode)
Page 103 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zur Signalgeschwindigkeit • es gilt: – Signalgeschwindigkeit stets kleiner gleich Autogeschwindigkeit – Straße leer: N → 0, V → Vmax , F → 0, FN = V, d.h. * Verursachte Störungen werden mit den Verursacherautos aus dem System entfernt, die anderen Teilnehmer merken nichts. – N und F steigen an: * Störungen entfernen sich langsamer als die Verursacher. * Die nachfolgenden Autos werden beeinträchtigt. – Maximalfluss erreicht:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* Signalgeschwindigkeit wird Null, Störung verharrt an der Quelle * katastrophal (und das gerade im Idealfall des maximalen Flusses!) – Dichte steigt weiter, Fluss sinkt wieder: * negative Signalgeschwindigkeit, Information läuft rückwärts , GAU * Stau (bei hoher Dichte besonders große negative Signalgeschwindigkeit)
Page 104 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zur Signalgeschwindigkeit (2) • grafische Veranschaulichung:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 105 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispielphänomene • Losfahren von Autos an der Ampel: – Signalgeschwindigkeit ist hier praktisch erleb-/erleidbar! – entspricht abruptem Dichteabfall, der aufgelöst wird (allerdings nicht sofort bzw. „auf einmal“, sondern im Laufe der Zeit) – immerhin: maximale Dichte an Ampel, also auch maximale negative Signalgeschwindigkeit! – Eingriffsmöglichkeiten: Position der Ampel, Synchronisierung aufeinander folgender Ampeln, Tempolimits etc. • Entstehung und Auflösung von Staus:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– eine Spur blockiert: plötzlicher Dichtesprung nach oben bzw. bei optimalem Einfädeln – Halbierung des Flusses – zweite Spur wird wieder frei: Dichtesprung nach unten
Page 106 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Graphische Veranschaulichung • Dichte und Geschwindigkeit am Anfang
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Verkehrsverhältnisse bewegen sich fort mit FN
Page 107 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.4.
Wärmeleitung: Modell und numerische Lösung
Auflösung des Raums • Differentialgleichungen bei Populationsmodellen: – nur eine unabhängige Variable, nämlich die Zeit – gewöhnliche Differentialgleichungen (ODE) • Verkehrsgleichung:
Populationsdynamik: . . . Regelungstechnik: . . .
– zwei unabhängige Variablen: Zeit und eine Raumrichtung
Verkehrsfluss: . . .
– partielle Ableitungen nach beiden Variablen, deshalb partielle Differentialgleichung (partial differential equation, PDE)
Wärmeleitung: Modell . . .
– wieder (wie bei Populationen) kontinuierliches Modell, obwohl diskrete Welt (Autos) • Einsatz der Modellwerkzeuge ODE und PDE durchaus auch alternativ: – höhere Genauigkeit bei PDE (räumliche Heterogenitäten) – höherer Aufwand bei PDE (Diskretisierung von Zeit und Raum) – Beispiel: Populationsdynamik
Page 108 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Raumauflösende Populationsmodelle • einerseits: rein zeitabhängige Modelle manchmal zu grob – Bevölkerungsentwicklung in den USA in den 1850ern (California gold rush): stark ausgeprägte räumliche Ost-West-Komponente – Migration der Weltbevölkerung im 21. Jahrhundert: auch hier räumlicher Aspekt, prognostizierte Süd-Nord-Wanderung – Vorhersage von Heuschreckenplagen in Afrika: flächige Ausbreitung, diffusive und andere Effekte • deshalb Zielgröße eher p(x, t) oder p(x, y, t) anstelle von p(t) • Erhöhung der Komplexität:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Modelle: Wie hängen räumliche und zeitliche Ableitungen zusammen (Wellengleichung wäre eine Möglichkeit)? – Numerik: Erhöhung von Speicher- und Rechenaufwand u.v.m. • andererseits: höhere Genauigkeit überhaupt erforderlich? – Europa: viele Auswanderer (egal, wo das Gold gefunden wird) – USA: wichtig zu wissen, wo Infrastruktur bereitzustellen Page 109 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Systemschicht vs. Prozessebene • Systemschicht: – Interesse nur an makroskopischen Phänomenen wie Wirkungsgrad, Kosteneffizienz etc. • Prozessschicht: – Verständnis der die Makroebene bestimmenden mikroskopischen Prozesse • Gretchenfragen: – Was aus der Mikrowelt brauche ich nur zum Zwecke der Erkenntnis? – Was aus der Mikrowelt darf ich nicht ohne signifikante Konsequenzen für die Makrowelt vernachlässigen?
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Beispiel: Bauelemente und Schaltkreise * für Standard-Schaltplan reicht Kirchhoff; keine „Elektronen“ * für EMV (elektromagnetische Verträglichkeiten, z.B. Wechselwirkungen zwischen zwei benachbarten Leiterbahnen) reicht Kirchhoff nicht Page 110 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.4.1.
Modellierung mit PDE
• Raumauflösung ist essentiell für viele Probleme bzw. Phänomene aus Physik oder Kontinuumsmechanik: – Strömungsmechanik/Thermodynamik: * Wo entsteht ein Tornado? * Ist die gegebene Karosserie aerodynamisch günstig? – Strukturmechanik? * Hält das Gebäude einer Belastung stand? * Wo sind Sollbruchstellen? – Verfahrenstechnik: * Wo wird es wie heiß in einem (nuklearen / chemischen / biologischen) Reaktor?
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Elektromagnetismus: * Wo im Transistor ist die Elektronendichte wie hoch? – Geologie: * Wann, wo und wie heftig wird es Erdbeben geben? Page 111 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Mit oder ohne Zeit? • stationäre Phänomene: – keine Zeit-, sondern nur Ortsabhängigkeit – bei eindimensionalem Raum: wieder ODE! – Beispiele: Gleichgewichtssituationen * Auflösung einer löslichen Substanz in Wasser ohne äußere Einflüsse (die Dinge ändern sich zwar, streben aber zu einer stationären Lösung) * Temperaturverteilung in konstant beheiztem Raum • instationäre Phänomene:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Zeit- und Ortsabhängigkeit – auf jeden Fall PDE! – Beispiele: * Oszillationen: mechanische Belastung einer Schiffschaukel * Turbulenz (vgl. Strömung über Wehr, Wirbel im Nachlauf startender Flugzeuge) * Regelkreis (ständig angepasste Randbedingungen)
Page 112 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Beispiel: Wärmeleitung • Kernproblem der Thermodynamik – Wärme beeinflusse die Oberfläche eines Objekts – Ziel: Ausbreitung bzw. Verteilung vorhersagen • Beispiele: – ein Hitzdraht – ein Kochtopf auf einer Herdplatte
Populationsdynamik: . . .
– Kühlwasser im Reaktor eines Kernkraftwerks
Regelungstechnik: . . .
– ein Zimmer im Winter: wo die Heizung platzieren?
Verkehrsfluss: . . .
– ein Zimmer im Sommer: Aufheizung an Fensterflächen
Wärmeleitung: Modell . . .
• Objekt der Begierde i.W.: Temperatur T T (x; t) oder
T (x, y; t) oder
T (x, y, z; t)
• entscheidend: – Rand- und Anfangsbedingungen – Materialeigenschaften (Wärmeleitfähigkeit etc.)
Page 113 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Szenario • Gebiet Ω ⊂ IR3 , Gebietsrand ∂Ω , Zeitintervall [0, τ ]
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Start: Wasser kalt, Platte kalt • dann: Einschalten der Platte – Platte schnell auf Zieltemperatur, Wasser erhitzt sich allmählich – beachte Kühlwirkung des Außenraums des Topfes
Page 114 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein Modell für die Wärmeleitung • Standardbeispiel für eine physikalische Herleitung • mehrteiliges Modell: – die PDE: * beschreibt die Wechselwirkungen von Temperaturänderungen in Bezug auf Raum (hier jetzt 3D) und Zeit * eine Gleichung für eine gesuchte Funktion ausreichend * Gleichung bestimmt Schar von Lösungen
– Anfangs- und Randbedingungen:
* zur eindeutigen Festlegung der gesuchten Lösung * Anfangsbedingungen: Temperaturfeld im gesamten Gebiet zum Startzeitpunkt (wie bei ODE) * Randbedingungen: Temperaturdaten (Werte oder Änderungen) am Rand des betrachteten Gebiets zu allen Zeiten
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• beide Teile sind herzuleiten
Page 115 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Herleitung der Wärmeleitungsgleichung • die Lösung vorneweg, die Wärmeleitungsgleichung: κ · (Txx + Tyy + Tzz ) = κ ·
∂2T ∂x2
+
∂2T ∂ 2 T ∂T = + = Tt ∂y 2 ∂z 2 ∂t
oder kurz κ · ∆T = Tt mit dem Laplace-Operator
Populationsdynamik: . . .
∆ = Txx + Tyy + Tzz
Regelungstechnik: . . . Verkehrsfluss: . . .
• kurze (!) Herleitung als kleiner Ausflug in die Physik:
Wärmeleitung: Modell . . .
– Start ist das grundlegende Prinzip der Energieerhaltung – zeitliche Temperaturänderungen in einem Teilgebiet D des Gebiets sind bedingt entweder durch Wärmefluss durch die Oberfläche von D oder durch Wärmequellen und -senken im Innern von D : ∂ ρcT dV = q dV + k∇T · n dS ∂t D D ∂D – Dichte ρ , spezifische Wärme c , Quellterm q , Wärmeleitfähigkeit k, äußere Normale n, Volumenelement dV bzw. Oberflächenelement dS, Gradient ∇
Page 116 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Veranschaulichung der Energieerhaltung • idealisiertes Volumenelement D : – Wärmeänderung in D durch * A: Wärmefluss ins Element hinein bzw. aus dem Element heraus * B: Wärmequellen oder -senken im Innern
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Erhaltung der Wärme oder Umwandlung der Wärmeenergie (z.B. in kinetische Energie) Page 117 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Herleitung (2) • Transformation mit dem Gaußschen Integralsatz in ein reines Volumenintegral: (ρcTt − q − k ∆T ) dV = 0 D
• Integral verschwindet für beliebiges Teilgebiet D , daher muss Integrand Null sein: q k Tt = κ ∆T + , κ = ρc ρc • κ > 0 wird thermischer Diffusionskoeffizient genannt (Laplace-Operator beschreibt Diffusions- oder Ausgleichsprozesse) • ohne äußeren Einfluss q = 0 ergibt sich eben die Wärmeleitungsgleichung:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Tt = κ∆T • stationärer Fall: Laplace-Gleichung ∆T = 0 Page 118 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Herleitung (3) • dies ist Beispiel für Standard-Vorgehensweise bei der Modellierung in der Physik (leicht verkürzt) – Start bei einem grundlegenden physikalischen Gesetz: * der berühmte Apfel * die Erhaltung von Masse * die Erhaltung von Energie oder ... – daraus Gewinnung einer Bilanzgleichung (typischerweise ein summatorischer Zusammenhang mit Integralen)
Populationsdynamik: . . . Regelungstechnik: . . .
– Zuhilfenahme der Analysis: Vereinfachung
Verkehrsfluss: . . .
* Gaußscher Integralsatz (ein üblicher Verdächtiger) * „Integral verschwindet stets“ zieht „Integrand Null“ nach sich – Vereinfachung durch physikalische Einschränkungen
Wärmeleitung: Modell . . .
* keine Quellterme etc. – am Ende dann die gewünschte (Differential-) Gleichung
Page 119 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Typen von Randbedingungen • nun zu (einigen) möglichen Randbedingungen: – Dirichlet Randbedingungen: T auf Rand vorgegeben T (t; x, y, z) = ψ(t; x, y, z) auf
∂Ω
* Temperatur T selbst wird auf dem Rand festgeschrieben * d.h.: definiertes Heizen oder Kühlen – Neumann Randbedingungen: Normalableitung (Ableitung in Richtung der äußeren Normalen) von T auf Rand vorgegeben ∂T (t; x, y, z) = ∇T (t; x, y, z) · n = ϕ(t; x, y, z) auf ∂n
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
∂Ω
Wärmeleitung: Modell . . .
* Wärmefluss durch Rand bzw. Teile des Rands festgeschrieben * Wert Null: keine Temperaturunterschiede, folglich kein Wärmetransport in das Gebiet hinein oder aus dem Gebiet hinaus, d.h. vollständige Isolation ϕ(t; x, y, z) = 0 auf ∂Ω • Frage: Was ist physikalisch sinnvoll (oft nicht trivial)? Page 120 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Rand- und Anfangsbedingungen • möglich sind: – reine Dirichlet-Randbedingungen * Temperaturvorgabe auf dem gesamten Rand – Mix aus Dirichlet- und Neumann-Randbedingungen * teils Vorgabe der Temperatur, teils Vorgabe des Flusses – reine Neumann-Randbedingungen dagegen nicht: * reichen nicht zur Gewährleistung einer eindeutigen Lösung * mit T ist auch jedes T + const. Lösung (Funktionen sind durch ihre Ableitungen alleine noch nicht eindeutig charakterisiert) * dies physikalisch auch vernünftig: Fluss legt Absolutbetrag nicht fest!
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Anfangsbedingungen: – Wert der Temperatur zur Startzeit im gesamten Gebiet, z.B. T (0; x, y, z) = φ(x, y, z) in Ω Page 121 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zurück zum Szenario • Gebiet Ω ⊂ IR3 , Zeitintervall [0, τ ]
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Anfangsbedingung: T = const. im ganzen Topf • Randbedingungen: – Platte: evtl. Hochfahren, dann konstant heiß (also Dirichlet) – Seitenwand und oben: konstanter Wärmeabtransport (Neumann)
Page 122 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Tut’s nicht auch eine Raumdimension? • Gebiet Ω ⊂ IR3 , Zeitintervall [0, τ ]
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• perfekte Platte:
Wärmeleitung: Modell . . .
– bei z = 0 gleichmäßig (flächig) beheizt – größere Unterschiede wohl nur in der Höhe – Eine Raumdimension (z-Koordinate) scheint auszureichen. • z > 0 : Störung der Homogenität der Temperatur in einer Ebene durch Kühlwirkung am Rand? Page 123 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Reichen zwei Raumdimensionen? • Gebiet Ω ⊂ IR3 , Zeitintervall [0, τ ]
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• reale Platte:
Wärmeleitung: Modell . . .
– Beheizung von der Plattenmitte aus – Unterschiede in der Höhe, im Abstand von Mitte – Zwei Raumdimensionen (z- und r- Koordinate)! • aber: Konkrete Position (Winkel) wirklich egal?
Page 124 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Also doch drei Raumdimensionen? • Gebiet Ω ⊂ IR3 , Zeitintervall [0, τ ]
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . .
• kaputte Platte:
Wärmeleitung: Modell . . .
– Beheizung von Plattenmitte aus, nicht rotationssymmetrisch – Unterschiede in der Höhe, im Abstand von Mitte, im Winkel – Drei Raumdimensionen (z-, r- und ϕ -Koordinate)! – gleich gute, aber adäquatere Beschreibung als (x, y, z)-System
Page 125 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zur Lösung des Modells • Wärmeleitung ist beliebtes Anschauungsbeispiel: – allgemeine Aussagen zu Existenz und Eindeutigkeit von Lösungen möglich – In einfachen eindimensionalen Konfigurationen kann die Lösung explizit bestimmt werden (Separation der Variablen, Fourier-Methode). – insofern Beispiel eines sehr einfachen Rand- bzw. Rand-Anfangswertproblems • allgemein – Existenz und Eindeutigkeit von Lösungen oft offen – nur numerische Lösung möglich
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– insbesondere im höherdimensionalen Fall (drei Raumdimensionen plus Zeit) oft extrem aufwändig • Numerik von PDE eines der wichtigsten Themen der numerischen Simulation!
Page 126 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Klassifizierung: Typen von PDE • Wärmeleitungsgleichung ist einfaches Beispiel einer linearen PDE 2. Ordnung in d Dimensionen: d
ai,j (x) · Txi ,xj (x) +
i,j=1
d
ai (x) · Txi (x) + a(x) · T (x) = f (x)
i=1
– linear: nur Linearkombination verschiedener Ableitungen, keine Produkte etc. – 2. Ordnung: nur Funktion selbst, erste und zweite Ableitungen
Populationsdynamik: . . . Regelungstechnik: . . .
– d = 4 : drei Raumdimensionen und die Zeit, also x = (x, y, z; t)
Verkehrsfluss: . . .
– Wärmeleitungsgleichung: alle Koeffizientenfunktionen konstant!
Wärmeleitung: Modell . . .
κ · ∆T − Tt = κ · Txx + κ · Tyy + κ · Tzz − Tt = 0 ⎞ ⎛ ⎛ ⎞ κ 0 0 0 0 ⎜ 0 κ 0 0⎟ ⎜0⎟ ⎟ ⎜ ⎟ A = (aij )i,j = ⎜ ⎝ 0 0 κ 0⎠ , b = (ai )i = ⎝ 0 ⎠ , a = f = 0 0 0 0 0 −1 Page 127 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Klassifizierung (2) • Man unterscheidet drei Typen: – elliptische PDE: * Matrix A der Koeffizienten ai,j ist positiv oder negativ definit – hyperbolische PDE: * Matrix A hat einen positiven und d − 1 negative Eigenwerte oder umgekehrt – parabolische PDE: * ein Eigenwert von A ist Null, die anderen haben gleiches Vorzeichen * Rang von A zusammen mit dem Vektor ai ist maximal bzw. voll, d.h. Wert d (wenn nach einer Variablen keine zweite Ableitung auftritt, dann wenigstens eine erste!)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• einfache Beispiele: – elliptisch: Laplace- oder Potenzialgleichung ∆u = 0 – parabolisch: Wärmeleitungsgleichung ∆u = ut – hyperbolisch: Wellengleichung
Page 128 of 165 Modellbildung und Simulation
∆u = utt
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
3.4.2.
Numerik von PDE
• Klassifizierung wichtig: – verschiedene Klassen erfordern komplett verschiedene Numerik – hier: Beschränkung auf den elliptischen Fall, genauer auf die einfache Laplace-Gleichung ∆u = 0 • verschiedene Diskretisierungsverfahren: – Finite Differenzen (FD): * direkte Approximation aller auftretenden Ableitungsterme (vgl. ODE) * nahe liegender, direkter Ansatz * leicht zu implementieren, allerdings wenig theoretischer Hintergrund – Finite Volumen (FV):
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* direkte Implementierung der kontinuumsmechanischen Erhaltungssätze auf kleinen Volumen, insb. bei Strömungen – Finite Elemente (FEM): * Variationsansatz, betrachte eine leicht abgeschwächte Variante der PDE anstatt ihrer selbst * komplizierter in der Implementierung, aber schöne Theorie
Page 129 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Finite Differenzen Verfahren • gegeben: Gebiet Ω ⊆ IRd , d ∈ {1, 2, 3} • darauf wird definiert: (regelmäßiges) Gitter Ωh • Gitterabstand oder Maschenweite h = (hx , hy , hz ) • ersetze Ableitungen durch Differenzenquotienten: – erste Ableitungen: Vorwärts-, Rückwärts- oder zentrale Differenzen
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
∂u . u(ξ + hx ) − u(ξ) u(ξ) − u(ξ − hx ) u(ξ + hx ) − u(ξ − hx ) , , (ξ) = ∂x hx hx 2hx – zweite Ableitungen: Standard 3-Punkt-Stern (vgl. Abschnitt 3.1.2) ∂ 2 u . u(ξ + hx ) − 2u(ξ) + u(ξ − hx ) = ∂x2 h2x – Laplace-Operator in 2D oder 3D: 5-Punkt- oder 7-Punkt-Stern – Es gibt auch breitere Sterne (Involvierung von mehr Nachbarn).
Page 130 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Finite Differenzen Verfahren (2) • in jedem inneren Gitterpunkt: – setze eine Differenzengleichung an – Unbekannte (Freiheitsgrade): diskrete Näherungswerte für die Funktionswerte in den Gitterpunkten – ein Freiheitsgrad pro Gitterpunkt und pro (skalarer) unbekannter Größe • in Punkten auf dem oder nahe am Rand:
Populationsdynamik: . . .
– konkrete Gestalt der Gleichung hängt von Randbedingungen ab: * Dirichlet: keine Differenzengleichung in Randpunkten (dort ist Funktionswert ja vorgegeben, somit keine Unbekannte) * Neumann: spezielle Gestalt der Differenzengleichung am Rand (Einbau der Randbedingungen, vgl. RWP bei ODE in Abschnitt 3.1)
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Diskretisierung führt auf lineares Gleichungssystem (LGS) – im Allgemeinen dünn besetzt – schnelle (iterative) Lösungsverfahren lebenswichtig Page 131 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Finite Differenzen Verfahren (3) • einfaches Beispiel: Poisson-Gleichung auf dem Einheitsquadrat −∆u = f
on
]0, 1[2
• äquidistantes quadratisches Gitter: h = hx = hy = 1/N • Anzahl der Freiheitsgrade (Unbekannten): M = (N − 1)2 • Standard-5-Punkt-Stern ergibt lineares System Ax = b mit – M × M -Matrix A (die punktweisen Differenzengleichungen) – M -Vektoren b (rechte Seite) und x (gewünschte Werte von u )
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• A ist eine dünn besetzte Matrix (sparse matrix) mit Bandstruktur; nach Multiplikation mit h2 : – Dirichlet Randbedingungen: * alle Diagonalelemente sind 4, pro Zeile 2 bis 4 Nicht-Nullen (−1) * Randwerte auf die rechte Seite – Neumann Randbedingungen (z.B., verschiedene Möglichkeiten): * Diagonalelemente sind 2 oder 3 nahe dem Rand und 4 sonst; dementsprechend viele – 1en (2 bis 4) in jeder Zeile * Paare (1, −1) entlang dem Rand auf die rechte Seite
Page 132 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskretisierung eines Dirichlet-Rands • Dirichlet:
• Punkt auf dem Rand, Wert gegeben
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Differenzengleichung: −u(x − h, y) + 2u(x, y) − u(x + h, y) −u(x, y − h) + 2u(x, y) − u(x, y + h) + h2 h2 = f (x, y) • Matrixzeile je nach Lage (nach Multiplikation mit h2 ): (. . . − 1 . . . − 1 . . . 4 . . . − 1 . . . − 1 . . .), (. . . − 1 . . . − 1 . . . 4 . . . − 1 . . .), (. . . − 1 . . . 4 . . . − 1 . . .) • Bekanntes (Randwerte): auf rechte Seite der Gleichung bringen
Page 133 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskretisierung eines Neumann-Rands • Neumann:
z.B.
u(x,y+h)−u(x,y) h
am oberen Rand bekannt
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Differenzengleichung: u(x − h, y) − 2u(x, y) + u(x + h, y) u(x, y − h) − 2u(x, y) + u(x, y + h) + = −f (x, y) h2 h2 • Matrixzeile je nach Lage (nach Multiplikation mit h2 ): (. . . − 1 . . . − 1 . . . 4 . . . − 1 . . . − 1 . . .), (. . . − 1 . . . − 1 . . . 3 . . . − 1 . . .), (. . . − 1 . . . 2 . . . − 1 . . .)
Page 134 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Finite Differenzen Verfahren (4) • Genauigkeit typischerweise quadratisch: uberechnet − uexakt = O(h2 ) = O(N −2 ) • Fluch der Dimension: – es werden O(N d ) Punkte bei d Dimensionen benötigt • Ansatzpunkte für Verbesserungen:
Populationsdynamik: . . .
– Sterne höherer Ordnung (kubisch, quartisch, ...):
Regelungstechnik: . . .
* mehr als zwei benachbarte Punkte pro Raumrichtung heranziehen * Problem: Matrix wird voller (weniger Nicht-Nullen) – größere Sparsamkeit mit Gitterpunkten:
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* lokal verfeinerte Gitter verwenden (adaptive Gitter) * Problem: Vorgehensweise an den Nahtstellen – welcher Wert?
Page 135 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Finite Elemente Methoden • keine direkte Diskretisierung von Ableitungen, sondern Transformation der PDE in eine schwache Form (Integralform, Skalarprodukt) • fünf wesentliche Schritte: – Substrukturierung oder Gittererzeugung: zerlege das Gebiet in einzelne Parzellen endlicher Ausdehnung (finite Elemente) – schwache Form: erfülle die PDE nicht punktweise überall, sondern nur noch abgeschwächt (in Form eines Skalarprodukts) bzw. gemittelt (als Integral) – endlich dimensionaler Ansatzraum: ersetze die kontinuierliche Lösung in der schwachen Form durch geeignete diskrete Approximationen
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– System linearer Gleichungen: stelle (wieder) das zugehörige System linearer Gleichungen auf (eine Gleichung für jeden Freiheitsgrad) – Lösung des linearen Systems: Einsatz geeigneter Iterationsverfahren zur Bestimmung der Approximation der Lösung
Page 136 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Substrukturierung, Gittererzeugung • unterteile das Problemgebiet in finite Elemente: – Vorgehen aus den Ingenieurwissenschaften (Statik): * zerlege komplexe Objekte in Standardkomponenten, deren Verhalten einfach zu beschreiben ist * leite daraus das Verhalten des Gesamtobjekts ab
– in 3D erhält man ein FE-Netz (geeignete Datenstruktur!) mit * * * *
Elementen: 3D Atome (Würfel, Tetraeder, ...) Flächen: 2D Oberflächenstrukturen (Dreiecke, Quadrate, ...) Kanten: 1D Randstrukturen der Elemente Gitterpunkten oder Knoten: hier leben die Unbekannten
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Zu jedem Knoten gehört eine Ansatzfunktion ϕk : – endlicher Träger: von Null verschieden nur in Nachbarelementen – alle Ansatzfunktionen zusammen spannen den linearen und endlichdimensionalen Ansatzraum auf und bilden eine Basis – in diesem Ansatzraum Vn wird die Näherung der Lösung gesucht • einfachste Basis von Ansatzfunktionen ϕk in 1D: stückweise lineare Stützpunktbasis
Page 137 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schwache Form der PDE • L bezeichne den Differentialoperator (z.B. ∆ ) • anstelle von Lu = f
auf Ω betrachte Lu · ψl dΩ = f · ψl dΩ ∀ψl Ω
Ω
für eine Menge von Testfunktionen ψl – Methode der gewichteten Residuen oder Galerkin-Ansatz
Populationsdynamik: . . .
– somit ein weiterer linearer Raum, jetzt von diesen Testfunktionen aufgespannt der Testraum Wn
Regelungstechnik: . . .
– falls Test- und Ansatzraum
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* identisch: Ritz-Galerkin-Ansatz * verschieden: Petrov-Galerkin-Ansatz – Linearität: nur Basisfunktionen müssen schwache Form erfüllen – Schreibweise: Bilinearform a(., .) und Linearform b(.): a(u, ψl ) = b(ψl ) ∀ψl ∈ Wn – das sind n lineare Gleichungen (n : Dimension des Testraums)
Page 138 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskrete Approximation und Gleichungen • ersetze (exakte/kontinuierliche) Lösung u in obigen n Gleichungen durch eine diskrete Approximation in Vn (damit n Unbekannte ak ): un = αk ϕk ∈ Vn k
• in der schwachen Form:
αk ϕk , ψl a(un , ψl ) = a =
Populationsdynamik: . . .
k
Regelungstechnik: . . .
αk · a(ϕk , ψl ) = b(ψl ) ∀ψl ∈ Wn
Verkehrsfluss: . . .
k
Wärmeleitung: Modell . . .
• alle a(ϕk , ψl ) und b(ψl ) : – hängen nicht ab von der jeweiligen Approximation für u, sondern nur vom Problem (vgl. ihre Definition!) – ein für alle Mal zu Beginn zu berechnen – führt auf System von n linear unabhängigen linearen Gleichungen in n Unbekannten: Ax = b mit der sogenannten Steifheitsmatrix A (also wieder LGS wie bei Finiten Differenzen)
Page 139 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Lösung des resultierenden LGS • Eigenschaften von A (und unseres LGS): – insb. bei Ritz-Galerkin Vn = Wn ist A oft symmetrisch positiv definit (SPD) – Traumkonstellation: A diagonal, d.h. alle Ansatz- bzw. Testfunktionen sind (bi-) orthogonal (selten, und wenn, dann schwer erreichbar) ai,j = a(ϕi , ϕj ) = Lϕi · ϕj dΩ = δi,j Ω
– immerhin: Ansatz- und Testfunktionen haben i.A. nur lokalen Träger – deshalb ist A typischerweise dünn besetzt
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– deshalb und aufgrund der oft hohen Zahl von Unbekannten werden iterative Löser benötigt • Strategie folglich: – wähle Ansatz- und Testräume mit guten Approximationseigenschaften – konstruiere für diese Räume Basen, die in „schönen“ Matrizen und damit in schnell zu lösenden LGS resultieren – Wesentlich besser als Stützpunktbasen (gleiche Basisfunktion in jedem Gitterpunkt) sind hier hierarchische Basen, wie sie in der Vorlesung „Hierarchie und Rekursion in numerischen Algorithmen“ betrachtet werden.
Page 140 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Standard-Iterationsverfahren für LGS • iterative Lösung großer (dünn besetzter) LGS: – eine der zentralen numerischen Aufgabenstellungen in der numerischen Simulation – treten bei der Diskretisierung von ODE (RWP) und PDE auf • direkte Löser oft nicht wettbewerbsfähig: – Zahl der Unbekannten zu groß (vgl. PDE in 3D)
Populationsdynamik: . . .
– klassische Elimination zerstört die (dünne) Struktur der Matrix
Regelungstechnik: . . .
– wozu exakte Lösung bei Approximationen? (gilt insb. im nichtlinearen Fall, wo ein LGS in jedem äußeren Iterationsschritt zur Behandlung der Nichtlinearität auftritt)
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Ziel: Performance nach der Art „für 3 Stellen braucht man 10 Schritte“ – unabhängig von der Zahl der Unbekannten • typisch jedoch für die klassischen Iterationsverfahren: – Konvergenzgeschwindigkeit verschlechtert sich mit wachsender Problemgröße! Page 141 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Grundlegendes zu Iterationsverfahren • betrachte Iterationsverfahren, starte bei x(0) ∈ IRn und ende (hoffentlich) nahe bei Lösung x von Ax = b: x(0) → x(1) → . . . → x(i+1) → . . . → lim x(i) = x i→∞
• Konvergenzgeschwindigkeit: x − x(i+1) < γ · x − x(i) s für ein 0 < γ < 1 , Konvergenzordnung s
Regelungstechnik: . . .
• typisches Verhalten einfacher Iterationsverfahren für LGS: s = 1,
γ = O(1 − n−k ),
Populationsdynamik: . . .
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
k ∈ {0, 1, 2, ...}
(n : Anzahl der Gitterpunkte) • Strategie: suche Verfahren mit – nur O(n) arithmetischen Operationen pro Iterationsschritt (Kosten; so viel Aufwand wohl mindestens nötig) – Konvergenzverhalten gemäß γ < 1 − const. (Nutzen) • zwei große Familien: Relaxations- und Krylov-Raum-Verfahren
Page 142 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Relaxationsverfahren • manchmal auch Glätter genannt: – Richardson-Iteration – Jacobi-Iteration – Gauß-Seidel-Iteration – Überrelaxation (SOR) oder gedämpfte Methoden • Ansatzpunkte:
Populationsdynamik: . . .
– Fehler e der momentanen Näherung (Ursache): unbekannt
Regelungstechnik: . . .
– Residuum r der momentanen Näherung (Wirkung): ermittelbar
Verkehrsfluss: . . .
r
(i)
(i)
= b − Ax
(i)
= Ax − Ax
(i)
Wärmeleitung: Modell . . .
(i)
= A(x − x ) = −Ae
(mangels Besserem als Fehlerindikator verwendet) • wie Residuum für verbesserte Approximation verwerten? – Richardson: Residuum direkt als Korrektur – Jacobi/Gauß-Seidel: eine Komponente von r zu Null machen – SOR/gedämpft: wie zuvor, aber Korrektur etwas zu hoch / niedrig
Page 143 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Wichtige Relaxationsverfahren • Richardson-Iteration: for i = 0,1,... for k = 1,...,n:
(i+1)
xk
(i)
(i)
:= xk + rk
Hier wird einfach das Residuum r(i) komponentenweise als Korrektur für die aktuelle Näherung x(i) herangezogen. • Jacobi-Iteration: for i = 0,1,... for k = 1,...,n: for k = 1,...,n:
Populationsdynamik: . . .
yk := (i+1)
xk
1 akk
:=
(i) · rk (i) xk +
Regelungstechnik: . . .
yk
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– In jedem Teilschritt k eines Schritts i wird eine Korrektur yk berechnet und gespeichert. – Sofort angewendet, würde diese zum (momentanen) Verschwinden der k-Komponente des Residuums r(i) führen (leicht durch Einsetzen zu verifizieren). – Gleichung k wäre mit dieser aktuellen Näherung für x somit exakt gelöst – ein Fortschritt, der im folgenden Teilschritt zu Gleichung k + 1 natürlich gleich wieder verloren ginge. – Allerdings werden diese Komponentenkorrekturen nicht sofort, sondern erst am Ende eines Iterationsschritts durchgeführt (zweite k-Schleife).
Page 144 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Wichtige Relaxationsverfahren (2) • Gauß-Seidel-Iteration: for i = 0,1,... for k = 1,...,n:
(i)
rk
k−1 (i+1) n (i) := bk − j=1 akj xj − j=k akj xj
yk :=
1 akk
(i)
(i+1)
· rk ,
xk
(i)
:= xk + yk
– Hier wird also dieselbe Korrektur wie beim Jacobi-Verfahren berechnet, der Update wird jetzt allerdings immer sofort und nicht erst am Ende des Iterationsschritts vollzogen.
Populationsdynamik: . . . Regelungstechnik: . . .
– Damit liegen beim Update von Komponente k für die Komponenten 1 bis k − 1 bereits die modifizierten neuen Werte vor.
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Manchmal führt in jeder der drei skizzierten Methoden ein Dämpfen (Multiplikation der Korrektur mit einem Faktor 0 < α < 1) bzw. eine Überrelaxation (Faktor 1 < α < 2) zu einem besseren Konvergenzverhalten: (i+1)
xk
(i)
:= xk + αyk .
– Im Gauß-Seidel-Fall ist vor allem die Version mit α > 1 gebräuchlich, man spricht hier von SOR-Verfahren (Successive Over Relaxation). Page 145 of 165
– Im Jacobi-Fall wird dagegen meistens gedämpft.
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskussion: Additive Zerlegung der Systemmatrix • Für eine kurze Konvergenzanalyse der obigen Verfahren benötigen wir eine algebraische Formulierung (anstelle der algorithmischen). • Alle gezeigten Ansätze basieren auf der einfachen Idee, die Matrix A als Summe A = M + (A − M ) zu schreiben, wobei M x = b sehr einfach zu lösen und der Unterschied A − M bzgl. einer Matrixnorm nicht zu groß sein sollte. • Mit Hilfe eines solchen geeigneten M werden sich Richardson-, Jacobi-, Gauß-Seidel- und SOR-Verfahren schreiben lassen als (i+1)
(i)
+ (A − M )x
Mx
Populationsdynamik: . . . Regelungstechnik: . . .
= b
Verkehrsfluss: . . .
bzw., nach x(i+1) aufgelöst, (i+1)
x
:= M
−1
b−M
−1
Wärmeleitung: Modell . . . (i)
(A−M )x
= M
−1
b−(M
−1
(i)
A−I)x
(i)
= x +M
−1 (i)
r
.
• Darüber hinaus zerlegen wir A additiv in seinen Diagonalteil DA , seinen strikten unteren Dreiecksteil LA sowie seinen strikten oberen Dreiecksteil UA : A =: LA + DA + UA . Damit können wir die folgenden Beziehungen zeigen: – – – –
Richardson: M Jacobi: M Gauß-Seidel: M SOR: M
:= I , := DA , := DA + LA , := α1 DA + LA .
Page 146 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskussion: Additive Zerlegung der Systemmatrix (2) • Bei Betrachtung der algorithmischen Formulierungen für das Richardsonsowie das Jacobi-Verfahren erweisen sich die ersten beiden Gleichungen als offensichtlich: – Bei Richardson wird das Residuum direkt als Korrektur genommen, als Vorfaktor ergibt sich somit die Identität I. – Bei Jacobi wird das Residuum durch das Diagonalelement dividiert, als Vorfaktor ergibt sich somit die Inverse des Diagonalanteils DA . • Weil die Gauß-Seidel-Iteration ein Spezialfall des SOR-Verfahrens ist (α = 1), reicht es aus, obige Formel für M für den allgemeinen SOR-Fall zu zeigen. Aus dem Algorithmus folgt unmittelbar ⎞ ⎛ k−1 n (i+1) (i) (i+1) (i) xk := xk + α ⎝bk − akj xj − akj xj ⎠ /akk
j=1
j=k
−1 b − LA x(i+1) − (DA + UA )x(i) ⇔ x(i+1) := x(i) + αDA
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
1 1 DA x(i+1) = DA x(i) + b − LA x(i+1) − (DA + UA )x(i) α α 1 1 (i+1) ⇔ + (1 − )DA + UA x(i) = b DA + LA x α α ⇔
⇔ M x(i+1) + (A − M )x(i) = b ,
Page 147 of 165 Modellbildung und Simulation
womit die Behauptung für das SOR-Verfahren bewiesen ist.
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskussion: Allgemeines Konvergenzverhalten • Was die Konvergenz angeht, so gibt es zwei unmittelbare Konsequenzen aus dem Ansatz M x(i+1) + (A − M )x(i) = b : – Falls die Folge x(i) konvergiert, dann ist der Grenzwert die exakte Lösung x unseres Systems Ax = b. – Für die Analyse werde angenommen, dass die Iterationsmatrix −M −1 (A− M ) (d.h. die Matrix, die auf e(i) angewandt wird, um e(i+1) zu erhalten; s.u.) symmetrisch sei. Dann ist der Spektralradius ρ (d.h. der betragsgrößte Eigenwert) die für das Konvergenzverhalten entscheidende Größe: ∀x(0) ∈ IRn : lim x(i) = x = A−1 b ⇔ ρ < 1.
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
i→∞
Um das zu sehen, subtrahiere man M x + (A − M )x = b von der Gleichung ganz oben: M e(i+1) + (A − M )e(i) = 0
⇔
e(i+1) = −M −1 (A − M )e(i) .
Wenn alle Eigenwerte betragsmäßig kleiner 1 sind und somit ρ < 1 gilt, werden alle Fehlerkomponenten in jedem Iterationsschritt reduziert. Im Falle ρ > 1 wird sich mindestens eine Fehlerkomponente aufschaukeln. – Ziel bei der Konstruktion iterativer Verfahren muss natürlich ein möglichst kleiner Spektralradius der Iterationsmatrix sein (möglichst nahe bei Null).
Page 148 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Diskussion: Konvergenzaussagen • Es gibt eine Reihe von Resultaten zur Konvergenz der verschiedenen Verfahren, von denen einige bedeutende erwähnt werden sollen: – Notwendig für die Konvergenz des SOR-Verfahrens ist 0 < α < 2. – Falls A positiv definit ist, dann konvergieren sowohl das SOR-Verfahren (für 0 < α < 2) als auch die Gauß-Seidel-Iteration. – Falls A und 2DA − A beide positiv definit sind, dann konvergiert das Jacobi-Verfahren. – Falls A strikt diagonal dominant ist (d. h. aii > j=i | aij | für alle i), dann konvergieren das Jacobi- und das Gauß-Seidel-Verfahren. – In bestimmten Fällen lässt sich der optimale Parameter α bestimmten (ρ minimal, so dass Fehlerreduktion pro Iterationsschritt maximal).
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• Die Gauß-Seidel-Iteration ist nicht generell besser als das Jacobi-Verfahren (wie man aufgrund des sofort vollzogenen Updates vermuten könnte). Es gibt Beispiele, in denen Erstere konvergiert und Letzteres divergiert, und umgekehrt. In vielen Fällen kommt das Gauß-Seidel-Verfahren jedoch mit der Hälfte der Iterationsschritte des Jacobi-Verfahrens aus. Page 149 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Zum Spektralradius typischer Iterationsmatrizen • Offensichtlich ist ρ nicht nur entscheidend für die Frage, ob die Iterationsvorschrift überhaupt konvergiert, sondern auch für deren Qualität, also ihre Konvergenzgeschwindigkeit: Je kleiner ρ ist, desto schneller werden alle Komponenten des Fehlers e(i) in jedem Iterationsschritt reduziert. • In der Praxis haben die obigen Resultate zur Konvergenz leider eher theoretischen Wert, da ρ oft so nahe bei 1 ist, dass – trotz Konvergenz – die Anzahl der erforderlichen Iterationsschritte, bis eine hinreichende Genauigkeit erreicht ist, viel zu groß ist. • Ein wichtiges Beispielszenario ist die Diskretisierung partieller Differentialgleichungen:
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Typisch ist, dass ρ von der Problemgröße n und somit von der Auflösung h des zugrunde liegenden Gitters abhängt, also beispielsweise 1 ρ = O(1 − h2l ) = O 1 − l 4 bei einer Maschenweite hl = 2−l . – Dies ist ein gewaltiger Nachteil: Je feiner und folglich auch genauer unser Gitter ist, umso erbärmlicher wird das Konvergenzverhalten unserer iterativen Verfahren. Bessere iterative Löser (z.B. Mehrgitterverfahren, die aber den Rahmen dieser Vorlesung sprengen würden) sind also ein Muss!
Page 150 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ganz anders: Steilster Abstieg • alternative Sehweise für positiv definites A : x löst Ax = b
⇔
x minimiert f (x) = 0.5 · xT Ax − bT x + c
(Eindeutigkeit des Minimums, da A positiv definit) • somit neue Strategie: – suche nach Minimum der Funktion f – möglicher Weg: Methode des steilsten Abstiegs T
repeat(i) :
αi
r(i) r(i) ; r(i)T Ar(i) = x(i) + αi r(i) ; = r(i) − αi Ar(i) ; =
x(i+1) r(i+1)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– sucht nach maximaler Verbesserung in Richtung des negativen Gradienten) −f (x(i) ) = r(i) – eindimensionale Suche nach Minimum: minαi f (x(i) + αi r(i) )! – vernünftige Vorgehensweise, aber trotzdem stark heuristisch
Page 151 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Steilster Abstieg (2) • Geht’s noch simpler? – 1D-Suche längs der Koordinatenachsen (statt Gradienten) – führt zu Gauß-Seidel-Iteration (man sieht: alles ist verwandt!) • Konvergenzgeschwindigkeit: – beliebig langsam – erzielter Fortschritt kann immer wieder verloren gehen (unabhängiges lokales Herumoptimieren in wechselnde und nicht zusammenhängende Richtungen)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• entscheidende Größe hier: – spektrale Konditionszahl von A (Eigenwerte betragsmäßig) κ(A) =
λmax (A) λmin (A)
– gut: möglichst klein, also begrenzte Breite des Spektrums (auch bei zunehmender Auflösung des Gitters) – ergo: algorithmische Entwicklung muss in Richtung einer Zähmung der Eigenwerte gehen (diesmal der Systemmatrix anstelle der Iterationsmatrix wie bei den Relaxationsverfahren)
Page 152 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Verbesserung: Konjugierte Richtungen • Verbesserung gegenüber steilstem Abstieg: – orthogonale Suchrichtungen, Fehler nach i Schritten soll orthogonal auf allen vorigen Suchrichtungen stehen – nichts zerstört, folglich: im Prinzip direktes Verfahren (nach n Schritten im Optimum) – allerdings sind n Iterationen in der Praxis zu viel, daher Einsatz als iterative Methode (semi-iteratives Verfahren)
Populationsdynamik: . . .
– neue Suchrichtungen: x(i+1) = x(i) + αi d(i)
Regelungstechnik: . . .
T
– optimale Orthogonalität wäre 0 = d(i) e(i+1) , aber Fehler ja unbekannt T
– daher Konjugation: 0 = d(i) Ae(i+1) (u und v A-orthogonal oder konjugiert , falls uT Av = 0 )
Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– Algorithmus: starte mit d(0) = r(0) und iteriere: T
repeat(i) :
αi x(i+1) r(i+1)
d(i) r(i) = (i)T (i) ; d Ad = x(i) + αi d(i) ; = r(i) − αi Ad(i) ; Page 153 of 165
– noch zu tun: Konstruktion der konjugierten Richtungen d(i)
Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schließlich: Konjugierte Gradienten • obiges Verfahren + Konstruktion der konjugierten Richtungen • Konstruktionsprinzip: Gram-Schmidt Konjugation der Residuen • keine Details oder Herleitung, nur der fertige Algorithmus: T
repeat(i) :
αi x(i+1) r(i+1)
d(i) r(i) = (i)T (i) ; d Ad = x(i) + αi d(i) ; = r(i) − αi Ad(i) ; T
βi+1 d(i+1)
r(i+1) r(i) ; r(i)T r(i) = r(i+1) + βi+1 d(i) ; =
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• schneller als steilster Abstieg, aber immer noch n-abhängig! • Suchräume bilden sogenannte Krylov Sequenz: span{d(0) , . . . , d(i−1) } = span{d(0) , Ad(0) , . . . , Ai−1 d(0) } = span{r(0) , Ar(0) , . . . , Ai−1 r(0) } • andere bekannte Krylov-Verfahren: GMRES, Bi-CGSTAB
Page 154 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Schnelle iterative Löser für LGS • entscheidender Nachteil der bisher diskutierten Lösungsverfahren: Konvergenz verlangsamt sich bei zunehmender Gitterfeinheit! • deshalb jetzt: Suche nach möglicher Abhilfe – Relaxation: explizite Anwendung des Mehrgitterprinzips – Krylov/cg: Vorkonditionierung (typischerweise auch mehrgitterartig) • zunächst Vorkonditionierung:
Populationsdynamik: . . .
– entscheidende Größe für die Konvergenz des cg-Verfahrens: Kondition der Matrix
Regelungstechnik: . . .
– PDE: diese Kondition wächst dramatisch mit n
Wärmeleitung: Modell . . .
Verkehrsfluss: . . .
– deshalb: modifiziere Matrix, um Kondition zu verbessern Ax = b ⇔ M −1 Ax = M −1 b ⇔ W −1 AW −T y = W −1 b, M s.p.d., W W T = M, y = W T x, M −1 A und W −1 AW −T
wobei ähnlich
– M bzw. W müssen nicht explizit konstruiert werden, sie müssen nur angewandt werden – Vorteil der W -Schreibweise: Symmetrie von A vererbt sich!
Page 155 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Strategien für Vorkonditionierer • einfachste Wahl: • beste Wahl:
M = I (billig, aber nutzlos)
M = A (perfekt, aber so teuer wie Ax = b )
• einige mögliche Kompromisse: – diagonaler oder Jacobi Vorkonditionierer: M = DA – GS oder SOR werden nicht benutzt (nicht symmetrisch) – SSOR Vorkonditionierer:
Populationsdynamik: . . .
M (1/2) := α−1 DA + LA ; M (1) := α−1 DA + UA ; α −1 M := M (1) (M (1/2) )−1 DA α−2
Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– unvollständige Faktorisierung, z.B. ILU: berechne näherungsweise Faktoren L und U anstelle der exakten bei direkten Lösern – sparse approximate inverse: konstruiere einfaches B mit min I − AB2 , B
M −1 = B
– Multilevel Vorkonditionierer: folgen dem Mehrgitterprinzip (folgt)
Page 156 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Das Mehrgitterprinzip • starting point: Fourier mode analysis of the errors – decompose the error e(i) = x(i) − x into its Fourier components (Fourier transform) – observe how they change/decrease under a standard relaxation like Jacobi or Gauß-Seidel (in a two-band sense): * The high frequency part (with respect to the underlying grid) is reduced quite quickly. * The low frequency part (w.r.t. the grid) decreases only very slowly; actually the slower, the finer the grid is. – This behaviour is annoying
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
* the low frequencies are not expected to make troubles, but we can hardly get rid of them on a fine grid; – but also encouraging * the low frequencies can be represented and, hopefully, tackled on a coarser grid – there is no need for the fine resolution.
Page 157 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein einfaches Beispiel • 1D Laplace equation, u(0) = u(1) = 0 (exact solution 0) • equidistant grid, 65 points, 3-point stencil, damped Jacobi method with damping parameter 0.5 • start with random values in [0, 1] for u in the grid points • After 100 (!) steps, there is still a maximum error bigger than 0.1 due to low-frequency components! • therefore the name smoothers for relaxation schemes: – They reduce the strongly oscillating parts of the error quite efficiently. – They, thus, produce a smooth error which is very resistent.
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• the idea: work on grids of different resolution and combine the effects in an appropriate way
Page 158 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Ein einfaches Beispiel (2)
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
Page 159 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Grobgitterkorrektur • sequence of equidistant grids on our domain: Ωl , l = 1, 2, . . . , L, with mesh width hl = 2−l • let Al , bl denote corresponding matrix, right-hand side, ... • combine work on two grids with a correction scheme:
smooth the current solution xl ; form the residual rl = bl − Al xl ; restrict rl to the coarse grid Ωl−1 ;
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
provide a solution to Al−1 el−1 = rl−1 ; prolongate el−1 to the fine grid Ωl ; add the resulting correction to xl ; if necessary, smooth again ;
Page 160 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Grobgitterkorrektur (2) • the different steps of this 2-grid algorithm : – the pre-smoothing: reduce high-frequency error components, smooth error, and prepare residual for transfer to coarse grid – the restriction: transfer from fine grid to coarse grid * injection : inherit the coarse grid values and forget the others * (full) weighting : apply some averaging process – the coarse grid correction: provide an (approximate) solution on the coarse grid (direct, if coarse enough; some smoothing steps otherwise)
Populationsdynamik: . . .
– the prolongation: transfer from coarse grid to fine grid
Wärmeleitung: Modell . . .
Regelungstechnik: . . . Verkehrsfluss: . . .
* usually some interpolation method
– the post-smoothing: sometimes reasonable to avoid new high-frequency error components
Page 161 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Der V-Zyklus • recursive application of 2-grid scheme leads to multigrid methods • there, the coarse grid equation is solved by coarse grid correction, too; the resulting algorithmic scheme is called V-cycle :
smooth the current solution xl ; form the residual rl = bl − Al xl ; restrict rl to the coarse grid Ωl−1 ; solve Al−1 el−1 = rl−1 by coarse grid correction ;
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
prolongate el−1 to the fine grid Ωl ; add the resulting correction to xl ; if necessary, smooth again ;
• on the coarsest grid: direct solution • number of smoothing steps: typically small (1 or 2)
Page 162 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Weitere Mehrgitterzyklen • the V-cycle is not the only multigrid scheme: – the W-cycle: after each prolongation, visit the coarse grid once more, before moving on to the next finer grid
Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
– W-cycle is sometimes advantageous with respect to speed of convergence
Page 163 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Nested Iteration und Full Multigrid • two more famous multigrid schemes: – the nested iteration: start on coarsest grid Ω1 , smooth, prolongate to Ω2 , smooth, prolongate to Ω3 , and so on, until finest grid is reached; now start V-cycle – full multigrid: replace ‘smooth´ steps above by ‘apply a V-cycle´; combination of improved start solution and multigrid solver Populationsdynamik: . . . Regelungstechnik: . . . Verkehrsfluss: . . . Wärmeleitung: Modell . . .
• multigrid idea is not limited to rectangular or structured grids: we just need a hierarchy of nested grids (works for triangles or tetrahedra, too) • also without underlying geometry: algebraic multigrid methods Page 164 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation Hans-Joachim Bungartz
Grundlegende Konvergenzresultate • Cost (storage and computing time): – 1D :
c · n + c · n/2 + c · n/4 + c · n/8 + . . . ≤ 2c · n = O(n)
– 2D :
c · n + c · n/4 + c · n/16 + c · n/64 + . . . ≤ 4/3c · n = O(n)
– 3D :
c · n + c · n/8 + c · n/64 + c · n/512 + . . . ≤ 8/7c · n = O(n)
– i.e.: work on coarse grids is negligible compared to finest grid • Benefit (speed of convergence):
Populationsdynamik: . . .
– always significant acceleration compared with pure use of smoother (relaxation method)
Regelungstechnik: . . .
– in most cases even ideal behaviour γ = O(1 − const.)
Wärmeleitung: Modell . . .
Verkehrsfluss: . . .
– effect: * constant number of multigrid steps to obtain a given number of digits * overall computational work increases only linearly with n
Page 165 of 165 Modellbildung und Simulation
3.Kontinuierliche Modellbildung und Simulation
Hans-Joachim Bungartz
