Numerische Methoden in der angewandten Thermodynamik ... · Numerische Methoden in der angewandten...

Numerische Methoden in der angewandten Thermodynamik

Numerische Grundlagen

Skriptum

Institut für Verbrennungskraftmaschinen und Thermodynamik

Technische Universität Graz

Dr. Günter OFFNER

++43 316 787 2103

[email protected]

SS 2015

2

Vorwort

Liebe Kolleginnen und Kollegen!

Die Durchführung dieser Lehrveranstaltung in den vergangenen Jahren zeigte, dass auf Grund

der Komplexität und auf Grund des Umfangs des Stoffes der alleinige Vortrag an der Tafel

oder via Powerpoint nicht ausreicht, um einen repräsentativen Einblick in dieses

Vorlesungsthema zu geben.

Es zeigte sich auch, dass dem eigentlichen Vorlesungsinhalt ein einführender Teil

vorangestellt werden muss, welcher auf „Numerische Grundlagen“ eingeht.

Als Konsequenz wurde dieses Skriptum ausgearbeitet. Das Skriptum umfasst neben dem Teil

„Numerische Grundlagen“ auch noch ein Kapitel „Einführung in Visual Basic“. In diesem ist

eine Kurzzusammenfassung der in der Programmiersprache Visual Basic (VB) gängigsten

Befehle enthalten. Zum Zweck eines besseren Verständnisses sind im Skriptum neben einer

Vielzahl an Beispielen auch diverse Verweise auf kleine, zum Skriptum gehörige,

Computerprogramme enthalten. Diese sind in der Programmiersprache VB geschrieben und

sollen den diskutierten Stoff mit weiteren anschaulichen Beispielen unterlegen. Die

Computerprogramme werden auch im Rahmen der Vorlesung erörtert, sollen aber auf jeden

Fall auch von jedem/r einzelnen ausprobiert und „erforscht“ werden.

Graz im Februar 2015 Günter Offner

3

Inhaltsverzeichnis

1 EINFÜHRUNG IN VISUAL BASIC ...................................................................... 5

1.1 Visual Basic Editor ...................................................................................................................................... 5

1.1.1 Starten des VB Editors ......................................................................................................................... 5

1.1.2 Programmieren von VB Code .............................................................................................................. 5

1.1.3 Ausführen von Programmen ................................................................................................................ 5

1.2 Technische Vorbemerkungen ..................................................................................................................... 5

1.3 Operatoren ................................................................................................................................................... 5

1.4 Zuweisungen ................................................................................................................................................ 6

1.5 Eingabe / Ausgabe Operationen ................................................................................................................. 6

1.6 Sprunganweisungen .................................................................................................................................... 6

1.7 Kontrollstrukturen ...................................................................................................................................... 7

1.7.1 Vollständige Alternative ...................................................................................................................... 7

1.7.2 Fallauswahl .......................................................................................................................................... 7

1.7.3 Zählergesteuerter Zyklus ...................................................................................................................... 7

1.7.4 Kopfgesteuerte Schleifen ..................................................................................................................... 7

1.7.5 Fußgesteuerte Schleifen ....................................................................................................................... 8

1.8 Datentypen ................................................................................................................................................... 8

1.8.1 Deklaration des Datentyps von Variablen ............................................................................................ 8

1.8.2 Deklaration von Arrays ........................................................................................................................ 9

1.8.3 Bemerkungen zu Datentypen ............................................................................................................... 9

1.9 Funktionen, Prozeduren und Makros ........................................................................................................ 9

1.10 Gültigkeitsbereich und Lebensdauer von Variablen ........................................................................... 9

2 NUMERISCHE GRUNDLAGEN ........................................................................ 11

2.1 Einleitung ................................................................................................................................................... 12

2.2 Numerischer Algorithmus ........................................................................................................................ 14

2.3 Numerischer Fehler – Residuum – Kondition ........................................................................................ 15

2.4 Lineare Gleichungssysteme ...................................................................................................................... 16

2.4.1 Eigenschaften von linearen Gleichungssystemen .............................................................................. 16

2.4.1.1 Normen und Skalarprodukt ...................................................................................................... 17

2.4.1.2 Symmetrie, Orthogonalität und Definitheit von Matrizen ........................................................ 17

2.4.1.3 Besetzungsstruktur ................................................................................................................... 18

2.4.2 Methoden zum Lösen linearer Gleichungssysteme ............................................................................ 19

2.4.2.1 Direkte Verfahren ..................................................................................................................... 19

2.4.2.1.1 Gaußsches Eliminationsverfahren ....................................................................................... 19

2.4.2.1.2 Cholesky – Verfahren .......................................................................................................... 20

2.4.2.2 Iterative Verfahren.................................................................................................................... 21

2.4.2.2.1 Abbruchkriterien .................................................................................................................. 22

2.4.2.2.2 Jacobi – Verfahren ............................................................................................................... 22

2.4.2.2.3 Gauß – Seidel – Verfahren................................................................................................... 23

2.4.2.2.4 Relaxationsverfahren ........................................................................................................... 23

2.4.2.2.5 Methode der konjugierten Gradienten ................................................................................. 24

2.4.2.2.6 Moderne iterative Verfahren ................................................................................................ 25

4

2.4.2.3 Vorkonditionierung .................................................................................................................. 26

2.4.2.4 Kompakte Speicherung ............................................................................................................ 27

2.5 Diskretisierungsverfahren ........................................................................................................................ 28

2.5.1 Diskretisierung des Problemgebiets ................................................................................................... 28

2.5.1.1 Netztypen.................................................................................................................................. 29

2.5.1.2 Netzstruktur .............................................................................................................................. 30

2.5.1.3 Erzeugung strukturierter und unstrukturierter Gitter ................................................................ 32

2.5.2 Diskretisierung von Ortsvariablen ..................................................................................................... 32

2.5.2.1 Finite – Differenzen – Methode ................................................................................................ 33

2.5.2.2 Finite – Volumen – Methode .................................................................................................... 34

2.5.2.3 Multigrid – Verfahren (Mehrgitter – Verfahren) ...................................................................... 35

2.5.3 Diskretisierung von zeitlichen Variablen (Zeitintegration) ................................................................ 37

2.5.3.1 Grundlagen ............................................................................................................................... 38

2.5.3.2 Explizite Einschrittverfahren .................................................................................................... 39

2.5.3.3 Implizite Einschrittverfahren .................................................................................................... 41

2.5.3.4 Kombinationen mehrerer Einschrittverfahren .......................................................................... 42

2.5.3.5 Explizite und implizite Mehrschrittverfahren ........................................................................... 42

2.6 Eigenschaften von Diskretisierungsverfahren ........................................................................................ 43

2.6.1 Konsistenz .......................................................................................................................................... 44

2.6.2 Stabilität ............................................................................................................................................. 45

2.6.3 Konvergenz ........................................................................................................................................ 46

2.7 Übungsbeispiele ......................................................................................................................................... 47

3 ABBILDUNGSVERZEICHNIS .......................................................................... 55

4 LINK- UND LITERATURVERZEICHNIS ........................................................... 56

4.1 Einführung in Visual Basic (VBA) ........................................................................................................... 56

4.2 Numerische Grundlagen ........................................................................................................................... 56

5

1 Einführung in Visual Basic

Die im Rahmen dieser Vorlesung erläuterten Algorithmen werden auch in Programmform zur Verfügung gestellt. Die einzelnen Programme wurden in der Programmiersprache Visual Basic (VBA oder VB) programmiert, da der Visual Basic Editor im Rahmen von Microsoft EXCEL Anwendungen zur Verfügung steht. Ziel dieser Programme ist es, die in der Vorlesung besprochenen Algorithmen auch in programmierter Form zu sehen und ausprobieren zu können.

1.1 Visual Basic Editor Mit Hilfe des Visual Basic Editors kann neben weiteren Funktionalitäten, Quellcode eingegeben und bearbeitet werden. Im Folgenden wird eine Kurzeinführung in diese Programmierumgebung sowie in die Programmiersprache Visual Basic gegeben.

1.1.1 Starten des VB Editors

Das Starten des VB Editors erfolgt über die Tastenkombination Alt+F11 oder über Extras – Makro – Visual Basic Editor.

1.1.2 Programmieren von VB Code

Das Programmieren erfolgt im VB Editor. Dabei stehen sämtliche Standardfunktionalitäten wie Kompilieren, Debuggern, … im Dialogmenü zur Verfügung.

1.1.3 Ausführen von Programmen

Das Ausführen eines kompilierten Programms erfolgt wiederum über Extras – Makro – Makros … . Da selbst geschriebene Programme in der Regel häufig aufgerufen werden, ist es dienlich, die Visual Basic Symbolleiste (und somit auch den darin enthaltenen Startbutton) über Ansicht – Symbolleisten zu aktivieren.

1.2 Technische Vorbemerkungen

• Es gibt keine Zeichen um das Ende eines Befehls zu markieren. Am Ende einer Zeile muss der Befehl fertig sein oder mit einem Zeilenfortsetzungszeichen (Leerzeichen) gefolgt von einem Unterstrich) fortgesetzt werden.

• Groß- und Kleinschreibung wird nicht unterschieden.

• Einzelne Zeilen werden als Kommentar erkannt, wenn sie mit einem Hochkomma beginnen.

• Da es für selbst geschriebene Programme keine „Undo“ – Funktion gibt, sollte in der Testphase stets mit Kopien gearbeitet werden.

1.3 Operatoren Folgende im Rahmen dieser Lehrveranstaltung benötigte Operatoren stehen in VBA zur Verfügung.

• Arithmetisch: ^, *, /, +, -, \ (Ganzzahldivision), mod

• Zeichenketten: +, & („15“ & 3 ergibt „153“)

• Vergleichsoperatoren: =, <, >, <=, >=, <>

6

• Logische Operatoren: and, or, not, xor, imp (Implikation), eqv (Äquivalenz)

1.4 Zuweisungen Die Zuweisung von Werten einer Variable kann auf zwei Arten erfolgen.

• Wertübergabe: =

• Parameterübergabe in Prozeduren: :=

1.5 Eingabe / Ausgabe Operationen Zur Ein- und Ausgabe stehen unterschiedliche Werkzeuge zur Verfügung

• InputBox(…): Abfragen einer Eingabe des Benutzers

• MsgBox: Benutzerinformation

• Debug.print: Ausgabe ins Direktfenster des Debuggers

• Worksheets("EXCEL Beispiel 3").Cells(i, j): Eingabe/Ausgabe aus/in EXCEL Sheet „EXCEL Beispiel 3“ Zelle mit Zeilenindex „i“ und Spaltenindex „j“

Beispiel 1.1:

Antwort = InputBox(„Fragetext“)

Beispiel 1.2:

MsgBox „Text“ & Variablenwert & „weiterer Text“

Beispiel 1.3:

Debug.pring „Text“

Beispiel 1.4:

iZeile = 2

iSpalte = 2

Worksheets("EXCEL Beispiel 3").Cells(iZeile, iSpalte) = 20.0

1.6 Sprunganweisungen Die im Folgenden aufgelisteten Sprunganweisungen stehen zur Verfügung.

• Goto label

• Exit do

• Exit for

• Exit function

• Exit property

• Exit sub

7

1.7 Kontrollstrukturen

1.7.1 Vollständige Alternative

Beispiel 1.5:

If (variable1 = 1) then

variable2 = 1000

Else if (variable1 = 2) then

Variable2 = 2000

Else

Variable2 = 3000

End if

1.7.2 Fallauswahl

Beispiel 1.6:

Case 1, 3, 5, 7, 8, 10, 12

ZahlDerTage = 31

Case 4, 6, 9, 11

ZahlDerTage = 30

Case 2 ´ Schaltjahre werden nicht berücksichtigt

ZahlDerTage = 28

Else

MsgBox “Falsch Montatseingabe”

End select

1.7.3 Zählergesteuerter Zyklus

Beispiel 1.7:

For iZaehler = 1 to 10

…

Next iZaehler

1.7.4 Kopfgesteuerte Schleifen

Beispiel 1.8:

Do While meineZahl > 10

meineZahl = meineZahl – 1

iZaehler = iZaehler + 1

Loop

Beispiel 1.9:

Do unitil meineZahl = 10


8


Loop

1.7.5 Fußgesteuerte Schleifen

Beispiel 1.10:

Do



Loop while meineZahl > 10

Beispiel 1.11:

Do



Loop unitil meineZahl = 10

1.8 Datentypen

1.8.1 Deklaration des Datentyps von Variablen

Wird für eine Variable kein Datentyp festgelegt, so wird automatisch der Datentyp Variant angenommen. Variablen vom Typ Variant können Zeichenfolgen, Datums-, Zeit-, boolesche oder numerische Werte enthalten und diese Werte automatisch umwandeln. Da keine explizite Deklaration zwingend erforderlich ist, werden Variablen, die vom Programmierer nur versehentlich falsch geschrieben wurden, als neue Variablen betrachtet. Es folgt dann eine oft aufwendige Fehlersuche. Mit der Option Explicit am Modulanfang kann die ausdrücklich empfohlene Deklaration von Variablen erzwungen werden.

Folgende Datentypen stehen zur Verfügung:

• Boolean: true, false

• Byte: ganze Zahl zwischen 0 und 255

• Integer: ganze Zahl zwischen - 32 767 und 32 767

• Long: ganze Zahl zwischen - 2 147 483 648 und 2 147 483 648

• Single: Gleitkommazahl ca. von – 3.4E38 bis 3.5E38 (8 Stellen Genauigkeit)

• Double: Gleitkommazahl ca. von – 1.8E308 bis 1.8E308 (16 Stellen Genauigkeit)

• Currency: Festkommazahl mit 15 Stellen vor und 4 Stellen nach dem Komma

• String (fester Länge): bis zu 65 400 Zeichen

• String (dynamischer Länge): bis ca. 2 Milliarden Zeichen

Die Deklaration des Datentyps erfolgt mit dem Schlüsselwort DIM.

Beispiel 1.12:

Die Zählvariable iZaehler wird als ganze Zahl deklariert.

DIM iZaehler as integer

9

1.8.2 Deklaration von Arrays

Aus den einfachen Datentypen können zusammengesetzte Datentypen gebildet werden. Für die Anwendungen im Rahmen dieser Vorlesung sind in diesem Zusammenhang aufeinander folgende, indizierte Elemente gleichen Datentyps (= Arrays) von Bedeutung.

Beispiel 1.13:

Das Datenfeld Temperatur, bestehend aus 100 Gleitkommazahlen wird deklariert. Anschließend werden in einer Schleife alle Felder Temperatur mit dem Wert 25.0 belegt.

DIM Temperatur(1 to 100) as double


For iZaehler = 1 to 100

Temperatur(iZaehler) = 25.0

Next iZaehler

1.8.3 Bemerkungen zu Datentypen

Das Datentyp - Konzept von VBA unterscheidet sich von dem anderer Programmiersprachen. Der Programmierer kann durch die Deklaration eines Datentyps für eine Variable den nötigen Speicherplatz und die Zugriffsgeschwindigkeit steuern. Die Deklaration erfüllt jedoch nicht den Zweck, den Programmierer vor Flüchtigkeitsfehlern zu schützen. Daher ist gerade bei der Deklaration auf Genauigkeit und Konsistenz zu achten.

Beispiel 1.14:


iZaehler = 2.5

Weist der Programmierer einer Integer-Variablen den Wert „2.5“ zu, so wird ohne Fehlermeldung gerundet und der Wert „2“ verwendet.

Beispiel 1.15:

DIM Temperatur, Dichte as double

Die Deklaration liefert Dichte als double und Temperatur als variant.

1.9 Funktionen, Prozeduren und Makros Kleinste selbständige Einheiten eines VBA - Programms sind Prozeduren und Funktionen. Funktionen geben einen einzelnen Wert zurück (vgl. EXCEL – Funktionen Summe oder Cos). Makros sind Prozeduren, die von EXCEL automatisch generiert werden, indem die Aktionen des Nutzers aufgezeichnet werden. Funktionen, Prozeduren und Makros werden als Teil des Projektes in so genannten Modulen gespeichert. Bevor man mit dem Programmieren beginnt, muss zunächst ein neues Modul erzeugt werden (Option: Einfügen Modul in der Visual Basic Entwicklungsumgebung). Verfügt man bereits über von EXCEL aufgezeichnete Makros, so wurde das Modul bereits automatisch eingefügt. Funktionen werden wie die Standard – EXCEL – Funktionen über die Option Einfügen/Funktion aufgerufen (Kategorie: „benutzerdefiniert“). Prozeduren werden genau wie die von EXCEL – automatisch generierten Makros über die Option Extras/Makro/Makros gestartet.

1.10 Gültigkeitsbereich und Lebensdauer von Variablen Variablen, die nur in einer Prozedur oder Funktion bekannt sein sollen, deklariert man innerhalb derselben mit dem Schlüsselwort DIM (mit dem Schlüsselwort Private kann

10

analoges bewirkt werden). Mit dem Schlüsselwort Static können lokale Variable deklariert werden, welche ihren Wert zwischen den einzelnen Prozeduraufrufen beibehalten. Diese Variablen werden erst wieder beim Zurücksetzen des Moduls neu initialisiert. Variablen, die in allen Prozeduren/Funktionen eines Moduls gültig sein sollen, deklariert man im Deklarationsteil eines Moduls (Abschnitt vor der 1. Prozedur/Funktion). Variablen, die darüber hinaus in allen Modulen Verwendung finden, muss man mit dem Schlüsselwort Public deklarieren. Auf Modulebene deklarierte Variablen erhalten ihren Wert, bis das Modul zurückgesetzt bzw. neu gestartet wird (Option: Ausführen/Zurücksetzen). Wird eine Variable nicht deklariert, so ist sie standardmäßig nur lokal gültig. Die verschiedenen Deklarationsarten unterscheiden sich auch in der Dauer, während der eine Variable ihren Wert behält.

