Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Simulation komplexer technischer Anlagen
Teil II: Elemente zum Bau virtueller Anlagenkomponenten
Kapitel 8: Algorithmen 2: Gewöhnliche Differentialgleichungen
Inhalt• Gewöhnliche Differentialgleichungen
• Euler- und Differenzen-Verfahren
• Verfahren höherer Ordnung
• Systeme gewöhnlicher Differentialgleichungen
Experimente:
Explizite und implizite Zeitdiskretisierung am Beispiel der Wärmeleitgleichung
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Das sollten Sie heute lernen
– Was ist eine gewöhnliche Differentialgleichung– Numerische Lösung gew. Dgl‘en– Differenzenverfahren implizit und explizit– Eulerverfahren– Diskretisierung gew. Dgl‘en– Struktur von Programmen zur Lösung von Dgl‘en– Ursachen für instabiles Verhalten
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Grundbegriffe aus der Theorie der Differentialgleichungen
Gleichungen zwischen Funktionen und ihren Ableitungen heißen Differentialgleichungen (Dgl). Man unterscheidet gewöhnliche Differentialgleichungen und partielle Differentialgleichungen. Gewöhnliche Differentialgleichungen enthalten nur gewöhnliche Ableitungen. Als Beispiel sei die Schwingungsgleichung angeführt:
Differentialgleichungen mit einer unabhängigen Variablen sind immer gewöhnliche Differentialgleichungen. Differentialgleichungen mit mehreren unabhängigen Variablen enthalten gewöhnlich partielle Ableitungen nach den einzelnen Variablen. Man spricht dann von partiellen Differentialgleichungen. Sie werden entweder direkt oder in Operatorform angeschrieben. Sei L der Operator
so lautet die Schwingungsgleichung L y = f(t)
L kann verschieden definiert werden.
Die Ordnung n einer Differentialgleichung gibt die höchste Ableitung in der Differentialgleichung an.
)(22 tfyyy
222
2L
dtd
dt
d
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Gewöhnliche DifferentialgleichungenZu lösen sei in a t b
mit y(a) = yo
Typischerweise hat bei solchen Problemen die unabhängige Variable die Bedeutung der Zeit. Yo ist dann ein Anfangswert.
Von den Problemen, die im Rahmen dieser Vorlesung behandelt werden, fordern wir:
a) Sie müssen eine eindeutige Lösung y(t) haben.
b) Die Lösung darf nur vom Anfangswert abhängen.
c) Sie darf sich nur wenig ändern, wenn yo oder f wenig geändert werden.
Eine Differentialgleichung heißt
– linear, wenn keine Produkte von Ableitungen auftreten und die Koeffizienten nicht von den abhängigen Variablen oder ihren Ableitungen abhängen;
– halb-linear, wenn in den Randbedingungen nichtlineare Funktionen der Abhängigen oder ihren Ableitungen vorkommen;
– quasi-linear, wenn die Differentialgleichung in ihrer höchsten Ableitung linear ist,– nicht-linear, wenn die Differentialgleichung Potenzen auch der höchsten
Ableitung enthält.
),( tyfydt
d
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Lösung gewöhnlicher Differentialgleichungen
Integriert man
so erhält man
Das Intervall tn bis tn+1 heißt Zeitschritt n+1 für einen Zeitschritt gilt:
Daraus ergeben sich folgende Möglichkeiten, yn+1 zu bestimmen:
1. Integration der rechten Seite nach dem Newton-Verfahren
Euler- und Runge-Verfahren.2. Entwicklung der rechten Seite nach Lagrange-Funktionen und anschließende Integration
Adams-Verfahren.3. Näherung der Ableitung (linke Seite) durch eine Approximation der Ordnung n
Gear-Verfahren. Die Verfahren werden im Folgenden kurz erläutert.
btamittyfydt
d ),(
2
11 2
1......),(),(
),()()(tt
ta
tt dttyfdttyf
dttyfbaayby
dttyfyy nt
ntnn ,11
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Euler-VerfahrenDer Integrand wird durch einen konstanten Wert genähert. Dazu gibt es drei Möglichkeiten:
a) f (y,t) = f (yn,tn)
b) f (y,t) = f (yn+1, tn+1)
c) f (y,t) = f (yn+, tn+ ) mit 0 1
Die rechte Seite wird mit h = tn+1 - tn wird
Daraus folgen 3 Bestimmungsgleichungen für
a) explizites Verfahren: yn+1 = yn + h f (yn, tn)
b) implizites Verfahren: yn+1 = yn + h f (yn+1, tn+1) (entspricht Iterationsvorschrift)
c) modifiziertes Euler-Verfahren: yn+1 = yn + h f (yn+ , tn+)
Setzt man = 0, 5, so folgt Prediktorschritt
Korrektorschritt
Euler-Verfahren entsprechen der Differenzennäherung
),(222/1 tyfhyy nnn
1
nhf fdt n
n
tt
) t,(y fh y y21/n21/nn1n
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Diskretisierung einer elliptischen Differentialgleichung mit der Methode der finiten Differenzen am Beispiel der eindim. stationären Wärmeleitgleichung mit inneren Wärmequellen und vorgegebener RandtemperaturModelliert wird ein isolierter Stab der Laenge 1 m mit der Wäremeleitfähigkeit Lambda = 15 W/mK. Die innere Wärmequelle sei konstant über die gesamte Stablänge. Die linke und rechte Randtemperatur, sowie die Stärke der inneren Wärmequelle können variiert werden.Das Lösungsgebiet wird durch n Lösungspunkte beschrieben. Berechnet werden insgesamt 5 Lösungen, wobei n jeweils um 2 erhöht wird.
1-dim. stationäre Wärmeleitgleichung
Diskretisierung der 1-dim stationären Wärmeleitgleichung
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Diskretisierung und Lösung der Bewegungsgleichung
Daraus ergeben sich für die Systemmatrix M und die rechte Seite R bei 5 Zeitschritten mit folgenden Bedingungen: Anfangshöhe , Anfangsgeschwin-digkeit , Zeitschrittweite t.Jetzt muss nur noch das lineare Gleichungssystem gelöst werden:My=R. Wobei y der Lösungsvektor der Fallhöhe zu dem diskreten Zeitpunkt darstellt.Aufgabe: Bestimme Zeitschrittweite, sodass exakte Lösung und Näherung in der Zeichengenauigkeit übereinstimmen.
Die Differentialgleichung der Bewegungsgleichung lautet . Ihre ana-lytische Lösung ist . Um die Dgl. zu lösen muss sie zunächst in eine Differenzengleichung umgewandelt werden. An der Stelle i gilt für konstante Zeitschrittweite t: . Die Diskretisierung erfolgt auf dem Maschenrand. Am linken Rand (i=0) sind als Anfangs-bedingungen und vorzugeben. Damit lautet die Gleichung für die einzelnen Punkte:
ddt y g
2
2
y0 vy y
t00 1
y = y
= g
y = g
=
00
0 1 0
0 1 22
y y t v t
y y t
²
²
y0
v0
Der Versuch wird durchKlick gestartet
ddt
yy y y
tgi i i
2
22 1
2
2
R
yt g v t
t gt gt gt g
02
02
2
2
2
M
1 0 0 0 0 01 1 0 0 0 0
1 2 1 0 0 00 1 2 1 0 00 0 1 2 1 00 0 0 1 2 1
y t gt v t y( ) 12
20 0
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Systeme gewöhnlicher DifferentialgleichungenEin System der Ordnung m wird durch m-Gleichungen definiert
Diskretisiert man dieses System, so kann man für jede Gleichung eine der beschriebenen Methoden verwenden. Die Auswahl muss nach physikalischen Gesichtspunkten geschehen.
Haben die Gleichung verschiedene Zeitkonstanten, wird das Gleichungssystem "steif". Dann breiten sich Fehler stark aus (implizite Lösungen). Hat man Nichtlinearitäten zu betrachten, so müssen Newton- oder Newton-Raphson-Methoden zur Lösung verwendet werden.
),,2
,1
(
),,2
,1
(2
2
),,2
,1
(1
1
tmyyy
mf
dtm
dy
tmyyyf
dt
dy
tmyyyf
dt
dy
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
21
21
xiTiTi
T
dtiT
2121/i
2
2
xiTiTiT
dx
T
Die Wärmeleitgleichung als BeispielDie Wärmeleitgleichung ist eigentlich eine partielle Differentialgleichung. Sie lautet
Diskretisiert man zunächst den x-Raum, so erhält man x xi, T Ti und
oder am Ortspunkt i
Das ist eine gewöhnliche Differentialgleichung, allerdings aus einem System von Differentialgleichungen für alle diskreten Punkte des Ortsraumes. Zur Lösung des Problems benötigen wir
• Die Länge des Stabes• Die Zahl der Punkte i• Werte für T am linken und rechten Rand• Einen Wert von • Die Dauer der Simulation• Die Zahl der Zeitschritte• Die Temperaturverteilung zum Zeitpunkt 0.
Die Diskretisierung ist konsistent, wenn Orts- und Zeitableitung am gleichen Raum-Zeit-Punkt erfolgen .
2
2
x
T
t
T
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Die Wärmeleitgleichung in diskreter Form
Explizites Verfahren
Implizites Verfahren
Gemischtes Verfahren (Zwischenschrittverfahren)
Fehlerfortpflanzung beim expliziten Verfahren
Für cond g = verringert sich Fehler
wegen
folgt, dass beschränkt ist, t und x hängen also voneinander ab.
)121(1 TniTniTniTniTni
)11
1211(1 TniTniTniTniTni
)121(1 TniTniTniTniTni
TnTnTngTn )(1
1)(
1
Tng
gTn
)(1 Tng 2)/( xt
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Die Wärmeleitgleichung lautet .Für den Versuch sollen die folgenden Bedingungen gelten.
Die Diskretisierung erfolgt nach dem Differenzenverfahren mit • Lösungspunkte auf dem Maschenrand • Konstanten Maschenweiten
• Ansatz für Lösung
• Ansatz für Differentiale
• Konsistenzbedingung
ddt
T ddx
T2
2
T t T x T x TL L R R( ) , ( ) , ( )0 0 T T
x x t t dtdxi n , 2
T x t Tin( , )
ddt
T T Tt
ddx
TT T T
x
n n n
i
i i i
1 2
21 1
2
2
,
ddt
T ddx
Tn n
ii
=
2
2
Präsentation
Lösung explizitLösung implizit
Präsentation
Diskretisierung der Wärmeleitgleichung
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Lösung der Wärmeleitgleichung nach dem expliziten Verfahren 1/2
Beim expliziten Verfahren erfolgt der Ansatz für die Diskretisierung am Zeitpunkt . Es gilt also .Die Wärmegleichung in diskreter Form lautet: Nun wird die Gleichung so umgestellt, dass die Temperaturen eines Zeitpunk-tes auf einer Seite stehen: , in Matrixschreibwei-se: . E ist die Einheitsmatrix und B eine Tridiagonalmatrix. Berücksichtigt man die Vorgabe Konstanter Randtemperaturen so gilt für B:
Die Anfangsbedingungen, d.h. im vorgegebenen Beispiel die Anfangstemper-aturen des Stabs müssen in den Startvektor eingehen. (Für sich mit der Zeit ändernde Randbedingungen müsste bei dieser Diskretisierung der Lösungs-vektor vor jedem Zeitschritt in den Werten, die den Rand betreffen - hier erste und letzte Komponente - modifiziert werden.)
T T T T Tin
in
in
in
in
11 12( )
E T B Tn n 1T T T Ti
nin
in
in
11 11 2 ( )
B
1 0 0 0 0 0 0 01 2 0 0 0 0 0
0 1 2 0 0 0 00 0 1 2 0 0 00 0 0 1 2 0 00 0 0 0 1 2 00 0 0 0 0 1 20 0 0 0 0 0 0 1
tn 0
T0
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Lösung der Wärmeleitgleichung nach dem expliziten Verfahren
Um das Gleichungssystem , zu lösen, kann die Einheits-matrix wegfallen und der jeweils neu berechnete Temperatur-vektor wird rekursiv als Ausgangsvektor für den folgenden Zeitschritt benützt. Das explizite Verfahren konvergiert aber nur für , bei größeren Werten zeigt sich, dass die Lösung instabil ist.
Bei diesem Versuch wird das Temperaturprofil eines idealisierten Stabs (ein- dimensional) bei fortlaufender Zeit visualisiert. Die Anfangstemperatur beträgt 20 Grad. Der Stabanfang wird konstant auf 100 Grad gehalten, das Stabende konstant auf 20 Grad. Für die Temperaturleitfähigkeit wird 1 angesetzt. Die Matrix B wird ebenfalls dargestellt und der erste berechnete Temperatur-vektor in schriftlicher Form gezeigt.
Aufgabe:Bestimmen Sie die Zahl der Zeitschritte bis zur stationären Lösung für einen Stab der Länge 5 und 10.Untersuchung der Stabilität undder Genauigkeit.
E T B Tn n 1
Der Versuch wird durchKlick gestartet
T B Tn n 1
Tn1 Tn
0 5.
2/2
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Lösung der Wärmeleitgleichung nach dem impliziten 1/1 Verfahren
Beim impliziten Verfahren erfolgt der Ansatz für die Diskretisierung am Zeitpunkt . Es gilt also für die Wegdiskretisierung.Die Wärmegleichung in diskreter Form lautet: Nun wird die Gleichung so umgestellt, dass die Temperaturen eines Zeitpunk-tes auf einer Seite stehen: , in Matrixschreibwei-se: . E ist die Einheitsmatrix und B eine Tridiagonalmatrix. Berücksichtigt man die Vorgabe Konstanter Randtemperaturen so gilt für B:
Die Anfangsbedingungen, d.h. im vorgegebenen Beispiel die Anfangstemper-aturen des Stabs müssen in den Startvektor eingehen.(Für sich mit der Zeit ändernde Randbedingungen müsste bei dieser Diskretisierung der Lösungs-vektor vor jedem Zeitschritt in den Werten, die den Rand betreffen - hier erste und letzte Komponente - modifiziert werden.)
1
T0
tn1
T T T T Tin
in
in
in
in
1
11 1
112( )
T T T Tin
in
in
in
1
1 11
12( ) ) B T E Tn n 1
B
1 0 0 0 0 0 0 01 2 0 0 0 0 0
0 1 2 0 0 0 00 0 1 2 0 0 00 0 0 1 2 0 00 0 0 0 1 2 00 0 0 0 0 1 20 0 0 0 0 0 0 1
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und EnergiesystemeUm das Gleichungssystem , zu lösen, kann die Einheits-matrix wegfallen und der jeweils neu berechnete Temperatur-vektor wird rekursiv als Ausgangsvektor für den folgenden Zeitschritt benützt. Das implizite Verfahren konvergiert für alle Werte von . Der Rechenaufwand beim impliziten Verfahren ist wegen der damit verbundenen Lösung größer als beim expliziten, dafür bleibt aber das Verfahren für alle Werte stabil.Bei diesem Versuch wird das Temperaturprofil eines idealisierten Stabs (ein- dimensional) bei fortlaufender Zeit visualisiert. Die Anfangstemperatur beträgt 20 Grad. Der Stabanfang wird konstant auf 100 Grad gehalten, das Stabende konstant auf 20 Grad. Für die Temperaturleitfähigkeit wird 1 angesetzt. Die Matrix B wird ebenfalls dargestellt und der erste berechnete Temperatur-vektor in schriftlicher Form gezeigt.Aufgabe:Bestimmen Sie die Zahl der Zeitschritte bis zur stationären Lösung für einen Stab der Länge 5 und 10.Untersuchung der Stabilität undder Genauigkeit.
Tn1 Tn
B T E Tn n 1
B T Tn n 1
Lösung der Wärmeleitgleichung nach dem impliziten Verfahren 1/2
Der Versuch wird durchKlick gestartet
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Aufbau eines Programms zur Lösung von Dgl‘en
Die numerische Lösung von Differentialgleichungen kann - zumindest für einfache Probleme - schnell und übersichtlich programmiert werden. Im Folgenden ist ein typisches Flussdiagramm gezeigt.
Zeitfortschaltung, Eigenwertiteration,
Nichtlinearität
Lösung des Gleichungssystems
Berechnung der rechten Seiten
Neue
Matrizen
Nicht linear
ja
nein
nein ja (bei Quellrechnung, Endzeitpunkt,
Endgenauigkeit)
Eingabe und ihre Verarbeitung Geometrie, Materialdaten,
Randbedingungen, Anfangswerte
Ausgabe
Nein linear
Erzeugung des Gleichungssystems
Ende
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Runge-Kutta-Verfahren
Verwendet man zur Integration der rechten Seite Verfahren höherer Ordnung, so erhält man die Klasse der Runge-Kutta-Verfahren zur Lösung von gewöhnlichen Differentialgleichungen.
a) Integration mit Trapez-Regel
Verfahren von Heun
Lösung iterativ mit Startwert
b) Iteration mit Simpson-Regel
Die Simpson-Regel verwendet die Punkte tn,tn+1/2 und t n+1 zur Integration, f (y,t) muss also an diesen Punkten genähert werden. Dies leistet gerade das Runge- Kutta-Verfahren der Ordnung 4.
Die Zwischenwerte werden wie folgt genähert:
Für f (y,t) = f (t) degeneriert das Verfahren zur Simpson-Formel.
)),(),((2
111 tytyyy nnnnnn ffh
),(1 tyyy nnnn
fho
),(),2
(2),2
(2),(6/14321
( ytytytytyy hfh
fh
ffhnnnnnnn
),(),2
(2/1
),2
(2/1
3412
231
yhtfhyyyh
tfhyy
yh
tfhyyyy
nnnn
nnn
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Beispiel zum Runge-Kutta-Verfahren
Gegeben sei das Anfangswertproblem: = y2 mit
y (0) = - 4 0 t 0,3 h = 0,1
Die exakte Lösung lautet
Für und
Für den Schritt n+1 folgt:
2)2)2)(201(
101()4,(
2)2)(201(
201
101
4
2)2)2)(201(
201()3,2(
22201
201
3
2)2)(201()2,2(2
201
2
2)1,n(tf
)41(/1
2
nynynyyhntfund
nynynynyy
nynynyyhntfundnynynyy
nynyyhntfundnynyy
wirdnyy
ty
y
nyy 1
222220/120/110/12
22220/120/122220/12260/1
1
ny
ny
ny
ny
ny
ny
ny
ny
ny
iy
ny
ny
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Adams-Verfahren
Eine weitere Verbesserung der Bestimmung der rechten Seite erhält man dadurch, dass man den Integranden über eine Funktionsentwicklung darstellt.
Als Entwicklungskoeffizienten können die Werte der diskreten Zeitpunkte verwendet werden.
Entwicklungsfunktionen sind dann wieder die Lagrange-Polynome .
Zum Zeitschritt n ist folgende Entwicklung möglich:
mit m n
Damit lässt sich f (y, t) integrieren.
Die mi sind tabelliert.
m
i
miin
tin
yftyf0
),(),(
n
i
mim
infdtjn
t
jnt
mijn
t
jnt
m
iin
fdttyf
0111 0
,
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Baskford-Adams-VerfahrenMan unterscheidet zwei Fälle:
a) j = 1 Verfahren nach Adams-Baskford
Bestimmung von yn+1:
Integration zwischen tn und tn+1
Entwicklung von f bis zur Stelle tn
Aus Entwicklung bis tn wird Verlauf extrapoliert - Prediktor-Schritt. Für die Integrale ni gilt:
I
ni
0 1 2 3 4
ni 1
2ni 1 1
12ni 5 1 -1
24ni 9 19 -5 1
720ni 251 646 -264 106 -19
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Adams-Moulton-Verfahrenb) j = 0 Verfahren nach Adams-Moulton
Bestimmung von yn+1
Integration zwischen tn und tn+1
Entwicklung von f bis zur Stelle tn+1 (ersetze n durch n+1 in allgemeiner Formel)
Entwicklung verwendet schon Endwert - Korrektor-Schritt oder implizit. Lösungen nur iterativ.
Für die Integrale ni gilt:
I
ni
0 1 2 3 4
ni 1
2ni 1 1
12ni 5 1 -1
24ni 9 19 -5 1
720ni 251 646 -264 106 -19
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Gear-Verfahren - 1
kn
r
kknn
kn
r
kknnn
nrkkn
r
nrk
r
kknn
rk
r
kkn
ybfa
yoder
ybayafydt
dmanerhält
kfürdt
d
abund
dt
damit
tdt
dyty
dt
dDaraus
tyty
11
1
11
11
110
10
11
01
1
01
)()(
)()(
Die Klasse der Gear-Verfahren erhält man, wenn man nicht den Integranden,sondern die Lösung nach Lagrange-Funktionen entwickelt und anschließend ableitet
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Gear-Verfahren - 2
nnn
knn
kk
r
kjj jnkn
jnnk
knrk
r
j jnn
yftyund
bt
a
tta
pb
rfürtt
ttp
rfürpdt
d
tta
1
1
11
1 11
11
1
1 11
11
konstantt und 1 rfür manerhält Daraus
1
11
1Die allgemeine Form der Koeffizienten lautet:
Das entspricht dem Ergebnis nach Euler
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Gear-Verfahren -3
nn
n
nnnnn tt
t
ttttta
1
1
1111
1111
man erhält 2 r Für
nn
nn
nn
nn
tt
ttp
tt
ttp
1
12
1
111
111
111
1
11
21k
3
1
3
4
3
2
2
43
3
1
3
4
3
2
3
1
3
4
3
22
3 a
nnnn
nnnn
n
nnnn
kk
yyfty
ft
yyyy
dt
d
yyftyund
bbk
p
tka
pb
t
konstanttfür manerhält Daraus
:1 giltoderimplizitifferenzenRückwärtsddieFür
Daraus folgt
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Stabilität -1
Drei Fehlerquellen können auftreten
a) Näherung von
b) Näherung des Integrals der rechten Seite,
c) Bestimmung von y.
Werden Fehler durch die Zeitfortschaltung verkleinert, so heißt ein Verfahren stabil. Analog der notwendigen Bedingung für die Verkleinerung von Fehlern bei der Iteration gilt auch bei der Zeitfortschaltung.
ydt
d
Bedingungnotwendigealsygtfo
ygyg
ygygcondwegen
gcondwennstabilistygyn
1)(lg
)()(
)(
1,1
Simulation technischer Systeme, WS 02/03 Kap. 8: Algorithmen-2
Universität Stuttgart
Institut für Kernenergetik und Energiesysteme
Stabilität -2
Mehrschrittverfahren kann man darstellen als
Dann gibt es zu dieser Gleichung ein charakteristisches Polynom p ()
Für p () = 0 gibt es m-Nullstellen und man kann zeigen:
Ein Verfahren ist stabil, wenn für alle Nullstellen gilt
und Nullstellen mit höchstens als einfache Nullstellen vorkommen.
Verfahren nach Adam und Gear sind stabil. Die Verfahren von Runge-Kutta sind schwach stabil. Explizite Euler-Verfahren sind für große Zeitschritte instabil.
0
11
22
11
aammam
mamp
1
1
),(10
....1211
tyFhmn
yany
ma
ny
ma
ny