Numerik von Dierentialgleichungen

Einführung und Überblick

Winfried Auzinger

Institut für Analysis and Scientic Computing

www.asc.tuwien.ac.at/∼winfried

`Anwendungsgebiete der Mathematik', SS 2017

1 / 35

Einleitung

In der VO `Anwendungsgebiete der Technischen Mathematik' wirdan vielen Stellen auf mathematische Modelle Bezug genommen, dieauf gewöhnliche oder partielle Dierentialgleichungen (DGL)führen.

Zugrunde liegt immer ein Modellierungsprozess; dieser führt(auÿer in Trivialfällen) nicht direkt zur Problemlösung, sondern`nur' auf eine zu lösende Gleichung.

Wenn z.B. die unabhängige Variable der Zeit t entspricht, soerhält man eine Dierentialgleichung immer dann, wenn die(zeitliche) Änderung eines Zustandes vom Zustand selbst abhängt,so dass Ableitungen (Geschwindigkeit, Beschleunigung, . . . ) in denGleichungen auftreten(im Gegensatz zu algebraischen Gleichungen).

2 / 35

Einleitung

Die mathematische Modellierung dient der Simulation eines realenVorganges auf einem Computer z.B. um Vorhersagen über dasVerhalten eines Systems zu gewinnen (wann ist die nächsteSonnennsternis?) oder aufwendige bzw. undurchführbareExperimente (1000 verschiedene Crashtests mit einemTestfahrzeug ?!) virtuell durchführen zu können.

Die Entwicklung auf dem Informatiksektor erlaubt immerkomplexere Simulationen und die Numerik vonDierentialgleichungen liefert eine der zentralen mathematischenGrundlagen, die diese erst ermöglichen.

In den Medien liest und hört man meist nur von Prozessoren undTaktfrequenzen und coolen Graken etc. die Mathematikdahinter wird kaum wahrgenommen.

3 / 35

EinleitungNotation

Einfachste Standardform einer gewöhnlichen DGL (ODE):

y ′(t) = f (t, y(t))

für eine gesuchte Funktion y(t) (y ′ = dy/dt).

• t . . . unabhängige Variable (in Anwendungen oft die Zeit,oder z.B. räumliche Koordinate)

• f (t, y) . . . gegebene Datenfunktion (ergibt sich aus zugrundeliegenden Modell)

Gleiche Schreibweise für Systeme von DGL:

• mehrere DGL für mehrere gesuchte Funktioneny1(t), . . . , yn(t) simultan gegeben

• Gleichungen untereinander verkoppelt

• Hier: y = (y1, . . . , yn) und f = (f1, . . . , fn) . . . vektoriell

4 / 35

Beispiel: PendelDas mathematische Modell

Mathematisches Modell für

• Pendelschwingung unter Einwirkung der Gravitation

• `Punktmasse' am Ende konzentriert (Modellvereinfachung!)

• Grundlage: 2. Newtonsches Gesetz( Kraft = Masse × Beschleunigung )

• Variablen: Zeit t (z.B. in s), Auslenkungswinkel ϕ(t)

• DGL. 2. Ordnung (Notation: ϕ = Ableitung nach t)

ϕ(t) = − ω2 sin(ϕ(t))

ω :=√g/l (g . . . Erdbeschleunigung, l . . . Pendellänge)

5 / 35

Beispiel: PendelDie Dierentialgleichung

• DGL ist nichtlinear in ϕ(t), Ordnung=2 (ϕ tritt auf)

• Konkrete Lösung durch zusätzliche Vorgabe eines konkretenAnfangszustandes (Auslenkung ϕ(0), Winkelgeschwindigkeitϕ(0) zum Zeitpunkt t = 0)

• Oft wird die Ableitung ψ(t) := ϕ(t) als eigenständigeZustandsvariable betrachtet. äquivalentes DGL-System 1.Ordnung:

ϕ(t) = ψ(t)

ψ(t) = −ω2 sin(ϕ(t)) ,

d.h. y = (ϕ,ψ) und f (t, y) = f (y) = (ψ,−ω2 sin(ϕ)).

• Exakte Lösung kann nicht explizit formelmäÿig angegebenwerden.

6 / 35

Beispiel: PendelAnmerkungen

Viele vereinfachende Annahmen:

