Universität Stuttgart Institut für Kernenergetik und Energiesysteme Simulation technischer...

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




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




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




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




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




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




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




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




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




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




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




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




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




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


T B Tn n 1

Tn1 Tn

0 5.

2/2




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



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





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




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




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




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

,




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




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




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




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




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




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




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

Universität Stuttgart Institut für Kernenergetik und Energiesysteme Simulation technischer...

Documents

Transcript of Universität Stuttgart Institut für Kernenergetik und Energiesysteme Simulation technischer...