Beispiel 1.16:

Die Variable x soll nur im aktuellen Modul bekannt sein. Die Variable y soll in allen Modulen bekannt sein.

DIM x as double

Public y as double

Private z as double

Beispiel 1.17:

Die Variable x soll eine lokale Variable sein und daher nur in der aktuellen Funktion/Prozedur bekannt sein. Die Variable behält ihren Wert nur bis zur Beendigung der Prozedur bzw. der Funktion. Die Variable soll zwar lokal sein, soll aber zwischen den einzelnen Prozeduraufrufen ihren Wert beibehalten. Die Variable y soll erst wieder beim Zurücksetzen des Moduls neu initialisiert werden.

Sub Prozedurname()

DIM x as double

Static y as double

End sub

11

2 Numerische Grundlagen

Die Funktionalität oder Effizienz von technischen oder naturwissenschaftlichen Systemen ist stets durch gewisse Eigenschaften der Systeme bestimmt. Kenntnisse über diese Eigenschaften sind meist der Schlüssel für das Verständnis oder ein Ansatzpunkt für die Optimierung der Systeme.

Im Bereich des Maschinenbaus spielen in diesem Zusammenhang insbesondere Festkörper- und Strömungseigenschaften, wie z. B.

• Deformationen oder Spannungen,

• Strömungsgeschwindigkeiten, Druck oder Temperaturverteilungen,

• Widerstands oder Auftriebskräfte,

• Druck- oder Energieverluste,

• Wärmeübergangs- oder Stofftransportraten,

• …

eine wichtige Rolle. Für Ingenieuraufgaben ist die Untersuchung solcher Eigenschaften in der Regel im Zuge von Neu- oder Weiterentwicklungen von Produkten und Prozessen von Bedeutung.

Um Kenntnisse über Eigenschaften von Systemen zu gewinnen, werden unterschiedliche Vorgehensweisen angewandt:

• Theoretische Methoden,

• Experimentelle Untersuchungen,

• Numerische Simulationen.

Theoretische Methoden, d. h. analytische Betrachtungen der das Problem beschreibenden Gleichungen, sind bei praxisrelevanten Problemstellungen nur sehr bedingt anwendbar. Die Gleichungen, die zur realistischen Beschreibung der Prozesse herangezogen werden müssen, sind meist so komplex (in der Regel Systeme von nichtlinearen partiellen Differentialgleichungen), dass sie nicht mehr analytisch lösbar sind. Vereinfachungen, die erforderlich wären, um eine analytische Lösung möglich zu machen, sind oft nicht zulässig und führen zu ungenauen Ergebnissen (und damit möglicherweise zu falschen Schlussfolgerungen). Allgemeiner gültige Überschlagsformeln, wie sie von Ingenieuren gerne benutzt werden, lassen sich für komplexere Systeme in der Regel nicht aus rein analytischen Betrachtungen ableiten.

Bei experimentellen Untersuchungen ist man bestrebt durch Versuche (an Modellen oder am realen Objekt) unter Einsatz von Geräten und Messinstrumenten an die notwendigen Systeminformationen heranzukommen. Dies bereitet in vielen Fällen Probleme:

• Messungen an realen Objekten sind oft schwierig oder gar unmöglich, weil beispielsweise die Dimension zu klein oder zu groß ist (z. B. Mikrosystemtechnik oder Erdatmosphäre), die Vorgänge zu langsam oder zu schnell ablaufen (z. B. Korrosionsvorgänge oder Explosionen, die Objekte nicht direkt zugänglich sind (z. B. menschlicher Körper) oder der zu untersuchende Vorgang im Verlauf der Messung gestört werden muss (z. B. Quantenmechanik).

• Schlussfolgerungen aus Modellversuchen zum realen Objekt sind, beispielsweise auf Grund unterschiedlicher Randbedingungen, oft nicht direkt vollziehbar (z. B. Flugzeug im Windkanal und im realen Flug).

12

• Experimente verbieten sich aus Sicherheits- bzw. Umweltgründen (z. B. Folgen eines Tankerunglücks oder eines Unfalls in einem Kernreaktor)

• Experimente sind oft sehr teuer und zeitraubend (z. B. Crashversuche, Windkanalkosten, Modellanfertigung, Parametervariationen, nicht alle interessierenden Größen können gleichzeitig gemessen werden).

Neben (oder besser zwischen) theoretischen und experimentellen Vorgehensweisen, hat sich in den letzten Jahren die numerische Simulation, in diesem Zusammenhang auch oft als Technisch – Wissenschaftliches Rechnen bezeichnet, als weitgehend eigenständige Wissenschaft etabliert, in deren Rahmen solche Untersuchungen mittels numerischer Berechnungsmethoden am Computer durchgeführt werden. Die Vorteile numerischer Simulation gegenüber rein experimentell durchgeführten Untersuchungen liegen auf der Hand:

• Numerische Ergebnisse können oft schneller und mit geringeren Kosten erhalten werden.

• Parametervariationen sind am Rechner meist sehr leicht durchführbar.

• Eine Simulation liefert meist umfassende Informationen durch globale und gleichzeitige Berechnung verschiedener problemrelevanter Größen (z. B. Temperatur, Druck, Luftfeuchtigkeit und Wind bei Wettervorhersage).

Eine wichtige Voraussetzung diese Vorteile auch zu nutzen, ist natürlich die Zuverlässigkeit der Berechnungen. Die Möglichkeiten hierzu haben sich auf Grund der Entwicklungen in den letzten Jahren signifikant verbessert (Abbildung 2.1), was sehr zum Aufschwung numerischer Simulationstechniken beigetragen hat. Dies bedeutet jedoch nicht, dass experimentelle Untersuchungen überflüssig sind oder werden. Es ist in erforderlich, beide Gebiete, Numerik und Experiment, weiterzuentwickeln und, am besten in komplementärer Weise einzusetzen, um zu optimalen Lösungen für die unterschiedlichen Aufgabenstellungen zu gelangen.

Abbildung 2.1: Beschleunigung im Bereich der Rechnergeschwindigkeit (grün) und der numerischen Verfahren (blau) in den letzten Dekaden

2.1 Einleitung Zur Illustration der unterschiedlichen Teilaspekte, die bei der Anwendung numerischer Simulationsmethoden zur Lösung ingenieurwissenschaftlicher Problemstellungen eine Rolle spielen, ist in Abbildung 2.2 die generelle Vorgehensweise schematisch dargestellt.

Der erste Schritt besteht in der geeigneten mathematischen Modellierung der zu untersuchenden Vorgänge bzw. bei Anwendungen eines fertigen Programmpakets, in der Auswahl des Modells, welches dem konkreten Problem am besten angepasst ist. Diesem Aspekt kommt eine entscheidende Bedeutung zu, da die Berechnung in der Regel keine nutzbringenden Ergebnisse bringen wird, wenn ihr kein adäquates Modell zu Grunde liegt.

1

10

100

1000

10000

100000

1970 1980 1990 2000

Numerische Verfahren

Rechengeschwindigkeit

13

Abbildung 2.2: Vorgehensweise bei der Anwendung numerischer Simulation zur Lösung von Ingenieurproblemen

Das aus der Modellierung resultierende kontinuierliche Problem, in der Regel im Rahmen der Kontinuums Mechanik abgeleitete Systeme von Differential- oder Integralgleichungen, muss anschließend durch ein geeignetes diskretes Problem approximiert werden, d. h. die zu berechnenden Größen müssen durch eine endliche Anzahl von Werten angenähert werden. Dieser Prozess wird als Diskretisierung bezeichnet und beinhaltet im Wesentlichen zwei Aufgaben:

• die Diskretisierung des Problemgebiets,

• die Diskretisierung der Gleichungen.

Die Diskretisierung des Problemgebiets approximiert das kontinuierliche Gebiet (in Raum und Zeit) durch eine endliche Anzahl von Teilgebieten, in denen dann numerische Werte der unbekannten Variablen bestimmt werden. Die Beziehungen zur Berechnung dieser Werte werden durch die Diskretisierung der Gleichungen gewonnen, welche die kontinuierlichen Systeme durch algebraische approximiert. Im Gegensatz zu einer analytischen Lösung, stellt eine numerische Lösung einen Satz von diskretisierten Problemgebiet zugeordneten Werten dar, aus dem der Verlauf der Variablen näherungsweise konstruiert werden kann. Als Diskretisierungsverfahren stehen im Wesentlichen drei unterschiedliche Ansätze zur Verfügung:

• die Finite – Differenzen – Methode (FD),

• die Finite – Volumen – Methode (FV),

• die Finite – Elemente – Methode (FE).

Der nächste Schritt im Ablauf der Simulation besteht in der Lösung der algebraischen Gleichungssysteme (die eigentliche Berechnung), wobei man es oftmals mit Gleichungen mit mehreren Millionen Unbekannten zu tun hat.

Die Berechnung liefert dann zunächst eine, meist riesige Menge an Zahlen, welche in der Regel nicht unmittelbar einer Interpretation zugänglich sind. Für die Auswertung der Berechnungsergebnisse ist daher insbesondere eine geeignete Visualisierung der Ergebnisse von Bedeutung, ohne die eine sinnvolle Aufbereitung der Datenflut meist nicht möglich ist. Nachdem Ergebnisse in interpretierbarer Form vorliegen, ist es notwendig diese hinsichtlich ihrer Qualität zu prüfen. Bei allen zuvor durchgeführten Schritten werden zwangsläufig Fehler gemacht, und man muss sich darüber Klarheit verschaffen, in welcher

TTeecchhnniisscchheess PPrroobblleemm PPrroobblleemmllöössuunngg ??????

DDiiffffeerreennttiiaallgglleeiicchhuunngg

MMaatthheemmaattiisscchhee MMooddeelllliieerruunngg

GGiitttteerreerrzzeeuugguunngg DDiisskkrreettiissiieerruunngg

BBiillddiinnffoorrmmaattiioonneenn

VViissuuaalliissiieerruunngg AAuusswweerrttuunngg

AAnnaallyyssee IInntteerrpprreettaattiioonn

VVaalliiddiieerruunngg

GGlleeiicchhuunnggssssyysstteemmee

NNuummeerriisscchhee LLöössuunngg AAllggoorriitthhmmeenn

14

Größenordnung diese sind (z. B. Referenzexperimente für Modellfehler, systematische Rechnungen für numerische Fehler). Oftmals ist es nach dieser Validierung notwendig, das Modell anzupassen oder die Rechnung mit einer verbesserten Diskretisierungsgenauigkeit erneut durchzuführen.

Zusammenfassend lässt sich festhalten, dass für die Entwicklung und Anwendung numerischer Berechnungsverfahren im Bereich der Ingenieurwissenschaften, insbesondere die folgenden Aufgabengebiete von besonderer Bedeutung sind:

• Mathematische Modellierung kontinuumsmechanischer Vorgänge,

• Entwicklung und Analyse numerischer Verfahren,

• Implementierung von Methoden in Computerprogramme,

• Anwendung der Berechnungsverfahren auf konkrete Problemstellungen,

• Auswertung und Interpretation der Berechnungsergebnisse.

Beispiel 2.1:

Eine typische, allgemeine Aufgabenstellung einer gewöhnlichen Differentialgleichung könnte lauten: das Anfangswertproblem der Ordnung n ,

( ) ( ) ( ) ( )( ) 0,,...,,, =′ tytytytFn , mit gegebenen Anfangswerten

( ) ( ) ( ) ( ) ( )1

00

1

0000 ,...,, −− =′=′= nnytyytyyty in einem gegebenen Definitionsgebiet ist

zu lösen (Existenz und Eindeutigkeit können als gesichert betrachtet werden.).

Beispiel 2.2:

Das im Rahmen dieser Vorlesung immer wieder verwendetes Demonstrationsbeispiel lautet: gegeben sei ein Stab der Länge L dessen Ortskoordinaten in x gemessen werden. Die Temperatur der beiden Stabenden sei über der Zeit t in der Form von

( ) ( ) ( ) ( )thtLtgt == ,,,0 φφ gegeben, und die Temperatur zum Zeitpunkt 0t durch

( ) ( ) Lxxfx <<= 0,0,φ bekannt. φ bezeichnet die Temperatur des Stabes, α die spezifische Wärme (welche als konstant angenommen wird), ρ die ebenfalls als konstant vorausgesetzte Massendichte und κ den konstanten Wärmeleitkoeffizienten. Das thermische Verhalten des Stabes sei durch die folgende Wärmeleitungsgleichung beschrieben:

2

2

2

2

xd

xt ∂

∂⋅=

∂

∂⋅

⋅=

∂

∂ φφ

αρ

κφ. (1)

2.2 Numerischer Algorithmus Allgemein formuliert, berechnet ein numerischer Algorithmus aus einem gegebenen Datensatz der Dimension m , ( )myyy ,...,, 21 unter Verwendung einer endlichen Folge von

Operationen f gesuchte Ergebnisse der Dimension n , ( )nxxx ~,...,~,~21 :

( ) ( )mn yyyfxxx ,...,,~,.,~,~2121 = . (2)

Im Allgemeinen approximiert die numerische Lösung ( )nxxx ~,.,~,~21 die exakte Lösung

( )nxxx ,.,, 21 , ist mit dieser aber nicht ident. Der Unterschied wird als Fehler bezeichnet und

kann unterschiedliche Ursachen haben:

• Fehler in Eingabedaten (z. B.: ungenaue Messdaten als Eingabedaten)

15

• Rundungsfehler (z. B. in Abhängigkeit der Maschinengenauigkeit des verwendeten Computers)

• Abbruchfehler (z. B. resultierend aus dem Abbruch eines iterativen Prozesses)

• Modellfehler (z. B. Fehler resultierend aus der Vernachlässigung von Größen im Zuge der Herleitung des mathematischen Modells)

2.3 Numerischer Fehler – Residuum – Kondition Grundsätzlich werden zwei Arten von numerischen Fehlern angegeben: der absolute Fehler und der relative Fehler. Der absolute numerische Fehler

absε berechnet sich aus der

Differenz zwischen der errechneten Lösung x~ und der exakten Lösung x

( ) xxxabs −= ~~ε . (3)

Der relative numerische Fehler errechnet sich indem der absolute Fehler durch die exakte Lösung x dividiert (skaliert) wird

( )x

xxxrel

−=

~~ε . (4)

Das Residuum kann durch Einsetzen der berechneten Lösung in die Ausgangsgleichung berechnet werden. Der numerische Fehler (absolut und relativ) kann im Allgemeinen nicht bestimmt werden, da die explizite Lösung nicht bekannt ist. Ob ein kleines Residuum einen kleinen Fehler impliziert, hängt sowohl vom Problem als auch vom verwendeten Verfahren ab.

Ein Problem heißt gut konditioniert, wenn kleine Änderungen in den Problemdaten zu kleinen Änderungen im Ergebnis führen. Kondition ist also eine Eigenschaft des Problems und nicht des Verfahrens. Eine genaue Definition der Kondition kann problemspezifisch vorgenommen werden.

Beispiel 2.3:

Durch iteratives Lösen des linearen Gleichungssystems 0bxA =+⋅ wird eine Lösung x~ errechnet. Der absolute Fehlervektor ist somit ( ) xxxεabs −= ~~ . Das

Residuum xr~ errechnet sich durch Einsetzen von x~ in das lineare Gleichungssystem

xrbxA ~~ =+⋅ .

Beispiel 2.4: EXCEL Beispiel 1

Gesucht ist die Lösung der Gleichung xex

−= . Diese wird mit dem Iterationsverfahren ix

i ex−

+ =1 für ,...3,2,1=i mit 11 =x approximiert. Die exakte Lösung ist 0,56714329.

Dem EXCEL Beispiel können die Verläufe des Residuums, des absoluten Fehlers und des relativen Fehlers entnommen werden. Die Gleichung x

ex−= kann als gut

konditioniert angesehen werden.


Gesucht ist die Lösung der Gleichung ( )15.0 2 +⋅= xx . Diese wird mit dem

Iterationsverfahren ( )15.0 2

1 +⋅=+ ii xx für ,...3,2,1=i mit 1.11 =x approximiert. Dem

EXCEL Beispiel ist der Verlauf des Residuums zu entnehmen. Das Problem ist ( )15.0 2 +⋅= xx ist schlecht konditioniert.

16

Beispiel 2.6:

Die Kondition eines linearen Gleichungssystems 0bxA =+⋅ wird durch die Konditionszahl

( )

∞

⋅=

sonst

regulärfalls A AAA

-1

κ (5)

der Matrix A beschrieben.

2.4 Lineare Gleichungssysteme Die Diskretisierung eines stationären Problems oder eines instationären Problems mit impliziter Zeitintegration, sei es mit der Finite – Volumen – oder der Finite – Elemente – Methode, führt auf große dünn besetzte Systeme algebraischer Gleichungen. Das Lösungsverfahren für diese Gleichungssysteme ist ein sehr wichtiger Bestandteil eines jeden numerischen Berechnungsverfahrens. Häufig wird mehr als 50% der Gesamtrechenzeit für die numerische Lösung dieser Systeme benötigt.

In diesem Kapitel werden Lösungsverfahren welche zur Lösung linearer Gleichungssysteme der Form

=

+

⋅

0

0

0

b

b

b

x

x

x

aa

aa

aaa

nnnnn

n

MMMOM

L

2

1

2

1

1

2221

11211

(6)

oder kurz

0bxA =+⋅ (7)

verwendet werden können. Es werden beispielhaft einige typische Lösungsverfahren behandelt, um charakteristische Eigenschaften solcher Methoden, insbesondere im Hinblick auf deren Effizienz, herauszuarbeiten. Es existiert eine Vielzahl weiterer, auf bestimmte Klassen von Systemen speziell zugeschnittene Methoden. In diesem Zusammenhang wird auf Spezialliteratur verwiesen.

2.4.1 Eigenschaften von linearen Gleichungssystemen

Jede lineare Abbildung zwischen endlich dimensionalen reellen Vektorräumen kann durch eine Matrix

nm

mnm

n

aa

aa

aaa

×ℜ∈

=

1

2221

11211

OM

L

A (8)

repräsentiert werden. Um die Menge aller Matrizen in unterschiedliche Klassen einteilen zu können, werden im Folgenden einige grundlegende Definitionen eingeführt.

17

2.4.1.1 Normen und Skalarprodukt In diesem Kapitel werden Matrix- und Vektornormen sowie das Skalarprodukt zusammengefasst ohne auf jeweilige elementare Eigenschaften näher Bezug zu nehmen. Die Matrixnorm einer beliebigen gegebenen Matrix nmA ×ℜ∈ ist eine Abbildung der Form

ℜ→ℜ= ×nm

...A . (9)

Beispiele für Matrixnormen sind die Zeilensummennorm

∑=

=∞=

n

j

ijmi

a1

,...,1maxA , (10)

die Spaltensummennorm

∑=

==

m

i

ijnj

a1

,...,11maxA (11)

und die Frobeniusnorm

∑∑= =

=m

i

n

j

ijFa

1 1

2A (12)

Ein weiteres, in der Praxis wesentliches Beispiel ist die Euklidische Matrixnorm

2

2

2sup

x

xAA

nx0

⋅=

ℜ∈≠

, (13)

welche durch die Euklidische Vektornorm für einen beliebigen Vektor nx ℜ∈ durch

∑=

=n

i

ix1

2

2x , (14)

induziert wird.

Üblicherweise wird auf dem nℜ das euklidische Skalarprodukt verwendet. Dies ist für zwei

beliebige Vektoren nx ℜ∈ und ny ℜ∈

( ) ∑=

⋅=n

i

ii yx1

yx, . (15)

2.4.1.2 Symmetrie, Orthogonalität und Definitheit von Matrizen Zu einer gegebene Matrix

nn

nnn

n

aa

aa

aaa

×ℜ∈

=

1

2221

11211

OM

L

A (16)

heißt

18

nn

nnn

n

aa

aa

aaa

×ℜ∈

=

1

2212

12111

OM

L

tA (17)

die zu A transponierte Matrix.

Die Matrix nnA ×ℜ∈ heißt symmetrisch, falls tAA = gilt.

Die Matrix nnA ×ℜ∈ heißt orthogonal, falls IAAAA t-1 =⋅=⋅ gilt.

Eine Matrix nnA ×ℜ∈ heißt positiv semidefinit, falls ( ) 0≥⋅ xx,A für alle nx ℜ∈ gilt.

Eine Matrix nnA ×ℜ∈ heißt positiv definit, falls ( ) 0>⋅ xx,A für alle 0x n \ℜ∈ gilt.

Eine Matrix nnA ×ℜ∈ heißt negativ semidefinit, falls A- positiv semidefinit ist.

Eine Matrix nnA ×ℜ∈ heißt negativ definit, falls A- positiv definit ist.

Beispiel 2.7:

Gegeben sei eine Matrix

=

987

654

321

A . Die zu A transponierte Matrix tA ist

=

963

852

741tA .

Beispiel 2.8:

Die Matrix

⋅=

21-0

1-21-

01-21

hA ist für ein beliebiges +ℜ∈h eine symmetrische und

positiv definite (alle Eigenwerte sind positiv) Matrix.

Beispiel 2.9:

Für jeden beliebigen Winkel [ [πα 2,0∈ ist die Matrix

=

αα

αα

cossin

sin - cosA

orthogonal.

2.4.1.3 Besetzungsstruktur

Eine Matrix nnA ×ℜ∈ heißt linke untere Dreiecksmatrix, falls ijaij >∀= 0 gilt.

Eine Matrix nnA ×ℜ∈ heißt rechte obere Dreiecksmatrix, falls ijaij <∀= 0 gilt.

Eine Matrix nnA ×ℜ∈ heißt Diagonalmatrix, falls ijaij ≠∀= 0 gilt.

Weist eine Matrix eine strukturierte Besetzung (z. B. Bandstruktur) auf, so wird die Matrix als strukturiert (sparse) bezeichnet. In diesem Zusammenhang sind speziell Matrizen mit Bandstruktur von Bedeutung. Dabei spielt der Begriff der Bandbreite eine wesentliche Rolle.

19

Die Bandbreite ist die Anzahl der unteren Nebendiagonalen + 1 + die Anzahl der oberen Nebendiagonalen. Ist eine Matrix voll besetzt, so wird ihre Besetzungsstruktur als dicht (dense) bezeichnet.

Beispiel 2.10:

Die Bandmatrix

−

−−

−−

−

=

2100

1210

0121

0012

A (18)

hat die Bandbreite 3.

2.4.2 Methoden zum Lösen linearer Gleichungssysteme

Grundsätzlich gibt es zwei Klassen von Verfahren, welche zum Lösen linearer Gleichungssysteme herangezogen werden können: Direkte Verfahren und Iterative Verfahren. Direkte Verfahren nützen in der Regel die besondere Gestalt von Gleichungssystemen nicht aus, wodurch vollbesetzte Zwischenmatrizen entstehen, die einerseits den verfügbaren Speicherplatz überschreiten und andererseits zu unakzeptablen Rechenzeiten führen. Des Weiteren entstehen solche Gleichungssysteme zumeist durch eine Diskretisierung der zugrunde liegenden Aufgabenstellung, wodurch auch die exakte Lösung des Gleichungssystems nur eine Approximation der gesuchten Lösung darstellt. Folglich erweist sich eine Näherungslösung für das Gleichungssystem mit einem Fehler in der Größenordnung des Diskretisierungsfehlers als ausreichend. Hierzu eignen sich iterative Methoden hervorragend.

2.4.2.1 Direkte Verfahren Unter einem direkten Verfahren zur Lösung eines linearen Gleichungssystems versteht man eine Rechenvorschrift, die unter Vernachlässigung von Rundungsfehlern die exakte Lösung in endlich vielen Schritten ermittelt. Direkte Verfahren werden nur selten zur unmittelbaren Lösung großer linearer Gleichungssysteme verwendet. Sie werden jedoch häufig in einer unvollständigen Form als Vorkonditionierer innerhalb iterativer Methoden genutzt und zur Lösung von Subproblemen eingesetzt. Weiters stellen direkte Verfahren bei einer kleinen Anzahl von Unbekannten oftmals eine effiziente Vorgehensweise dar. Im Rahmen dieser Vorlesung werden exemplarisch nur das Gaußsche Eliminationsverfahren sowie das Cholesky - Verfahren behandelt. Details zu weiteren Verfahren (z. B. QR – Zerlegung) können der einschlägigen Fachliteratur entnommen werden.

2.4.2.1.1 Gaußsches Eliminationsverfahren

Das Gaußsche Eliminationsverfahren formt die Matrix nnA ×ℜ∈ in ein äquivalentes Produkt

RLA ⋅= (19)

mit der unteren Dreiecksmatrix

=

− 1

0

1

001

1,1

21

nnn ll

l

L

OOM

MO

L

L (20)

20

und der oberen Dreiecksmatrix

=−

nn

nn

n

r

r

r

rrr

00

0

,1

22

11211

L

OOM

MO

L

R (21)

um. Dieser Vorgang wird als LR – Zerlegung (oder auch als LU – Zerlegung) bezeichnet. Unter Verwendung von L und R wird das lineare Gleichungssystem durch Vorwärtseinsetzen

0bcL =+⋅ (22)

und durch anschließendes Rückwärtseinsetzen

0cxR =−⋅ (23)

gelöst.

Zur Analyse des Rechenaufwandes werden lediglich die zeitaufwendigen Multiplikationen

und Divisionen berücksichtigt. Es ergibt sich eine Anzahl von nnnZGauß ⋅++⋅=3

2

3

1 23

wesentlichen Operationen.

Beispiel 2.11:

Beträgt die Rechenzeit für eine wesentliche Operation 1 Mikrosekunde, so benötigt das Gaußsche Eliminationsverfahren für eine reelle Matrix nnA ×ℜ∈ die in Tabelle 1 angeführten Rechenzeiten.

n Zeit

1000

10000

100000

ca. 4 min 30 sec

ca. 4 Tage

ca. 10 Jahre 7 Monate

Tabelle 1: Rechenzeiten des Gaußschen Eliminationsverfahrens


Das lineare Gleichungssystem

−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das Gaußsche

Eliminationsverfahren gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die Gauß Methode mit OK zu bestätigen. Das Residuum in der Größenordnung 4E-16 resultiert aus der Gleitkommaarithmetik.

2.4.2.1.2 Cholesky – Verfahren Für symmetrische, positiv definite Matrizen kann der beim Gaußschen Eliminationsverfahren benötigte Aufwand zur Berechnung einer LR – Zerlegung verringert werden. Dabei wird die Matrix nnA ×ℜ∈ in ein Produkt

21

tLLA ⋅= (24)

mit der unteren Dreiecksmatrix

=

nnnn lll

ll

l

L

OM

MO

L

21

2221

11

0

00

L (25)

zerlegt. Diese Zerlegung wird als Cholesky – Zerlegung bezeichnet. Ausgehend von dieser Zerlegung, wird die Lösung des linearen Gleichungssystems wiederum durch Vorwärtseinsetzen

0bcL =+⋅ (26)

und anschließendes Rückwärtseinsetzen

0cxLt =−⋅ (27)

bestimmt.

Zur Analyse des Rechenaufwandes werden wiederum die zeitaufwendigen Multiplikationen

und Divisionen berücksichtigt. Es ergibt sich eine Anzahl von nnnZCholesky ⋅+⋅+⋅=3

1

2

3

6

1 23

wesentlichen Operationen. Für große n benötigt die Cholesky – Zerlegung somit nur etwa die Hälfte der Operationen des Gaußschen Eliminationsverfahrens.



−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das Cholesky –

Verfahren gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die Cholesky – Methode mit OK zu bestätigen. Das Residuum in der Größenordnung 4E-16 resultiert aus der Gleitkommaarithmetik.

2.4.2.2 Iterative Verfahren In diesem Kapitel werden klassische Iterationsverfahren zur Lösung linearer Gleichungssysteme der Form

0bxA =+⋅ (28)

mit einer regulären, reellen und quadratischen Matrix nn×∈ℜA betrachtet.

Ein iteratives Verfahren ermittelt sukzessive Näherungen kx an die exakte Lösung bA- -1 ⋅ durch wiederholtes Ausführen einer festgelegten Rechenvorschrift

( )k1k xPx =+ (29)

für

0,1,2,...k = (30)

22

bei gewähltem Startvektor n∈ℜ0x .

2.4.2.2.1 Abbruchkriterien

Bei einem iterativen Verfahren muss man entscheiden, wann das Verfahren abgebrochen werden soll. Ideal wäre ein Abbruch, wenn der Fehler hinreichend klein ist. Da der numerische Fehler im Lauf des iterativen Verfahrens nicht bekannt ist, wird in der Praxis das Residuum zur Konvergenz- und Abbruchkontrolle verwendet. Im Allgemeinen steht man allerdings vor dem Problem, dass der Zusammenhang zwischen Residuum und Fehler unbekannt ist.

Ein mögliches Abbruchkriterium ist

( ) ε≤=

ini

xr~,...,1

max . (31)

Dabei wird der maximale Betrag der Einträge des Residuenvektors xr~ in Relation zu einem

vorgegebenen ε gesetzt. Eine weitere Variante summiert alle Beträge der Einträge des Residuenvektors xr~ vergleicht die Summe mit dem vorgegebenen ε .

ε≤∑=

n

ii

1

~xr (32)

Beide genannten Kriterien haben den Nachteil, dass sie jeweils das absolute Residuum betrachten. Eine mögliche Verbesserung wäre, die jeweils ausgewerteten Komponenten des Residuenvektors mit dessen betragsgrößter Komponente zu skalieren.

2.4.2.2.2 Jacobi – Verfahren Das Jacobi – Verfahren zur iterativen Lösung eines linearen Gleichungssystems setzt Nichtverschwindende Diagonalelemente der Matrix nnA ×ℜ∈ voraus und ist in Komponentenschreibweise für die Iteration k durch

nibxaa

xn

jij

i

k

jij

ii

k

i ,,2,11

1

1 K=

+⋅−= ∑≠=

+ (33)

gegeben. Beim Jacobi – Verfahren wird die neue Iterierte 1+k

ix ausschließlich mittels der

alten Iterierten k

ix ermittelt. Die Methode wird aus diesem Grund auch als

Gesamtschrittverfahren bezeichnet und ist folglich unabhängig von der gewählten Nummerierung der Unbekannten x .



−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das Jacobi

Verfahren gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die Jacobi Methode mit OK

zu bestätigen. Als Startwert wird der angegebene Vektor ( )t111 verwendet. 37

Iterationen werden benötigt, bis das eingestellte Abbruchkriterium erfüllt ist.

23

2.4.2.2.3 Gauß – Seidel – Verfahren Analog wie beim Jacobi – Verfahren wird auch beim Gauß – Seidel – Verfahren ein lineares Gleichungssystem mit einer regulären Matrix nnA ×ℜ∈ betrachtet, die ebenfalls keine Nullelemente auf der Hauptdiagonale aufweist. In Komponentenschreibweise lautet das Verfahren für die Iteration k

nibxaxaa

xn

ij

i

k

jij

i

j

k

jij

ii

k

i ,,2,11

1

1

1

11 K=

+⋅+⋅−= ∑∑

+=

−

=

++ . (34)

Aus dieser Darstellung wird deutlich, dass beim Gauß – Seidel – Verfahren zur Berechnung der i -ten Komponente der 1+k -ten Iterierten neben den Komponenten der alten k -ten

Iterierten kx die bereits bekannten ersten 1−i Komponenten der 1+k -ten Iterierten 1kx + verwendet werden. Das Verfahren wird daher auch als Einzelschrittverfahren bezeichnet und ist abhängig von der gewählten Nummerierung der Unbekannten x . Verglichen mit dem Jacobi – Verfahren kann beim Gauß – Seidel – Verfahren auf Grund des kleineren Spektralradius der Iterationsmatrix eine schnellere Konvergenz erwartet werden.



−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das Gauß –

Seidel Verfahren gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die Gauß – Seidel Methode mit OK zu bestätigen. Als Startwert wird der

angegebene Vektor ( )t111 verwendet. 11 Iterationen werden benötigt, bis das

eingestellte Abbruchkriterium erfüllt ist.

2.4.2.2.4 Relaxationsverfahren

Das Ziel von Relaxationsverfahren liegt in einer durch geeignete Gewichtung bedingte Verbesserung der Konvergenzgeschwindigkeit der Methode. Zur Gewichtung wird ein Relaxationsparameter ω mit +ℜ∈ω verwendet. 10 << ω dämpft zu Schwingungen neigenden Lösungsverlauf (Unterrelaxiation) während ω<1 einen langsam konvergierenden Lösungsverlauf beschleunigt (Überrelaxation). Der Ansatz von Relaxationsverfahren besteht in einer Gewichtung der Korrektur der i -ten Komponente der 1+k -ten Iterierten

nixxxk

i

k

i

k

i ,,2,111 K=∆⋅+= ++ ω . (35)

Dieser Relaxationsansatz ergibt im Fall des Gauß – Seidel – Verfahrens eine Iterationsvorschrift in Komponentenschreibweise

( ) nibxaxaa

xxn

ij

i

k

jij

i

j

k

jij

ii

k

i

k

i ,,2,111

1

1

11 K=

+⋅+⋅−⋅−= ∑∑

+=

−

=

++ ωω . (36)

Dieses Verfahren wird auch als sukzessives Überrelaxationsverfahren (SOR) bezeichnet. Das Verfahren konvergiert höchstens für Relaxationsparameter ( )2,0∈ω .


24


−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das SOR

Verfahren gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die SOR Methode mit OK

zu bestätigen. Als Startwert wird der angegebene Vektor ( )t111 verwendet. 8

Iterationen werden benötigt, bis das eingestellte Abbruchkriterium erfüllt ist.

2.4.2.2.5 Methode der konjugierten Gradienten Die Methode der konjugierten Gradienten (CG - Verfahren) kombiniert das Gradientenverfahren mit dem Verfahren der konjugierten Richtungen, Abbildung 2.3.

Abbildung 2.3: Herleitung des CG – Verfahrens

Das Verfahren nutzt die Residuenvektoren zur Definition der konjugierten Suchrichtungen, wodurch ein problemangepasstes Vorgehen in Bezug auf die Auswahl der Suchrichtungen vorliegt. Zudem wird die Optimalität des Verfahrens der konjugierten Richtungen erhalten.

Die Idee dieses Verfahrens ist, dass die Lösung eines symmetrischen, positiv definiten Systems

0bxA =+⋅ (37)

zugleich das Funktional

( ) vbvAv tt ⋅+⋅⋅⋅=2

1vF (38)

minimiert. Man versucht also das Minimum des quadratischen Funktionals F zu finden. Im Laufe des Verfahrens der konjugierten Gradienten wird dieses Minimum iterativ bestimmt, wobei man zur schrittweisen Verbesserung der aktuellen Lösung den Gradienten von F

( ) RbvA =+⋅=vFgrad (39)

verwendet. Es kann gezeigt werden, dass für ein System mit n Unbekannten spätestens nach n Iterationen die Lösung des Gleichungssystems mit hinreichender Genauigkeit berechnet werden kann.

Der Gradient einer Funktion weist immer in die Richtung der lokal stärksten Zunahme des Funktionswertes. Weil das Funktional zu minimieren ist, zeigt die Auswertung der Beziehung (39) daher den denkbar ungünstigsten Weg zur Verbesserung der aktuellen Lösung an.

VVeerrffaahhrreenn ddeess sstteeiillsstteenn AAbbssttiieeggss

BBaassiiss:: GGrraaddiieenntteenn aallss SSuucchhrriicchhttuunngg VVoorrtteeiill:: PPrroobblleemmoorriieennttiieerrtthheeiitt NNaacchhtteeiill:: KKoonnvveerrggeennzzvveerrhhaalltteenn

VVeerrffaahhrreenn ddeerr kkoonnjjuuggiieerrtteenn RRiicchhttuunnggeenn

BBaassiiss:: KKoonnjjuuggiieerrtthheeiitt ddeerr SSuucchhrriicchhttuunngg VVoorrtteeiill:: OOppttiimmaalliittäätt NNaacchhtteeiill:: FFeehhlleerrvveerrllaauuff

VVeerrffaahhrreenn ddeerr kkoonnjjuuggiieerrtteenn GGrraaddiieenntteenn

BBaassiiss:: GGrraaddiieenntteenn zzuurr BBeerreecchhnnuunngg kkoonnjjuuggiieerrtteerr SSuucchhrriicchhttuunngg VVoorrtteeiillee:: PPrroobblleemmoorriieennttiieerrtthheeiitt,, OOppttiimmaalliittäätt

25

Daher wird die entgegen gesetzte Richtung als Relaxationsrichtung kp verwendet. Zur

Bestimmung der Relaxationsrichtung verwendet man das zuletzt berechnete Residuum 1−kr

( )( )

1

22

111

,

, −

−−

−−− ⋅+−= k

kk

kkkk p

rr

rrrp . (40)

Dabei sind jeweils zwei aufeinander folgende Relaxationsrichtungen kp und 1+kp konjugiert. Die Schrittweite ergibt sich aus

( )( )kk

kk

kpAp

rr

⋅=

−−

,

,q

11

(41)

bestimmt. Der vollständige CG Algorithmus für symmetrische, positiv definite System ist in Abbildung 2.4 zusammengefasst. Eine Erweiterung des beschriebenen Algorithmus für unsymmetrische, nicht positiv definite Systeme ist möglich und wird in der entsprechenden Fachliteratur erörtert.

( )( )

( )( )

testKonvergenz

k

kk

k

k

kk

k

kk

k

k

k

k

kk

kk

kk

k

zrr

pxx

zp

rr

pAz

prp

rr

rr

rpb,xArx

⋅+=

⋅+=

=

⋅=

⋅+−=

=>

=

−=+⋅=

−

−

−−

−−

−

−−

−−

−

q

q

,

,q

e

,

,e

:1k falls

1,2,3,....k :Iteration

, von Wahl:Start

1

1

11

1

1

1

22

11

1

01000

(42)

Abbildung 2.4: Iterationsvorschrift des konjugierten Gradientenverfahrens



−

=

⋅

1

2

3

621

2114

145

3

2

1

x

x

x

soll durch das CG Verfahren

gelöst werden. Dazu ist das Makro SolveLinearEquationSystemDriver auszuführen und beim Abfragen der einzelnen Methoden die CG Methode mit OK zu bestätigen.

Als Startwert wird der angegebene Vektor ( )t111 verwendet. 3 Iterationen werden

benötigt, bis das eingestellte Abbruchkriterium erfüllt ist.

2.4.2.2.6 Moderne iterative Verfahren Die Entwicklung moderner numerischer Algorithmen zur Simulation praxisrelevanter Problemstellungen hat zu einer großen Nachfrage hinsichtlich effizienter, schneller und robuster iterativer Gleichungssystemlöser geführt. Derzeit stehen eine Vielzahl von unterschiedlichen Verfahren zur Verfügung. Beispiele für moderne Verfahren sind GMRES, GCS, BiCGSTAB, TFQMR und QMRCGSTAB. Details zu den einzelnen Verfahren sind der

26

Fachliteratur zu entnehmen. Die historische Entwicklung sowie die Zusammenhänge dieser Verfahren sind in Abbildung 2.5 dargestellt.

Abbildung 2.5: Zusammenhänge zwischen unterschiedlichen iterativen Gleichungssystemlösern

2.4.2.3 Vorkonditionierung Auf den Zusammenhang zwischen Fehler- und Residuenvektor, sowie auf den damit einher gehenden Begriff der Konditionszahl des Gleichungssystems wurde bereits eingegangen. So wurde auf den Vorteil einer kleinen Konditionszahl der Matrix nnA ×ℜ∈ des linearen Gleichungssystems 0bxA =+⋅ bereits verwiesen.

Motiviert durch diese Zusammenhänge wird in der Praxis versucht, das Gleichungssystem in ein äquivalentes System dessen Matrix eine niedrigere Konditionszahl hat, umzuformen. Solche Techniken werden als Vorkonditionierung (Präkonditionierung) des Systems bezeichnet und haben sich in der praktischen Anwendung als effizientes Mittel zur Beschleunigung und Stabilisierung der darauf angewandten Lösungsverfahren erwiesen. Neben der Verwendung eines geeigneten Gleichungssystemlösers ist die Wahl des Vorkonditionierers von entscheidender Bedeutung für das resultierende Gesamtverfahren. Beispiele für die unterschiedlichen Varianten sind Skalierungen, polynomiale Vorkonditionierer und splitting – assoziierte Vorkonditionierer.

Die Vorkonditionierung wird durch Multiplikation des zu lösenden Gleichungssystems 0bxA =+⋅ mit geeignet gewählten regulären Matrizen nn

LP ×ℜ∈ und nn

RP ×ℜ∈ erreicht. Dann heißt

PR

LPRL

xPx

0bPxPAP

⋅=

=⋅+⋅⋅⋅ (43)

das zum Gleichungssystem 0bxA =+⋅ gehörende vorkonditionierte System. Gilt IPL ≠ ,

so heißt LP Linksvorkonditionierer und das System linkspräkonditioniert. Ist IPR ≠ , so heißt

bAx = ( ) bAxxF −=bxAt =

CCGG –– VVeerrffaahhrreenn,, 11995522

BBiiCCGG –– VVeerrffaahhrreenn,, 11997755

CCGGSS –– VVeerrffaahhrreenn,, 11998899

BBiiCCGGSSTTAABB –– VVeerrffaahhrreenn,, 11999922

BBiiCCGGSSTTAABB((ll)) –– VVeerrffaahhrreenn,, 11999933

GGMMRREESS –– VVeerrffaahhrreenn,, 11998866

QQMMRR –– VVeerrffaahhrreenn,, 11999911

TTFFQQMMRR –– VVeerrffaahhrreenn,, 11999933

QQMMRRCCGGSSTTAABB –– VVeerrffaahhrreenn,, 11999944

SSiimmuullttaannee BBeettrraacchhttuunngg von:

MMiinniimmiieerruunngg:

SSppeeiicchheerrppllaattzzrreedduukkttiioonn dduurrcchh QQuuaassiimmiinniimmiieerruunngg

VVeerrmmeeiidduunngg vvoonn MMuullttiipplliikkaattiioonn mmiitt

OOsszziillllaattiioonnssmmiinniimmiieerruunngg dduurrcchh eeiinnddiimmeennssiioonnaallee RReessiidduueennmmiinniimmiieerruunngg

VVeerrmmeeiidduunngg vvoonn MMuullttiipplliikkaattiioonn mmiitt t

A

tA

ÜÜbbeerreeiinnssttiimmmmeennddee GGrruunnddiiddeeee EErrwweeiitteerruunngg aauuff ll--ddiimmeennssiioonnaallee RReessiidduueennmmiinniimmiieerruunngg

27

RP Rechtsvorkonditionierer und das System rechtspräkonditioniert. Ein links- und rechtspräkonditioniertes System heißt beidseitig präkonditioniert.

2.4.2.4 Kompakte Speicherung Praxisrelevante Problemstellungen führen häufig auf große schwach besetzte Gleichungssysteme. Die Speicherung derartiger Gleichungssysteme wird gewöhnlich erst durch die Vernachlässigung der Nullelemente der Matrix, die teilweise über 99% der Matrixkoeffizienten darstellen, ermöglicht. Aus diesem Grund wird im folgenden Kapitel kurz auf einige Aspekte des Speicherbedarfs in Abhängigkeit von der Besetzungsstruktur eingegangen.

Im Fall des Cholesky – Verfahrens ist eine nähere Betrachtung des Speicherbedarfs sinnvoll. Es zeigt sich, dass auf Grund der Struktur der unteren Dreiecksmatrix L die Speicherung

von lediglich ( )12

1+⋅⋅= nnS Werten ungleich Null erforderlich ist (anstatt von 2

nS = , wie

dies im Fall einer vollbesetzten Matrix der Fall ist). Zu diesem Zweck werden die Elemente der unteren Dreiecksmatrix L zeilenweise in einem Vektor

( )t

nnlllllll L333231222111=L (44)

gespeichert. Diese Speicherstrategie wird als kompakte Speicherung bezeichnet und wird im Fall von schwach besetzten Matrizen angewandt.

Beispiel 2.18:

Eine linke untere Dreiecksmatrix

=

54.003.0002.0

043.002.0

0032.0

0002

L (45)

kann kompakt in der Form

( )t54.003.0002.043.002.032.02=L (46)

gespeichert werden.

Eine weitere kompakte Speicherungsvariante findet im Fall von Bandmatrizen ihre Anwendung. Besitzt die Matrix nnA ×ℜ∈ Bandstruktur, z. B.

=

nna

aa

aaa

aa

00

00

00

3332

232221

1211

K

OOOM

O

MO

L

A , (47)

so kann diese ebenfalls in Vektorform

( )taaaaaa L233222122111=A (48)

gespeichert werden. In diesem Vektor werden die Bandwerte spaltenweise gespeichert. Um die Werte im Vektor den einzelnen Spalten der Matrix eindeutig zuordnen zu können, ist ein

28

ganzzahliger Zusatzvektor mit n Einträgen (entsprechend der Spaltenanzahl) erforderlich. In diesem Zusatzvektor werden die Startindizes der einzelnen Spalten im Vektor (48) gespeichert. Für das angeführte Beispiel hat der Zusatzvektor die Einträge

( )tL9631=Sky . (49)

Diese Speichervariante wird als Skyline – Speicherung bezeichnet.

Beispiel 2.19:

Eine Bandmatrix

−

−−

−−

−

=

2100

1210

0121

0012

A (50)

kann kompakt in der Form

( )t2112112112 −−−−−−=A (51)

mit dem Skyline – Vektor

( )t9631=Sky . (52)

gespeichert werden.

2.5 Diskretisierungsverfahren Der Einsatz digitaler Rechenautomaten (Computer) beschränkt sich auf algebraische Probleme. Um eine Differentialgleichung mit numerischen Methoden lösen zu können, müssen daher alle auftretenden Ableitungen (in der Zeit und im Ort) durch algebraische Ausdrücke approximiert werden. Solche Ausdrücke heißen endlich bzw. finit. Der dazugehörige Vorgang wird als Diskretisierung der Gleichung(en) bezeichnet. Weiters muss aber auch eine Approximation des kontinuierlichen Raums und/oder der Zeit durch eine endliche Anzahl von Teilgebieten vorgenommen werden. Dieser Vorgang heißt Diskretisierung des Problemgebiets.

2.5.1 Diskretisierung des Problemgebiets

Die Diskretisierung des Problemgebiets approximiert das kontinuierliche Gebiet (in Raum und Zeit) durch eine endliche Anzahl von Teilgebieten, in denen dann numerische Werte der unbekannten Variablen bestimmt werden. Diese erfolgt in der Regel durch die Definition einer geeigneten Gitterstruktur, welche das Gebiet überdeckt. In der Praxis ist diese Gittererzeugung für komplexe Problemgeometrien oftmals der aufwendigste Teil einer numerischen Untersuchung. Es gilt hierbei einerseits die Geometrie möglichst exakt zu modellieren und andererseits ein „gutes“ Gitter im Hinblick auf eine effiziente Berechnung erzeugen. Hierbei ist zu berücksichtigen, dass eine enge Wechselwirkung zwischen der Diskretisierung des Problemgebiets, der Diskretisierung der Gleichungen und dem verwendeten Lösungsverfahren besteht. Grundsätzlich kann der Zusammenhang zwischen dem numerischen Lösungsverfahren und der Gitterstruktur folgendermaßen charakterisiert werden: je regelmäßiger ein Gitter ist, desto effizienter sind die Lösungsalgorithmen für die Berechnung, desto unflexibler ist es aber im Hinblick auf die Modellierung komplexer Geometrien.

29

Beispiel 2.20:

Bei Verwendung der Finite – Elemente – Methode ist der „Aspect ratio“ ein wesentliches Beurteilungskriterium für die Qualität eines Berechnungsnetzes. In diesem Zusammenhang wird auf den FE – Teil der Vorlesung verwiesen.

2.5.1.1 Netztypen Der Typ eines Gitters steht in enger Wechselwirkung mit dem für die Berechnung verwendeten Diskretisierungs- und Lösungsverfahren. Es lassen sich eine Reihe verschiedener Unterscheidungsmerkmale der Gitter ausmachen, denen in diesem Zusammenhang eine Bedeutung zukommt. Neben der Unterscheidung nach der Topologie (Abbildung 2.6) der Gitterzellen (meist Dreiecke oder Vierecke bzw. Tetraeder, Hexaeder, Prismen, seltener allgemeine Polygone bzw. Polyeder), welche auch krummlinige Ränder aufweisen können, kann man bei numerischen Gittern zunächst generell die folgenden Typen unterscheiden:

• Randapproximierende Netze: Rand des Problemgebiets wird durch Gitterlinien approximiert.

• Kartesische Netze: Gebiet wird mit regulärem Gitter überdeckt, am Rand können irreguläre Gitterzellen auftreten, die einer speziellen Behandlung bedürfen.

• Überlappende Netze: Unterschiedliche Gebietsbereiche werden unabhängig voneinander diskretisiert. An den Schnittstellen der Bereiche werden Überlappungszonen zugelassen, die speziell behandelt werden müssen.

Beispiel 2.21:

Abbildung 2.6: Beispiele zwei- und dreidimensionaler Gitterelemente

Beispiel 2.22:

Abbildung 2.7: Beispiel eines randapproximierenden Netzes (links) eines kartesischen Netzes (Mitte) und eines überlappenden Netzes (rechts)

30

2.5.1.2 Netzstruktur Ein für die praktische Anwendung sehr wichtiges Unterscheidungsmerkmal numerischer Gitter ist die (logische) Anordnung der Gitterzellen. Generell lassen sich in der Praxis verwendete Gitter diesbezüglich in zwei Klassen unterteilen:

• Strukturierte Netze

• Unstrukturierte Netze

Strukturierte Netze sind durch eine regelmäßige Anordnung der Gitterpunkte charakterisiert. Diese Gitter können dabei durchaus krummlinig sein, logisch sind sie jedoch rechteckig (bzw. quaderförmig im dreidimensionalen Fall). Dies bedeutet, dass es Richtungen gibt, entlang derer die Anzahl der Gitterpunkte gleich ist. Die regelmäßige Anordnung hat zur Folge, dass die Nachbarschaftsbeziehungen zwischen den Gitterpunkten stets einem festen Muster folgen und damit programmtechnisch relativ einfach umgesetzt werden können. Es müssen nur die Koordinaten der Gitterpunkte gespeichert werden, die Identität der Nachbarpunkte, welche für die Diskretisierung benötigt werden, ist durch die Gitterstruktur eindeutig festgelegt.

Bei unstrukturierten Gittern existiert im Allgemeinen keinerlei Regelmäßigkeit in der Anordnung der Gitterpunkte. Die Möglichkeit die Gitterzellen unregelmäßig über das Problemgebiet zu verteilen, ermöglicht die größtmögliche Flexibilität zur Modellierung der Problemgeometrie, da die Gitterzellen optimal an die Ränder des Problemgebiets angepasst werden können. Neben den Koordinaten der Gitterpunkte, müssen für derartige Gitter auch die Nachbarschaftsbeziehungen zu benachbarten Gitterzellen mittels aufwendiger Datenstrukturen gespeichert werden.

Beispiel 2.23:

Abbildung 2.8: Beispiel eines strukturierten Netzes (links) und eines unstrukturierten Netzes (rechts)

In Tabelle 2 sind die wichtigsten Vor- und Nachteile strukturierter und unstrukturierter Gitter zusammengestellt. Eine einfache Schlussfolgerung für die praktische Anwendung hieraus ist, dass man zunächst versuchen sollte, das Gitter möglichst strukturiert zu wählen. Abweichungen von dieser Struktur sollten nur dann eingeführt werden, wenn dies aufgrund der Problemgeometrie erforderlich ist, um die Gitterqualität nicht zu schlecht werden zu lassen. Letzteres muss eventuell auch im Zusammenspiel mit einer aus Genauigkeitsgründen notwendigen lokalen Gitterverfeinerung gesehen werden, welche sich für unstrukturierte Gitter vergleichsweise einfach realisieren lässt (bei strukturierten Gittern ist man hier stark eingeschränkt, will man die Struktur beibehalten).

Oftmals wird die Gitterstruktur mit der Form der Gitterzellen in Verbindung gebracht. So werden Dreiecks- bzw. Tetraedergitter als unstrukturiert und Vierecks- bzw. Hexaedergitter als strukturiert bezeichnet. Die Form der Gitterzellen hat aber zunächst nichts mit der

31

Struktur des Gitters zu tun. Es lassen sich ohne weiteres strukturierte Dreiecksgitter und unstrukturierte Vierecksgitter realisieren. Auch die Art der Diskretisierung, welche für ein vorgegebenes Gitter zur Approximation der Gleichungen verwendet wird, wird häufig mit der Gitterstruktur in Verbindung gebracht: Finite – Differenzen – und Finite – Volumen – Verfahren mit strukturierten Gittern und Finite – Elemente – Verfahren mit unstrukturierten Gittern. Auch dies ist nicht gerechtfertigt. Jede der Methoden kann für strukturierte oder unstrukturierte Gitter formuliert werden, wobei natürlich die spezifischen Eigenschaften der Diskretisierungsmethoden die eine oder andere Vorgehensweise als vorteilhafter erscheinen lassen.

Tabelle 2: Übersicht über Vor- und Nachteile strukturierter und unstrukturierter Netze (jeweils relativ zueinander gesehen)

Es gibt auch Mischformen zwischen strukturierten und unstrukturierten Gittern, deren Verwendung sich für viele Problemstellungen als sehr vorteilhaft erweist, da es damit teilweise gelingt, die Vorteile beider Ansätze zu kombinieren. Die wichtigsten Varianten solcher Gitter sind:

• Blockstrukturierte Netze

• Hierarchisch strukturierte Netze.

Die Abbildung 2.9 zeigt Beispiele für diese Gittertypen. Blockstrukturierte Gitter sind lokal innerhalb eines jeden Blocks strukturiert, aber im Allgemeinen global unstrukturiert (irreguläre Blockanordnung). Sie stellen daher einen Kompromiss zwischen den geometrisch unflexiblen strukturierten Gittern und den numerisch aufwendigen unstrukturierten Gittern dar. Mit blockstrukturierten Gittern ist eine adäquate Modellierung komplexer Geometrien uneingeschränkt möglich. Innerhalb der einzelnen Blöcke können effiziente Lösungsverfahren eingesetzt werden, wobei jedoch besonderes Augenmerk auf die Behandlung der Blockübergänge gerichtet werden muss. Lokale Gitterverfeinerungen sind nur blockweise möglich, wobei in diesem Fall zusätzlich die Problematik so genannter hängender Knoten auftreten kann (d. h. keine stetigen Gitterlinien über die Blockgrenzen), welche im Lösungsverfahren besonderer Beachtung bedarf. Eine Blockstruktur bietet ferner eine natürliche Basis für die Parallelisierung eines Berechnungsverfahrens.

Hierarchisch strukturierte Gitter sind dadurch gekennzeichnet, dass ausgehend von einem (block-) strukturierten Gitter gewisse Bereiche des Problemgebiets wiederum in strukturierter Weise lokal verfeinert sind. Mit derartigen Gittern hat man nahezu alle Freiheiten für eine lokale adaptive Gitterverfeinerung, kann aber dennoch innerhalb der einzelnen Teilbereiche auf effiziente strukturierte Löser zurückgreifen. Solche Gitter stellen hinsichtlich dieser Aspekte einen Kompromiss zwischen strukturierten und unstrukturierten Gittern dar. Auch hier sind die Übergangszonen zwischen unterschiedlich fein diskretisierten Bereichen, die besonderer Beachtung bedürfen.

Beispiel 2.24:

Eigenschaft Strukturiert UnstrukturiertModellierung komplexer Geometrien - +Lokale Netzadaptivität - +Automatische Vernetzung - +CPU-Zeit für Vernetzung + -Zeitaufwand für Programmierung + -Datenmanagement + -Lösung des linearen Gleichungssystems + -Parallelisierung und Vektorisierung des Lösers + -

32

Abbildung 2.9: Beispiel eines blockstrukturierten Gitters (links) und eines hierarchisch strukturierten Gitters (rechts)

Für Finite – Elemente Berechnungen im Bereich der Strukturmechanik werden in der Regel unstrukturierte Gitter verwendet. Im Bereich der Strömungsmechanik, wo die Anzahl der Gitterpunkte in der Regel wesentlich größer ist (und damit die Effizienz des Lösers bedeutsamer wird), dominieren (block-) strukturierte Gitter.

2.5.1.3 Erzeugung strukturierter und unstrukturierter Gitter Die beiden am häufigsten verwendeten Methoden zur Erzeugung strukturierter Gitter sind die algebraische Gittererzeugung und die elliptische Gittererzeugung.

Bei der algebraischen Gittererzeugung werden die Punkte im Inneren eines Berechnungsgebiets durch eine Interpolationsvorschrift aus den Randpunkten ermittelt. Algebraische Methoden sind sehr einfach zu implementieren, erfordern wenig Rechenaufwand und besitzen für einfachere Geometrien eine ausreichende Flexibilität zur schnellen Erzeugung „vernünftiger“ Gitter. Für komplexere Problemgebiete sind diese Methoden jedoch weniger gut geeignet. Irregularitäten (z. B. Knicke) an den Rändern des Problemgebiets pflanzen sich ins Innere fort und die Kontrolle der Glattheit und der Verzerrung des Gitters ist vergleichsweise schwierig.

Eine alternative Vorgehensweise zur Erzeugung strukturierter Gitter bieten Techniken, die auf der Lösung von geeigneten partiellen Differentialgleichungen basieren. In diesem Zusammenhang ist die elliptische Gittererzeugung in der Praxis am weitesten verbreitet. Ein Vorteil dieser Methode ist, dass auch im Fall von Randirregularitäten glatte Gitter im Inneren resultieren. Auch die Metrik des Gitters muss in der Regel numerisch bestimmt werden, was aber aufgrund der Glattheit des Gitters im Allgemeinen unkritisch ist. Im Vergleich zu algebraischen Methoden ist sind diese Methoden allerdings aufgrund der Lösung einer partiellen Differentialgleichung aufwendiger.

Eine wesentliche Motivation zur Erzeugung unstrukturierter Gitter ist, damit den Gittergenerierungsprozess möglichst zu automatisieren. Der Idealfall wäre, dass, ausgehend von der Beschreibung der Problemgeometrie durch Randkurven bzw. –flächen, ohne weiteren Eingriff des Anwenders ein „brauchbares“ Gitter erzeugt wird. Da man diesem Idealziel mit Dreiecken bzw. Tetraedern am nächsten kommt, werden meist diese Zellentypen für unstrukturierte Gitter verwendet. Ein gängiges Beispiel einer Methode zur Erzeugung unstrukturierter Gitter ist die Delaunay – Triangulierung.

2.5.2 Diskretisierung von Ortsvariablen

Im Rahmen dieses Kapitels werden die Verfahren und deren Anwendung zur Durchführung der Ortsdiskretisierung im Überblick dargestellt. Nähere Informationen können der Fachliteratur entnommen werden.

33

2.5.2.1 Finite – Differenzen – Methode Bei der Finite – Differenzen – Methode werden Differentialausdrücke durch eine entsprechende Differenzenformel ersetzt. Die entstehenden Gleichungen, welche die Lösungen an diskreten räumlichen und zeitlichen Punkten verknüpfen, werden Differenzengleichungen genannt. Die Bildung der Differenzengleichung basiert auf der Definition der Ableitung

( ) ( ) ( ) ( )

x

xfxxf

x

xfxxf

x

f

xx ∆

−∆+≈

∆

−∆+=

∂

∂→∆ 0

lim . (53)

Der Fehler, der dabei gemacht wird, hängt unmittelbar von x∆ ab, wird also von der Feinheit des Gitters bestimmt.

Für kartesische äquidistante Gitter lassen sich die Differenzenformeln mittels Taylorreichenentwicklung herleiten. Die Taylorreihe um einen beliebigen Punkt x mit der Schrittweite hx =∆ für eine Funktion f lautet

( ) ( ) ...!3!2!1 3

33

2

22

+∂

∂+

∂

∂+

∂

∂+=+

xxx x

fh

x

fh

x

fhxfhxf (54)

Für 0>h ist (54) eine Vorwärtsentwicklung, für 0<h ist (54) eine Rückwärtsentwicklung. Zur Approximation der ersten Ableitung durch den vorwärtigen Differenzenquotionen wird die Reihe nach dem zweiten Term abgebrochen. Es ergibt sich die finite Differenzenformulierung für die erste Ableitung mit einem Abbruchfehler der Ordnung ( )xO ∆ :

( ) ( ) ( )xO

x

xfxxf

x

f

x

∆+∆

−∆+=

∂

∂. (55)

Die Differenzengleichung der ersten Ableitung einer Funktion φ im Punkt ix ist

( )xOxx

ii

xi

∆+∆

−≈

∂

∂ + φφφ 1 . (56)

Beispiel 2.25:

Im Fall der Verwendung eines rückwärtigen Differenzenquotionen resultiert die Differenzengleichung

( )xOxx

ii

xi

∆+∆

−≈

∂

∂ −1φφφ. (57)

Beispiel 2.26:

Die Kombination des vorwärtigen und des rückwärtigen Differenzenquotionen liefert den zentralen Differenzenquotionen. Die dazugehörige Differenzengleichung lautet

( )xOxx

ii

xi

∆+∆⋅

−≈

∂

∂ −+

2

11 φφφ. (58)

Beispiel 2.27:

34

Die Differenzengleichung mit einem Abbruchfehler der Ordnung ( )3xO ∆ zur

Approximation der ersten Ableitung einer Funktion φ im Punkt ix lautet

( )3321

6

291811xO

xx

iiii

xi

∆+∆⋅

⋅−⋅+⋅−⋅≈

∂

∂ −−− φφφφφ. (59)

Die Herleitung der Differenzengleichungen für höhere Ableitungen müssen die Taylorreihen um mehrere Punkte entwickelt und kombiniert werden.

Beispiel 2.28:

Die zentrale Differenzengleichung zur Approximation der zweiten Ableitung einer Funktion φ im Punkt

ix lautet

( )2

2

11

2

2 2xO

xx

iii

xi

∆+∆

+⋅−≈

∂

∂ −+ φφφφ. (60)

Die Formel hat einen Abbruchfehler der Ordnung ( )2xO ∆ .

2.5.2.2 Finite – Volumen – Methode Finite – Volumen – Methoden (FV) werden hauptsächlich für die numerische Lösung von Problemen der Strömungsmechanik eingesetzt. Die Anwendungsmöglichkeiten der Methode sind jedoch keineswegs auf Problemstellungen aus diesem Bereich beschränkt. Eine der wichtigsten Eigenschaften von Finite – Volumen – Methoden ist, dass die Erhaltungsprinzipien, die den mathematischen Modellen kontinuumsmechanischer Problemstellungen zu Grunde liegen, per Definition auch für die diskretisierten Gleichungen erfüllt sind.

In diesem Kapitel werden lediglich die wichtigsten Grundlagen der Finite – Volumen Diskretisierung zur Anwendung auf kontinuumsmechanische Problemstellungen erörtert. Zur klaren Darstellung der wesentlichen Prinzipien wird nur der zweidimensionale Fall diskutiert.

Ausgangspunkt für eine Finite – Volumen Diskretisierung ist eine Zerlegung des Problemgebiets in eine endliche Anzahl von finiten Volumen (Kontrollvolumen), deren Vereinigung das ganze Problemgebiet überdeckt (Abbildung 2.10).

Weiters muss bei den zu berechnenden unbekannten Variablen festgelegt werden, an welcher Stelle diese ausgewertet werden müssen. Man kann hierfür beispielsweise eckenorientiert oder zellenorientiert vorgehen. Abbildung 2.11 zeigt ein zweidimensionales Volumen mit einer zellenorientierten Anordnung der Variablen.

Für jedes dieser Kontrollvolumen werden die Erhaltungsgleichungen in Integralform formuliert. Diese stehen normalerweise direkte aus den entsprechenden kontinuumsmechanischen Erhaltungssätzen zur Verfügung, können aber auch durch Integration der entsprechenden Differentialgleichung gewonnen werden. Die resultierende integrale Bilanzgleichung bildet den Ausgangspunkt für die weitere Diskretisierung des betrachteten Problems mit dem Finite – Volumen – Verfahren. Die eigentliche Diskretisierung besteht im Wesentlichen darin, die Integrale der Bilanzgleichung mit geeigneten Mittelwerten der jeweiligen Integranden an den Seiten des jeweiligen Kontrollvolumens zu approximieren und diese in Beziehung zu den unbekannten Funktionswerten in den Volumenmittelpunkten zu setzen.

Die für die einzelnen Volumen resultierenden Bilanzgleichungen werden zu einem Gesamtgleichungssystem assembliert, welches in weiterer Folge unter Verwendung geeigneter Methoden zu lösen ist.

35

Abbildung 2.10: Diskretisierung eines ebenen Gebietes Ω (beschränkt durch Γ ) mittels rechteckiger finiter Volumina

Abbildung 2.11: Zweidimensionales Beispiel: Zellenorientierte Anordnung der Variablen für das Finite – Volumen – Verfahren

2.5.2.3 Multigrid – Verfahren (Mehrgitter – Verfahren) Die Idee der Mehrgitter – Verfahren basiert auf der Tatsache, dass ein iterativer Lösungsalgorithmus gerade solche Fehlerkomponenten einer Näherungslösung sehr effizient eliminiert, deren Wellenlängen der Gittermaschenweite entsprechen. Die langwelligen Fehler können hingegen mit einer solchen Methode nur sehr langsam abgebaut werden. Der Grund hierfür ist, dass mittels der Diskretisierungsvorschrift für jeden Gitterpunkt jeweils nur lokale Nachbarschaftsbeziehungen aufgestellt werden, was zur Folge hat, dass der globale Informationsaustausch (z. B. die Fortpflanzung von Randwerten in das Innere des Lösungsgebiets) bei Iterationsmethoden nur sehr langsam erfolgt. Zur Erläuterung dieses Sachverhalts dient folgendes Beispiel.

Beispiel 2.29:

ΓΩ

36

Das eindimensionale Diffusionsproblem 02

2

=∂

∂

x

φ für 10 << x und ( ) ( ) 010 == φφ soll

unter Verwendung einer zentralen Differenzendiskretisierung 2. Ordnung auf einem äquidistanten Gitter mit 1−N inneren Gitterpunkten gelöst werden. Die diskreten Gleichungen lauten 02 11 =+⋅− −+ iii φφφ für 1,...,1 −= Ni . Die Iterationsvorschrift für

das Jacobi – Verfahren zur Lösung dieses Gleichungssystems lautet ( ,...1,0=k ):

2

111k

i

k

ik

i−++ +

=φφ

φ für 1,...,1 −= Ni und 011

0 == ++ k

N

k φφ . (61)

Angenommen es sind die im oberen Teil von Abbildung 2.12 gezeigten Näherungslösungen vor. Man erkennt das unterschiedliche Fehlerreduktionsverhalten für die verschiedenen Startwerte:

• Im linken Fall ist der iterative Algorithmus sehr effizient. Bereits nach einer Iteration hat man die richtige Lösung.

• Im rechten Fall ist die Verbesserung gering. Es werden viele Informationen benötigt, um die exakte Lösung zu erhalten.

Abbildung 2.12: Fehlerreduktion bei Anwendung des Jacobi – Verfahrens zur Lösung des eindimensionalen Diffusionsproblems für

unterschiedliche Startwerte

Abbildung 2.13: Zusammenhang zwischen Informationsaustausch und Gittergröße auf Grund von Nachbarschaftsbeziehungen von

Diskretisierungsvorschriften

Im Normalfall sind in einem Startwert verschiedene Fehlerkomponenten überlagert vorhanden. Die hochfrequenten Anteile werden sehr schnell, die niederfrequenten Anteile sehr langsam reduziert. Die Idee eines Mehrgitter – Verfahrens besteht darin, die niederfrequenten Komponenten durch einen Iterationsprozess auf einer Hierarchie von sukzessive vergröberten Netzen zu reduzieren. Ein Mehrgitteralgorithmus verlagert die Berechnung nach einigen Iterationen (oft Glättungsiterationen genannt, da die Fehlerfunktion danach glatt, d. h. frei von hochfrequenten Anteilen, ist) von einem feinen auf ein grobes Gitter, in welchem z. B. nur jeder zweite Gitterpunkt auftritt. Glatte Funktionen können hierbei

37

ohne wesentlichen Informationsverlust auf gröberen Gittern dargestellt werden. Durch die Übertragung des Lösungsprozesses auf gröbere Gitter sehen die Fehler relativ zur Gitterweite dort hochfrequenter aus und man erreicht damit eine effiziente Reduktion der auf dem feinen Gitter niederfrequenten Fehlerkomponenten. Anschaulich lässt sich dies durch den schnelleren Informationsaustausch auf gröberen Gittern deuten (Abbildung 2.17). Die Effizienz der Mehrgittermethode liegt darin begründet, dass man einerseits eine deutlich effizientere Fehlerreduktion erreicht und andererseits der zusätzliche Aufwand für die Berechnungen auf gröberen Gittern auf Grund der geringeren Anzahl von Gitterpunkten relativ gering ist.

Beispiel 2.30:

Algebraische Multigrid – Verfahren (AMG) sind moderne Verfahren, welche beispielsweise in der Fluiddynamik im Einsatz sind. Der Verfeinerungsschritt (Prolongation) sowie der Vergröberungsschritt (Restriktion) sind bei AMG - Verfahren charakteristisch. Dabei wird auf ein Netz mit höherem Feinheitsgrad (Prolongation) oder mit niederem Feinheitsgrad (Restriktion) übergegangen und die zu lösenden Gleichungen werden darauf jeweils berechnet. Dies ist in folgender Abbildung dargestellt.

Abbildung 2.14: Schematische Darstellung der Restriktion / Prolongation bei AMG Verfahren

2.5.3 Diskretisierung von zeitlichen Variablen (Zeitintegration)

Bei vielen praktischen Anwendungen sind die zu untersuchenden Vorgänge instationär, so dass für deren numerische Simulation die Lösung zeitabhängiger Modellgleichungen notwendig ist. Die Zeit nimmt in den Differentialgleichungen eine gewisse Sonderstellung ein, da, anders als bei den Ortskoordinaten, aufgrund des Kausalitätsprinzips eine ausgezeichnete Richtung existiert. Dieser Tatsache müssen auch die Diskretisierungstechniken für die Zeit Rechnung tragen. In diesem Kapitel werden die wichtigsten Techniken zur Zeitdiskretisierung behandelt.

38

2.5.3.1 Grundlagen Bei instationären Vorgängen sind die physikalischen Größen, zusätzlich zur Ortsabhängigkeit, auch von der Zeit t abhängig. Typische zeitlich abhängige Probleme sind Transportvorgänge und Schwingungsvorgänge. Während die Gleichungen für instationäre Transportvorgänge nur die 1. Ableitung nach der Zeit enthalten, tritt bei Schwingungsvorgängen auch die 2. Zeitableitung auf. Im ersten Fall nennt man das Problem parabolisch, im zweiten Fall hyperbolisch.

Um zeitabhängige Probleme vollständig zu definieren, benötigt man, zusätzlich zu den Randbedingungen (die ebenfalls zeitabhängig sein können), auch Anfangsbedingungen. Für Transportvorgänge benötigt man nur eine Anfangsverteilung der unbekannten Funktion, während bei Schwingungsproblemen zusätzlich auch noch eine Anfangsverteilung für die erste Zeitableitung vorgegeben werden muss.

Zur numerischen Lösung zeitabhängiger Probleme wird üblicherweise zunächst eine Ortsdiskretisierung mit einer der in den vorangegangenen Kapiteln beschriebenen Techniken durchgeführt. Dies resultiert dann in einem System von gewöhnlichen Differentialgleichungen (in der Zeit). Im Folgenden wird zur Vereinfachung der Schreibweise die rechte Seite der aus der Ortsdiskretisierung resultierenden Gleichung (egal, ob diese mit einer Finite Differenzen - Methode, einer Finite Volumen – Methode, einer Finite Elemente Methode oder einer anderen Methode erhalten wird) durch den Operator L ausgedrückt:

( )nL

tφ

φ=

∂

∂. (62)

( )tφφ = bezeichnet den Vektor der unbekannten Funktion. Zur Zeitdiskretisierung, d. h. zur Diskretisierung des Systems gewöhnlicher Differentialgleichungen, können analoge Techniken, wie für die Ortskoordinaten benützt werden (Finite Differenzen, Finite Volumen, Finite Elemente). Im Rahmen dieser Vorlesung wird die Betrachtung auf die am häufigsten verwendeten Finite Differenzen Ansätze beschränkt.

Zunächst wird das zu betrachtende Zeitintervall [ ]Tt ,0 in einzelne, im Allgemeinen nicht –

äquidistante Teilintervalle nt∆ zerlegt:

,...2,1,01 =∆+=+ nttt nnn. (63)

Zur weiteren Vereinfachung der Schreibweise wird ein Variablenwert zum Zeitpunkt nt mit

einem Index n bezeichnet, z. B.:

( )( ) ( )n

n LtL φφ = . (64)

Nach dem Kausalitätsprinzip kann die Lösung zur Zeit 1+nt nur von davorliegenden

Zeitpunkten nt , 1−nt , … abhängen. Da die Zeit damit sozusagen eine „Einbahnkoordinate“ ist,

muss die Lösung für 1+nt als Funktion der Randbedingungen und der Lösungen zu früheren

Zeiten bestimmt werden. Bei der Zeitdiskretisierung handelt es sich daher um eine Extrapolation. Die gesuchte Variable wird, startend von den vorgegebenen Anfangsbedingungen bei 0t , sukzessive zu den Zeitpunkten 1t , 2t , … berechnet. Der

zeitliche Verlauf von φ kann dann als Folge verschiedener räumlicher Werte zu den diskreten Zeitpunkten dargestellt werden (siehe Abbildung 2.15 für ein räumlich zweidimensionales Problem). Es sei an dieser Stelle bemerkt, dass für Transportprobleme häufig auch stationäre Lösungen mit einem Zeitdiskretisierungsverfahren aus den zeitabhängigen Gleichungen als Grenzwert für ∞→t berechnet werden, was jedoch in der Regel ein wenig effizientes Verfahren liefert.

39

Abbildung 2.15: Zusammenhang zwischen räumlicher und zeitlicher Diskretisierung

Es ist immer mindestens eine bereits bekannte Zeitebene notwendig, um die Zeitableitung zu diskretisieren. Benützt man nur diese eine bekannte Zeitebene, d. h. die Werte zum Zeitpunkt nt , spricht man von Einschrittverfahren, werden mehrere bekannte Zeitebenen

verwendet, d. h. Werte zu den Zeitpunkten nt , 1−nt , … von Mehrschrittverfahren. Weiters

werden die Methoden zur Zeitdiskretisierung, je nach Wahl der Zeitpunkte, an denen die rechte Seite ausgewertet wird, generell in zwei Klassen eingeteilt:

• Explizite Verfahren: Diskretisierung der rechten Seite nur zu vorherigen (bereits bekannten) Zeitpunkten ( ),..., 11 −+ = nnn

F φφφ ,

• Implizite Verfahren: Diskretisierung der rechten Seite auch zum neuen (unbekannten) Zeitpunkt ( ),...,, 111 −++ = nnnn

F φφφφ .

Mit F ist jeweils eine beliebige Diskretisierungsvorschrift bezeichnet, für deren Wahl im Folgenden noch Beispiele angegeben werden. Die Unterscheidung in explizite und implizite Verfahren ist ein sehr wichtiges Merkmal, da weit reichende Unterschiede bezüglich der Eigenschaften des numerischen Verfahrens auftreten. Einzelheiten der beiden Verfahrensklassen werden in den nächsten Abschnitten behandelt.

2.5.3.2 Explizite Einschrittverfahren Das einfachste Beispiel eines expliziten Zeitdiskretisierungsverfahrens ist das explizite Euler Verfahren. Dieses erhält man durch eine Approximation der Zeitableitung zum Zeitpunkt

nt

mittels einer Vorwärtsdifferenzenformel

( ) ( )n

n

nn

n Ltt

tφ

φφφ=

∆

−≈

∂

∂ +1

. (65)

Dies entspricht der Approximation der Zeitableitung der Komponenten iφ von φ an der Stelle

nt mittels der Steigung der Geraden durch die Punkte n

iφ und 1+n

iφ . Das Verfahren besitzt

einen Abbruchfehler 1. Ordnung (bezüglich der Zeit) und ist auch unter dem Namen Eulersche Polygonzugmethode bekannt. Gleichung (65) kann explizit nach 1+nφ aufgelöst werden:

2φ

1φ

0φ

40

( )n

n

nnLt φφφ ⋅∆+=+1 . (66)

Auf der rechten Seite stehen nur Werte der bereits bekannten Zeitebene, so dass die Gleichungen für die Werte zum Zeitpunkt 1+nt an den verschiedenen Gitterpunkten

vollständig entkoppelt sind und unabhängig voneinander berechnet werden können. Dies ist charakteristisch für explizite Verfahren.

Beispiel 2.31:

Gegeben sei die partielle Differentialgleichung

( )43421

φ

φ

αρ

κφ

L

xt 2

2

∂

∂⋅

⋅=

∂

∂ mit den

Randbedingungen

( ) ( ) ( ) ( )thtLtgt == ,,,0 φφ und den Anfangswerten ( ) ( ) Lxxfx <<= 0,0,φ . Zur numerischen Lösung der Gleichung wird diese mittels der Methode der finiten Differenzen diskretisiert. Die Diskretisierung des Ableitungsterms auf der rechten

Seite 2

11

2

2 2

xx

iii

∆

+⋅−≈

∂

∂ −+ φφφφ wird mit einem zentralen Ansatz 2. Ordnung

durchgeführt. Zur Diskretisierung der Zeitableitung auf der linken Seite wird der explizite Euleransatz (65) angewandt. Die resultierende Differenzengleichung lautet

2

1

112

xt

nnnn

i

n

i iii

∆

+⋅−⋅

⋅=

∆

−−+

+ φφφ

αρ

κφφ.


Die Auswertung von 2

1

112

xd

t

nnnn

i

n

i iii

∆

+⋅−⋅=

∆

−−+

+ φφφφφ mit gegebenem

αρ

κ

⋅=d sowie

gegebenen Anfangs- und Randwerten wird in dem EXCEL Beispiel veranschaulicht. Dabei ist die Ortsdiskretisierung horizontal aufgetragen, während die Zeitdiskretisierung vertikal verläuft.

In Abbildung 2.16 ist die Vorgehensweise bei der Anwendung des expliziten Euler – Verfahrens auf die eindimensionale Diffusionsgleichung veranschaulicht. Man erkennt, dass eine Reihe von Zeitschritten benötigt wird, bis sich eine Änderung der Randbedingungen im Inneren des Rechengebietes bemerkbar macht. Je feiner das Gitter ist, desto langsamer ist die Ausbreitung der Information (bei gleichem Zeitschritt). Dieser Umstand führt zu einer Begrenzung der Zeitschrittweite (Stabilitätsbedingung), die quadratisch von der räumlichen Auflösung abhängt und mit feiner werdendem Ortsgitter immer restriktiver wird. Diese Begrenzung ist ausschließlich numerisch bedingt und völlig unabhängig vom tatsächlichen zeitlichen Verlauf der Lösung.

41

Abbildung 2.16: Vorgehensweise und Informationsfluss bei Anwendung des expliziten Euler - Verfahrens auf die eindimensionale Diffusionsgleichung

Es gibt eine Vielzahl weiterer expliziter Einschrittverfahren, die sich vom expliziten Euler – Verfahren in der Approximation der Flüsse und rechten Seiten unterscheiden (z. B.: Runge – Kutta - Verfahren). Diese werden an dieser Stelle nicht näher behandelt.

2.5.3.3 Implizite Einschrittverfahren

Approximiert man die Zeitableitung zum Zeitpunkt 1+nt durch eine Rückwärtsdifferenz 1.

Ordnung, so ergibt sich das implizite Euler – Verfahren:

( ) ( )1

1

1 ++

+ =∆

−≈

∂

∂ n

n

nn

n Ltt

tφ

φφφ. (67)

Diese Beziehung unterscheidet sich vom expliziten Euler – Verfahren nur durch die Auswertung der rechten Seite, die in diesem Fall für die neue (unbekannte) Zeitebene berechnet wird. Eine Konsequenz ist, dass damit eine explizite Auflösung nach 1+n

iφ nicht

mehr möglich ist, da alle Variablen der neuen Zeitebene miteinander gekoppelt sind. Diese Eigenschaft ist charakteristisch für implizite Verfahren. Zur Berechnung jeder neuen Zeitebene ist, wie im stationären Fall, die Lösung eines linearen Gleichungssystems erforderlich. Im impliziten Fall wirken sich Änderungen der Randbedingungen bereits im gleichen Zeitschritt im ganzen Lösungsgebiet aus, so dass die beim explizite Euler – Verfahren auftretenden Stabilitätsprobleme nicht auftreten. Unabhängig vom Feinheitsgrad der Ortsdiskretisierung und von t∆ erweist sich das Verfahren als stabil. Das implizite Euler – Verfahren hat, ebenso wie das explizite, nur einen Abbruchfehler 1. Ordnung in der Zeit. Es ist aufwendiger als die explizite Methode, da mehr Speicher (für die Koeffizienten der Quellterme) und mehr Rechenzeit (zum Lösen des linearen Gleichungssystems) pro Zeitschritt benötigt werden. Es liegt jedoch keine stabilitätsbedingte Begrenzung für die Größe des Zeitschritts vor. Der höhere Aufwand des Verfahrens wird in der Regel durch die Möglichkeit, sehr große Zeitschritte zu wählen, mehr als kompensiert, so dass das Verfahren in den meisten Fällen insgesamt wesentlich effizienter ist als die explizite Variante.

Beispiel 2.33:

Gegeben sei die partielle Differentialgleichung

( )43421

φ

φ

αρ

κφ

L

xt 2

2

∂

∂⋅

⋅=

∂

∂ mit den

Randbedingungen

( ) ( ) ( ) ( )thtLtgt == ,,,0 φφ und den Anfangswerten ( ) ( ) Lxxfx <<= 0,0,φ . Zur numerischen Lösung der Gleichung wird diese mittels der Methode der finiten Differenzen diskretisiert. Die Diskretisierung des Ableitungsterms auf der rechten

Seite 2

11

2

2 2

xx

iii

∆

+⋅−≈

∂

∂ −+ φφφφ wird mit einem zentralen Ansatz 2. Ordnung

durchgeführt. Zur Diskretisierung der Zeitableitung auf der linken Seite wird der implizite Euleransatz (67) angewandt. Die resultierende Differenzengleichung lautet

2

1111

112

xt

nnnn

i

n

i iii

∆

+⋅−⋅

⋅=

∆

−++++

−+φφφ

αρ

κφφ.


42

Die implizite Differenzengleichung 2

1111

112

xt

nnnn

i

n

i iii

∆

+⋅−⋅

⋅=

∆

−++++

−+φφφ

αρ

κφφ soll für das

gesamte Berechnungsgebiet bei gegebenen Anfangs- und Randwerten zeitintegriert werden. Dabei ist die Ortsdiskretisierung horizontal aufgetragen, während die Zeitdiskretisierung vertikal verläuft. Dazu wird die Differenzengleichung zu

( ) 02 21

1

121

1 =⋅∆+⋅⋅∆+⋅⋅∆⋅+∆−⋅⋅∆ ++

++− 4342144444444 344444444 21

ii b

n

i

A

n

i

n

i

n

i xctdtxdt φφφφ umgeformt. Das

resultierende lineare Gleichungssystem wird mit dem Cholesky – Verfahren gelöst. Dazu ist das Makro EulerImplicitDriver auszuführen.

Weitere implizite Einschrittverfahren (z. B.: das Crank – Nicolson – Verfahren) werden im Rahmen dieser Vorlesung nicht näher behandelt.

2.5.3.4 Kombinationen mehrerer Einschrittverfahren Die Grundidee von Prädiktor – Korrektor – Verfahren ist eine kombinierte Anwendung von expliziten und impliziten Formeln. Mit dem expliziten Verfahren (dem sogenannten Prädiktor

– Verfahren) wird in einem ersten Schritt ein guter Startwert 1~ +nφ ermittelt. Dies ist bei Verwendung eines expliziten Einschrittverfahrens

( )nnnLt φφφ ⋅∆+=+1~

. (68)

Dieser Startwert wird in einem zweiten Schritt zur Berechnung der neuen Lösung 1+nφ unter Verwendung eines impliziten Verfahrens (hier ein Einschrittverfahren) angewendet

( )11 ~ ++ ⋅∆+= nnnLt φφφ . (69)

Der Fehler eines derartig kombinierten Verfahrens ist gleich dem des impliziten Verfahrens, der stets der betragsmäßig kleinere ist. Der Vorteil von Prädiktor – Korrektor – Verfahren besteht darin, dass alle Komponenten direkt (explizit) ausgewertet werden können und dass die aufwändigere Lösung eines Gleichungssystems auf diese Weise vermieden wird.

2.5.3.5 Explizite und implizite Mehrschrittverfahren Mehrschrittverfahren benützen mehr als zwei Zeitebenen zur Approximation der zeitlichen Ableitungen. Das entsprechende Diskretisierungsschema kann beispielsweise durch Annahme eines polynomialen Verlaufs der Variablen in der Zeit (z. B. quadratisch bei drei Zeitebenen) oder durch eine geeignete Taylor – Reihenentwicklung definiert werden. Das Berechnen einer Lösung mit einem Mehrschrittverfahren muss immer mit einem Einschrittverfahren gestartet werden, da als Anfangsbedingung nur die Lösung bei 0t

bekannt ist. Erst nachdem die Lösungen bei 1t , …, 2−pt berechnet ist, kann mit einer p -

Zeitebenen Methode weitergerechnet werden. Im Verlauf der Rechnung müssen alle Variablenwerte aus den verwendeten (bekannten) Zeitebenen gleichzeitig mit den neuen Werten gespeichert werden, was bei großen Systemen und vielen Zeitebenen einen großen Speicherplatzbedarf zur Folge haben kann. Je nach Anzahl der verwendeten Zeitebenen, der Approximation der Zeitableitung und der Auswertung der rechten Seite lassen sich die unterschiedlichsten Mehrschrittverfahren definieren.

Eine wichtige Klasse expliziter Mehrschrittverfahren, die in der Praxis häufig verwendet werden, sind die Adams – Bashforth – Verfahren. Diese lassen sich durch Polynominterpolation für beliebige Ordnungen ableiten. In der Praxis werden jedoch nur die Verfahren bis zur Ordnung 4 benützt.

Eine wichtige Klasse von impliziten Methoden bilden die so genannten BDF – Verfahren (Backward – Differencing - Formula), die sich für beliebige Ordnung, durch Approximation

43

der Zeitableitung bei 1+nt mit Rückwärtsdifferenzenformeln unter Einbeziehung einer

entsprechenden Anzahl vorheriger Zeitpunkte ableiten lassen. Das BDF – Verfahren 1. Ordnung entspricht dem impliziten Euler – Verfahren. Insbesondere das BDF – Verfahren 2. Ordnung wird in der Praxis häufig verwendet. Bei diesem wird die unbekannte Funktion mittels einer durch die Funktionswerte zu den Zeitpunkten 1−nt , nt und 1+nt festgelegten

Parabel approximiert. Bei vergleichbar guten Stabilitätseigenschaften ist das Verfahren sowohl bezüglich des Rechenaufwands als auch bezüglich der Implementierung nur unwesentlich aufwendiger als das implizite Euler – Verfahren. Lediglich die Werte 1−nφ müssen zusätzlich gespeichert werden. Das Verfahren ist jedoch für äquidistante Gitter von 2. Ordnung, besitzt also in der Regel bei etwa gleichem Aufwand eine wesentlich bessere Genauigkeit. Ab der Ordnung 3 verschlechtern sich die Stabilitätseigenschaften der Verfahren mit wachsender Ordnung, so dass die Verwendung von BDF – Verfahren höherer Ordnung als 4 nicht ratsam ist.

2.6 Eigenschaften von Diskretisierungsverfahren In diesem Kapitel wird zusammenfassend auf Eigenschaften von numerischen Berechnungsmethoden eingegangen, denen für die praktische Anwendung eine besondere Bedeutung zukommt, da sie entscheidend für die Zuverlässigkeit und Funktionalität der darauf basierenden Verfahren sowie die „richtige“ Interpretation der damit erzielten Ergebnisse wird.

Bei der Anwendung von Diskretisierungsverfahren hat man für eine unbekannte Funktion ( )t,xφ für einen beliebigen Gitterpunkt Px und eine beliebige Zeit

nt die folgenden drei

Werte zu unterscheiden:

• die exakte Lösung der Differentialgleichung ( )nt,Pxφ ,

• die exakte Lösung der diskreten Gleichung n

Pφ und

• die tatsächlich berechnete Lösung n

Pφ~

.

Diese Werte stimmen im allgemeinen nicht überein, da bei den verschiedenen eingeführten Approximationen und der numerischen Lösung Fehler gemacht werden. Die Eigenschaften, die nachfolgend diskutiert werden, betreffen im Wesentlichen die Beziehung zwischen diesen Werten und die damit verbundenen Fehler.

Neben Modellfehlern, die an dieser Stelle außer Acht gelassen werden, hat man es entsprechend den unterschiedlichen Näherungswerten, bei der Anwendung eines numerischen Verfahrens im Wesentlichen mit zwei Arten von numerischen Fehlern zu tun:

• Der Diskretisierungsfehler ( ) n

Pnt φφ −,Px ist die Differenz der exakten Lösung der

Differentialgleichung und der exakten Lösung der diskreten Gleichung.

• Der Lösungsfehler n

P

n

P φφ −~

ist die Differenz der exakten Lösung der diskreten

Gleichungen und der tatsächlich berechneten Lösung.

Der für die praktische Anwendung wichtige Gesamtnumerikfehler ( ) n

Pnt φφ~

, −Px ist die

Differenz der exakten Lösung der Differentialgleichung und der tatsächlich berechneten Lösung und setzt sich aus dem Diskretisierungs- und dem Lösungsfehler zusammen. Der Lösungsfehler beinhaltet auch Anteile, die sich auf Grund einer eventuell nur näherungsweisen Lösung der diskreten Gleichungen ergeben. Dieser Anteil, der sich auf anhand von Residuen sehr gut kontrollieren lässt, ist jedoch in der Regel sehr klein, so dass dieser in diesem Kapitel nicht weiter betrachtet wird.

44

Für die Beziehungen zwischen den verschiedenen Lösungen und Fehlern, die in Abbildung 2.17 schematisch dargestellt sind, spielen insbesondere die Begriffe Konsistenz, Stabilität und Konvergenz eine wichtige Rolle. Diese werden nachfolgend kurz erläutert.

Abbildung 2.17: Zusammenhang zwischen Lösungen, Fehlern und Eigenschaften

2.6.1 Konsistenz

Man nennt eine Diskretisierungsmethode konsistent, wenn die diskretisierten Gleichungen beim Übergang 0, →∆∆ tx in die ursprünglichen Differentialgleichungen übergehen. Die Konsistenz definiert also einen Zusammenhang zwischen der exakten Lösung der Differentialgleichung ( )nt,Pxφ und der exakten Lösung der diskretisierten Gleichung n

Pφ . Die

Prüfung der Konsistenz kann durch eine Analyse des Abbruchfehlers n

Pτ erfolgen. Allgemein lässt sich der Abbruchfehler durch die Differenz zwischen der Differentialgleichung (interpretiert durch die Taylor – Reihenentwicklung von Ableitungen) und der diskreten Gleichung, wenn darin die exakten Variablenwerte eingesetzt werden, definieren. Strebt der Abbruchfehler für 0, →∆∆ tx gegen Null, so ist das Verfahren konsistent.

Beispiel 2.35:

Die Ableitungen der Diffusionsgleichung 2

2

xt ∂

∂⋅

⋅=

∂

∂ φ

αρ

κφ seien mit einem expliziten

Euleransatz in der Zeit und mit einem zentralen Differenzenansatz 2. Ordnung im Ort diskretisiert. Auf Grund der Taylor – Reihenentwicklung der Ableitungen

( )

( )4

4

42

2

11

2

2

2

2

21

12

2

2

xOx

x

xx

tOt

t

ttn

i

n

i

n

i

n

i

n

i

n

i

n

i

n

i

n

i

∆+

∂

∂⋅

∆+

∆

+⋅−=

∂

∂

∆+

∂

∂⋅

∆+

∆

−=

∂

∂

−+

+

φφφφφ

φφφφ

(70)

EExxaakkttee LLöössuunngg ddiisskkrreettee

GGlleeiicchhuunngg n

Pφ

EExxaakkttee LLöössuunngg DDiiffffeerreennttiiaall--gglleeiicchhuunngg

( )nt,Pxφ

BBeerreecchhnneettee LLöössuunngg

n

Pφ~

GGeessaammttnnuummeerriikkffeehhlleerr

KKoonnvveerrggeennzz

BBeerreecchhnnuunngg LLöössuunnggssffeehhlleerr

SSttaabbiilliittäätt

DDiisskkrreettiissiieerruunngg DDiisskkrreettiissiieerruunnggssffeehhlleerr

KKoonnssiisstteennzz

45

und der Betrachtung deren Abbruchfehler ( ) ( )2xOtO

n

P ∆+∆=τ , kann erkannt

werden, dass 0→n

Pτ für 0, →∆∆ tx gilt. Das Verfahren ist konsistent mit Konsistenzordnung 1 bezüglich der Zeit und Konsistenzordnung 2 bezüglich des Raumes.

Bei genügend kleinem x∆ und t∆ bedeutet höhere Ordnung des Abbruchfehlers auch eine höhere Genauigkeit der Lösung. Im Allgemeinen sagt die Ordnung des Abbruchfehlers aber nur, wie schnell die Fehler abklingen, wenn das Gitter verfeinert wird. Über die absoluten Werte des Lösungsfehlers auf einem gegebenen Gitter, die auch von der Lösung selbst abhängen, sagt die Ordnung des Verfahrens nichts aus.

2.6.2 Stabilität

Das Konzept der Stabilität dient zur Herstellung eines Zusammenhangs zwischen der

tatsächlich berechneten Lösung n

Pφ~

und der exakten Lösung n

Pφ der diskreten Gleichung. Stabilität ist also eine Eigenschaft des Verfahrens und nicht des Problems. Es existieren in der Literatur eine ganze Reihe unterschiedlicher Definitionen der Stabilität, die für verschiedene Untersuchungen sinnvoll sind. Im Rahmen dieser Vorlesung wird eine einfache, den erörterten Anwendungen angepasste Definition verwendet. Eine

Diskretisierungsmethode heißt stabil, falls der Lösungsfehler n

P

n

P φφ −~

im ganzen

Lösungsgebiet und für alle Zeitschritte begrenzt bleibt.


Das Verfahren ix

i ex−

+ =1 zur Lösung der Gleichung x

ex−= ist stabil.


Durch Umformen von xex

−= erhält man die Gleichung xx −=ln . Diese wird mit

dem Iterationsverfahren ii xx ln1 −=+ für ,...3,2,1=i mit 567.01 =x approximiert. Trotz

der Wahl des Startwertes „nahe“ der exakten Lösung 0,56714329 konvergiert das Verfahren nicht. Das Verfahren ist instabil.

Die generelle Vorgehensweise bei der Überprüfung von Stabilität eines Verfahrens besteht darin, zu untersuchen, wie sich kleine Störungen (verursacht z. B. durch Rundungsfehler) auf die späteren Zeitebenen auswirken. Die entscheidende Frage dabei ist, ob diese Störungen durch die Diskretisierungsvorschrift gedämpft (dann ist das Verfahren stabil) oder verstärkt (dann ist das Verfahren instabil) werden. Eine derartige Stabilitätsanalyse kann mit verschiedenen Methoden durchgeführt werden. Die gängigsten Methoden (für Details wird auf Spezialliteratur verwiesen) dazu sind die von Neumann – Analyse, die Matrix – Methode, und die Methode der kleinsten Störungen (Perturbationsmethode). Für allgemeine Probleme sind solche Analysen, wenn sie überhaupt auf analytischem Weg durchgeführt werden können, in der Regel sehr schwierig und aufwendig.

Beispiel 2.38:


2

xd

t ∂

∂⋅=

∂

∂ φφ mit

αρ

κ

⋅=d seien mit einem

expliziten Euleransatz in der Zeit und mit einem zentralen Differenzenansatz 2. Ordnung im Ort diskretisiert. Die resultierende diskrete Gleichung lautet

( ) n

i

n

i

n

i

n

ix

td

x

tdφφφφ ⋅

∆

∆⋅⋅−++⋅

∆

∆⋅= −+

+

2112

121 . (71)

46

Eine einfache heuristische Betrachtung des Problems verlangt eine Zunahme von 1+n

iφ bei einer Erhöhung von n

iφ . Eine solche gleichsinnige Änderung von n

iφ und 1+n

iφ ist im Allgemeinen nur gesichert, wenn alle Koeffizienten der diskreten

Gleichung (71) positiv sind. Da 02

>∆

∆⋅

x

td gilt, braucht nur der Koeffizient

0212

>∆

∆⋅⋅−

x

td (72)

gefordert werden. Eine einfache Umformung ergibt dafür die für Stabilität zu

erfüllende Bedingung d

xt

⋅

∆<∆

2

2

. Das Verfahren ist also nur bedingt stabil.


Die iterative Auswertung von 2

1

112

xd

t

nnnn

i

n

i iii

∆

+⋅−⋅=

∆

−−+

+ φφφφφ mit gegebenem

αρ

κ

⋅=d sowie gegebenen Anfangs- und Randwerten wird in dem EXCEL Beispiel

veranschaulicht. Dabei ist die Ortsdiskretisierung horizontal aufgetragen, während die Zeitdiskretisierung vertikal verläuft. Bei unterschiedlicher Wahl von t∆ und x∆ kann

die für Stabilität zu erfüllende Bedingung d

xt

⋅

∆<∆

2

2

überprüft werden.

Beispiel 2.40:


2

xd

t ∂

∂⋅=

∂

∂ φφ mit

αρ

κ

⋅=d seien mit einem

impliziten Euleransatz in der Zeit und mit einem zentralen Differenzenansatz 2. Ordnung im Ort diskretisiert. Die resultierende diskrete Gleichung lautet

( ) n

i

n

i

n

i

n

ix

td

x

tdφφφφ ++⋅

∆

∆⋅=

∆

∆⋅⋅+ +

−+

++ 1

1

1

12

1

221 . (73)

Alle Koeffizienten der diskreten Gleichung (73) sind positiv, so dass keine Probleme hinsichtlich einer nicht gleichgesinnten Änderung von φ zu erwarten sind. Das Verfahren ist unabhängig von der Wahl der Zeitschrittweite und somit unbedingt stabil.

2.6.3 Konvergenz

Eine entscheidende Forderung an ein Diskretisierungsverfahren ist, dass die numerisch berechnete Lösung bei immer feiner werdenden Gittern immer besser mit der exakten Lösung der Differentialgleichung übereinstimmt. Diese Eigenschaft wird als Konvergenz des Verfahrens bezeichnet.

Die Begriffe Konsistenz, Stabilität und Konvergenz stehen in enger Beziehung zueinander. Für lineare Probleme stellt das fundamentale Äquivalenztheorem von Lax den Zusammenhang her, das unter gewissen Voraussetzungen an das kontinuierliche Problem (diese werden hier nicht näher erläutert), lautet: Für ein konsistentes Diskretisierungsschema ist die Stabilität eine notwendige und hinreichende Bedingung für dessen Konvergenz. Auf

47

Grund des Lax – Theorems kann man bei der Analyse eines Diskretisierungsschemas wie folgt vorgehen:

• Analyse der Konsistenz: man erhält die Ordnung des Schemas und den Abbruchfehler.

• Analyse der Stabilität: man erhält Informationen über das Fehlerverhalten.

Hieraus erhält man dann, ohne weitere Untersuchungen, Aussagen über die Konvergenz des Verfahrens, welche die entscheidende Eigenschaft ist. Insbesondere erhält man aus der Stabilitätsanalyse Informationen über die Wahl der Gitter- bzw. Zeitschrittweite, was für die praktische Anwendung von besonderer Bedeutung ist.

Für nichtlineare Probleme gilt das Lax – Theorem zwar nicht, aber Stabilität und Konsistenz sind auch hier wichtige Anhaltspunkte für ein vernünftig funktionierendes Verfahren.

2.7 Übungsbeispiele

2.7.1. Wie lautet eine allgemeine Formulierung eines Anfangswertproblems der Ordnung n mit gegebenen Anfangswerten?

2.7.2. Wie lautet die allgemeine Wärmeleitungsgleichung zur Beschreibung des thermischen Verhaltens eines Stabs der Länge L . Geben Sie Anfangs- und Randwerte an.

2.7.3. Erklären Sie die Begriffe absoluter Fehler, relativer Fehler und Residuum.

2.7.4. Listen Sie Ursachen auf, die numerische Fehler zur Folge haben können.

2.7.5. Gegeben sei ein lineares Gleichungssystem bxA =⋅ mit

=

21

12A und

−=

2

2b sowie ein numerisch ermittelter Lösungsvektor

−=

004.2

01.2~x .

Bestimmen Sie den absoluten Fehlervektor sowie den Residuenvektor.

2.7.6. Folgt aus einem kleinen Residuum ein kleiner absoluter Fehler? Welche Eigenschaft des Problems spielt in diesem Zusammenhang eine Rolle?

2.7.7. Die Lösung der Gleichung xex

−= wird mit einem Iterationsverfahren ix

i ex−

+ =1 für ,...3,2,1=i mit 11 =x approximiert (exakte Lösung ist

0.56714329). Der Verlauf der Iterationslösung ist in der folgenden Abbildung dargestellt. Welche Aussage kann über die Kondition des Problems getroffen werden?

0

0.2

0.4

0.6

0.8

1

1.2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Exakte Lösung

Iterationslösung

48

2.7.8. Die Lösung der Gleichung ( )15.0 2 +⋅= xx wird mit einem Iterationsverfahren

( )15.0 2

1 +⋅=+ ii xx für ,...3,2,1=i mit 1.11 =x approximiert (exakte Lösung ist

1.0). Der Verlauf des Residuums ist in der folgenden Abbildung dargestellt. Welche Aussage kann über die Kondition des Problems getroffen werden?

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Iterationslösung

Residuum

2.7.9. Gegeben seien zwei lineare Gleichungssysteme 11 bxA =⋅ und 22 bxA =⋅

mit

=

4000230001

800600A ,

=

10000

2001b und

=

10001

2002b . Bestimmen Sie die

beiden Lösungen 1x und 2x . Welche Aussage kann man über die Kondition dieses Problems machen. Begründen Sie diese Aussage. Wie kann man die Kondition beschreiben?

2.7.10. Gegeben sei ein lineares Gleichungssystem 11 bxA =⋅ mit

=

4000230001

800600A und

=

10000

2001b . Bestimmen Sie die Konditionszahl

der Matrix A .

2.7.11. Welche besonderen Eigenschaften oder Strukturen von Systemmatrizen die bei der Lösung von linearen Gleichungssystemen zu berücksichtigt sind, können Sie nennen?

2.7.12. Bestimmen Sie die Zeilensummennorm, die Spaltensummennorm sowie die

Frobeniusnorm für die Matrix

=

21

12A .

2.7.13. Bestimmen Sie die Euklidische Vektornorm für den Vektor

−=

2

2x .

Bestimmen Sie für

−=

2

2x und die Matrix

=

21

12A die Euklidische

Matrixnorm.

2.7.14. Geben sei eine allgemeine Matrix nnA ×ℜ∈ . Wie lautet die Definition folgender Matrizeneigenschaften:

• Symmetrisch

• Orthogonal

49

• Positiv semidefinit

• Positiv definit

• Negativ semidefinit

• Negativ definit

2.7.15. Die Systemmatrix einer eindimensionalen Diskretisierung hat für ein

beliebiges +ℜ∈h die Form

⋅=

21-0

1-21-

01-21

hA . Bestimmen Sie die

Eigenwerte der Matrix. Was kann über die Eigenschaften der Matrix A

gesagt werden.

2.7.16. Zeigen Sie, dass die Matrix

=

αα

αα

cossin

sin - cosA für beliebige Winkel

[ [πα 2,0∈ orthogonal ist.

2.7.17. Geben sei eine allgemeine Matrix nnA ×ℜ∈ . Wie lautet die Definition folgender Matrizeneigenschaften:

• Linke untere Dreiecksmatrix

• Rechte obere Dreiecksmatrix

• Diagonalmatrix

• Strukturierte Matrix

• Dichte Matrix

2.7.18. Wie ist die Bandbreite einer Matrix definiert? Bestimmen Sie die Bandbreite der folgenden Matrix A .

−

−

−

−

−

−

−

−

−

−

−

−

=

410100000000

141010000000

014101000000

101410100000

010141010000

001014101000

000101410100

000010141010

000001014101

000000101410

000000010141

000000001014

A

2.7.19. Diskutieren Sie den Gauß-Algorithmus zur Berechnung von 0bxA =+⋅ . Welche Verfahrensschritte müssen durchgeführt werden? Geben Sie den Rechenaufwand dieses Verfahrens an.

50

2.7.20. Diskutieren Sie das Cholesky - Verfahren zur Berechnung von 0bxA =+⋅ . Welche Verfahrensschritte müssen durchgeführt werden? Welche Voraussetzungen müssen für die Anwendung dieses Verfahrens erfüllt sein? Geben Sie den Rechenaufwand dieses Verfahrens an.

2.7.21. Welche charakteristischen Lösungsschritte passieren bei der iterativen Lösung linearer Gleichungssysteme? Nennen Sie einige iterative Verfahren zur Lösung von 0bxA =+⋅ . Diskutieren Sie unterschiedliche Abbruchkriterien zum Beenden des iterativen Prozess.


=

21

12A und

−=

2

2b sowie ein mittels Iterationsverfahren in der Iteration i ermittelter

Lösungsvektor

−=

004.2

01.2~x . Weiters sei 05.0=ε gegeben. Überprüfen Sie,

ob das in der Iteration berechnete Ergebnis das Abbruchkriterium ( ) ε≤

=i

nixr~

,...,1max erfüllt.


=

21

12A und

−=

2

2b sowie ein mittels Iterationsverfahren in der Iteration i ermittelter

Lösungsvektor

−=

004.2

01.2~x . Weiters sei 05.0=ε gegeben. Überprüfen Sie,

ob das in der Iteration berechnete Ergebnis das Abbruchkriterium ε≤∑=

n

ii

1

~xr

erfüllt.

2.7.24. Welche Aussagen können betreffs des Konvergenzverhaltens unterschiedlicher iterativer Verfahren zur Lösung von linearen Gleichungssystemen getroffen werden?

2.7.25. Diskutieren Sie das Jacobi - Verfahren zur Lösung von 0bxA =+⋅ . Geben Sie die verwendete Iterationsvorschrift an.

2.7.26. Diskutieren Sie das Gauß-Seidel - Verfahren zur Lösung von 0bxA =+⋅ . Geben Sie die verwendete Iterationsvorschrift an. Warum hat dieses Verfahren im Vergleich zum Jacobi - Verfahren bessere Konvergenzeigenschaften?

2.7.27. Diskutieren Sie das SOR - Verfahren zur Lösung von 0bxA =+⋅ . Geben Sie die verwendete Iterationsvorschrift an. Wann und wie wird Überrelaxiation, wann und wie wird Unterrelaxation angewendet?

2.7.28. Gegeben ist das lineare Gleichungssystem

=

−

−

−

+

⋅

0

0

0

1

2

2

101

020

102

3

2

1

x

x

x

.

Rechnen Sie eine Jacobi-Iteration sowie eine Iteration unter Verwendung

51

des Gauß-Seidel – Verfahrens. Verwenden Sie den Startvektor ( )t

100=0x .


=

−

−

−

+

⋅

0

0

0

1

2

2

101

020

102

3

2

1

x

x

x

.

Rechnen Sie jeweils eine Iteration unter Verwendung des SOR – Verfahrens mit 9.0=ω bzw. 1.1=ω . Verwenden Sie den Startvektor ( )t

100=0x .

2.7.30. Das Verfahren der konjugierten Gradienten kombiniert das Verfahren des steilsten Abstiegs mit dem Verfahren der konjugierten Richtungen. Treffen Sie jeweils Aussagen über die verwendete Basis, die Vor- und Nachteile.

2.7.31. Welches Prinzip liegt dem konjugierten Gradienten Verfahren zur Lösung von 0bxA =+⋅ zu Grunde? Erklären Sie die Bestimmung der Relaxationsrichtung sowie die Bestimmung der Schrittweite.


=

−

−

−

+

⋅

0

0

0

1

2

2

101

020

102

3

2

1

x

x

x

.

Rechnen Sie eine Iteration unter Verwendung des konjugierten Gradienten Verfahrens. Verwenden Sie den Startvektor ( )t

100=0x .

2.7.33. Geben Sie 5 Beispiele für moderne iterative Verfahren an.

2.7.34. Erklären Sie an einem allgemeinen Beispiel 0bxA =+⋅ das Prinzip der Vorkonditionierung. Nennen Sie unterschiedliche Vorkonditionierungsvarianten.

2.7.35. Bei Anwendung des Cholesky – Verfahrens resultiert eine untere

Dreiecksmatrix

=

15.21

035.1

002

L . Schreiben sie die Elemente der

Dreiecksmatrix in kompakter Speicherung. Welche Aussage kann betreffs der Anzahl der in kompakter Form zu speichernden Elemente für eine allgemeine Matrix nn×ℜ=L gemacht werden?

2.7.36. Gegeben Sei eine Systemmatrix A . Die Matrix resultierend aus der

Diskretisierung der 2D Poissongleichung ( )2

2

2

2

,yx

yxf∂

∂+

∂

∂=

φφ mittels zentraler

Differenzenquotienten ji

jijijijijijif

yx,2

1,,1,

2

,1,,1 22=

∆

+−+

∆

+− +−+− φφφφφφ. Schreiben

Sie die Elemente der Matrix unter Verwendung einer kompakten Speichermethode (Skyline – Speicherung). Analysieren Sie den Speicherbedarf der kompakten Methode in Relation zur Speicherung aller Matrixelemente. Bestimmen Sie die Bandbreite sowie die Semibandbreite der Matrix. Welches Verfahren kann zur Lösung des resultierenden linearen Gleichungssystems 0bxA =+⋅ verwendet werden? Geben Sie eine Aufwandsabschätzung an.

52

−

−

−

−

−

−

−

−

−

−

−

−

=

410100000000

141010000000

014101000000

101410100000

010141010000

001014101000

000101410100

000010141010

000001014101

000000101410

000000010141

000000001014

A

2.7.37. Erläutern Sie den Begriff Diskretisierung.

2.7.38. Welche Methoden zur Ortsdiskretisierung wurden in der Vorlesung besprochen?

2.7.39. Was ist der Unterschied zwischen strukturierten und unstrukturierten Netzen? Vergleichen Sie die beiden Netztypen bezüglich „Modellierung komplexer Geometrien“, „Lokale Netzadaptivität“ und „Automatisierte Vernetzung“. Welche Möglichkeiten der gemischten Anwendung der beiden Netztypen gibt es?

2.7.40. Welche Vernetzungsstrategien für strukturierte und unstrukturierte Netze gibt es?

2.7.41. Leiten Sie die Rückwärtsdifferenzenformel 1. Ordnung zur Approximation

von x∂

∂φ unter Verwendung einer Taylor – Reihe her.

2.7.42. Leiten Sie die Differenzenformel xx

ii

∆⋅

−≈

∂

∂ −+

2

11 φφφ durch Kombination der

Rückwärts- sowie der Vorwärtsdifferenzenformel 1. Ordnung her.

2.7.43. Diskutieren Sie die Methode der Finiten Volumen.

2.7.44. Erklären Sie die Strategie von Multigrid – Verfahren.

2.7.45. Erklären Sie die Begriffe Zeitdiskretisierung sowie Ortsdiskretisierung am

Beispiel des 1D Konvektions-Diffusionsproblems 2

2

21x

cx

ct ∂

∂⋅+

∂

∂⋅=

∂

∂ φφφ.

Welche unterschiedlichen Methoden wurden in der Vorlesung erörtert? Wie müssen Anfangswerte sowie Randwerte berücksichtigt werden?

2.7.46. Für das 1D Konvektions-Diffusionsproblem 02

2

21 =∂

∂⋅+

∂

∂⋅

xc

xc

φφ mit den

Randbedingungen ( ) 000

=∂

∂==

=xxx

φφ und ( ) 11

1

=∂

∂==

=xxx

φφ und der

Schrittweite 1.0=∆x soll mit den Zentraldifferenzenformeln für erste und

53

zweite Ableitungen ein so genanntes Differenzensystem aufgestellt werden. Leiten Sie die Differenzengleichung für äquidistanten Punktabstand her. Aus wie vielen Gleichungen besteht das Gleichungssystem? Stellen Sie das komplette Gleichungssystem in der Form „Koeffizientenmatrix“ * „Vektor der Unbekannten“ = „Vektor der rechten Seite“ auf.

2.7.47. Erklären Sie die Begriffe Kausalität der Zeitdiskretisierungsmethoden, Einschrittverfahren, Mehrschrittverfahren, explizite Verfahren, implizite Verfahren.

2.7.48. Wenden Sie das explizite sowie das implizite Euler Verfahren zur Zeitdiskretisierung der eindimensionalen Wärmeleitungsgleichung

2

2

xt ∂

∂⋅

⋅=

∂

∂ φ

αρ

κφ an. Die Diskretisierung der Zeitableitung soll mittels des

zentralen Differenzenansatzes 2. Ordnung erfolgen. Diskutieren Sie die beiden Verfahren bezüglich Lösungsschema und Abbruchfehlerordnung. Was versteht man unter Prädiktor – Korrektor Verfahren? Kombinieren Sie am Beispiel der Wärmeleitungsgleichung das explizite und das implizite Euler Verfahren zu einem Prädiktor – Korrektor Verfahren.

2.7.49. Erklären Sie am Beispiel der eindimensionalen Wärmeleitungsgleichung die Begriffe “Exakte Lösung der Differentialgleichung”, ”Exakte Lösung der diskreten Gleichung” und “Numerische Lösung der diskreten Gleichung”. Über welche Fehlertypen hängen diese Begriffe zusammen? Wie sind die Begriffe “Konsistenz”, “Konvergenz” und “Stabilität” in diesem Zusammenhang definiert?

2.7.50. Was versteht man unter dem Begriff Konsistenz einer diskreten Gleichung? Wie kann die Konsistenz überprüft werden? Bestimmen Sie die Konsistenz am Beispiel der eindimensionalen Wärmeleitungsgleichung unter Verwendung von Orts- und Zeitdiskretisierungsmethoden eigener Wahl.

2.7.51. Die Ableitungen der Diffusionsgleichung 2

2

xd

t ∂

∂⋅=

∂

∂ φφ mit constd = seien mit

einem Euleransatz in der Zeit und mit dem zentralen Differenzenansatz 2. Ordnung im Ort diskretisiert. Erklären Sie den Begriff Stabilität an den Beispielen Explizites Euler Verfahren und Implizites Euler Verfahren. Leiten Sie die Stabilitätsbedingung für das explizite Euler Verfahren her.

2.7.52. Das Verfahren ix

i ex−

+ =1 wird zur Lösung der Gleichung xex

−= verwendet. Die folgende Abbildung zeigt den Verlauf der Lösungen der ersten 20 Iterationen. Was kann über die Stabilität des Verfahrens gesagt werden?

54

0

0.2

0.4

0.6

0.8

1

1.2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Exakte Lösung

Iterationslösung

2.7.53. Das Verfahren ii xx ln1 −=+ wird zur Lösung der Gleichung xex

−= verwendet.

Die folgende Abbildung zeigt den Verlauf der Lösungen der ersten 20 Iterationen. Was kann über die Stabilität des Verfahrens gesagt werden?

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1 2 3 4 5 6 7 8

Exakte Lösung

Iterationslösung

2.7.54. Erklären Sie den Begriff Konvergenz eines numerischen Verfahrens.

55

3 Abbildungsverzeichnis

Abbildung 2.1: Beschleunigung im Bereich der Rechnergeschwindigkeit (grün) und der

numerischen Verfahren (blau) in den letzten Dekaden .................................................. 12

Abbildung 2.2: Vorgehensweise bei der Anwendung numerischer Simulation zur Lösung von

Ingenieurproblemen ...................................................................................................... 13

Abbildung 2.3: Herleitung des CG – Verfahrens ................................................................... 24

Abbildung 2.4: Iterationsvorschrift des konjugierten Gradientenverfahrens .......................... 25

Abbildung 2.5: Zusammenhänge zwischen unterschiedlichen iterativen

Gleichungssystemlösern ............................................................................................... 26

Abbildung 2.6: Beispiele zwei- und dreidimensionaler Gitterelemente ................................. 29

Abbildung 2.7: Beispiel eines randapproximierenden Netzes (links) eines kartesischen Netzes

(Mitte) und eines überlappenden Netzes (rechts) ........................................................... 29

Abbildung 2.8: Beispiel eines strukturierten Netzes (links) und eines unstrukturierten Netzes

(rechts) ......................................................................................................................... 30

Abbildung 2.9: Beispiel eines blockstrukturierten Gitters (links) und eines hierarchisch

strukturierten Gitters (rechts) ........................................................................................ 32

Abbildung 2.10: Diskretisierung eines ebenen Gebietes Ω (beschränkt durch Γ ) mittels

rechteckiger finiter Volumina ....................................................................................... 35

Abbildung 2.11: Zweidimensionales Beispiel: Zellenorientierte Anordnung der Variablen für

das Finite – Volumen – Verfahren ................................................................................ 35

Abbildung 2.12: Fehlerreduktion bei Anwendung des Jacobi – Verfahrens zur Lösung des

eindimensionalen Diffusionsproblems für unterschiedliche Startwerte .......................... 36

Abbildung 2.13: Zusammenhang zwischen Informationsaustausch und Gittergröße auf Grund

von Nachbarschaftsbeziehungen von Diskretisierungsvorschriften ................................ 36

Abbildung 2.14: Schematische Darstellung der Restriktion / Prolongation bei AMG Verfahren

..................................................................................................................................... 37

Abbildung 2.15: Zusammenhang zwischen räumlicher und zeitlicher Diskretisierung .......... 39

Abbildung 2.16: Vorgehensweise und Informationsfluss bei Anwendung des expliziten Euler

- Verfahrens auf die eindimensionale Diffusionsgleichung ............................................ 41

Abbildung 2.17: Zusammenhang zwischen Lösungen, Fehlern und Eigenschaften ............... 44

56

4 Link- und Literaturverzeichnis

4.1 Einführung in Visual Basic (VBA)

• http://www.ti5.tu-harburg.de/manual/vba5/httoc.htm

• http://mypage.bluewin.ch/reprobst/WordFAQ/VBA.htm

• http://www.jumper.ch

• http://people.freenet.de/a.seeberger/index.html

• http://www.excel-cd.de

• http://vb-fun.de/vb/index.htm

• http://excel-inside.de

• http://www.fh-potsdam.de/~Bauing/fachgebiete/vba/startseite.htm

• http://www.arstechnica.de/computer/msoffice/vbaindex.html

• http://www.excel-vba.de

• http://www.excel-center.de/index.php

• http://www.office-loesung.de

• http://www.excel-center.de/forum

• http://www.herber.de/forum

• Hilfetexte von VBA

• Doberenz W., Kowalski T.: Microsoft Access – Programmierung. Microsoft Press, 2001

• Schumann A.: Programmieren mit Visual Basic for Applications (VBA), Seminar für Diplomanden und Doktoranden, FR Geoinformatik, Skriptum, 2004

4.2 Numerische Grundlagen

• Bathe K. – J.: Finite – Elemente – Methode, Springer, 1986

• Björck A., Dahlquist G.: Numerische Methoden, Oldenburg, 1979

• Bronstein I. N., Semendjajew K. A.: Taschenbuch der Mathematik, Teubner, 1991

• Fischer G.: Lineare Algebra, Vieweg, 1978

• Hackbusch W.: Iterative Lösung großer schwach besetzter Gleichungssysteme, Teubner, 1993

• Heuberger C., Wallner H.: Technische Numerik, Skriptum, 2001

• Hirsch C.: Numerical Computation of Internal and External Flows, Wiley, Chichester, 1988

• Meister A.: Numerik linearer Gleichungssysteme, Vieweg, 2005

• Rung T., Xue L., Yan J., Schatz M., Thiele F.: Numerische Methoden in der Thermo- und Fluiddynamik, Skriptum, 2002

• Schäfer M.: Numerik im Maschinenbau, Springer, 1999

• Schwarz H. R.: Numerische Mathematik, Teubner, 1997

57

• Smith G. D.: Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford Applied Mathematics and Computing Science Series, 1978

• Steinbach O.: Lösungsverfahren für lineare Gleichungssysteme, Teubner, 2005

Numerische Methoden in der angewandten Thermodynamik ... · Numerische Methoden in der angewandten...

Documents

Transcript of Numerische Methoden in der angewandten Thermodynamik ... · Numerische Methoden in der angewandten...