• Gesamte Masse im Endpunkt konzentriert

• Räumliche Ausdehnung klein, so dass Annahme g = const.gerechtfertigtm

• Starrer Körper (keine Verformungen)

• Reibungsfreie Aufhängung; kein Luftwiderstand

• . . . . . .

• Derartige Modellfehler hat man so gut wie immer;

eine exakte formelmäÿige Lösung wäre daher ohnehin nicht`wirklich exakt'.

• Numerische Approximation erforderlich in der Praxis der Normalfall.

7 / 35

Beispiel: PendelModellvereinfachung (Linearisierung)

Für hinreichend kleine Auslenkungen ϕ(t):• sinϕ ≈ ϕ

• vereinfachte, linearisierte DGL

ϕ(t) = − ω2 ϕ(t)

mit exakter Lösung

ϕ(t) = C1 sin(ωt) + C2 cos(ωt)

• . . . keine numerische Approximation erforderlich.• Beachte jedoch: Auch sin, cos müssen am Computerapproximiert werden!

• Grenze zwischen `exakter' und `numerischer' Lösung ist nichtganz scharf deniert.

• Exakte Lösungsdarstellungen (sofern verfügbar) dienen eherzum theoretischen Verständnis; sie eignen sich auch oft nichtfür die numerische Umsetzung.

8 / 35

Ein `Urproblem': (Numerische) IntegrationWas heiÿt `integrierbar' ?

• Einfacher Spezialfall einer DGL:

Löse y ′(t) = f (t) ( f (t, y) ≡ f (t) )

• Lösung zu Anfangswert y(0) = y0 :

y(t) = y0 +

∫ t

0

f (τ) dτ . . . bestimmtes Integral von f

• Konkrete Darstellung mittels Stammfunktion F (t) von f (t)(d.h. F ′(t) ≡ f (t)) :

y(t) = y0 + (F (t)− F (0)))

. . . setzt voraus, dass F (t) explizit bekannt ist und Implementierung/Approximation von F am Rechnerverfügbar ist (z.B. sin, cos, . . .)

• Jede (z.B.) stetige Funktion f (t) ist integrierbar, aber dieStammfunktion ist i.all.g. nicht explizit angebbar.

9 / 35

Ein `Urproblem': (Numerische) IntegrationNumerische Integration

• Integration erfordert numerische Approximation für

(i) das Integral∫ t

0

f (τ) dτ oder

(ii) die Stammfunktion F (t).

• (i) ist die allgemeinere Strategie: Man verwendetApproximationsmethoden, die sich (im Prinzip) für beliebigeIntegranden einsetzen lassen.

• Approximationsprinzip: Diskretisierung, z.B.∫ t

0

f (τ) dτ ≈∑i

f (ti ) ∆ti

über einem geeignet gewählten Gitter ti.• Die nachfolgende betrachteten Verfahren für DGL basierendauf der Idee der Diskretisierung angewendet auf die DGL (odereine äquivalente Integralgleichung).

10 / 35

Ein `Urproblem': (Numerische) IntegrationLineare Dierentialgleichungen

• Die Lösung y(t) einer linearen DGL, z.B.

y ′(t) = a y(t) + b(t)

kann ebenfalls direkt als Integral angegeben werden:

y(t) = etay0 +

∫ t

0

e(t−τ)ab(τ) dτ

• Dennoch ist direkte Diskretisierung der DGL (siehe unten)meist sinnvoller.

• Für Systeme linearer DGL, z.B. (y(t) ∈ Rn)

y ′(t) = A · y(t) + b(t)

gilt analog

y(t) = etA · y0 +

∫ t

0

e(t−τ)A · b(τ) dτ

Auswertung der Matrix-Exponentialfunktion etA (Lösungdes homogenen Problems) ist i. allg. nichttrivial.

11 / 35

Numerischer Zugang: DiskretisierungEuler-Verfahren

• Basierend auf der Approximation

y ′(t) ≈ y(t + h)− y(t)

h,

ersetzt man die DGL über einem `Gitter' ti = i h durch dieDierenzengleichung mit Schrittweite h,

yi − yi−1h

= f (ti−1, yi−1), i = 1, 2, . . .

• Rekursion für die Näherungswerte yi ≈ y(ti ) :

yi = yi−1 + h f (ti−1, yi−1), i = 1, 2, . . .

• Approximationsfehler (`lokaler Diskretisierungsfehler') in jedemeinzelnen Schritt:

y(t + h)− y(t)

h− y ′(t) = O(h) für h→ 0

12 / 35

Numerischer Zugang: DiskretisierungFehlerbetrachtung

• Lokaler Diskretisierungsfehler: Genauere Darstellung(a priori - Schätzung) :

y(t + h)− y(t)

h− y ′(t) = h y ′′(t) +O(h 2) für h→ 0

• ` A priori ' : Zeigt qualitativ, wie lokaler Fehler von h und der(unbekannten) Lösung y(t) abhängt.

Dies ist jedoch (klarerweise) keine berechenbare Gröÿe.

• A priori Schätzungen dienen dem theoretischen Verständnisund bilden die theoretische Basis für die Konstruktionpraktikabler a posteriori Schätzer.

• Der `wirkliche' (globale) Fehler ist

yi − y(ti ) = ?, i = 1, 2, . . .

13 / 35

Numerischer Zugang: DiskretisierungStabilität und Konvergenz

• Die Frage nach der Stabilität eines Verfahrens ist omnipräsentin der Numerik:

Diverse unvermeidliche Fehler (lokale Fehler, auch:Modellungenauigkeiten, Approximationsfehler bei Auswertungder Datenfunktionen, . . . ) sollen sich nicht unnatürlich starkauf auf das Endergebnis auswirken.

• Satz: Das explizite Euler-Verfahren ist stabil.

(Eigentlich: stabil für nicht allzu groÿe Schrittweiten h.)

Beweis: Studium der Dierenzengleichung für den globalenFehler und a priori Abschätzung seiner Lösung.

• Folgerung: Für den globalen Fehler yi − y(ti ) gilt

yi − y(ti ) = O(h) für h→ 0

d.h. das Verfahren ist konvergent von der Ordnung 1.14 / 35

Numerischer Zugang: DiskretisierungFAQ, (i): Relevante Fragen zu Diskretisierungsverfahren

• Welches Diskretisierungsverfahren eignet sich für welche Klassevon Dierentialgleichungen?

Das `Super-Universalverfahren' gibt es nicht.

• ∃ beliebig hoher Ordnung p > 1 . Fehlerverhalten:

yi − y(ti ) = O(hp) .

Bei einem Verfahren z.B. 4. Ordnung (p = 4, z.B. klassisches`Runge-Kutta'-Verfahren) geht der Fehler auf das (ca.)2−p = 2−4 - fache, also um einen Faktor 1/16 zurück, wenn dieSchrittweite h halbiert wird.

• höherer Rechenaufwand pro Schritt lohnt sich.

• Aber: Wird inezienter für `schwierigere' Probleme mitGlattheitsdefekten.

• Weitere Problematik: Abklingverhalten mit stark gedämpfterDynamik kann Instabilität verursachen.

15 / 35

Numerischer Zugang: DiskretisierungBeispiel: Flug einer Rakete im inhomogenen Gravitationsfeld

• Dierentialgleichung für die Flughöhe h = h(t) :

h′′(t) = − g r2

(r + h(t))2,

g = Erdbeschleunigungr = Erdradius

• Grak zeigt exakte Lösung, Euler-Approximation undRunge-Kutta Approximation angewendet auf äquivalentesSystem 1.Ordnung, mit Geschwindigkeit h′(t) = v(t) alseigener Variable.

16 / 35

Numerischer Zugang: DiskretisierungBeispiel: Gedämpfter Oszillator (Pendel, kleine Auslenkung, mit Reibung)

• Dierentialgleichung für den Winkel ϕ(t) :

ϕ(t) + 2ρ ϕ(t) + ω2ϕ(t) = 0,ρ = Dämpfungsparameterω = Eigenfrequenz

• Grak zeigt exakte Lösung, Euler-Approximation undRunge-Kutta Approximation, angewendet auf äquivalentesSystem 1.Ordnung, für den Fall ρ = ω 0.

• Beide Verfahren verhalten sich instabil, obwohl exakte Lösungsehr stabil (glatter Verlauf, leicht abklingend).

17 / 35

Numerischer Zugang: DiskretisierungFAQ, (ii): Relevante Fragen zu Diskretisierungsverfahren

• . . . Beispiel für steifes Verhalten:

Originalproblem ist `super-stabil', Störungen werden sehrschnell weggedämpft

Aber: Verfahren reagieren empndlich auf Störungen (fallsnicht h extrem klein)

• Abhilfe: Z.B. modiziertes, `implizites' Euler-Verfahren:

yi − yi−1h

= f (ti , yi ), i = 1, 2, . . .

Dierenzengleichung ist implizit, d.h. die Bestimmung von yifür gegebenes yi−1 erfordert die Lösung einer Gleichung (i. allg.eines Systems):

Löse z − h f (ti , z) = yi−1 ⇒ yi = z .

• Grak zeigt

18 / 35

Numerischer Zugang: DiskretisierungFAQ, (iii): Relevante Fragen zu Diskretisierungsverfahren

• Wie kann man sicherstellen, dass wichtige probleminhärenteEigenschaften auch in der numerischen Lösung möglichstgetreu wiedergegeben werden?

• Beispiel Reibungsfrei schwingendes Pendel:

Beispiel für ein konservatives mechanisches System, dessengesamte (potentielle + kinetische) Energie über die Zeithinweg konstant ist.

• Zur verlässlichen numerischen Integration über längereZeiten: Benötige Verfahren, die ebenfalls energieerhaltendsind.

• Die Euler-artigen Verfahren sind dafür ungeeignet. Es wurdenspezielle Verfahren für diesen Zweck entwickelt.

• Anwendung: Z.B. in der Astronomie, exakte langfristigeVorhersage von Planetenbahnen.

19 / 35

Numerischer Zugang: DiskretisierungBeispiel: Die äuÿeren Planeten im Sonnensystem

• Beispiel: c© E.Hairer, C. Lubich, E.Wanner

• Grak zeigt verschiedene (qualitativ falsche und richtige)Simulationen der Planetenbahnen (Langzeitverhalten).

(J=Jupiter, S=Saturn, U=Uranus, N=Neptun, P=Pluto)20 / 35

Numerischer Zugang: DiskretisierungFAQ, (iv): Relevante Fragen zu Diskretisierungsverfahren

• Wie kann man den Fehler a posteriori, d.h. für eine konkretberechnete Näherungslösung verlässlich berechenbarschätzen ?

• Eine Möglichkeit (Grundidee Rückwärtsanalyse):

Interpretiere die berechnete Folge von Näherungswerten yi(bzw. eine die yi interpolierende Funktion y(t)) als Lösungeiner `gestörten' DGL.

Der Störungsterm lässt sich mit Hilfe der yi und der gegebenenDGL berechnen bzw. genau schätzen.

• ∃ verschiedene Techniken, um aus dieser Störung auf ihrenEekt (= globaler Fehler!) zurückzuschlieÿen.

Allerdings auch hier nicht 100% Sicherheit auchFehlerschätzer operiert diskret und unterliegt einem (allerdingsgeringeren) Diskretisierungsfehler.

21 / 35

Numerischer Zugang: DiskretisierungFAQ, (v): Relevante Fragen zu Diskretisierungsverfahren

• Wie ist die (im allgemeinen variable) Schrittweite h (= hi ) zuwählen, damit man eine vorgegebene Genauigkeitsforderungmöglicht verlässlich und ohne unnötig hohen Rechenaufwanderfüllt?

• Die Grundlage dafür besteht in einer verlässlichena-posteriori Fehlerschätzung (siehe oben).

Darauf basierend wurden verschiedene Algorithmen entwickelt,bei denen versucht wird, die Schrittweite so zu steuern, dassder Fehler möglichst gleichmäÿig kontrolliert wird.

adaptive Methoden

22 / 35

RandwertproblemeProblemtyp; Beispiel

• Randbedingungen (an mehreren Stellen) ergeben globalesProblem

• Beispiel: Poisson-Gleichung (1D, linear)

u′′(x) = f (x), x ∈ (a, b)

+ Randbedingungen am Rand des Intervalles [a, b],z.B. u(a) = α, u(b) = β

• Anwendung: Z.B. in Mechanik Lösung u(x) . . . Verformung einer gespannten, elastischen

Saite, wobei

f (x) = gegebene Dichtefunktion, beschreibt `load'

• Allgemeiner: Nichtlineares Problem(z.B. nichtlineares Spannungs-Verformungs-Verhalten)

u′′(x) = f (x , u(x), u′(x)), x ∈ (a, b)

23 / 35

RandwertproblemeNumerische Behandlung: Z.B. Dierenzenverfahren für u′′(x) = f (x , u(x), u′(x))

• Approximiere Ableitung über diskretem Gitter xi,Schrittweite h

• Z.B. mit xi = a + i h :

u′′(xi ) ≈u(xi−1)− 2u(xi ) + u(xi+1)

h 2

• . . . Ergibt groÿes (algebraisches) Gleichungssystem für dieui ≈ u(xi ):

ui−1 − 2ui + ui+1

h 2= f (xi , ui ,

ui−ui−1h ), i = 1, 2, . . .

+ Randbedingungen

System Au = F (u), u = (ui ) ∈ Rn; A ∈ Rn×n

• Allgemeiner: Schrittweite h = hi variabel → Adaptivität

24 / 35

RandwertproblemeGroÿes, schwach besetztes Gleichungssystem

• Systemmatrix A (entspricht Approximation von u′′):

A =1

h 2

−2 1

1 −2 1

. . .. . .

. . .

. . .. . .

. . .

1 −2 1

1 −2

• Gleichungssystem ist `schwach besetzt' (tridiagonal),leicht lösbar falls n nicht extrem groÿ.

• Beachte: n→∞ für h→ 0. Benötige Stabilität für numerische Umsetzbarkeit am Digitalrechner Konvergenz der Approximation für h→ 0

• Stabilität: Untersuchung der Inversen A−1 für h→ 0 das ist Lineare Algebra.

25 / 35

Partielle Dierentialgleichungen (PDE)Stationäre Probleme

• Beispiel: Poisson-Gleichung (2D, linear)

∆u(x , y) = f (x , y), (x , y) ∈ Ω ⊂ R2

(∆u := (∇·∇) u = uxx + uyy . . . `Laplace-Operator' ),

+ Randbedingungen am Rand des Gebietes Ω

• Ein Anwendungsbeispiel: Lösung u(x , y) . . . Verformung einer gespannten, elastischen

Membran, wobei

f (x , y) = gegebene Dichtefunktion, beschreibt `load'

• Poisson Gleichung tritt auch auf in Strömungsmechanik,Elektrostatik, . . .

• Analog in 3D

26 / 35

Partielle Dierentialgleichungen (PDE)Stationäre Probleme: FD-Verfahren (i)

• Approximiere die PDE durch eine partielle Dierenzen-gleichung über einem diskreten Gitter xi , yj ⊂ Ω

r r r r r r rr r r r r r rr r r r r r rr r r r r r rr r r r r r rr r r r r r rr r r r r r r

b b b b b b b b bbbbbbbbbb bbbbbbbbb bbbbbbbbbh

h

• . . . führt auf groÿes lineares Gleichungssystemfür die ui ,j ≈ u(xi , yj); `sehr groÿ' im 3D-Fall

27 / 35

Partielle Dierentialgleichungen (PDE)Stationäre Probleme: FD-Verfahren (ii)

System-Matrix mit spezieller Struktur:

A =

T I

I T I

. . .. . .

. . .

. . .. . .

. . .

I T I

I T

∈ Rn2×n2

mit

T =1

h 2

−4 1

1 −4 1

. . .. . .

. . .

. . .. . .

. . .

1 −4 1

1 −4

∈ Rn×n

28 / 35

Partielle Dierentialgleichungen (PDE)Stationäre Probleme: Finite-Elemente -Methoden (FEM, (i)

• Approximiere die Lösung in einem endlich-dimensionalenRaum von approximierenden Funktionen (z.b. stückweiselinear, polynomial über kleinen Teilgebieten, den `FinitenElementen'), wobei

• verschiedene Kriterien für die Approximation, z.B. Minimierungeines gewissen Fehlerfunktionals

( Variationsrechnung)

• Flexibel anpassbar an allgemeine Geometrien

• groÿe, schwach besetzte Gleichungssysteme

• Adaptivität

• Wesentlich: Geeignete Datenstrukturen für Umsetzung amDigitalrechner

29 / 35

Partielle Dierentialgleichungen (PDE)Stationäre Probleme: Finite-Elemente -Methoden (FEM, (ii)

• Beispiel: Poisson-Gleichung mit `unglatter' Lösung

• Welche Mathematik brauchen wir dafür:Höhere Analysis (insbesondere Variationsrechnung),Funktionalanalysis, Approximationstheorie, Lineare Algebra . . .

30 / 35

Partielle Dierentialgleichungen (PDE)Instationäre Probleme Problemtypen

• Wärmeleitungsgleichung (Diusionsgleichung):

∂∂t u(x , t) = ∆x u(x , t)

• Wellengleichung:

∂2

∂t2u(x , t) = ∆xu(x , t)

• + Verallgemeinerungen

• Diskretisierung in Ort x ∈ Rd (d = 1, 2, 3) und t ∈ R(simultan oder `nacheinander')

• Je nach Typ verschiedene Zugänge erforderlich

• Weiters: Nichtlineare Eekte, z.B. Schockwellen (Überschall) Systeme von PDEs,

z.B. Navier-Stokes, Maxwell-Gleichungen . . .

31 / 35

Partielle Dierentialgleichungen (PDE)Instationäre Probleme Beispiel: Schrödinger-Gleichung (i)

• Grundgleichung der Quantenmechanik (n Elektronen):

i∂∂tψ(x , t) = H(ψ(x , t))

für die komplexwertige Wellenfunktion ψ. Dabei ist:

t . . . Zeit x = (x1, . . . , xn), xj ∈ R3 . . . Position der Elektronen

|ψ(x , t) | beschreibt die Wahrscheinlichkeits-dichte dafür, dass sich zum Zeitpunkt t die Teilchen(Elektronen) an den betreenden Positionen xj benden.

• Numerische Behandlung durch `direkte Diskretisierung'honungslos, weil extrem hochdimensional: x ∈ R3n !

• `State of the art': Spezielle approximierende Ansatzfunktionen,dann erst Diskretisierung.

32 / 35

Partielle Dierentialgleichungen (PDE)Instationäre Probleme Beispiel: Schrödinger-Gleichung (ii)

• Eine Lösung der homogenen Schrödinger-Gleichung( ein `free particle', 1D, Reψ(x , t) ):

• Zeitintegration: aufwendig; für die Integrationsschritte sindeziente Approximationsverfahren erforderlich.

33 / 35

Partielle Dierentialgleichungen (PDE)Lösung groÿer Gleichungssysteme

• Fast überall treten sehr groÿe Gleichungssysteme auf(nach Linearisierung): Au = b

• Direkte Lösung oft langsam und/oder speicherintensiv• Alternative: Iterative Verfahren

• Z.B. iterative Minimierung einer Residualnorm ‖Ax − b‖• Beispiel: cg-Verfahren für positiv denite Systeme:

# cg for A.x=b, starting from x[0]

r[0] := b-A.x[0]; d[0] := r[0];

for i=0,1,... do

alpha[i] := (r[i],r[i])/(A.d[i],d[i])

x[i+1] := x[i]+alpha[i]*d[i]

r[i+1] := r[i]-alpha[i]*A.d[i]

beta[i] := (r[i+1],r[i+1])/(r[i],r[i])

d[i+1] := r[i+1] + beta[i]*d[i]

end

• Aber: Muss problemabängig angepasst werden fürvernünftige Konvergenzrate.

34 / 35

Denition der Numerik?`Numerische Mathematik'; `Scientic Computing'

• . . . nicht nur Rundungsfehler und so . . .

• Die Numerik ist der konstruktive Zweig der Mathematik, dersich mit der konkreten rechnerischen Lösung von Problemenbefasst, insbesondere mittels Algorithmen, die sich aufDigitalrechnern implementieren lassen.

• Konstruktivität, Ezienz und stabile Umsetzung erfordertdurchwegs Approximation.

• Die Numerik ist die Freundin der Angewandten Mathematik.Sie gilt aber auch ein eigenständiges mathematisches Teilgebiet(Numerische Mathematik). Entwicklung eng verbunden mitder Existenz und Weiterentwicklung des Digitalcomputers.

• `Scientic Computing': der zur Implementierung aufHochleistungsrechnern hin orientierte, praktische undangewandte Aspekt der Numerik.

35 / 35