ifiCompu
SC
nt cie
S
tin
g
Statistische Auswertung Numerischer Simulationen Diplomarbeit
Martin Krosche
Institut f¨ur Wissenschaftliches Rechnen Technische Universit¨at Braunschweig Braunschweig, Germany
Erkl¨arung: Ich versichere, dass ich die Arbeit ohne fremde Hilfe und ohne Benutzung anderer als der angegebenen Quellen angefertigt habe, und dass die Arbeit in gleicher oder a¨ hnlicher Form noch keiner anderen Pr¨ufungsbeh¨orde vorgelegen hat und von dieser als Teil einer Pr¨ufungsleistung angenommen wurde. Alle Ausf¨uhrungen, die w¨ortlich oder sinngem¨aß u¨ bernommen wurden, sind als solche gekennzeichnet.
Braunschweig, den
Martin Krosche
Inhaltsverzeichnis 1
Einleitung
2
hochdimensionale Integration 2.1 Monte Carlo Methode . . . . . . . . . . . . 2.2 Quasi-Monte Carlo Methode . . . . . . . . 2.2.1 Halton-Folge . . . . . . . . . . . . 2.2.2 Hammersley-Folge . . . . . . . . . 2.2.3 arithmetische Folge . . . . . . . . . 2.3 Smolyak Algorithmus . . . . . . . . . . . . 2.3.1 Gauss-Legendre Quadraturformel . 2.3.2 Clenshaw-Curtis Quadraturformeln 2.3.3 Kronrod-Patterson Quadraturformel 2.3.4 Symmetrien . . . . . . . . . . . . . 2.3.5 Vorteile der Symmetrieausnutzung . 2.4 Experimente . . . . . . . . . . . . . . . . .
3
4
5
1
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
System zur statistischen Auswertung numerischer Simulationen 3.1 Teilsystem Quadraturalgorithmen . . . . . . . . . . . . . . . 3.2 Teilsystem objektive Funktion . . . . . . . . . . . . . . . . . 3.3 Interaktion zwischen den Teilsystemen . . . . . . . . . . . . . 3.4 Symmetrieausnutzung . . . . . . . . . . . . . . . . . . . . . . Beispielmodell und dessen Auswertung 4.1 Semantik . . . . . . . . . . . . . . 4.1.1 Skalierungsfunktionen . . . 4.2 Simulation und Bewertung . . . . . 4.3 Ellipsen als Einschl¨usse . . . . . . . 4.3.1 Semantik . . . . . . . . . . 4.3.2 Simulation und Bewertung .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . . . . . . . .
. . . .
. . . . . .
Design und Implementierung 5.1 grundlegende Klassen und Projekte . . . . . . . . . . . . . . . . 5.2 System zur statistischen Auswertung numerischer Simulationen 5.3 Werkzeuge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Component Template Library (CTL) . . . . . . . . . . . 5.3.2 Parallel Finite Element Program (ParaFEP) . . . . . . . 5.3.3 smolpack und finpack . . . . . . . . . . . . . . . . . .
i
. . . . . . . . . . . .
2 3 5 6 6 7 8 11 11 14 14 17 22
. . . .
25 29 30 32 32
. . . . . .
33 33 35 38 38 38 42
. . . . . .
44 46 50 52 52 52 52
6
Beispiele 6.1 Motivationsbeispiel . . . . . . . . 6.2 konkrete Beispiele . . . . . . . . 6.2.1 Semantik . . . . . . . . . 6.2.2 Simulation und Bewertung 6.2.3 Quader . . . . . . . . . . 6.2.4 Stoßlasche . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
53 53 55 55 56 57 61
7
Zusammenfassung
76
8
Ausblick
77
9
Danksagung
77
ii
Abbildungsverzeichnis 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
MC Gitter, 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . QMC Gitter, 2D Halton-Punktfolge . . . . . . . . . . . . . . . . QMC Gitter, 2D Hammersley-Punktfolge . . . . . . . . . . . . . QMC Gitter, 2D arithmetische Punktfolge . . . . . . . . . . . . . Gitter des Smolyak Algorithmus mit Gauss-Legendre Quadraturformel, skaliert auf [0, 1]2 . . . . . . . . . . . . . . . . . . . . . . Gitter des Smolyak Algorithmus mit Clenshaw-Curtis Quadraturformel, skaliert auf [0, 1]2 . . . . . . . . . . . . . . . . . . . . . . Gitter des Smolyak Algorithmus mit Kronrod-Patterson Quadraturformel, skaliert auf [0, 1]2 . . . . . . . . . . . . . . . . . . . . . Gitter des Smolyak Algorithmus mit verz¨ogerter KronrodPatterson Quadraturformel, skaliert auf [0, 1]2 . . . . . . . . . . . . Tensorprodukt der eindimensionalen Basispolynome einer Basis des Polynomraumes P≤p 2 . . . . . . . . . . . . . . . . . . . . . . . ¯ ≤p f¨ur den Polynomraum P ¯ ≤p . . . . . . . Konstruktion der Basis B 2 d ≤p ¯ ≤p ≤p ˜ u¨ ber Grad p. . . Dimension der Polynomr¨aume P50 , P50 und P 50 Funktion f mit n = 5, m = 100, ε = 0.05 und d = 1. . . . . . . . Funktion f mit n = 5, m = 100, ε = 0.05 und d = 2. . . . . . . . relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 10. . . . . . . . . . . . . . . . . . . . . . . . . . . relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 20. . . . . . . . . . . . . . . . . . . . . . . . . . . relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 10 mit Symmetrieausnutzung. . . . . . . . . . . . Gesamtsystem, Perspektive 1. . . . . . . . . . . . . . . . . . . . Gesamtsystem, Perspektive 2. . . . . . . . . . . . . . . . . . . . Semantik-Funktion. . . . . . . . . . . . . . . . . . . . . . . . . . Inverse Verteilungsfunktion zur Generierung von Einschl¨ussen. . . Semantik-Funktion zur Generierung von Material mit Einschluss. Ellipse. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inverse Verteilungsfunktion zur Generierung von Ellipsen als Einschl¨ussen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . gleiche Ellipsen trotz unterschiedlicher Drehwinkel. . . . . . . . . eingeschr¨ankter Drehwinkel f¨ur Ellipsen. . . . . . . . . . . . . . Generierung einer Ellipse. . . . . . . . . . . . . . . . . . . . . . Semantik-Funktion mit Kegelschnittgleichung. . . . . . . . . . . Material mit 3 Ellipsen. . . . . . . . . . . . . . . . . . . . . . . . Simulation eines 2D FE Modells. . . . . . . . . . . . . . . . . . . Verzeichnisbaum. . . . . . . . . . . . . . . . . . . . . . . . . . . iii
4 7 8 9 12 13 15 16 20 21 23 24 25 26 27 28 30 31 32 35 36 39 40 40 41 42 43 43 44 45
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
UML Klassendiagramm, Auszug aus den Klassen Monte carlo und Point gen. . . . . . . . . . . . . . . . . . UML Klassendiagramm, Auszug aus den Klassen Smolyak und One dim quadrature. . . . . . . . . . . . . . . . . . . . . . statistische Auswertung einer konkreten objektiven Funktion. . . . Gasbetonstein Nr.1. . . . . . . . . . . . . . . . . . . . . . . . . . Gasbetonstein Nr.2. . . . . . . . . . . . . . . . . . . . . . . . . . simuliertes Quader Modell. . . . . . . . . . . . . . . . . . . . . . relativer Fehler bei 3 Kreisen. . . . . . . . . . . . . . . . . . . . . R¨aumliche Materialverteilungen, generiert aus Integrationspunkten der Monte Carlo Methode und der Quasi-Monte Carlo Formeln. Verteilung der maximalen Verschiebung bei 3 Kreisen. . . . . . . R¨aumliche Materialverteilungen, generiert aus Integrationspunkten des Smolyak Algorithmus. . . . . . . . . . . . . . . . . . . . Verteilungsfunktion zur Skalierung der Halbachsen bei 10 Ellipsen. simulierte Stoßlasche. . . . . . . . . . . . . . . . . . . . . . . . . relativer Fehler bei 10 Ellipsen. . . . . . . . . . . . . . . . . . . . R¨aumliche Materialverteilungen, generiert aus aufeinanderfolgenden 50-dimensionalen Halton-Punkten. . . . . . . . . . . . . . R¨aumliche Materialverteilungen, generiert aus den ersten 50dimensionalen Halton-Punkten. . . . . . . . . . . . . . . . . . . . R¨aumliche Materialverteilungen, generiert aus aufeinanderfolgenden 50-dimensionalen Punkten der arithmetischen Folge. . . . R¨aumliche Materialverteilungen, generiert aus aufeinanderfolgenden 50-dimensionalen Punkten der Monte Carlo Methode. . . Vergleich zweier Zufallsgeneratoren. . . . . . . . . . . . . . . . . relativer Fehler bei 20 Ellipsen. . . . . . . . . . . . . . . . . . . . relativer Fehler bei 30 Ellipsen. . . . . . . . . . . . . . . . . . . . R¨aumliche Materialverteilungen, generiert aus den ersten 150dimensionalen Halton-Punkten. . . . . . . . . . . . . . . . . . . . Verteilungen der maximalen Verschiebung bei 10, 20 und 30 Ellipsen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . relativer Fehler bei Verwendung der Monte Carlo Methode. . . . . relativer Fehler bei Verwendung der auf der arithmetischen Folge basierenden Quasi-Monte Carlo Methode. . . . . . . . . . . . . . relativer Fehler bei Verwendung der auf der Halton-Punktfolge basierenden Quasi-Monte Carlo Methode. . . . . . . . . . . . . .
iv
47 48 51 54 54 58 59 59 60 61 65 65 67 68 68 68 69 69 70 71 72 72 73 74 75
Tabellenverzeichnis 1 2 3 4 5 6 7
unterschiedliche Symmetriegruppen bei vorgegebener Dimension. Dimensionsunabh¨angigkeit bis zur jeweiligen Smolyak Ordnung f¨ur verschiedene Symmetriegruppen. . . . . . . . . . . . . . . . . Anzahl der Basispolynome des Polynomraums d-dimensionaler symmetrischer Polynome beschr¨ankten totalen Grades. . . . . . . Spezifikationen des Quader Modells. . . . . . . . . . . . . . . . . Anzahlen der Kreise mit unterschiedlich vielen verschiedenen Parametern zu unterschiedlichen Smolyak Ordnungen. . . . . . . . . Anzahlen der Ellipsen mit unterschiedlich vielen verschiedenen Parametern zu unterschiedlichen Smolyak Ordnungen. . . . . . . dimensionsunabh¨angige Spezifikationen der Stoßlasche. . . . . .
v
18 19 21 57 62 63 64
1
Einleitung
Die numerische Simulation von Systemen der realen Welt unterliegt h¨aufig Unsicherheiten. So k¨onnen Modelle solcher Systeme mit Zufall behaftet sein. Als n¨utzliches Mittel bieten sich Statistikl¨osungen an, bei denen die Simulation mehrfach durchgef¨uhrt wird, und auf Basis der erhaltenen Ergebnisse beispielsweise Momente wie der Erwartungswert oder die Varianz berechnet werden. Dies f¨uhrt auf die L¨osung von Integralen. Da oftmals der analytische L¨osungsweg f¨ur Integrale zu aufwendig oder gar unm¨oglich ist, bieten sich numerische Quadraturverfahren zur Approximation von Integralen an. Allerdings haben Systeme der realen Welt mitunter eine Vielzahl an Freiheitsgraden. Da klassische numerische Quadraturverfahren dem Fluch der Dimension unterliegen, also mit anwachsender Dimensionszahl der Rechenaufwand exponentiell steigt, sind diese ungeeignet. Monte Carlo und Quasi-Monte Carlo Methoden sowie der Smolyak Algorithmus sind hingegen unabh¨angig oder fast unabh¨angig von der Dimensionszahl. Die Wahl einer der zuletzt genannten Quadraturverfahren h¨angt zum Einen von den Eigenschaften des Integranden, also der Funktionenklasse, zum Anderen von der Dimensionszahl ab. In dieser Diplomarbeit wird ein System zur statistischen Auswertung numerischer Simulationen implementiert. Das System besteht aus zwei Teilsystemen, den Quadraturalgorithmen und der objektiven Funktion. Die Teilsysteme sind voneinander weitestgehend entkoppelt, so dass eines dieser leicht ausgetauscht werden kann ohne das andere g¨anzlich oder auch nur teilweise neu implementieren zu m¨ussen. Die Quadraturalgorithmen liefern das Fundament zur Approximation von Integralen und der statistischen Auswertung, die objektive Funktion, verkn¨upft mit einer Statistikfunktion, den Integranden. Die objektive Funktion geht aus einem Experten-gest¨utzten Modell hervor. Genauer setzt sich die objektive Funktion aus der Modellbildung, der anschließenden Simulation des Modells und einer Bewertung des Simulationsergebnisses zusammen. Im Umfang dieser Arbeit ist die Implementierung der klassischen (naiven) Monte Carlo Methode, der Quasi-Monte Carlo Methode und des Smolyak Algorithmus eingeschlossen. Desweiteren wird ein modifizierter Smolyak Algorithmus implementiert. Dieser nutzt die Symmetrie einer symmetrischen objektiven Funktion aus. Dadurch m¨ussen mit unter bei weitem weniger Funktionswerte ausgewertet werden. Das Resultat ist ein erheblich geringerer Rechenaufwand. Voraussetzung hierf¨ur ist neben der symmetrischen Funktion ein genaues Wissen des Experten u¨ ber die Symmetrie. Als Beispiel einer objektiven Funktion wird die Simulation eines FiniteElemente Modells durch das Finite-Elemente Programm ParaFEP (Parallel Finite Element Program) betrachtet und als bewertendes Merkmal die maximale Verschiebung extrahiert. Hierbei gilt es die implementierten hochdimensionalen Quadraturverfahren und insbesondere den modifizierten Smolyak Algorithmus 1
auf Konvergenzverhalten zu untersuchen.
2
hochdimensionale Integration
Bei der numerischen Quadratur wird versucht den Wert eines Integrals n¨aherungsweise zu bestimmen. Der Rechenaufwand vieler klassischer numerischer Verfahren zur Approximation multivariater Integrale w¨achst exponentiell mit der Zahl der Dimension des Problems bei einer fest vorgegebenen Genauigkeit. Dies wird als Fluch der Dimension (engl.: curse of dimension) bezeichnet [5]. H¨aufig k¨onnen schon Probleme mit scheinbar geringer Dimension von beispielsweise 5 oder 10 als hochdimensionale Probleme bezeichnet werden, da diese mit einem hohen Rechenaufwand verbunden sind [12]. Numerische Quadraturverfahren wie Monte Carlo und Quasi-Monte Carlo Methoden sowie der Smolyak Algorithmus sind hingegen unabh¨angig bzw. bis auf einen logarithmischen Faktor unabh¨angig von der Dimensionszahl und wirken so dem Fluch der Dimension entgegen [6, 3]. Jede dieser Verfahren ist zum Einen f¨ur eine bestimmte Klasse von Funktionen optimiert [20, 3]. Zum Anderen eignen sie sich f¨ur unterschiedliche Dimensionszahlen [9, 20]. Experimente in [21] ergaben, dass sich bei geringer Dimension Tensorprodukt Quadratur, f¨ur Dimensionen von einigen Hundert Smolyak Quadratur, f¨ur sehr hohe Dimensionen Quasi-Monte Carlo Methoden und f¨ur extrem hohe Dimensionen Monte Carlo Methoden anbieten [9]. Das Integral einer Funktion f (x) der Funktionenklasse F u¨ ber dem ddimensionalen Raum R ⊆ Rd Z (1) Id f := f (~x)d~x R
wird u¨ ber eine numerische Quadraturformel (2)
Ql,d f :=
nl,d X
wl,i f (~xl,i)
i=1
approximiert, wobei der Funktionswert eines Integrationspunktes ~xl,i (mitunter auch als Sampling-Punkt, St¨utzstelle oder Knoten bezeichnet) mit dem Gewicht wl,i belegt wird. l ∈ N gibt die Ordnung (engl.: level) der Quadraturformel an und wird durch die Anzahl der verwendeten Integrationspunkte nl,d spezifiziert, so dass mit h¨oherer Ordnung auch eine h¨ohere Anzahl an Integrationspunkten verwendet wird, formeller ausgedr¨uckt nl,d < nl+1,d [5]. Eine Quadraturformel ist somit die Summe gewichteter Funktionswerte, ausgewertet an den Positionen ~xl,i mit i = 1, . . . , nl,d . 2
Die Menge aller nl,d Integrationspunkte Γl,d := {~xl,i |1 ≤ i ≤ nl,d } ⊂ R einer Quadraturformel wird als Gitter (engl.: grid) bezeichnet. Gilt f¨ur eine Quadraturformel, dass die Integrationspunkte der Ordnung l auch in l +1 vorkommen, so ist von geschachtelten Gittern (engl.: nested grids) die Rede. Eine Quadraturformel mit einem geschachtelten Gitter wird als geschachtelte Quadraturformel (engl.: nested quadrature formula) bezeichnet. Formeller ausgedr¨uckt, sind eine Quadraturformel und das zugrundeliegende Gitter geschachtelt, wenn Γl,d ⊂ Γl+1,d . [5]. Bei nicht-geschachtelten Quadraturformeln ist dies nicht notwendigerweise der Fall. Im Allgemeinen erwartet man, dass eine h¨ohere Anzahl von Knoten einen besseren Sch¨atzwert f¨ur Id f liefert, das heißt eine h¨ohere Genauigkeit als eine geringere Anzahl von Knoten. Dies ist allerdings h¨aufig nicht der Fall wie das Monte Carlo Beispiel 1 in [19] zeigt. Zur Absch¨atzung der Genauigkeit ist der Quadraturfehler (Integrationsfehler) [5] gegeben durch El,d f :=| Id f − Ql,d f | . Da nun mal typischerweise Id f gerade nicht bekannt ist (denn dies gilt es ja n¨aherungsweise zu bestimmen), muss auch hierf¨ur eine Absch¨atzung gefunden werden. Als eine M¨oglichkeit kann die Differenz zweier Quadraturformeln als Maß f¨ur den Fehler aufgefaßt werden | Id f − Ql,d f |≈| Qm,d f − Ql,d f | mit nm,d > nl,d [10].
2.1 Monte Carlo Methode Bei der klassischen (naiven) Monte Carlo (MC) Methode werden nd Integrationspunkte ~xi (pseudo-)zuf¨allig und unabh¨angig voneinander bez¨uglich der Gleichverteilung u¨ ber dem d-dimensionalen Raum R gew¨ahlt und der Durchschnittswert der zugeh¨origen Funktionswerte bestimmt. Der Durchschnittswert wird mit dem Raumvolumen Vd multipliziert, um so eine Approximation f¨ur Id f zu erhalten [16, 19]. So ergibt sich (3)
nd Vd X Qd := f (~xi ). nd i=1
3
F¨ur eine vorgegebene Genauigkeit Ed f bedarf es des Aufwandes −1/2
Ed f = O(nd
)
f¨ur Integranden beschr¨ankter Varianz [6]. Der bezeichnende Vorteil gegen¨uber der Quasi-Monte Carlo Methode (siehe Abschnitt 2.2) und dem Smolyak Algorithmus (siehe Abschnitt 2.3) ist die Unabh¨angigkeit des Fehlers von der Dimension. Zur Reduzierung des Integrationsfehlers um eine Dezimale muss allerdings die Anzahl der Integrationspunkte um den Faktor 100 erh¨oht werden [19]. Da (Pseudo)Zufallszahlen dazu neigen sich in Gebieten anzuh¨aufen bzw. genau dies nicht zu tun, so dass L¨ucken entstehen, eignet sich die Monte Carlo Methode f¨ur Funktionen hoher Varianz nicht. Allerdings bieten sich daher wiederum Monte Carlo Methoden f¨ur nicht glatte Integranden an [3]. Abbildung 1 zeigt ein Beispiel eines 2-dimensionalen Gitters der Monte Carlo Methode mit 64 Integrationspunkten. Um die langsame Konvergenz der Monte Carlo Methode zu beschleunigen, entstanden eine Reihe effizienter Techniken, wie unter anderem Importance Sampling oder Stratified Sampling [23, 19]. 0.9 ’MC Gitter’ 0.8
0.7
Koordinate y
0.6
0.5
0.4
0.3
0.2
0.1
0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
Abbildung 1: Beispiel eines 2D MC Gitters mit 64 Integrationspunkten.
4
1
2.2 Quasi-Monte Carlo Methode Quasi-Monte Carlo (QMC) Methoden unterscheiden sich von Monte Carlo Methoden lediglich in der Generierung der Integrationspunkte. Im Gegensatz zu (pseudo-)zuf¨alligen Punkten werden bei Quasi-Monte Carlo Methoden Folgen niedriger Diskrepanz (engl.: low discrepancy sequence) verwendet, um so eine m¨oglichst gleichm¨aßige Punktverteilung zu erreichen. Zur Kl¨arung dieses Begriffs gilt es weitere einzuf¨uhren und zu definieren. Gegeben seien die Integrationspunkte ~x1 , . . . , ~xn ∈ [0, 1]d einer QMC-Formel und der Quader Wt = {(u1 , . . . , ud) | 0 ≤ ui ≤ ti , i = 1, . . . , d} ⊆ [0, 1]d f¨ur ~t = (t1 , . . . , td ) ∈ [0, 1]d . So ergibt sich die Stern-Diskrepanz zu Dn∗ := D(~x1 , . . . , ~xn ) := sup | µ(Wt ) − t∈[0,1]d
1 X 1| n ~ xν ∈Wt
mit µ(Wt ) = t1 t2 · · · td
Eine Punktfolge im Rd ist von niedriger Diskrepanz, falls Dn∗ = O(n−1 lnd n) siehe [16]. Unter Hinzunahme der zweidimensionalen Halton-Punktfolge aus Abbildung 2 kann die Stern-Diskrepanz veranschaulicht werden. Man stelle sich auf jeder der beiden Achsen einen Schieberegler vor, u¨ ber die das Aufziehen jedes denkbaren Rechtecks mit linkem unteren Eckpunkt (0, 0) und rechtem oberen Eckpunkt von (0, 0) bis (1, 1) m¨oglich ist. F¨ur jedes dieser Rechtecke bestimmt man nun das Verh¨altnis zwischen der Anzahl der im Rechteck befindlichen Punkte und der Anzahl aller Punkte im Rechteck mit rechtem oberen Eckpunkt (1, 1). Dieses Verh¨altnis dient so als Maß f¨ur die Gleichm¨aßigkeit der Punkteverteilung. Die Stern-Diskrepanz liefert den schlechtesten Fall aus allen Verh¨altnissen. Somit gilt, je gleichm¨aßiger eine Zahlenfolge, desto geringer die Diskrepanz. Oder mit Bezug auf die MC Methode, je gleichm¨aßiger eine Punktverteilung, desto geringer die Wahrscheinlichkeit, dass L¨ucken oder Anh¨aufungen von Punkten in Teilgebieten des Raums, u¨ ber den integriert wird, entstehen. So eignet sich die QMC Methode auch f¨ur Funktionen hoher Varianz. F¨ur eine vorgegebene Genauigkeit Ed f bedarf es des Aufwandes d Ed f = O(n−1 d (log nd ) )
5
f¨ur Integranden mit beschr¨ankter Variation. Quasi-Monte Carlo Methoden sind nahezu unabh¨angig von der Dimension d bis auf einen logarithmischen Faktor, verbessern allerdings die Konvergenz [6]. Im Folgenden werden drei Punktfolgen niedriger Diskrepanz n¨aher beschrieben. 2.2.1 Halton-Folge Die Halton-Punktfolge ist durch ~xi = (φb1 (i), . . . , φbd (i)) ∈ [0, 1]d,
i = 1, 2, . . .
gegeben, wobei b1 , . . . , bd ∈ N>1 paarweise teilerfremde Zahlen sind. H¨aufig wird f¨ur bj die j-te Primzahl verwendet. Die Radix-inverse Funktion φ ist definiert als φb (i) = (0.d0 d1 . . . dk )b =
k X
dj b−1−j
j=0
mit i − 1 = (dk dk−1 . . . d1 d0 )b =
k X j=0
dj bj und dj ∈ 0, 1, . . . , b − 1.
F¨ur die Stern-Diskrepanz der Halton-Punktfolge gilt d
Dn∗
d 1Y ≤ + n n j=1
bj + 1 bj − 1 ln n + 2 ln bj 2
siehe [16]. Ein entscheidender Vorteil der Halton-Punktfolge gegen¨uber der Hammersley-Punktfolge (siehe Abschnitt 2.2.2) liegt in der Unabh¨angigkeit von der Anzahl der zu generierenden Punkte. So kann eine QMC-Approximation mittels der ersten n Halton-Punkte nachtr¨aglich um die n¨achsten m Halton-Punkte erweitert werden, falls das Ergebnis nicht zufriedenstellend ist. Abbildung 2 zeigt die ersten 64 zweidimensionalen Halton-Punkte mit den ersten zwei Primzahlen als Basis. 2.2.2 Hammersley-Folge Die Hammersley-Punktfolge ist durch i ~xi = ( , φb1 (i), . . . , φbd−1 (i)) ∈ [0, 1]d , n
6
i = 1, 2, . . . , n
1 QMC Gitter aus Halton-Punkten 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 2: QMC Gitter, die ersten 64 2D Halton-Punkte. gegeben, wobei n die Anzahl der zu generierenden Punkte und alles andere wie bei der Halton-Folge (siehe Abschnitt 2.2.1) ist. F¨ur die Stern-Diskrepanz der Hammersley-Punktfolge gilt Dn∗
d−1
d 1Y ≤ + n n j=1
bj − 1 bj + 1 ln n + 2 ln bj 2
siehe [7]. Die Hammersley-Punkte haben wie in Abschnitt 2.2.1 schon erw¨ahnt den Nachteil, dass sie von der Anzahl der zu generierenden Punkte abh¨angig sind. F¨uhrt die Anzahl n nicht zu einem zufriedenstellenden Ergebnis, so kann dieses Ergebnis nicht f¨ur eine h¨ohere Anzahl an Punkten weiter verwendet werden. Abbildung 3 zeigt die ersten 64 zweidimensionalen Hammersley-Punkte mit der ersten Primzahl als Basis. 2.2.3 arithmetische Folge Die arithmetische Punktfolge modulo 1 sei definiert durch ~xi = (αp1 (i), . . . , αpd (i)) ∈ [0, 1]d , 7
i = 1, 2, . . . .
1 QMC Gitter aus Hammersley-Punkten 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 3: QMC Gitter, die ersten 64 2D Hammersley-Punkte. Die Funktion αp : N → R ist gegeben durch αp (i) = i ∗ rp
mod 1,
wobei rp eine Nullstelle eines quadratischen Polynoms p mit ganzzahligen Koeffi¨ zienten ist. Uber die Kettenbruchentwicklung kann gezeigt werden, dass die Folge von niedriger Diskrepanz ist. Wie die Halton-Punktfolge ist auch diese Folge unabh¨angig von der Anzahl der zu generierenden Punkte. In dieser Diplomarbeit werden die Polynome pj (x) = x2 − bj ,
j ∈ {1, . . . , d} p mit der j-ten Primzahl bj und der Nullstelle rpj = bj verwendet. Abbildung 4 zeigt die ersten 64 zweidimensionalen Punkte dieser Folge.
2.3 Smolyak Algorithmus Bei der vollen Tensorprodukt Quadratur wird zur Approximation des multivariaten Integrals Id f eine eindimensionale Quadraturformel Ql,1 mit ein und derselben
8
1 ’QMC Gitter mit arithmetischer Punktfolge’ 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 4: QMC Gitter, die ersten 64 2D Punkte der arithmetischen Folge. Ordnung l in jeder Dimension verwendet Id f ≈ (Ql,1 ⊗ · · · ⊗ Ql,1 )f =
nl,1 X
i1 =1
···
nl,1 X
wl,i1 . . . wl,id f (xl,i1 , . . . , xl,id ).
id =1
Da in jeder Dimension nl,1 Knoten gew¨ahlt werden, ergeben sich insgesamt (nl,1 )d n¨otige Auswertungen der Funktion f . Dies unterliegt dem Fluch der Dimension [20]. Der Smolyak Algorithmus zur Approximation von Id f geht aus Tensorprodukten eindimensionaler Quadraturformeln der Ordnungen ki ∈ {1, . . . , lS } hervor. Ausgehend von einer solchen Reihe eindimensionaler Quadraturformeln f¨ur eine univariate Funktion f nki ,1
Qki ,1 :=
X
wki ,j f (xki ,j )
j=1
ergibt sich der Smolyak Algorithmus [17, 5] der Ordnung lS zu X d−1 lS +d−|~k|−1 (Qk1 ,1 ⊗ · · · ⊗ Qkd ,1 )f (4) QlS ,d f := (−1) |~k| − l lS ≤|~k|≤lS +d−1
9
mit dem Tensorprodukt der d eindimensionalen Quadraturformeln nkd ,1
nk1 ,1
(5)
(Qk1 ,1 ⊗ · · · ⊗ Qkd ,1 )f =
X
i1 =1
···
X
id =1
wk1 ,i1 · · · wkd ,id f (xk1 ,i1 , . . . , xkd ,id )
und dem Multi-Index ~k = [k1 , . . . , kd ] ∈ N+ zur Referenzierung der eindimensionalen Quadraturformeln unterschiedlicher Ordnungen ki ∈ {1, . . . , lS }, i = 1, . . . , d und d X ~ |k| := ki . i=1
Diese Schreibweise des Smolyak Algorithmus ist auch als Kombinationstechnik bekannt [5], eine alternative Schreibweise [5] lautet X QlS ,d f := (∆k1 ,1 ⊗ · · · ⊗ ∆kd ,1 )f |~k|≤lS +d−1
wobei ∆k,1 f := (Qk,1 − Qk−1,1 )f mit Q0,1 := 0. Der Vorteil des Smolyak Algorithmus ergibt sich aus der Nutzung von eindimensionalen Quadraturformeln h¨oherer Ordnungen in nur wenigen Dimensionen und niedriger Ordnungen in den restlichen Dimensionen [20]. Es wird nahezu Unabh¨angigkeit von der Dimensionszahl bis auf einen logarithmischen Faktor f¨ur r-fach stetig ableitbare Funktionen erreicht [20, 3] (d−1)(r+1) Ed f = O(n−r ). lS ,d (log nlS ,d )
Als eindimensionale Quadraturformeln finden geschachtelte Quadraturformeln wie die Clenshaw-Curtis Formel und die Kronrod-Patterson Formel [5] oder auch nicht-geschachtelte Quadraturformeln wie die Gauss(-Legendre) Formel und die Gauss-Hermite Formel [20] Anwendung. Geschachtelte eindimensionale Quadraturformeln haben u¨ blicherweise eine etwas schlechtere Integrationsordnung (engl.: polynomial degree of exactness), das heißt der Grad der Polynome, die exakt integriert werden, ist geringer als bei nicht-geschachtelten eindimensionalen Quadraturformeln [3], f¨ur n¨aheres siehe Abschnitte 2.3.1 bis 2.3.3. Allerdings m¨ussen nicht so viele Integrationspunkte tats¨achlich f¨ur die Approximation von Id f ausgewertet werden, da aufgrund der Schachtelung viele Punkte zusammen fallen. Geschachtelte eindimensionale Quadraturformeln k¨onnen daher als o¨ konomischer bezeichnet werden [3]. Da die eindimensionalen Integrationspunkte geschachtelt sind, sind es auch die d-dimensionalen Integrationspunkte. Anders gesprochen, das d-dimensionale Gitter ist ebenfalls geschachtelt. 10
Zusammenfassend kann also gesagt werden, dass die Qualit¨at des Smolyak Algorithmus und die Funktionenklasse, f¨ur die dieser optimiert ist, von der Wahl der eindimensionalen Quadraturformeln bzw. des verwendeten Gitters abh¨angt [3]. In den Abschnitten 2.3.1 bis 2.3.3 werden einige eindimensionale Quadraturformeln n¨aher erl¨autert. Symmetrien als auch die Vorteile deren Ausnutzung werden in den Abschnitten 2.3.4 und 2.3.5 beschrieben. 2.3.1 Gauss-Legendre Quadraturformel Die Gauss(-Legendre) Quadraturformel der Ordnung l, spezifiziert durch l Integrationspunkte, ist gegeben durch GLl :=
l X
wi f (xi ).
i=1
Die Integrationspunkte sind die Nullstellen des Legendre Polynoms Ll (x), die achsensymmetrisch zur Ordinate sind und im Intervall [−1, 1] liegen. Ll+1 (x) ist durch die Rekursionsformel (l + 1)Ll+1 (x) = (2l + 1)xLl (x) − lLl−1 (x) definiert mit L0 (x) = 1 und L1 (x) = x. Das Gewicht wi des Integrationspunktes xi ist 2 wi = , i = 1, . . . , l. 2 (1 − xi )L′l (xi )2 Die Integrationsordnung der Gauss-Legendre Quadraturformel betr¨agt 2l − 1, siehe [22, 23, 5]. Abbildung 5 zeigt das Gitter des Smolyak Algorithmus der Ordnung 9 f¨ur die Dimension 2 mit Gauss-Legendre Quadraturformel, skaliert auf [0, 1]2 . 2.3.2 Clenshaw-Curtis Quadraturformeln Die Clenshaw-Curtis Quadraturformeln der Ordnung l, spezifiziert durch l Integrationspunkte sind definiert durch CCl :=
l X
wi f (xi ).
i=1
Die Integrationspunkte sind die l Extremstellen des Tschebyschew Polynoms Tl−1 der ersten Art, liegen in [−1, 1] und sind durch π(i − 1) , i = 1, ..., l xi = − cos l−1 11
1 ’Smolyak Algorithmus mit Gauss-Legendre Quadraturformel’ 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 5: Gitter aus 281 Integrationspunkten des Smolyak Algorithmus der Ordnung 9 f¨ur die Dimension 2 mit Gauss-Legendre Quadraturformel, skaliert auf [0, 1]2 .
12
gegeben. Das Gewicht wi des Integrationspunktes xi berechnet sich zu 1 f¨ur i ∈ {1, l}, l gerade (l−1)2 1 f¨ur i ∈ {1, l}, l ungerade wi = l(l−1) l−1 X ∗[ 2 ] 2 2 l−1 1− cos 2k i−1 π sonst, 4k 2 −1 l−1 k=1
X∗
wobei bedeutet, dass der letzte Summand halbiert werden muss, falls l ungerade ist, siehe [1, 4, 16]. 1 ’Smolyak Algorithmus mit Clenshaw-Curtis Quadraturformel’
Koordinate y
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Koordinate x
Abbildung 6: Gitter aus 321 Integrationspunkten des Smolyak Algorithmus der Ordnung 7 f¨ur die Dimension 2 mit Clenshaw-Curtis Quadraturformel, skaliert auf [0, 1]2 .
Die in dieser Diplomarbeit verwendete Clenshaw-Curtis Quadraturformel ist geschachtelt und ergibt sich aus den ungeraden Ordnungen l = 2n−1 + 1, n ∈ N>1 mit der Quadraturformel CC1 := 2f (0) und einer Integrationsordnung von l − 1, siehe [16, 5]. Abbildung 6 zeigt das Gitter des Smolyak Algorithmus der Ordnung 7 f¨ur die Dimension 2 mit geschachtelter Clenshaw-Curtis Quadraturformel, skaliert auf [0, 1]2 .
13
2.3.3 Kronrod-Patterson Quadraturformel Eine Gauss-Kronrod Quadraturformel erweitert eine Gauss Quadraturformel mit n Integrationspunkten um weitere n + 1 Integrationspunkte, die miteinander koinzidieren, so dass eine Integrationsordnung von 2n + n ¯ + 2 mit ( n f¨ur n ungerade n ¯ := n − 1 sonst. erhalten wird. Die Kronrod-Patterson Quadraturformel ist eine geschachtelte Quadraturformel dieser Form, siehe [10, 5]. Die verz¨ogerte Kronrod-Patterson Quadraturformel aus [14] ist gegeben durch del KP(i) = KPm ,
3m + 3 3m + 4 ≤i≤ , 8 4
m = 1, 3, 7, 15, 31, . . . ,
wobei die Ordnung (i) ∈ N+ nicht mit der Anzahl der Integrationspunkte identisch ist. F¨ur die ersten Ordnungen ergibt sich so del KP(1) = KP1 ,
del del KP(2) = KP(3) = KP3 ,
del del del KP(4) = KP(5) = KP(6) = KP7 , . . .
Die Abbildungen 7 und 8 zeigen Gitter des Smolyak Algorithmus mit KronrodPatterson und verz¨ogerter Kronrod-Patterson Quadraturformel f¨ur die Dimension 2, skaliert auf [0, 1]2 . 2.3.4 Symmetrien Eine Funktion f : Rd → R ist symmetrisch in den Argumenten, falls (6)
f (x1 , . . . , xd ) = f (xp(1) , . . . , xp(d) )
f¨ur alle bijektiven Funktionen (Permutationen) p : {1, . . . , d} → {1, . . . , d}. Implementiert sind Symmetrien der Form f (~x1 , . . . , ~xm ) = f (~xp(1) , . . . , ~xp(m) ) f¨ur alle bijektiven Funktionen p : {1, . . . , m} → {1, . . . , m}, wobei die Komponenten von ~x ∈ Rd zu m ∈ N+ Bl¨ocken ~xi gleicher Gr¨oße s ∈ N+ mit m ∗ s = d zusammengefaßt sind, so dass f (~x) = f (x1 , . . . , xd ) = f (~x1 , . . . , ~xm )
14
1 ’Smolyak Algorithmus mit Kronrod-Patterson Quadraturformel’ 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 7: Gitter aus 321 Integrationspunkten des Smolyak Algorithmus der Ordnung 6 f¨ur die Dimension 2 mit Kronrod-Patterson Quadraturformel, skaliert auf [0, 1]2 .
15
1 ’Smolyak Algorithmus mit verz. Kronrod-Patterson Quadraturformel’ 0.9 0.8
Koordinate y
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1
0.2
0.3
0.4
0.5 Koordinate x
0.6
0.7
0.8
0.9
1
Abbildung 8: Gitter aus 321 Integrationspunkten des Smolyak Algorithmus der Ordnung 15 f¨ur die Dimension 2 mit verz¨ogerter Kronrod-Patterson Quadraturformel, skaliert auf [0, 1]2 .
16
mit ~xi = (xs∗(i−1)+1 , . . . , xs∗(i−1)+s ) f¨ur i ∈ {1, . . . , m}. Im Falle solch einer symmetrischen Funktion werden in dieser Arbeit die Bl¨ocke ~xi als Symmetriebl¨ocke bezeichnet. Dabei spezifiziert die Anzahl r der Komponenten eines Symmetrieblockes die vorliegende Symmetriegruppe Sr . F¨ur r = 1 ergibt sich die in Gleichung 6 beschriebene vollst¨andige Symmetriegruppe. Zwei Symmetriengruppen Sr1 und Sr2 sind unterschiedlich, falls r1 6= r2 ist. Ist allerdings r1 = k ∗ r2 mit k ∈ N+ , so ist Sr1 ⊆ Sr2 . 2.3.5 Vorteile der Symmetrieausnutzung Durch Symmetrieausnutzung m¨ussen je nach Symmetriegruppe (siehe Abschnitt 2.3.4) durch den Smolyak Algorithmus wesentlich weniger Integrationspunkte ausgewertet werden. Tabelle 1 zeigt dies anhand der Clenshaw-Curtis Quadraturformel bei vorgegebener Dimension f¨ur unterschiedliche Symmetriegruppen. F¨ur eine vorgegebene Symmetriegruppe Sr , stellt sich sogar f¨ur gen¨ugend große Dimensionen eine Dimensionsunabh¨angigkeit bis zu einer jeweiligen Ordnung des Smolyak Algorithmus ein, siehe Tabelle 2. F¨ur eine Dimension di liefert eine markierte Ordnung in der Tabelle die gleiche Anzahl an Integrationspunkten wie die gleiche Ordnung f¨ur die in den Tabellen vorangegangene Dimension di−1 . Desweiteren sind die einzelnen Komponenten der di−1 -dimensionalen Integrationspunkte mit den di−1 ersten bzw. letzten Komponenten (je nach Sichtweise) der di -dimensionalen Integrationspunkte identisch. Dieser Effekt soll nun anhand der vollst¨andigen Symmetriegruppe und symmetrischer Polynome mit beschr¨anktem totalen Grad mathematisch fassbar gemacht werden. Eine Basis des Polynomraumes P≤p d-dimensionaler Polynome mit bed schr¨anktem Grad ≤ p ergibt sich aus dem Tensorprodukt d ∗ (p + 1) eindimensionaler Basispolynome zu ( d ) Y ≤p ei Bd := xi | ei ≤ p, ei ∈ N . i=1
≤p Somit spannen die (p + 1)d Basispolynome aus B≤p d den Polynomraum Pd auf. Anders gesagt, gehen die Polynome aus P≤p d durch Linearkombinationen der Ba≤p d sispolynome hervor. | Bd |= (p + 1) bezeichne die Dimension von P≤p d . Die Basispolynome und daher auch die Polynome aus P≤p k¨ o nnen durch die volle d Tensorprodukt Quadratur mittels eindimensionaler Quadraturformel der Integrationsordnung p exakt integriert werden. Abbildung 9 zeigt das Tensorprodukt der eindimensionalen Basispolynome einer Basis des Polynomraumes P≤p 2 . ≤p ¯ Desweiteren ist eine Basis des Polynomraumes Pd d-dimensionaler Polyno-
17
Dimension 10 keine Symmetrie Ordnung #Auswertungen 1 1 2 21 3 221 4 1581 5 8801 6 41265 Symmetriegruppe S5 Ordnung #Auswertungen 1 1 2 11 3 116 4 796 5 4431 6 20663 Symmetriegruppe S2 Ordnung #Auswertungen 1 1 2 5 3 23 4 91 5 342 6 1214 Symmetriegruppe S1 Ordnung #Auswertungen 1 1 2 3 3 8 4 20 5 50 6 122 Tabelle 1: Anzahl der Funktionsauswertungen f¨ur unterschiedliche Smolyak Ordnungen und Symmetriegruppen bei vorgegebener Dimension 10 basierend auf der Clenshaw-Curtis Quadraturformel.
18
Symmetriegruppe S5 Dimension 10 Ordnung #Auswertungen 1 1 2 11 3 116 4 796 5 4431 6 20663 Dimension 20 Ordnung #Auswertungen 1 1 2 11 3 116 4 1016 5 8116 6 57998 Dimension 30 Ordnung #Auswertungen 1 1 2 11 3 116 4 1016 5 8116 6 60000
Symmetriegruppe S2 Dimension 10 Ordnung #Auswertungen 1 1 2 5 3 23 4 91 5 342 6 1214 7 4086 8 13222 9 41446 Dimension 20 Ordnung #Auswertungen 1 1 2 5 3 23 4 91 5 342 6 1214 7 4086 8 13874 9 45091 (b) Symmetriegruppe S2
(a) Symmetriegruppe S5
Tabelle 2: Dimensionsunabh¨angigkeit bis zur jeweiligen Smolyak Ordnung f¨ur (a) die Symmetriegruppe S5 und (b) die Symmetriegruppe S2 .
19
.. .
.. .
.
.. .
x1 xp2 x21 xp2 · · · xp1 xp2
..
xp2
.. .
x22
x1 x22 x21 x22 · · · xp1 x22
x2
x1 x2 x21 x2 · · · xp1 x2
1
x21
x1
···
xp1
Abbildung 9: Tensorprodukt der eindimensionalen Basispolynome einer Basis des Polynomraumes P≤p 2 .
me mit beschr¨anktem totalen Grad ≤ p gegeben durch ( d ) d Y X ¯ ≤p := B xei | ei ≤ p, ei ∈ N i
d
i=1
i=1
p+d
¯ ≤p ⊆ B≤p . Der Smolyak Algorithmus der Ordmit d Basispolynomen und B d d ¯ ≤2l+1 exakt, falls f¨ur die i-te eindimensionung l integriert die Polynome aus P d nale Quadraturformel, der durch den Smolyak Algorithmus genutzten Folge l eindimensionaler Quadraturformeln, die Integrationsordnung ≥ 2i − 1 ist mit ¯ ≤p f¨ur den i ∈ {1, . . . l} [14]. Abbildung 10 zeigt die Konstruktion der Basis B d ¯ ≤p . Polynomraum P 2 ˜ ≤p d-dimensionaler symmetrischer PolynoEine Basis des Polynomraumes P d me mit beschr¨anktem totalen Grad ≤ p ist definiert durch ( ) d d XY e φi X ≤p ˜ := B ei ≤ p, ed ≤ · · · ≤ e1 x | i
d
φ∈Pd i=1
i=1
˜ ≤p ⊆ B ¯ ≤p ⊆ B≤p . mit der Menge Pd der Permutationen auf {1, 2, . . . , d} und B d d d Nun gilt es die Anzahl der Basispolynome zu bestimmen, also die Dimension ˜ ≤p . F¨ur d ≥ p ist die Anzahl des durch diese aufgespannten Polynomraumes P d der Basispolynome wegen der Beschr¨anktheits- und Symmetrieeigenschaften unabh¨angig von der Dimension d. Dies soll durch Ausz¨ahlen der Basispolynome anhand des beschr¨ankten totalen Grades p ≤ 4 veranschaulicht werden. Hierzu 20
xp2 xp−1 x1 xp−1 2 2
x2
x1 x2
1
x1
..
.
.
.. .
..
.. .
· · · xp−1 1 x2 ···
xp−1 xp1 1
¯ ≤p . Abbildung 10: Konstruktion der Basis B¯ ≤p ur den Polynomraum P 2 d f¨
p 10 20 30 40 50 60 70 80 90 100
˜ ≤p | |B d 138 2713 28628 215307 1295970 6639348 30053953 123223638 465672548 1642992567
Tabelle 3: Anzahl der Basispolynome des Polynomraums d-dimensionaler symmetrischer Polynome beschr¨ankten totalen Grades mit d ≥ p.
21
werden alle den beiden Eigenschaften gen¨ugenden Folgen von d Exponenten f¨ur die totalen Grade pt ∈ {1, 2, 3, 4} betrachtet. So ergeben sich aus pt = 1 = |1 + 0 + 0 + 0{z+ 0 + · · · + 0} d Exponenten
pt = 2 = 2 + 0 + 0 + 0 + 0 + · · · + 0 = 1+1+0+0+0+···+0
pt = 3 = 3 + 0 + 0 + 0 + 0 + · · · + 0 = 2+1+0+0+0+···+0
= 1+1+1+0+0+···+0
pt = 4 = 4 + 0 + 0 + 0 + 0 + · · · + 0 = 3+1+0+0+0+···+0
= 2+2+0+0+0+···+0
= 2+1+1+0+0+···+0
= 1+1+1+1+0+···+0 ˜ ≤4 . Eine Dimension d ≥ p = 4 a¨ ußert sich 1 + 2 + 3 + 5 = 11 Basispolynome in B d lediglich in weiteren angeh¨angten Nullen, so dass die Anzahl der Basispolynome gleich bleibt. Tabelle 3 zeigt die Anzahl der Basispolynome f¨ur ausgew¨ahlte totale Grade ≤ p. ¯ ≤p ˜ ≤p In Abbildung 11 sind die Dimensionen der Polynomr¨aume P≤p 50 , P50 und P50 u¨ ber dem Grad p doppelt logarithmisch aufgetragen.
2.4 Experimente Die vorgestellten Quadraturverfahren sollen nun hinsichtlich einer ausgew¨ahlten glatten und vollst¨andig symmetrischen Funktion (siehe Abschnitt 2.3.4) auf Konvergenz gepr¨uft werden. Der Integrand f ist eine auf der Ordinate verschobene Sinus-basierte Funktion variierender Frequenz, die St¨orungen unterliegt. f ist definiert als f (~x) = f (x1 , . . . , xd ) =
d Y
(t(xi ) + s(xi ))
i=1
mit der Tr¨agerfunktion 3 x2 1 + (1 − x) sin nπ x − t(x) = 2 2 2 22
, n ∈ N+
1e+90 Polynome beschraenkten Grades Polynome beschraenkten totalen Grades symmetrische Polynome beschraenkten totalen Grades
1e+80
Dimension des Polynomraumes
1e+70
1e+60
1e+50
1e+40
1e+30
1e+20
1e+10
1 1
10 Grad p
100
≤p ¯ ≤p ˜ ≤p u¨ ber Grad p. Abbildung 11: Dimension der Polynomr¨aume P50 , P50 und P 50
und der St¨orfunktion mπ s(x) = ε sin x 2
, m ∈ N+ , ε ∈ R+ .
Wegen der Unabh¨angigkeit der d Variablen xi voneinander, gilt f¨ur die Stammfunktion F von f Z Z F (~x) = f (~x) d~x = f (x1 , . . . , xd ) dx1 . . .dxd =
Z Y d
(t(xi ) + s(xi )) dx1 . . .dxd =
i=1
d Z Y
(t(xi ) + s(xi )) dxi
i=1
=
d Z Y
t(xi ) dxi +
i=1
Z
s(xi ) dxi .
Die Stammfunktion T von t ergibt sich nach Anwendung der Substitutionsregel zu Z 3 x2 1 + x. cos nπ x − T (x) = t(x) dx = − 2nπ 2 2 23
2 f 1.9
1.8
1.7
1.6
1.5
1.4
1.3
1.2
1.1 0
0.2
0.4
0.6
0.8
1
Abbildung 12: Funktion f mit n = 5, m = 100, ε = 0.05 und d = 1. Desweiteren ist die Stammfunktion S von s Z mπ 2ε cos x . S(x) = s(x) dx = − mπ 2 Das Integral von t u¨ ber dem Intervall [0, 1] lautet Z 1 nπ 3 1 1 I[0,1] t = t(x) dx = [T (x)]0 = 1 − cos + . 2nπ 2 2 0 Zudem ist das Integral von s u¨ ber dem Intervall [0, 1] Z 1 mπ 2ε I[0,1] s = s(x) dx = [S(x)]10 = 1 − cos . mπ 2 0 Abbildungen 12 und 13 zeigen die Funktion f mit n = 5, m = 100 und ε = 0.05 f¨ur d = 1 und d = 2. In Abbildungen 14 und 15 ist der relative Fehler der Integralsapproximation u¨ ber dem Raum [0, 1]d durch das jeweilige Quadraturverfahren u¨ ber der Anzahl der Integrationspunkte f¨ur diese Funktion in den Dimensionen 10 und 20 dargestellt. Der Smolyak Algorithmus konvergiert unabh¨angig von der verwendeten eindimensionalen Quadraturformel am 24
f
4 3.5 3 2.5 2 1.5 1 0 0.2 0
0.4 0.2
0.4
0.6 0.6
0.8 0.8
1 1
Abbildung 13: Funktion f mit n = 5, m = 100, ε = 0.05 und d = 2. schlechtesten. Die Konvergenzeigenschaften verschlechtern sich bei Dimensionszuwachs zunehmend. Die Monte Carlo Methode zeigt gutes Konvergenzverhalten. Die Quasi-Monte Carlo Methoden konvergieren am schnellsten. Im Folgenden soll das Konvergenzverhalten des Smolyak Algorithmus bei unterschiedlicher Ausnutzung der Symmetrie betrachtet werden. Als eindimensionale Quadraturformel wird die Clenshaw-Curtis Quadraturformel genutzt. Abbildung 16 zeigt den erheblichen Vorteil der Symmetrieausnutzung, insbesondere bei vollst¨andiger Symmetriegruppe. Das Konvergenzverhalten des Smolyak Algorithmus hat sich entschieden verbessert, kann f¨ur die Funktion f die Quasi-Monte Carlo Formeln allerdings nicht u¨ berbieten.
3
System zur statistischen Auswertung numerischer Simulationen
Die statistische Auswertung numerischer Simulationen f¨uhrt auf das L¨osen von Integralen. In dieser Diplomarbeit werden Integrale u¨ ber dem Raum Ω = [0, 1]d
25
1 Smolyak mit Clenshaw-Curtis Smolyak mit Gauss-Legendre Smolyak mit Kronrod-Patterson Smolyak mit verz. Kronrod-Patterson MC QMC mit Halton-Folge QMC mit Hammersley-Folge QMC mit arithmetischer Folge
0.1
0.01
relativer Fehler
0.001
0.0001
1e-05
1e-06
1e-07
1e-08 10
100
1000 10000 #Funktionsauswertungen
100000
1e+06
Abbildung 14: relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 10.
26
10 Smolyak mit Clenshaw-Curtis Smolyak mit Gauss-Legendre Smolyak mit Kronrod-Patterson Smolyak mit verz. Kronrod-Patterson MC QMC mit Halton-Folge QMC mit Hammersley-Folge QMC mit arithmetischer Folge
1
0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06
1e-07 10
100
1000
10000 100000 #Funktionsauswertungen
1e+06
1e+07
Abbildung 15: relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 20.
27
1 Smolyak mit Clenshaw-Curtis ohne Symmetrie Smolyak mit Clenshaw-Curtis und S_5 Smolyak mit Clenshaw-Curtis und S_2 Smolyak mit Clenshaw-Curtis und S_1 MC QMC mit Halton-Folge
0.1
0.01
relativer Fehler
0.001
0.0001
1e-05
1e-06
1e-07
1e-08 1
10
100
1000
10000
100000
1e+06
#Funktionsauswertungen
Abbildung 16: relativer Fehler u¨ ber der Anzahl der Funktionsauswertungen f¨ur die Dimension 20 mit Symmetrieausnutzung durch den Smolyak Algorithmus mit den Symmetriegruppen S5 , S2 und S1 .
28
der Form m(f ) Id
(7)
:=
Z
m(f (~ω ))d~ω
Ω
bez¨uglich des Lebesgue-Maßes betrachtet. Die Funktion f : Ω → R wird als objektive Funktion bezeichnet, m : R → R als Statistikfunktion. Zur Berechnung des Erwartungswertes ist m(x) = x wegen der Wahrscheinlichkeitsdichtefunktion fDichte = 1. F¨ur das 2. Moment gilt m(x) = x2 . Beliebige Verteilungen k¨onnen durch Skalierung der Komponenten der Integrationspunkte u¨ ber inverse Verteilungsfunktionen erreicht werden. Hierzu im Verlauf dieses Kapitels n¨aheres. Mit Gleichung 7 als Basis besteht das in dieser Diplomarbeit zu entwickelnde System zur statistischen Auswertung numerischer Simulationen aus zwei Teilsystemen, den Quadraturalgorithmen und der objektiven Funktion, siehe Abbildung 17. Die Quadraturalgorithmen liefern das mathematische Fundament zur Approximation von Integralen und zur statistischen Auswertung. Die objektive Funktion f , verkn¨upft mit der Statistikfunktion m, bildet den Integranden. f geht aus einem Experten-gest¨utzten Modell hervor. Genauer umfasst f die Modellbildung, die anschließende Simulation des Modells und die Bewertung des Simulationsergebnisses. Die Bewertung extrahiert die interessierende Kenngr¨oße aus dem Simulationsergebnis, bricht das Ergebnis somit auf eine Dimension herunter. Es sei erw¨ahnt, dass auch n-wertige objektive Funktionen generell unterst¨utzt werden. Desweiteren ist die effiziente Auswertung mehrerer statistischer Funktionen m¨oglich, so dass auch Integrale der Form Z m( ~ f~) ~ ω ))d~ω := m( ~ f(~ ∀m ~ ∈ {m ~ 1, · · · , m ~ k} Id Ω
durch das in dieser Diplomarbeit entwickelte System approximiert werden k¨onnen. Dabei ist f~ : Ω → Rn . m ~ ist eine n-wertige Statistikfunktion, die in jeder der n Dimensionen die gleiche eindimensionale und einwertige Statistikfunktion m auswertet. Zum einfacheren Verst¨andnis wird allerdings im Folgenden lediglich der Fall einer einwertigen objektiven Funktion und genau einer Statistikfunktion betrachtet. Die beiden Teilsysteme sind weitestgehend voneinander entkoppelt, so dass eines dieser leicht ausgetaucht werden kann ohne das andere g¨anzlich oder auch ¨ nur teilweise neu implementieren zu m¨ussen. Ublicherweise wird die objektive Funktion ausgetauscht. Im Folgenden sollen die beiden Teilsysteme und deren Zusammenspiel n¨aher erl¨autert werden.
3.1 Teilsystem Quadraturalgorithmen Wie schon erw¨ahnt, liefert das Teilsystem Quadraturalgorithmen das Fundament zur Approximation von Integralen und zur statistischen Auswertung. 29
Gesamtsystem Quadraturalgorithmen Statistikauswertung
numerische Quadratur
f (~ω )
~ω
objektive Funktion
Abbildung 17: Gesamtsystem, Perspektive 1.
Als numerische Quadraturverfahren sind die klassische (naive) Monte Carlo Methode, Quasi-Monte Carlo Formeln, der klassische Smolyak Algorithmus und ein modifizierter Smolyak Algorithmus zur Ausnutzung von Symmetrien implementiert. Die Quasi-Monte Carlo Formeln basieren auf der Halton-, der Hammersley- und der arithmetischen Punktfolge aus Abschnitt 2.2.3. Als eindimensionale Quadraturformeln f¨ur die beiden Smolyak Algorithmen sind die Gauss-Legendre und die Clenshaw-Curtis Quadraturformel implementiert. Die Kronrod-Patterson Quadraturformel wird lediglich bis zur Ordnung 6 und die verz¨ogerte Kronrod-Patterson Quadraturformel bis zur Ordnung 48 angeboten. So stehen neben nicht-geschachtelten eindimensionalen Quadraturformeln auch geschachtelte zur Verf¨ugung. Eine Erweiterung oder Anpassung dieses Teilsystems wie beispielsweise das Einbinden neuer Punktgeneratoren oder eindimensionaler Quadraturverfahren ist komfortabel gestaltet worden, siehe hierzu Abschnitt 5. Die mathematischen Grundlagen dieser Diplomarbeit sind in Abschnitt 2 zu finden.
3.2 Teilsystem objektive Funktion Wie schon angef¨uhrt geht die objektive Funktion f aus einem Experten-gest¨utzten Modell hervor. f umfasst genauer die Modellbildung, die Simulation des Modells 30
Semantik
~ω
~e(~ω )
Algorithmus
objektive Funktion M(~x, ~e) Simulation
f (~ω)
u(~x, M) Bewertung
Abbildung 18: Gesamtsystem, Perspektive 2.
und die Bewertung des Simulationsergebnisses. Eine andere Perspektive auf das Gesamtsystem zur statistischen Auswertung numerischer Simulationen hebt die mathematische Sicht auf das System und insbesondere die objektive Funktion st¨arker hervor, siehe Abbildung 18. Der Algorithmus ist ein Synonym f¨ur das Teilsystem Quadraturalgorithmen. Ihm obliegt sowohl die numerische Quadratur als auch die statistische Auswertung. Die objektive Funktion f wird u¨ ber die Funktionenkomposition der Funktion Semantik, der Simulation und der Bewertung gebildet, also f = Bewertung(Simulation(Semantik)). Verkn¨upft mit der Statistikfunktion m ergibt sich der Integrand m(f ) aus Gleichung 7. Die Semantik-Funktion transformiert die durch den Algorithmus generierten Punkte aus Ω in ein f¨ur die Simulation verst¨andliches Format. Abstrakt gesprochen, verleiht sie den Integrationspunkten Bedeutung. Genauer erzeugt oder vervollst¨andigt die Semantik-Funktion aus einem Integrationspunkt ~ω das zu simulierende Modell M. Dazu bedarf es typischerweise Skalierungsfunktionen, die ~ω in ~e ∈ E innerhalb der Semantik-Funktion transformieren. E ist dabei ein vom Experten gew¨unschter Raum. Diese Skalierungsfunktionen sind die entsprechenden inversen Verteilungsfunktionen, siehe hierzu Abbildung 19. Jede Komponente von ~ω kann durch eine andere inverse Verteilungsfunktion skaliert werden. Somit k¨onnen auch die Intervalle, in denen die einzelnen skalierten Komponenten liegen, unterschiedlich sein. Bei Ausnutzung von Symmetrien m¨ussen allerdings die inversen Verteilungsfunktionen f¨ur vertauschbare Komponenten gleich sein, siehe hierzu Abschnitt 3.4. Das durch die Semantik generierte vollst¨andige Modell M wird durch die Simulation berechnet. Aus dem Ergebnis der Simulation extrahiert die Bewertung die f¨ur den Experten interessante Kenngr¨oße u und kann daher auch als function of interest bezeichnet werden. 31
Semantik
Ω ∋ ~ω =
ω1
. . . ωd
Skalierungsfunktion 1
. . .
Modellbildung/ -vervollst¨andigung
Skalierungsfunktion d
Abbildung 19: Semantik-Funktion.
3.3 Interaktion zwischen den Teilsystemen Das Teilsystem Quadraturalgorithmen generiert Integrationspunkte aus Ω. Diese werden zur Auswertung an die objektive Funktion f u¨ bergeben. Der Funktionswert f (~ω ) zu einem Integrationspunkt ~ω wird von der objektiven Funktion an den Quadraturalgorithmus zur¨uckgegeben und intern der auszuwertenden Statistikfunktion m u¨ bergeben. Wurden alle Integrationspunkte auf diese Weise ausgewertet, ist die Approximation des gew¨unschten Integrals abgeschlossen, siehe Abbildung 17. Genauer erh¨alt die Semantik-Funktion der objektiven Funktion die durch das Teilsystem Quadraturalgorithmen generierten Integrationspunkte. Die SemantikFunktion skaliert einen Intergrationspunkt ~ω in ~e. Aus ~e erzeugt oder vervollst¨andigt sie das zu simulierende Modell M. Im Anschluss an die Simulation von M extrahiert die Bewertung die f¨ur den Experten interessante Kenngr¨oße u aus dem Simulationsergebnis und gibt diesen als Funktionswert f (~ω ) an den Quadraturalgorithmus zur¨uck. Somit stellen die Semantik und die Bewertung die Schnittstellen zwischen den beiden Teilsystemen dar, siehe Abbildung 18.
3.4 Symmetrieausnutzung Zur Ausnutzung von Symmetrien bedarf es zweier Voraussetzungen. Zum Einen muss die objektive Funktion symmetrisch sein, zum Anderen muss der Experte 32
dem Quadraturalgorithmus die Symmetrie bekannt geben, siehe Abschnitt 2.3.4. Das in dieser Diplomarbeit implementierte System zur statistischen Auswertung numerischer Simulationen unterst¨utzt die Symmetrienausnutzung u¨ ber einen modifizierten Smolyak Algorithmus, siehe Abschnitt 2.3.5. Dabei muss zur Ansteuerung der Symmetrie dem modifizierten Smolyak Algorithmus die Symmetriegruppe bekannt gegeben werden. Zudem muss, wie oben erw¨ahnt, sichergestellt sein, dass die objektive Funktion von dieser Symmetrie ist, indem der Experte die Semantik-Funktion entsprechend definiert.
4
Beispielmodell und dessen Auswertung
Anhand des in Abschnitt 3 vorgestellten Systems zur statistischen Auswertung numerischer Simulationen soll nun ein Beispiel f¨ur eine objektive Funktion betrachtet werden. Die objektive Funktion umfasst dabei die Modellvervollst¨andigung, Simulation und Bewertung eines Finite-Elemente (FE) Modells. Nach Abbildung 18 ergeben sich folgende Spezifikationen. Das durch den Experten ausgew¨ahlte Quadraturverfahren des Teilsystems Quadraturalgorithmen (in der Abbildung als Algorithmus bezeichnet) generiert Integrationspunkte aus Ω. Die Semantik-Funktion der objektiven Funktion f transformiert aus einem Integrationspunkt ~ω die r¨aumliche Materialverteilung des FE Modells. Nach der Simulation mittels eines FE Programms wird die interessierende Kenngr¨oße des Simulationsergebnisses extrahiert und als Funktionswert f (~ω ) an das Quadraturverfahren zur¨uckgeliefert. Das Quadraturverfahren u¨ bergibt f (~ω ) intern der Statistikfunktion. In den Abschnitten 4.1 und 4.2 sollen auf die Spezifikationen im Einzelnen eingegangen und kurz einige Alternativen vorgestellt werden. Im Anschluss soll in Abschnitt 4.3 dieses Modell weiter konkretisiert werden.
4.1 Semantik Das aus einem Integrationspunkt ~ω durch die Semantik-Funktion generierte Material soll durch eine vorgegebene Anzahl von Einschl¨ussen a¨ hnlicher Form (beispielsweise Ellipsen oder Rechtecke bei einem 2D FE Model) definiert sein. Somit wird das Material durch die unterschiedlichen Materialien von Matrix und Einschluss sowie der Lage und Orientierung der in diesem zu findenen Einschl¨usse beschrieben. Um n Einschl¨usse mit jeweils m Parametern unabh¨angig voneinander generieren zu k¨onnen, werden n ∗ m Parameter ben¨otigt. Somit ergibt sich die Dimension der Integrationspunkte zu d = n ∗ m. 33
Durch die inversen Verteilungsfunktionen werden nun die einzelnen Komponenten von ~ω = (ω1 , . . . , ωd )T in die Parameter der n Einschl¨usse skaliert. Dies geschieht, indem jeweils m aufeinanderfolgende Komponenten von ~ω stets in der gleichen Reihenfolge auf die m Parameter eines Einschlusses abgebildet werden. So liefern die Komponenten ωi, . . . , ωi+m−1 mit i ∈ N+ und (i − 1) mod m = 0 die Parameter pei ,1 , . . . , pei,m des i-ten Einschlusses ei . Jede j-te Komponente eines solchen Blockes von m Komponenten wird dabei mittels derselben inversen Verteilungsfunktion skaliert, siehe Abbildung 20. Auf diese Weise ist sichergestellt, dass die Einschl¨usse unabh¨angig voneinander generiert werden. Allerdings muss auch erw¨ahnt werden, dass sich dadurch Einschl¨usse u¨ berlappen k¨onnen. Dies sei aber in dieser Diplomarbeit zugelassen. Zudem m¨ussen die inversen Verteilungsfunktionen so gew¨ahlt werden, dass die Parameter innerhalb f¨ur das FE Modell sinnvoller Intervalle liegen, denn an den Positionen des FE Modells, an denen das FE Programm das Material abfragt, muss das Material auch definiert sein. Die generierten Einschl¨usse und die vorgegebenen Materialparameter von Matrix und Einschluss bilden so das Material des FE Modells, siehe Abbildung 21. F¨ur die sp¨atere Simulation bedarf es zudem einer Materialfunktion (8)
M : Rk × Ω → Rl ,
die zu einem Punkt ~x = (x1 , · · · , xk ) des k-dimensionalen FE Modells mit k ∈ {2, 3} das Material M(~x, ~ω ), bestehend aus l Parametern, liefert. Genauer liefert M das Material f¨ur Einschluss, falls ~x mindestens in einem Einschluss ¨ liegt (durch die erlaubte Uberlappung von Einschl¨ussen kann ~x auch in mehreren Einschl¨ussen liegen), ansonsten das Material der Matrix. Da f¨ur jeden Punkt ~x der Integrationspunkt ~ω w¨ahrend der Simulation derselbe ist, und somit auch das Material dasselbe ist, wird ~ω lediglich zur Initialisierung der Funktion M u¨ bergeben. Somit kann man auch M(~x, ~ω ) = Mω~ (~x) schreiben. Da es keine Rolle spielt in welcher Reihenfolge die einzelnen Einschl¨usse generiert werden, liefern Integrationspunkte deren Bl¨ocke von m Komponenten lediglich vertauscht sind das gleiche Material, siehe auch Abbildung 28 in Abschnitt 4.3.1. Diese Tatsache resultiert in einer symmetrischen objektiven Funktion mit Symmetriegruppe Sm (siehe Abschnitt 2.3.4). F¨ur d = m ∗ n bedeutet dies genauer, dass zu einem gegebenen Integrationspunkt anstelle von insgesamt n! Integrationspunkten, die durch Permutationen der Symmetriebl¨ocke zustande kommen, durch Symmetrieausnutzung lediglich ein Integrationspunkt ausgewertet werden muss, der dann zus¨atzlich mit n! gewichtet wird. Dies ist nat¨urlich nur m¨oglich sofern das verwendete Quadraturverfahren eine solche Symmetrieausnutzung erlaubt. ¨ Eine M¨oglichkeit Symmetrie zu erhalten und dennoch Uberlappungen zu vermeiden, k¨onnte darin bestehen, nach der Generierung die Einschl¨usse nachzube34
Semantik
=
~ω ∈ Ω
ω1 .. . ωm . . . ω(n−1)∗m+1 .. . ωn∗m
Skalierung 1 .. . Skalierung m . . . Skalierung 1 .. . Skalierung m
pe1 ,1 .. . pe1 ,m . . . pen ,1 .. . pen ,m
···
Abbildung 20: Inverse Verteilungsfunktion zur Generierung von Einschl¨ussen.
handeln, indem u¨ berlappende Einschl¨usse so lange verkleinert werden bis diese ¨ sich nicht mehr u¨ berlappen. Auf weitere M¨oglichkeiten um das Uberlappungsproblem in den Griff zu bekommen (sofern u¨ berhaupt erw¨unscht) soll nicht einge¨ gangen werden. Wie schon angef¨uhrt soll in dieser Diplomarbeit die Uberlappung zugelassen werden. Alternativ zur Materialgenerierung k¨onnte auch das Material vorgegeben und die Geometrie des FE Modells durch einen Integrationspunkt erzeugt werden. Im Fokus dieser Diplomarbeit steht allerdings nur die Materialgenerierung. Im Folgenden sollen nun Verteilungen vorgestellt werden, deren Inverse bei dem soeben beschriebenen Modell Verwendung finden k¨onnten. 4.1.1 Skalierungsfunktionen Zur m¨oglichst genauen Definition einer Skalierungsfunktion u¨ ber eine inverse Verteilungsfunktion ist nicht nur die Wahl dieser entscheidend, sondern u¨ blicherweise auch eine Skalierung in ein entsprechendes Intervall erforderlich. Zu einer gegebenen Verteilungsfunktion F ist F˜ die durch x ∈ [xmin , xmax ] bedingte Verteilungsfunktion, das heißt, der Tr¨ager der bedingten Wahrscheinlichkeitsdichtefunktion liegt in [xmin , xmax ]. Demnach sind also F˜ (xmin ) = 0 und F˜ (xmax ) = 1. Somit liefert die zugeh¨orige inverse Verteilungsfunktion von F˜ zu ω ∈ [0, 1] ein x ∈ [xmin , xmax ]. Zu der Wahrscheinlichkeitsdichtefunktion f 35
Semantik
~ω ∈ Ω
=
Modellbildung/ -vervollst¨andigung
ω1 .. . ωm . . . ω(n−1)∗m+1 .. . ωn∗m
Skalierung 1 .. . Skalierung m . . . Skalierung 1 .. . Skalierung m
Material
FE Modell
Abbildung 21: Semantik-Funktion zur Generierung von Material mit Einschluss.
ergibt sich die bedingte Wahrscheinlichkeitsdichtefunktion f˜ =
f . F (xmax ) − F (xmin )
So ist die bedingte Verteilungsfunktion Z x Z x 1 ˜ ˜ f (t)dt F = f (t)dt = F (xmax ) − F (xmin ) xmin xmin 1 (F (x) − F (xmin )) mit x ∈ [xmin , xmax ]. F (xmax ) − F (xmin ) Der Erwartungswert f¨ur die bedingte Verteilung ergibt sich folglich zu Z xmax Z xmax 1 ˜ t f (t)dt. t f(t)dt = EF˜ = F (xmax ) − F (xmin ) xmin xmin =
Nun werden zwei Verteilungen n¨aher betrachtet, zum Einen die Rechteckverteilung, zum Anderen die im Folgenden beschriebene exponentielle Verteilung. In Abschnitt 6.2 werden die bedingten Formen dieser Verteilungen genutzt. Die Wahrscheinlichkeitsdichtefunktion der Gleichverteilung [18], auch Rechteckverteilung genannt, f¨ur ein Intervall [xmin , xmax ] ist gegeben durch ( 1 f¨ur xmin ≤ x ≤ xmax fR (x) = xmax −xmin 0 sonst. 36
Die Verteilungsfunktion lautet 0 x−xmin FR (x) = xmax −xmin 1
f¨ur x < xmin f¨ur xmin ≤ x ≤ xmax f¨ur x > xmax ,
und der Erwartungswert
xmin + xmax . 2 Der Tr¨ager der Wahrscheinlichkeitsdichtefunktion liegt in [xmin , xmax ]. F¨ur die Wahrscheinlichkeitsdichtefunktion der betrachteten exponentiellen Verteilung gilt ( λ2 xe−λx f¨ur x ≥ 0, λ > 0 fE (x) = 0 sonst ER =
Die Verteilungsfunktion ist gegeben durch ( 1 − (λx + 1)e−λx FE (x) = 0
f¨ur x ≥ 0, λ > 0 sonst
und der Erwartungswert nach zweifacher partieller Integration durch 2 ∞ Z ∞ Z ∞ x 2x 2 2 2 −λx 2 −λx EE = x fE (x)dx = − − 2− 3 λ x e dx = λ e λ λ λ −∞ 0 0 2 . λ Daraus ergeben sich f¨ur die bedingte Wahrscheinlichkeitsdichtefunktion ( 1 λ2 xe−λx f¨ur 0 ≤ xmin ≤ x ≤ xmax , λ > 0 f˜E (x) = F (xmax )−F (xmin ) 0 sonst , =
die bedingte Verteilungsfunktion 1 ((λxmin + 1)e−λxmin − (λx + 1)e−λx ) F (xmax )−F (xmin ) ˜ FE (x) = 0 37
f¨ur 0 ≤ xmin ≤ x, x ≤ xmax , λ>0 sonst
und den Erwartungswert Z xmax x f˜E (x)dx = E˜ = E
xmin
1 = F (xmax ) − F (xmin )
1 F (xmax ) − F (xmin )
Z
xmax
x fE (x)dx xmin
2 xmax ! x 2x 2 λ2 e−λx − − 2 − 3 λ λ λ xmin
Die inverse Verteilungsfunktion der Gleichverteilung ist analytisch, die der exponentiellen Verteilung durch Diskretisierung implementiert, siehe Abschnitt 5.1.
4.2 Simulation und Bewertung Nach Generierung der r¨aumlichen Materialverteilung aus einem Integrationspunkt ~ω durch die Semantik-Funktion kann das FE Modell simuliert werden. An diskreten Punkten innerhalb der Finiten-Elemente des FE Modells fragt das FE Programm w¨ahrend der Simulation das Material u¨ ber die Materialfunktion M aus Gleichung 8 ab. Das Simulationsergebnis wird durch die Bewertung auf die interessierende Kenngr¨oße heruntergebrochen, die dann als Funktionswert der objektiven Funktion dem Quadraturverfahren zur¨uckgegeben wird. Interessierende Kenngr¨oßen k¨onnen beispielsweise maximale Spannung oder maximale Verschiebung sein. Wie schon zu Beginn von Abschnitt 3 erw¨ahnt, k¨onnen auch m Kenngr¨oßen als m-wertiger Funktionswert dem Quadraturverfahren u¨ bergeben werden. Hier soll allerdings der Einfachheit halber eine einwertige objektive Funktion betrachtet werden.
4.3 Ellipsen als Einschlusse ¨ Nun wird das soeben beschriebene Modell zur Veranschaulichung weiter konkretisiert. Es soll ein zweidimensionales FE Modell mit Ellipsen als Einschl¨usse im Material betrachtet werden. Im Folgenden werden wieder die einzelnen Spezifikationen f¨ur die in Abbildung 18 aufgef¨uhrten Komponenten angegeben. 4.3.1 Semantik Aus einem Integrationspunkt ~ω ∈ Ω wird die r¨aumliche Materialverteilung eines zweidimensionalen FE Modells generiert. Dabei sind die Einschl¨usse im Material u¨ ber Ellipsen definiert. Somit ist das Material durch die vorgegebenen Materialparameter der Matrix und der Einschl¨usse sowie die Lage und Orientierung der Ellipsen gegeben. 38
hy
×
hx
φ
~xM
Abbildung 22: Ellipse.
Eine um den Ursprung rotierte und anschließend verschobene Ellipse kann durch f¨unf Parameter eindeutig beschrieben werden, den beiden Koordinaten xM und yM des Mittelpunktes ~xM , den beiden Halbachsen hx und hy (entlang der Koordinatenachsen vor Ellipsenrotation) und dem Drehwinkel φ, siehe Abbildung 22. F¨ur n Ellipsen wird zur unabh¨angigen Generierung dieser demnach ein Inte¨ grationspunkt ~ω der Dimension d = n ∗ 5 ben¨otigt. Uber inverse Verteilungsfunktionen werden die einzelnen Komponenten von ~ω , wie in Abschnitt 4.1 beschrieben (mit m = 5), entsprechend der oben genannten Parameterreihenfolge skaliert, siehe Abbildung 23. Die inversen Verteilungsfunktionen m¨ussen so gew¨ahlt werden, dass die Parameter der Ellipsen in f¨ur das FE Modell sinnvolle Intervalle abgebildet werden. F¨ur den Drehwinkel bietet sich das Intervall [0, π2 ] an, da so ausgeschlossen werden kann, dass zwei unterschiedliche Bl¨ocke von f¨unf Komponenten des Integrationspunktes ~ω auf gleiche Ellipsen abgebildet werden k¨onnen, siehe Abbildung 24. Damit jede denkbare Orientierung einer Ellipse f¨ur gegebene Halbachsen erreicht werden kann, m¨ussen allerdings die beiden Intervalle f¨ur die Halbachsen gleich sein, denn eine zu einer gegebenen Ellipse um π2 rotierte Ellipse erh¨alt man durch Vertauschen der Halbachsen, siehe Abbildung 25. Es gilt nun noch die Materialfunktion M aus Gleichung 8 aufzustellen. Dazu wird f¨ur die Parameter einer jeden Ellipse die Kegelschnittgleichnung aufgestellt, mittels welcher festgestellt werden kann, ob ein gegebener Punkt ~xp des FE Mo¨ dells in mindestens einer der generierten Ellipse liegt (durch die erlaubte Uberlappung kann ~xp auch in mehreren Ellipsen liegen). Ist dies so, liefert M das Material der Einschl¨usse, ansonsten der Matrix. Zur Ableitung der Kegelschnittgleichung einer um den Ursprung rotierten und anschließend verschobenen Ellipse (siehe Abbildung 26), muss ausgehend von der Kegelschnittgleichung f¨ur die zugeh¨orige unrotierte und nicht verschobene Ellipse die Transformationsreihenfolge umgekehrt und die inversen Transformationen verwendet werden. Die Kegelschnittgleichung einer unrotierten Ellipse mit dem 39
Semantik
=
~ω ∈ Ω
ω1 .. . ω5 . . . ω(n−1)∗5+1 .. . ωn∗5
Skalierung 1 .. . Skalierung m . . . Skalierung 1 .. . Skalierung m
xM1 .. . φ1 . . . xMn .. . φn
···
Abbildung 23: Inverse Verteilungsfunktion zur Generierung von Ellipsen als Einschl¨ussen.
gleiche
×
φ+π
φ
×
Ellipsen
Abbildung 24: gleiche Ellipsen trotz unterschiedlicher Drehwinkel.
40
φ ∈ [0◦ , 90◦ )
φ ∈ [0◦ , 90◦) vertauschte
×
φ
×
φ
Halbachsen
Abbildung 25: eingeschr¨ankter Drehwinkel von [0, π2 ] f¨ur Ellipsen; drei Ellipsen mit vertauschten Halbachsen.
Mittelpunkt im Ursprung lautet x2 y2 + =1 h2x h2y mit den Achsenabschnitten (Halbachsen) hx und hy . Alternativ kann diese auch folgend geschrieben werden ! 1 0 x h2x . und D := ~xT D~x = 1 mit ~x = 0 h12 y y
F¨ur die Generierung der Kegelschnittgleichung einer um den Ursprung rotierten und anschließend verschobenen Ellipse mit den Parametern ~xM = (xM , yM ), hx , hy und φ muss zun¨achst die zugeh¨orige nicht transformierte Ursprungsellipse invers verschoben werden. Daraus ergibt sich (~x − ~xM )T D(~x − ~xM ) = 1 Zum Schluss muss die invers verschobene Ursprungsellipse invers rotiert werden, demnach gilt (~x − ~xM )T Rot(−φ)T DRot(−φ)(~x − ~xM ) = 1 cos(α) − sin(α) . mit der Rotationsmatrix Rot(α) := sin(α) cos(α) 41
~pm
×
×
φ
Abbildung 26: Generierung einer Ellipse, Drehung mit anschließender Verschiebung.
Durch Ausschreiben dieser Gleichung erh¨alt man die gew¨unschte Kegelschnittgleichung der Form ax2 + by 2 + cxy + dx + ey + f = 1
mit a, b, c, d, e, f ∈ R.
Somit ist ein Punkt ~xp = (xp , yp )T in der rotierten und verschobenen Ellipse, falls gilt ax2p + byp2 + cxp yp + dxp + eyp + f ≤ 1
Abbildung 27 zeigt die Semantik-Funktion unter Einbeziehung der Koeffizientenbildung der Kegelschnittgleichung. Wie in Abschnitt 4.1 angef¨uhrt, besteht durch die unabh¨angige Generierung der Ellipsen und in diesem Zusammenhang die Zulassung von Ellipsen¨uberlappungen die Symmetriegruppe S5 f¨ur die objektive Funktion. F¨ur drei generierte Ellipsen zeigt Abbildung 28 alle Integrationspunkte, die zu ein und demselben Material f¨uhren. Jeder dieser Integrationspunkte ist u¨ ber die Symmetriebl¨ocke ~ωei mit i ∈ {1, 2, 3} dargestellt (siehe auch Abschnitt 2.3.4). Dabei wird aus Symmetrieblock ~ωei die Ellipse ei generiert. Von den sechs Integrationspunkten muss bei Symmetrieausnutzung lediglich einer ausgewertet, aber f¨ur sechs Materialien gewichtet werden. 4.3.2 Simulation und Bewertung Das durch die Semantik-Funktion generierte 2D FE Modell wird durch das FE Programm simuliert. W¨ahrend dieser Simulation wird das Material an diskreten Punkten innerhalb der Finiten-Elemente des 2D FE Modells abgefragt, siehe Abbildung 29. 42
Semantik
ω1 .. . ω5 . . . ω(n−1)∗5+1 .. . ωn∗5
Skalierung
xM1 .. . φ1 . . . xMn .. . φn
Ableitung der Kegelschnittgleichung
=
~ω ∈ Ω
a1 . . . f1 an . . . fn
···
Abbildung 27: Semantik-Funktion unter Einbeziehung der Koeffizientenbildung der Kegelschnittgleichung.
~ωe1 ~ωe1 ~ωe2 , ~ωe3 , ~ωe2 ~ωe3
e1 e2
~ωe2 ~ωe2 ~ωe1 , ~ωe3 , ~ωe1 ~ωe3
~ωe3 ~ωe3 ~ωe1 , ~ωe2 ~ωe1 ~ωe2
e3
Abbildung 28: Material mit drei Ellipsen, Symmetriegruppe S5 .
43
diskretisiertes 2D FE Modell
bb
b
bb
b
Material
Material? Material? Einschluss! Matrix!
bb
b
bb
b
Abbildung 29: Simulation eines 2D FE Beispielmodells durch FE Programm.
Die Bewertung soll hier aus dem Ergebnis der Simulation die maximale Verschiebung extrahieren und als Funktionswert der objektiven Funktion an das Quadraturverfahren des Teilsystems Quadraturalgorithmen zur¨uckgeben.
5
Design und Implementierung
Das in Abschnitt 3 vorgestellte System zur statistischen Auswertung numerischer Simulationen ist in der objektorientierten Programmiersprache C++ unter dem Betriebssystem Linux implementiert. Sofern sinnvoll, wurden die Klassen zur gr¨oßtm¨oglichen Flexibilit¨at als Templates umgesetzt. Da eine ganze Reihe von Klassen auch zur L¨osung anderer Probleme herangezogen werden k¨onnen, entstanden viele mehr oder weniger kleine Projekte, deren ¨ Uberf¨ uhrung in CTL-Komponenten direkt m¨oglich ist, siehe auch Abschnitt 5.3.1. So ergibt sich statt einer tiefen Verzeichnisstruktur ein breiter Verzeichnisbaum, siehe Abbildung 30. Da die einzelnen Projekte gr¨oßtenteils unabh¨angig voneinander implementiert wurden und somit zumeist einen allgmeing¨ultigeren Standpunkt vertreten, sollen diese vorerst nicht in das System aus Abschnitt 3 gepresst werden. In Abschnitt 5.1 werden grundlegende Klassen und Projekte vorgestellt. Im Anschluss folgt 44
Abbildung 30: Verzeichnisbaum.
45
in Abschnitt 5.2 der Zusammenhang zu dem System aus Abschnitt 3 anhand der implementierten objektiven Funktionen und deren statistischen Auswertung aus Abschnitt 6.
5.1 grundlegende Klassen und Projekte Im Folgenden werden Klassen bzw. Projekte dieser Diplomarbeit n¨aher beschrieben. Generell werden nicht alle Methoden einer Klasse betrachtet, sondern zu ¨ Gunsten der Ubersichtlichkeit vielmehr einige wenige extrahiert. Die Klasse Function im Verzeichnis functions ist eine Basisklasse f¨ur eine mathematische Funktion. Von dieser sind eindimensionale Polynome, Statistikfunktionen als auch die in Abschnitt 6 vorgestellten konkreten objektiven Funktionen abgeleitet. Die Basisklasse Distribution des gleichnamigen Projekts im Verzeichnis distribs spezifiziert die u¨ bergreifende Beschreibung einer Verteilung. Sie beinhaltet Methoden f¨ur die Wahrscheinlichkeitsdichtefunktion, die Verteilungsfunktion und die inverse Verteilungsfunktion. Die zur Wahrscheinlichkeitsdichtefunktion geh¨orende Methode ist neben dem Destruktor die einzige virtuelle Methode dieser Klasse. Zur Spezifizierung einer konkreten Verteilung muss eine abgeleitete Klasse lediglich diese Methode implementieren. Die Verteilungsfunktion als auch die inverse Verteilungsfunktion werden durch Diskretisierung und st¨uckweiser Integration mittels der in der Basisklasse zugeh¨origen implementierten Methoden gewonnen. Nach den Regeln der Vererbung ist es nat¨urlich auch m¨oglich, die Methoden der Verteilungsfunktion und inversen Verteilungsfunktion in der abgeleiteten Klasse zu u¨ berschreiben. Ganz a¨ hnlich verh¨alt es sich mit der Basisklasse Point gen des Projekts Point Generator im Verzeichnis pt gens zur Beschreibung eines allgemeing¨ultigen Punktgenerators, siehe Abbildung 31. In einer abgeleiteten Klasse muss zur Konkretisierung eines Punktgenerators nur die virtuelle Methode get point der Basisklasse implementiert werden. Der Konstruktor erh¨alt als Parameter die Dimension dim als auch die untere und obere Schranke des Intervalls [min, max], aus dem die einzelnen Komponenten der dim-dimensionalen Punkte generiert werden sollen. Hierzu stellt die Basisklasse eine Methode zur Punktskalierung bereit. Abgeleitet wurden ein auf der Methode ranf der Bibliothek RANDLIB [2] basierender Zufallsgenerator und Generatoren f¨ur die Halton- und HammersleyPunktfolge sowie die arithmetische Folge aus Abschnitt 2.2.3. Die Klasse Monte carlo des Projekts (Quasi)-Monte Carlo im Verzeichnis mc implementiert das Rahmenwerk f¨ur die Monte Carlo und Quasi-Monte Carlo Methode, siehe Abbildung 31. Spezifiziert wird eine dieser Methoden durch den im Konstruktor u¨ bergebenen Pointer pt gen auf einen von der Klasse Point gen abgeleiteten Punktgenerator. Durch diese Vorgehensweise reduziert 46
Monte_carlo Monte_carlo( f: const F*, fs_stat: const vector&, ints: vector&, pt_gen: Point_gen*) integrate( n: unsigned int) : void
Point_gen Point_gen( dim: unsigned int, min: real, max: real) get_point(pt: vector&) : void get_dim() : unsigned int get_min() : real get_max() : real
Abbildung 31: UML Klassendiagramm, Auszug aus den Klassen Monte carlo und Point gen.
sich der Aufwand des Generatoraustauschs lediglich auf den Konstruktoraufruf. Zudem muss dem Konstruktor die auszuwertende dim-dimensionale Funktion f und die mit dieser separat zu verkn¨upfenden eindimensionalen Statistikfunktionen fs stat u¨ bergeben werden. Die Funktionen sind als Template-Parameter F und F stat in der Klasse Monte carlo spezifiziert, um so dem Anwender keine Ableitung von der Klasse Function aufzuzwingen. Die approximierten Integralwerte werden in ints abgelegt. Die Dimension dim als auch der Raum [min, max]dim u¨ ber den integriert werden soll, sind durch pt gen gegeben. Demnach muss zur Approximation der Integrale lediglich der Methode integrate die Anzahl n der auszuwertenden Integrationspunkte mitgeteilt werden. integrate kann wiederholt aufgerufen werden, um weitere m Integrationspunkte auszuwerten. Allerdings ist diese Funktionsweise nur dann gew¨ahrleistet, wenn der genutzte Punktgenerator dies unterst¨utzt. Von den implementierten Punktgeneratoren ist lediglich der Hammersley-Punktgenerator dazu nicht in der Lage, da zur Generierung der Punkte die Anzahl der zu generierenden Punkte einbezogen wird, siehe Abschnitt 2.2.2. Die Klasse One dim quadrature des Projekts One dim quadrature im Verzeichnis one dim quad stellt die Basisklasse f¨ur eine eindimensionale Quadraturformel dar, siehe Abbildung 32. Eine abgeleitete Klasse zur Konkretisierung einer eindimensionalen Quadraturformel muss nur die abgeleitete Methode get spts wgts not scaled zur Generierung der unskalierten Punkte und Gewichte implementieren. Als unskaliert gelten Punkte und Gewichte, wenn sie f¨ur das urspr¨ungliche Intervall, typischerweise [−1, 1], der Quadraturformel generiert werden. Dem Konstruktor werden neben der Ord47
nung level die beiden Schranken des Intervalls [min, max] u¨ bergeben. Die Skalierung der Punkte und Gewichte in dieses Intervall erfolgt automatisch u¨ ber die in der Basisklasse implementierten Methoden get spts wgts und get spts wgts till level. Wegen weiteren implementierten Methoden in der Basisklasse zur komfortablen und flexiblen Nutzung einer eindimensionalen Quadraturformel muss get spts wgts not scaled auf der ihr u¨ bergebenen Ordnung level operieren. Als eindimensionale Quadraturformeln sind die Gauss-Legendre und die Clenshaw-Curtis Quadraturformel implementiert. Desweiteren werden die Kronrod-Patterson Quadraturformel bis zur Ordnung 6 und die verz¨ogerte Kronrod-Patterson Quadraturformel bis zur Ordnung 48 angeboten. Die Punkte und Gewichte der beiden letzten Quadraturformeln wurden der Software smolpack entnommen, siehe Abschnitt 5.3.3. Smolyak Smolyak( dim: unsigned int, sym: unsigned int, one_dim_quad: const One_dim_quadrature*) integrate( f: const F*, fs_stat: const vector&, ints: vector&) : void
Multi_index_gen Tensor_product
One_dim_quadrature One_dim_quadrature( level: unsigned int, min: real, max: real) get_spts_wgts( spts: vector&, wgts: vector&) : void get_spts_wgts_till_level( spts: vector&, wgts: vector&) : void integrate( f: const Function*) : real get_level() : unsigned int get_min() : real get_max() : real get_spts_wgts_not_scaled( level: unsigned int, spts: vector&, wgts: vector&) : void
Abbildung 32: UML Klassendiagramm, Auszug aus den Klassen Smolyak und One dim quadrature.
Die Klasse Smolyak des gleichnamigen Projekts im Verzeichnis smolyak umfasst den Smolyak Algorithmus als auch den modifizierten Smolyak Algorithmus zur Symmetrieausnutzung, siehe Abbildung 32. Zudem wird eine Implementierung f¨ur anisotrope Dimensionen angeboten. Bei anisotropen Dimensionen k¨onnen verschiedene Ordnungsfolgen der verwendeten eindimensionalen Quadraturformel f¨ur unterschiedliche Dimensionen genutzt werden, um so auch in h¨oher48
en Dimensionen von beispielsweise 360 mitunter akzeptable Ergebnisse erhalten zu k¨onnen, siehe [14]. Durch verschiedene Konstruktoren kann zwischen isotroper und anisotroper Anwendung gew¨ahlt werden. Alle Konstruktoren erm¨oglichen Symmetrieausnutzung. Die Spezifizierung der Symmetriegruppe erfolgt u¨ ber die Anzahl sym der in einem Symmetrieblock (siehe Abschnitt 2.3.4) zu findenen Komponenten. Sollen keine Symmetrien ausgenutzt werden ist sym gleich der Dimension dim. Desweiteren wird dem in Abbildung 32 gezeigten Konstruktor neben dim und sym lediglich ein Pointer one dim quad auf eine von der Basisklasse One dim quadrature abgeleitete eindimensionale Quadraturformel u¨ bergeben. Die Smolyak Ordnung level und die Schranken des Raumes [min, max]dim , u¨ ber den integriert werden soll, sind durch one dim quad definiert. Es stehen eine ganze Reihe weiterer Konstruktoren zur komfortablen und flexiblen Nutzung zur Verf¨ugung. Beispielsweise kann mitunter auch eine Menge von eindimensionalen Integrationspunkten und zugeh¨origen Gewichten anstelle eines Pointers one dim quad u¨ bergeben werden. So oder so, die eindimensionalen Quadraturformeln werden aus allgemeing¨ultiger Sicht betrachtet, so dass beliebige genutzt werden k¨onnen. Alle Integrationspunkte und zugeh¨origen Gewichte des Smolyak Algorithmus werden durch den Konstruktoraufruf generiert und bis zur Zerst¨orung der Instanz im Speicher gehalten. Dabei werden identische Integrationspunkte lediglich einmal gespeichert und die zugeh¨origen Gewichte aufaddiert. Die Speicherung hat Vor- und Nachteile. Der Nachteil des mitunter ben¨otigten hohen Speicheraufwands liegt klar auf der Hand. Allerdings k¨onnen die Vorteile einer geschachtelten eindimensionalen Quadraturformel genutzt werden, ohne genaue Kenntis u¨ ber die Schachtelung zu haben. So werden die ben¨otigten Funktionsauswertungen auf ein Minimum beschr¨ankt. Dies gilt auch f¨ur nicht-geschachtelte Quadraturformeln, die gleiche Punkte in verschiedenen Ordnungen aufweisen. Dadurch ergibt sich ein weiterer Nachteil. Die Generierung der Integrationspunkte und Gewichte des Smolyak Algorithmus ist nicht effizient. Allerdings ist der Programmcode und die Nutzung u¨ bersichtlicher, da zwischen konkreten eindimensionalen Quadraturformeln nicht unterschieden werden muss. Außerdem k¨onnte die Implementierung um die in [15] vorgestellte Baumstruktur zur effizienten Berechnung der Gewichte erweitert werden. Zur Symmetrieausnutzung ist die Speicherung der Integrationspunkte und Gewichte ohnehin n¨otig. Durch Ausnutzung von Symmetrien wird eine weitere wesentliche Reduzierung der sp¨ateren Auswertungen vorgenommen. In Anbetracht eines zeitlich aufwendigen Integranden kann so der mitunter hohe Speicheraufwand gerechtfertigt werden. Zur Ausnutzung der Symmetrien werden die Integrationspunkte zweifach gefiltert. Aus Effizienzgr¨unden muss eine Filterung so fr¨uh wie m¨oglich einsetzen, und zwar in der Multi-Index Generierung. Hierzu werden nur die Multi-Indizes generiert, dessen Symmetriebl¨ocke lediglich lexikographisch geordnet sind und die die vorgegebenen Einschr¨ankungen des Smolyak Algo49
rithmus erf¨ullen. Eine weitere Filterung mittels lexikographischer Ordnung muss auf den durch die Multi-Indizes resultierenden Integrationspunkten aufsetzen. Die Multi-Index Filterung wird durch eine Instanz der Klasse Multi index gen durchgef¨uhrt, die Filterung der Integrationspunkte durch eine Instanz der Klasse Tensor product im Verzeichnis smolyak, siehe Abbildung 32. Das Projekt Two dim display im Verzeichnis view ist eine mittels der C++-Bibliothek Qt der Firma Troll Tech implementierte grafische Oberfl¨ache zur Darstellung von Ellipsen. Abbildung 44 in Abschnitt 6.2.4 zeigt ein Anwendungsbeispiel. Zum komfortablen Plotten von diskretisierten mathematischen Funktionen wurde die Klasse Gnuplot support im Verzeichnis util implementiert. Diese stellt f¨ur zwei- und dreidimensionale Plots eine ganze Reihe von Methoden bereit und generiert Daten- und Skriptdateien f¨ur das kommandozeilenbasierte interaktive Plotprogramm GNUPLOT. Eine Skriptdatei kann in GNUPLOT mittels load Anweisung geladen werden. Es entstanden weitere Klassen bzw. Projekte auf die nicht weiter eingegangen werden soll.
5.2 System zur statistischen Auswertung numerischer Simulationen Das System zur statistischen Auswertung numerischer Simulationen aus Abschnitt 3 nutzt die implementierten Quadraturprojekte aus Abschnitt 5.1 zur Integralapproximation u¨ ber dem Raum Ω. Das Projekt Distribution aus Abschnitt 5.1 kann zur Skalierung der Integrationspunkte in der Semantik-Funktion der objektiven Funktion genutzt werden. Im Folgenden wird die Implementierung der statistischen Auswertung aus Abschnitt 6.2.4 beschrieben. In Anlehnung an Abbildung 18 aus Abschnitt 3.2 erzeugt der Quadraturalgorithmus Integrationspunkte aus Ω. Aus einem Integrationspunkt ~ω generiert die Semantik-Funktion der objektiven Funktion f die r¨aumliche Materialverteilung Mω~ (~x) des vorgegebenen 2D FE Modells. Anschließend wird das Modell durch das FE Programm ParaFEP (siehe Abschnitt 5.3.2) simuliert. Aus dem Simulationsergebnis wird die maximale Verschiebung durch die Bewertung extrahiert und als Funktionswert f (~ω ) an den Quadraturalgorithmus zur¨uckgegeben. Zur Simulation der FE Modelle laufen siebzehn ParaFEP Instanzen auf einem Cluster von achtzehn Prozessoren. Die ParaFEP Instanzen laufen parallel zueinander, befinden sich allerdings jeweils im seriell-optimierten Betrieb. Das heißt, eine ParaFEP Instanz l¨ost ein aus einem Integrationspunkt ~ω resultierendes FE Modell vollst¨andig alleine. Jede ParaFEP Instanz greift auf eine Material-Instanz zu. Die Instanzen sind als CTL-Komponenten umgesetzt, siehe
50
objektive Funktion f Rechenstation 2 Rechenstation 1 (Quadratur-) Algorithmus
Material
Semantik
~ω Parafepcall
~ω
f (~ω ) f (~ω) Bewertung
~ω
~x
Mω~ (~x)
ParaFEP Simulation
Abbildung 33: statistische Auswertung einer konkreten objektiven Funktion.
Abschnitt 5.3.1. Daher werden diese im Folgenden als ParaFEP- bzw. MaterialKomponenten bezeichnet. Die ParaFEP-Komponenten werden u¨ ber Prozessorgrenzen hinweg angesprochen. Hingegen erfolgt die Kommunikation zwischen einer ParaFEP-Komponente und der zugeh¨origen Material-Komponente aus Effizienzgr¨unden lokal. Die Hauptroutine l¨auft auf dem ersten Prozessor. Auf diesem findet keine FE Simulation statt. Abbildung 33 veranschaulicht die Auswertung anhand zweier Prozessoren. Auf dem ersten Prozessor l¨auft die Hauptroutine. Der Quadraturalgorithmus generiert einen Integrationspunkt ~ω und gibt diesen an das Frontend Parafepcall der objektiven Funktion f . In der Implementierung sind Parafepcall und f identisch. Parafepcall leitet ~ω an die ParaFEPKomponente auf dem zweiten Prozessor weiter. Die ParaFEP-Komponente u¨ bergibt ~ω an die Material-Komponente, welche die r¨aumliche Materialverteilung Mω~ (~x) des vorgegebenen Modells generiert. W¨ahrend der Simulation fragt die ParaFEP-Komponente das Material an den Integrationspunkten ~x ∈ R2 der FiniteElemente-Methode u¨ ber die Material-Komponente ab. Im Anschluss an die Simulation wird die maximale Verschiebung aus dem Ergebnis extrahiert und als Funktionswert f (~ω ) an den Algorithmus zur¨uckgeliefert. F¨ur die achtzehn Prozessoren u¨ bernimmt Parafepcall die Funktion des Schedulers und teilt die vom Quadraturalgorithmus erzeugten Integrationspunkte den siebzehn ParaFEPKomponenten zu.
51
5.3 Werkzeuge F¨ur das in dieser Diplomarbeit zu entwickelnde System aus Abschnitt 3 und zur Analyse hochdimensionaler Quadraturverfahren wurden eine Reihe Werkzeuge verwendet. Ausgew¨ahlte sind in den Abschnitten 5.3.1 bis 5.3.3 beschrieben. Zudem wurden, wie in Abschnitt 5.1 erw¨ahnt, die C++-Bibliothek Qt und das Plotprogramm GNUPLOT verwendet. Die UML-Diagramme aus Abschnitt 5.1 wurden mittels ArgoUML erstellt. 5.3.1 Component Template Library (CTL) Die folgenden Passagen wurden dem CTL Manual for Linux/Unix entnommen und in das Deutsche u¨ bersetzt. Die Component Template Library (CTL) ist eine Implementierung der auf der generischen C++ Template Programmierung basierenden Komponententechno¨ logie. Ahnlich wie CORBA kann sie zur Realisierung verteilter Komponentenbasierter Software Systeme benutzt werden, wobei eine Komponente ein St¨uck Software ist, das aus einer wohldefinierten Schnittstelle und einer Implementierung besteht. Schnittstelle und Implementierung sind durch einen Kommunikationskanal u¨ ber beispielsweise TCP/IP oder MPI miteinander verbunden. Auf diese Weise ist die Benutzung einer Komponente unabh¨angig vom Ort, also Ortstransparent. 5.3.2 Parallel Finite Element Program (ParaFEP) In ParaFEP wurden verschiedene Probleml¨osungen in den Bereichen Informatik, Numerik und Mechanik zu einem h-, p- und d-adaptiven dynamisch lastbalancierten durchgehend parallelen FEM-Code zusammengef¨uhrt. In ParaFEP sind eine parallele knotenregul¨are, lokale Hexaeder-Verfeinerung - objektorientierte Darstellung der eigentlichen Finiten-Elemente mit Parametrisierung in p-Ansatz, Integrationsordnung und Material - Kopplung von Schalen und Kontinuum - automatische Abarbeitung allgemeiner, asynchroner, paralleler Algorithmen - Repartitionierung verteilter Graphen - effektive parallele Matrixoperationen und Vorkonditionierer - Parallelisierung der Visualisierung - fomuliert und in Form von Algorithmen und C++ Klassen implementiert, siehe [11]. 5.3.3 smolpack und finpack Die C-Programme smolpack und finpack [13] implementieren den Smolyak Algorithmus zur Approximation von Integralen u¨ ber dem Raum [0, 1]d . smolpack unterst¨utzt isotrope Dimensionen, also dieselbe Ordnungsfolge der genutzten eindimensionalen Quadraturformel in jeder Dimension, finpack hingegen auch ani52
sotrope Dimensionen. Beide Programme bieten die verz¨ogerte Kronrod-Patterson Quadraturformel an, smolpack zus¨atzlich die Clenshaw-Curtis Quadraturformel. Die eindimensionalen Quadraturformeln sind bis zu einer bestimmten Ordnung als Menge an Punkten und zugeh¨origen Gewichten im Programmcode abgelegt. Die Programme implementieren Algorithmen zur effizienten Berechnung von Integrationspunkten und insbesondere von Gewichten. smolpack und finpack konnten zur Korrekheitspr¨ufung des in dieser Diplomarbeit implementierten Smolyak Algorithmus wesentlich beitragen. Zudem wurden die Integrationspunkte und Gewichte der Kronrod-Patterson Quadraturformel f¨ur das Interval [0, 1] aus smolpack entnommen und zur Implementierung der Kronrod-Patterson und verz¨ogerten Kronrod-Patterson Quadraturformel verwendet.
6
Beispiele
S¨amtliche Beispiele k¨onnen von dem in Abschnitt 4.3 beschriebenen Modell abgeleitet werden. Zum Einen wird ein m¨ogliches Beispiel der realen Welt vorgestellt. Zum Anderen werden die Simulationen von Beispielmodellen mittels dem in dieser Diplomarbeit implementierten Systems statistisch ausgewertet und die Ergebnisse analysiert.
6.1 Motivationsbeispiel Als Material des aus der realen Welt zu u¨ berf¨uhrenden Modells soll Gasbeton (Porenbeton) betrachtet werden. Dabei wird weder auf die Geometrie des Modells noch auf konkrete Parameter eingegangen, sondern lediglich ein anschauliches m¨ogliches Anwendungsbeispiel und kurz die Vorgehensweise zur Bildung der objektiven Funktion vorgestellt. Die Abbildungen 34 und 35 zeigen Nahaufnahmen aufgebrochener Gasbetonsteine. Wie hier zu sehen, sind die Lufteinschl¨usse Ellipsen a¨ hnlich und k¨onnen sich u¨ berlappen. So kann das in Abschnitt 4.3 beschriebene Modell als Basis f¨ur dieses Anwendungsbeispiel verwendet werden. Zur Bildung der objektiven Funktion gilt es die drei Funktionen Semantik, Simulation und Bewertung wie gehabt zu definieren. Aufgrund der ausf¨uhrlichen Erl¨auterung der Spezifikationen in Abschnitt 4.3 soll lediglich die SemantikFunktion n¨aher beleuchtet werden. Durch Aufbrechen einer Vielzahl von Gasbetonsteinen kann das Verh¨altnis zwischen den Lufteinschlussausmaßen und der zugeh¨origen H¨aufigkeit des Auftretens abgesch¨atzt werden. Bei Betrachtung der Einschl¨usse als Ellipsen kann so m¨oglicherweise festgestellt werden, dass Ellipsen mit Halbachsen aus [1.5mm, 2.0mm] a¨ ußerst selten, aus [0.8mm, 1.0mm] recht selten aber dennoch 53
Abbildung 34: aufgebrochener Gasbetonstein Nr.1.
Abbildung 35: aufgebrochener Gasbetonstein Nr.2.
54
h¨aufiger und aus [0.2mm, 0.6mm] typisch sind. Die inverse Verteilungsfunktion zur Skalierung der f¨ur die Halbachsen zur Verf¨ugung stehenden Komponenten eines durch das Quadraturverfahren generierten Integrationspunktes ω ~ kann nun so gew¨ahlt werden, dass die resultierenden Halbachsen mit den Messungen und H¨aufigkeitsanalysen u¨ bereinkommen. Beispielsweise k¨onnte die Inverse der in Abschnitt 4.1.1 vorgestellten bedingten exponentiellen Verteilungsfunktion f¨ur die Halbachsen Verwendung finden. Nach Festlegung der inversen Verteilungsfunktionen aller f¨unf Parameter einer Ellipse, kann v¨ollig analog zu Abschnitt 4.3 vorgegangen werden.
6.2 konkrete Beispiele Die konkreten Beispiele dieser Diplomarbeit k¨onnen wie schon erw¨ahnt von dem in Abschnitt 4.3 definierten 2D FE Modell abgeleitet werden. Dabei werden zwei Modelle mit unterschiedlichen Geometrien betrachtet, dessen zugeh¨orige objektive Funktionen sich zudem in der Symmetriegruppe unterscheiden. Die Parameter beider Modelle gen¨ugen nicht den Anforderungen eines aus der realen Welt u¨ berf¨uhrten Modells. Zu Beginn sollen in den Abschnitten 6.2.1 und 6.2.2 weitestgehend gemeinsame Spezifikationen der Beispielmodelle betrachtet werden. Im Anschluss folgen in den Abschnitten 6.2.3 und 6.2.4 neben den unterschiedlichen Spezifikationen bzw. Einschr¨ankungen die Ergebnisse und Analysen. 6.2.1 Semantik Da die generelle Aufgabe der Semantik-Funktion in Abschnitt 4.3 ausf¨uhrlich erl¨autert ist, sollen hier lediglich die Parameter und Skalierungsfunktionen f¨ur die f¨unf Parameter einer Ellipsen konkretisiert werden. Die Materialparameter der Matrix und der Einschl¨usse (Ellipsen) beschr¨anken sich jeweils auf den Elastizit¨atsmodul (E-Modul) und die Querkontraktion. F¨ur N Einschluss betr¨agt der E-Modul 1.0e8 mm 2 , die Querkontraktion 0.3. Der E-Modul N der Matrix ist hingegen 1.0e4 mm und die Querkontraktion 0.4. 2 Die inverse Verteilungsfunktion der Gleichverteilung, siehe Abschnitt 4.1.1, wird zur Skalierung der Koordinaten xM und yM des Ellipsenmittelpunktes ~xM mit den minimalen und maximalen Koordinaten der Modellgeometrie und zur Skalierung des Drehwinkels im Intervall [0, π2 ] (siehe hierzu auch Abschnitt 4.3) verwendet. F¨ur die beiden Halbachsen wird hingegen die gleiche bedingte exponentielle Verteilung aus Abschnitt 4.1.1 im gleichen Intervall genutzt, siehe auch Abschnitt 4.3. Die Bestimmung der Intervallgrenzen und des Verteilungsparameters λ soll einem weiteren Parameter unterliegen, dem Mischverh¨altnis. Das Mischverh¨altnis rm ist als Verh¨altnis zwischen der Fl¨ache An der n Ellipsen (Einschl¨usse im 55
Material) und der Fl¨ache A des Modells definiert. F¨ur den Fall, dass das Material nur genau eine Ellipse e1 enth¨alt, ist die Wahrscheinlichkeit, dass ein Punkt des 2D FE Modells genau in dieser liegt, πEh2x πEhx Ehy Ae1 pe 1 = = = A A A mit den Erwartungswerten Ehx und Ehy f¨ur die beiden Halbachsen der Ellipse. Ehx und Ehy sind identisch, da f¨ur beide Halbachsen wie oben erw¨ahnt die gleichen Voraussetzungen gelten. Somit ist die Wahrscheinlichkeit, dass ein Punkt genau nicht in dieser Ellipse liegt πEh2x . A F¨ur n Ellipsen ist somit die Wahrscheinlichkeit, dass ein Punkt in keiner dieser liegt n πEh2x n qn = qe1 ∗ · · · ∗ qen = (qe1 ) = 1 − , A da die Ellipsen unabh¨angig voneinander generiert werden. F¨ur ein vorgegebenes Mischverh¨altnis muss demnach Gleichung n πEh2x = rm pn = 1 − qn = 1 − 1 − A qe1 = 1 − pe1 = 1 −
erf¨ullt sein. Nach Ehx umgeformt, ergibt sich r √ (1 − n 1 − rm )A (9) Ehx = . π Mit dem Erwartungswert EE˜ der verwendeten bedingten exponentiellen Verteilung aus 4.1.1 gilt wegen EE˜ = Ehx 2 xmax ! x 2x 2 1 = 0. λ2 e−λx − − 2 − 3 Ehx − F (xmax ) − F (xmin ) λ λ λ xmin Da diese Gleichung nicht nach λ aufgel¨ost werden kann, wird mittels Bisektionsverfahren zur Nullstellenbestimmung ein passendes λ ∈ [ 5E1h , E5h ] f¨ur das Halbachsenintervall [
Ehx , 3Ehx ] 3
x
x
bestimmt.
6.2.2 Simulation und Bewertung Die Simulation eines Beispielmodells wird durch das FE Programm ParaFEP durchgef¨uhrt, siehe Abschnitt 5.3.2. Die Bewertung extrahiert aus dem Simulationsergebnis die maximale Verschiebung. 56
Geometrie L¨ange 10mm Breite 10mm Fl¨ache 10mm ∗ 10mm = 100mm2 Vernetzung #Elemente 1024 Material #Kreise 3 #Parameter pro Kreis 3, davon 2 variabel Mischverh¨altnis 0.5 xM gleichverteilt in [0mm, 10mm] yM gleichverteilt in [0mm, 10mm] rK 2.56256mm Tabelle 4: Spezifikationen des Quader Modells. 6.2.3 Quader Das von dem Modell aus Abschnitt 4.3 abgeleitete Modell eines Quaders unterliegt grundlegender Einschr¨ankungen. Die r¨aumliche Materialverteilung wird durch drei Kreise mit vorgegebenem einheitlichem Radius bestimmt. So ist ein Einschluss u¨ ber drei Parameter definiert, den Koordinaten xM und yM des Kreismittelpunktes ~xM sowie dem Radius rK . Lediglich xM und yM sind variabel, wodurch die Symmetriegruppe S2 vorliegt. Somit ergibt sich eine Dimension von d = 3 ∗ 2 = 6. Die inversen Verteilungsfunktionen zur Skalierung der Mittelpunktskoordinaten werden, wie in Abschnitt 6.2.1 beschrieben, gew¨ahlt. Der Radius rK ist der Erwartungswert aus Gleichung 9 im Abschnitt 6.2.1 mit einem Mischverh¨altnis rm = 0.5. Tabelle 4 faßt die Spezifikationen des Quader Modells zusammen. Im Folgenden soll die statistische Auswertung der maximalen Verschiebung des durch ParaFEP simulierten Quader Modells betrachtet werden. Abbildung 36 zeigt das simulierte Modell des ersten Integrationspunktes der Monte Carlo Methode mit dem Betrag der Verschiebung und u¨ berlagerter r¨aumlicher Materialverteilung. In Abbildung 37 ist der relative Fehler des Erwartungswertes f¨ur die maximale Verschiebung u¨ ber der Anzahl der durch das jeweilige hochdimensionale Quadraturverfahren ausgewerteten Integrationspunkte aufgetragen. Zur Bestimmung des relativen Fehlers wurden durch die Monte Carlo Methode 1962000 Integrationspunkte ausgewertet. Der Referenzwert ergab sich aus der Mittelung u¨ ber die sich aus den letzten 10000 Integrationspunkten ergebenen Integralswerte. Relative Fehler unterhalb einer Schranke von 10−3 sind aufgrund der Mittelung als 57
unsicher zu betrachten. Die Monte-Carlo Methode und beide Quasi-Monte Carlo Formeln zeigen a¨ hnlich gute Konvergenzverhalten. Dies a¨ ndert sich allerdings in h¨oheren Dimensionen, siehe Abschnitt 6.2.4. Abbildung 38 zeigt einige aus Integrationspunkten resultierende r¨aumliche Materialverteilungen. Der Smolyak Algorithmus konvergiert hingegen trotz Symmetrieausnutzung schlecht, hierzu sp¨ater mehr. F¨ur den Smolyak Algorithmus der Ordnungen 6, 8 und 12 mit GaussLegendre Quadraturformel ergaben sich sogar negative Integralwerte. Daher ist dieser nicht in der Abbildung 37 aufgef¨uhrt. Abbildung 39 zeigt die Verteilung der maximalen Verschiebung, generiert aus 1962000 ausgewerteten Integrationspunkten durch die Monte Carlo Methode.
Abbildung 36: simuliertes Quader Modell mit dem Betrag der Verschiebung, vergr¨oßert um den Faktor 100. Umrandung zeigt unskaliertes Modell mit r¨aumlicher Materialverteilung, generiert aus dem ersten Integrationspunkt der Monte Carlo Methode.
Das schlechte Konvergenzverhalten des Smolyak Algorithmus wird im Folgenden zuerst f¨ur den allgemeineren Fall der Ellipsen erl¨autert. Ist das Verh¨altnis l zwischen der Ordnung l des Smolyak Algorithmus und der Dimension d klein, d so sind die meisten Komponenten eines Integrationspunktes ~ω identisch. Dadurch resultieren viele identische Ellipsen oder Ellipsen, die sich nur in wenigen Parametern unterscheiden. So kann das Mischverh¨altnis rm im Durchschnitt nicht erreicht werden. Je gr¨oßer das Verh¨altnis dl ist, umso mehr Integrationspunkte ~ω haben mehr unterschiedliche Komponenten. Dennoch sind auch die ung¨unstigen Integrationspunkte vorhanden und gehen mitunter in großer Anzahl in die Integrals58
1 Monte Carlo Methode Quasi-Monte Carlo Methode mit arithmetischer Punktfolge Quasi-Monte Carlo Methode mit Halton-Punktfolge Smolyak Algorithmus mit Clenshaw-Curtis Quadraturformel Smolyak Algorithmus mit verzoegerter Kronrod-Patterson Quadraturformel
0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06 0
20000
40000
60000 80000 100000 120000 #(ausgewerteter) Integrationspunkte
140000
160000
180000
Abbildung 37: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte bei 3 Kreisen.
Abbildung 38: R¨aumliche Materialverteilungen, generiert aus dem 153. Integrationspunkt der Monte Carlo Methode, der Quasi-Monte Carlo Methode mit arithmetischer Punktfolge und zu guter letzt mit Halton-Punktfolge (von links nach rechts).
59
1 Verteilung der maximalen Verschiebung 0.9 0.8 0.7
Verteilung
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.005
0.01
0.015 0.02 maximale Verschiebung u
0.025
0.03
0.035
Abbildung 39: Verteilung der maximalen Verschiebung bei 3 Kreisen. approximation ein. Zudem kann l nicht beliebig hoch gew¨ahlt werden, da mit l die Anzahl an Integrationspunkten stark ansteigt. Desweiteren steigt mit wachsendem d bei gleicher Ordnung l und gleicher Folge eindimensionaler Quadraturformeln die Anzahl an Integrationspunkten polynomial an, sofern die erste eindimensionale Quadraturformel der durch den Smolyak Algorithmus genutzten Folge eindimensionaler Quadraturformeln nur genau einen Punkt umfasst. Ist dies nicht so, steigt die Anzahl der Integrationspunkte mit wachsendem d sogar exponentiell [14]. Um viele, in der Anzahl dominierende, Integrationspunkte mit vielen unterschiedlichen Komponenten zu erhalten, m¨ussten die Smolyak Ordnung l und die Anzahl der Punkte der eindimensionalen Quadraturformeln hoch gew¨ahlt werden. Doch dies ist aus den vorher genannten Gr¨unden auch unter Symmetrieausnutzung hier nicht umsetzbar. Da der Smolyak Algorithmus f¨ur sehr kleine Verh¨altnisse dl sogar instabiles Verhalten zeigt, ist dieser f¨ur die hohen Dimensionen und die Symmetriegruppe S5 aus Abschnitt 6.2.4 ungeeignet. Abbildung 40 zeigt einige aus dem Smolyak Algorithmus der Ordnung 7 mit Clenshaw-Curtis Quadraturformel resultierende r¨aumliche Materialverteilungen f¨ur das oben beschriebene Quader Modell. Die anschließenden Tabellen basieren auf dem Smolyak Algorithmus mit Clenshaw-Curtis Quadraturformel. F¨ur das oben beschriebene Modell der Symmetriegruppe S2 zeigt Tabelle 5 f¨ur zwei Kreise, also der Dimension 60
Abbildung 40: R¨aumliche Materialverteilungen, generiert aus dem 2., 11., 1211. und 1212. Integrationspunkt (von links nach rechts) des Smolyak Algorithmus der Ordnung 7 mit ClenshawCurtis Quadraturformel.
d = 4, die Anzahlen unterschiedlicher Ellipsen zu verschiedenen Ordnungen l des Smolyak Algorithmus. F¨ur die Ordnung 4 resultiert aus 5 auszuwertenden Integrationspunkten jeweils ein Kreis. Aus 66 auszuwertenden Integrationspunkten ergeben sich jeweils zwei unterschiedliche Kreise. Von den 66 Integrationspunkten liefern 34 zwei Kreise, die sich lediglich in einer Koordinate unterscheiden. Die restlichen 32 Integrationspunkte ergeben zwei Kreise, die sich in beiden Koordinaten unterscheiden. Im Ausblick auf die nachfolgenden Beispielmodelle der Symmetriegruppe S5 im Abschnitt 6.2.4 sind in Tabelle 6 zu unterschiedlichen Smolyak Ordnungen l und der Dimension 10 die Anzahlen der Ellipsen mit unterschiedlich vielen verschiedenen Parametern dargestellt. Diese Tabelle ist al¨ lerdings noch zu optimistisch, da sie vollst¨andige Uberdeckungen von Ellipsen nicht einbezieht. Denn unterscheiden sich zwei Ellipsen beispielsweise lediglich in der Halbachse, so u¨ berdeckt die Ellipse mit der gr¨oßeren Halbachse die zweite vollst¨andig. Bei zwei Ellipsen die sich in zwei Parametern unterscheiden ist ¨ nur dann eine vollst¨andige Uberdeckung der einen durch die andere ausgeschlossen, wenn diese sich in beiden Koordinaten des Mittelpunktes unterscheiden. In Testl¨aufen f¨ur die Dimension 50 und die Symmetriegruppe S5 ergab sich bei Verwendung der Clenshaw-Curtis Quadraturformel bei einer Smolyak Ordnung von 6 und 60000 Auswertungen ein negatives Integral bei positivem Integranden. Generell konnten bis zur Ordnung 7 mit 418550 Auswertungen nur stark vom exakten Wert abweichende Ergebnisse erzielt werden. 6.2.4 Stoßlasche Die maximale Verschiebung eines durch ParaFEP simulierten Modells einer Stoßlasche soll nun, basierend auf dem abstrakteren Modell aus Abschnitt 4.3 statistisch ausgewertet werden. 61
Ordnung 1 (#Auswertungen = 1) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 1 0 0 2 0 0 0 Ordnung 2 (#Auswertungen = 5) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 1 0 0 2 4 4 0 Ordnung 3 (#Auswertungen = 23) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 5 0 0 2 18 10 8 Ordnung 4 (#Auswertungen = 71) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 5 0 0 2 66 34 32 Ordnung 5 (#Auswertungen = 207) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 13 0 0 2 194 80 114 Ordnung 6 (#Auswertungen = 559) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 13 0 0 2 546 200 346 Ordnung 7 (#Auswertungen = 1479) #untersch. #ausgewert. untersch. Parameter Kreise Punkte 1 2 1 29 0 0 2 1450 452 998 Tabelle 5: Anzahlen der Kreise mit unterschiedlich vielen verschiedenen Parametern zu unterschiedlichen Smolyak Ordnungen bei Dimension 6 und Symmetriegruppe S2 . 62
#untersch. Kreise 1 2 #untersch. Kreise 1 2 #untersch. Kreise 1 2 #untersch. Kreise 1 2 #untersch. Kreise 1 2 #untersch. Kreise 1 2 #untersch. Kreise 1 2
Ordnung 1 (#Auswertungen = 1) #ausgewert. untersch. Parameter Punkte 1 2 3 4 1 0 0 0 0 0 0 0 0 0 Ordnung 2 (#Auswertungen = 11) #ausgewert. untersch. Parameter Punkte 1 2 3 4 1 0 0 0 0 10 10 0 0 0 Ordnung 3 (#Auswertungen = 116) #ausgewert. untersch. Parameter Punkte 1 2 3 4 11 0 0 0 0 105 25 80 0 0 Ordnung 4 (#Auswertungen = 796) #ausgewert. untersch. Parameter Punkte 1 2 3 4 11 0 0 0 0 785 145 320 320 0 Ordnung 5 (#Auswertungen = 4431) #ausgewert. untersch. Parameter Punkte 1 2 3 4 61 0 0 0 0 4370 350 1620 1760 640 Ordnung 6 (#Auswertungen = 20663) #ausgewert. untersch. Parameter Punkte 1 2 3 4 61 0 0 0 0 20602 1190 5380 9040 4480 Ordnung 7 (#Auswertungen = 85833) #ausgewert. untersch. Parameter Punkte 1 2 3 4 241 0 0 0 0 85592 2780 18260 35560 24640
5 0 0
.
5 0 0
.
5 0 0
.
5 0 0
.
5 0 0
.
5 0 512
.
5 0 4352
.
Tabelle 6: Anzahlen der Ellipsen mit unterschiedlich vielen verschiedenen Parametern zu unterschiedlichen Smolyak Ordnungen bei Dimension 10 und Symmetriegruppe S5 .
63
Geometrie L¨ange 405mm Breite 175mm Fl¨ache 405mm ∗ 175mm = 70875mm2 Vernetzung #Elemente 870 Material #Parameter pro Ellipse 5 Mischverh¨altnis 0.5 xM gleichverteilt in [0mm, 405mm] yM gleichverteilt in [0mm, 175mm] φ gleichverteilt in [0, π2 ] Tabelle 7: dimensionsunabh¨angige Spezifikationen der Stoßlasche. Die statistische Auswertung erfolgt f¨ur die drei Dimensionen 50, 100 und 150. Mit der Symmetriegruppe S5 , also f¨unf variablen Parametern pro Ellipse, werden so die Materialien entsprechend der Dimensionen durch 10, 20 und 30 Ellipsen definiert. Zur Bestimmung des Verteilungsparameters λ und der Intervallgrenzen der in Abschnitt 4.1.1 vorgestellten bedingten exponentiellen Verteilung zur Skalierung der Halbachsen wird in jeder Dimension das Mischverh¨altnis 0.5 betrachtet, siehe auch Abschnitt 6.2.1. Unabh¨angig von der durch die Dimensionszahl vorgegebene Anzahl an Ellipsen sollen diese somit im Durchschnitt 50 Prozent der Modellfl¨ache u¨ berdecken. Genauer wird allerdings als Fl¨ache der Stoßlasche die einh¨ullende Rechtecksfl¨ache betrachtet, da die Ellipsenmittelpunkte ohne Ber¨ucksichtigung der nicht definierten Bereiche des 2D FE Modells innerhalb des umgebenden Rechtecks skaliert werden. Zu Beginn werden die verschiedenen Dimensionen einzeln betrachtet. F¨ur eine gegebene Dimension wird die Konvergenz der unterschiedlichen hochdimensionalen Verfahren gepr¨uft. Dabei werden die Monte Carlo Methode und zwei Quasi-Monte Carlo Formeln, zum Einen basierend auf der Halton-Punktfolge, zum Anderen auf der arithmetischen Punktfolge aus Abschnitt 2.2.3, verwendet. Die Hammersley-Punktfolge findet aufgrund ihrer Abh¨angigkeit von der Anzahl der Integrationspunkte keine Anwendung. Der Smolyak Algorithmus wird aus den in Abschnitt 6.2.3 ausf¨uhrlich beschriebenen Gr¨unden nicht in Betracht gezogen. Im Anschluss werden die Ergebnisse der unterschiedlichen Dimensionen miteinander verglichen. Tabelle 7 faßt die dimensionsunabh¨angigen Spezifikationen des 2D FE Modells zusammen. Die verbleibenden dimensionsabh¨angigen Spezifikationen beschr¨anken sich haupts¨achlich auf die Skalierung der Halbachsen. F¨ur die Dimension 50 beschreiben 10 Ellipsen das zu generierende Materi64
1 bedingte exponentielle Verteilung
0.8
0.6
0.4
0.2
0 20
40
60
80
100
Abbildung 41: bedingte exponentielle Verteilungsfunktion zur Skalierung der Halbachsen bei 10 Ellipsen.
Abbildung 42: simuliertes Modell einer Stoßlasche mit dem Betrag der Verschiebung und u¨ berlagerter r¨aumlicher Materialverteilung, generiert aus dem ersten Integrationspunkt der Monte Carlo Methode.
65
al. Die Halbachsen hx und hy liegen im Intervall [12.9563mm, 116.607mm], verteilt mit der inversen bedingten exponentiellen Verteilungsfunktion aus Abschnitt 4.1.1 mit dem Verteilungsparameter λ = 0.0582708 und dem zugeh¨origen Erwartungswert 38.8689mm, siehe Abbildung 41. In Abbildung 42 ist das simulierte FE Modell mit dem Betrag der Verschiebung und zugeh¨origer u¨ berlagerter r¨aumlicher Materialverteilung, generiert aus dem ersten Integrationspunkt der Monte Carlo Methode, zu sehen. Abbildung 43 zeigt den relativen Fehler des Erwartungswertes f¨ur die maximale Verschiebung, aufgetragen u¨ ber der Anzahl der durch das jeweilige hochdimensionale Quadraturverfahren ausgewerteten Integrationspunkte. F¨ur den relativen Fehler wurde zur Bestimmung einer Referenz u¨ ber die letzten 10000 Auswertungen aus insgesamt 4772000 Auswertungen durch die Monte Carlo Methode gemittelt. Aufgrund der Mittelung sind die relativen Fehler unterhalb einer Schranke von 10−3 als unsicher zu betrachten. Die auf der Halton-Punktfolge basierende Quasi-Monte Carlo Formel zeigt das schlechteste Konvergenzverhalten. Dies erk¨art sich folgend. Zur Generierung der i-ten Komponente eines d-dimensionalen Halton-Punktes wird als Radix die i-te Primzahl pi verwendet. F¨ur jeweils pi Halton-Punkte in Folge generiert das Halton-Z¨ahlwerk f¨ur diese Komponente monoton steigende Werte und setzt dann wieder auf einen kleinen Wert zur¨uck. Innerhalb dieser Folge von pi Halton-Punkten unterscheiden sich zwei aufeinanderfolgende Halton-Punkte somit in der Komponente i umso weniger, je gr¨oßer i und somit pi ist. Folglich setzt zudem das HaltonZ¨ahlwerk bei vorgegebener Anzahl an Integrationspunkten umso seltener zur¨uck, je gr¨oßer pi ist. Somit a¨ ndern sich bei Dimensionszuwachs immer mehr Ellipsen u¨ ber immer l¨angere Punktfolgen von Punkt zu Punkt nur unmerklich, wodurch die Konvergenz weiter nachl¨aßt bzw. wom¨oglich ein anderer Integralwert angestrebt wird. Abbildung 44 veranschaulicht die Arbeitsweise des Halton-Z¨ahlwerks anhand der aus aufeinanderfolgenden Halton-Punkten resultierenden r¨aumlichen Materialverteilungen. Da mit der Generierung des ersten Halton-Punktes f¨ur jede Komponente die monoton steigende Wertefolge beginnt, entsteht die in Abbildung 43 zu erkennende Anfangsschwierigkeit, die mit Dimensionszuwachs aus den vorher genannten Gr¨unden zunimmt. Abbildung 45 veranschaulicht das Problem des Z¨ahlbeginns anhand der sich aus den ersten Halton-Punkten ergebenden r¨aumlichen Materialverteilungen. Die Problematik der Halton-Punktfolge besteht aufgrund der a¨ hnlichen Punktgenerierung auch f¨ur die Hammersley-Punktfolge. Die auf der arithmetischen Punktfolge basierende Quasi-Monte Carlo Formel als auch die Monte Carlo Methode zeigen gutes Konvergenzverhalten. Tendenzielle Aussagen u¨ ber eine Bevorzugung eines dieser beiden Methoden k¨onnen nicht gemacht werden. Die gute Konvergenzeigenschaft der Monte Carlo Methode kann darauf begr¨undet sein, dass der Integrand Unstetigkeiten an den Stellen aufweist, an den ein Integrationspunkt der FE Methode auf dem Rand einer Ellipse liegt. Die Abbildungen 46 und 47 zeigen die aus ausgew¨ahlten Integrationspunkten resultie66
renden r¨aumlichen Materialverteilungen. Die Verteilung der maximalen Verschiebung, generiert aus 4772000 ausgewerteten Integrationspunkten durch die Monte Carlo Methode, ist in Abbildung 52 gezeigt. 1 Monte Carlo Quasi-Monte Carlo mit arithmetischer Punktfolge Quasi-Monte Carlo mit Halton-Punktfolge 0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 43: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte bei 10 Ellipsen.
Nebenbei sei erw¨ahnt, dass die Monte Carlo Methode mit einem auf der Methode rand aus der C++-Standardbibliothek basierenden Punktgenerators gegen einen ab der f¨unften Stelle fehlerhaften Wert zu konvergieren schien, siehe Abbildung 48. Dies wurde zudem mit der Quasi-Monte Carlo Formel mit arithmetischer Punktfolge u¨ berpr¨uft. Daher wurde auf die Methode ranf der Bibliothek RANDLIB [2] zur¨uckgegriffen. F¨ur die Dimension 100 wird das zu generierende Material durch 20 Ellipsen beschrieben. Dabei liegen die Halbachsen hx und hy im Intervall [9.24051mm, 83.1646mm], verteilt mit der inversen bedingten exponentiellen Verteilungsfunktion aus Abschnitt 4.1.1 mit dem Verteilungsparameter λ = 0.0817025 und dem zugeh¨origen Erwartungswert 27.7215mm. In Abbildung 49 ist der relative Fehler des Erwartungswertes f¨ur die maximale Verschiebung, aufgetragen u¨ ber der Anzahl der durch das jeweilige hochdimensionale Quadraturverfahren ausgewerteten Integrationspunkte dargestellt. Eine Referenz zur Bestimmung des relativen Fehlers ergab sich aus der Mittelung u¨ ber die letzten 10000 67
Abbildung 44: R¨aumliche Materialverteilungen, generiert aus dem 1985., 1986. und 1987. Halton-Punkt (von links nach rechts) der Dimension 50.
Abbildung 45: R¨aumliche Materialverteilungen, generiert aus dem 4., 10. und 153. HaltonPunkt (von links nach rechts) der Dimension 50.
Abbildung 46: R¨aumliche Materialverteilungen, generiert aus dem 10., 1986. und 1987. Punkt der arithmetischen Folge (von links nach rechts) mit Dimension 50.
68
Abbildung 47: R¨aumliche Materialverteilungen, generiert aus dem 10., 1986. und 1987. Punkt der Monte Carlo Methode (von links nach rechts) mit Dimension 50.
0.1 ’Monte Carlo mit rand()’ ’Monte Carlo mit ranf()’ 0.01
relativer Fehler
0.001
0.0001
1e-05
1e-06
1e-07
1e-08 0
200000
400000
600000 800000 #Integrationspunkte
1e+06
1.2e+06
1.4e+06
Abbildung 48: scheinbar unterschiedliche Konvergenz der Monte Carlo Methode bei unterschiedlichen Punktgeneratoren.
69
Auswertungen von insgesamt 3592000 Auswertungen durch die Monte Carlo Methode. Dadurch sind relative Fehler unterhalb einer Schranke von 10−3 als unsicher zu betrachten. Die Monte Carlo Methode als auch die auf der arithmetischen Folge basierende Quasi-Monte Carlo Formel zeigen gutes Konvergenzverhalten. Die Konvergenz der auf der Halton-Punktfolge basierenden Quasi-Monte Carlo Formel hat sich aufgrund der oben beschriebenen Zusammenh¨ange weiter verschlechtert. Abbildung 52 zeigt die Verteilung der maximalen Verschiebung, generiert aus 3592000 Auswertungen durch die Monte Carlo Methode. 1 Monte Carlo Quasi-Monte Carlo mit arithmetischer Punktfolge Quasi-Monte Carlo mit Halton-Punktfolge 0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06
1e-07 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 49: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte bei 20 Ellipsen.
F¨ur die Dimension 150 ist die r¨aumliche Materialverteilung durch 30 Ellipsen definiert. Die Halbachsen hx und hy sind im Intervall [7.56656mm, 68.0991mm] mit der inversen bedingten exponentiellen Verteilungsfunktion aus Abschnitt 4.1.1 mit dem Verteilungsparameter λ = 0.0997775 und dem zugeh¨origen Erwartungswert 22.6997mm verteilt. In Abbildung 50 ist der relative Fehler des Erwartungswertes f¨ur die maximale Verschiebung u¨ ber der Anzahl der Integrationspunkte aufgetragen. Zur Bestimmung des relativen Fehlers ergab sich eine Referenz durch Mittelung u¨ ber die letzten 10000 Auswertungen der insgesamt 3434000 Auswertungen durch die Monte Carlo Methode. Dadurch gelten relative Fehler unterhalb einer Schranke von 10−3 als unsicher. Auch in dieser Dimension zei70
gen die Monte Carlo Methode und die auf der arithmetischen Folge basierende Quasi-Monte Carlo Formel gute Konvergenzeigenschaften. Hingegen hat sich die auf der Halton-Punktfolge basierende Quasi-Monte Carlo Methode aus den oben genannten Gr¨unden zunehmend verschlechtert. Abbildung 51 zeigt die aus den ersten Halton-Punkten resultierenden r¨aumlichen Materialverteilungen. In Abbildung 52 ist die Verteilung der maximalen Verschiebung, generiert aus 3434000 Auswertungen durch die Monte Carlo Methode aufgezeigt. 1 Monte Carlo Quasi-Monte Carlo mit arithmetischer Punktfolge Quasi-Monte Carlo mit Halton-Punktfolge 0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 50: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte bei 30 Ellipsen.
Abbildung 52 zeigt die Verteilungen der maximalen Verschiebung bei 10, 20 und 30 Ellipsen. Abschließend sollen nun die verwendeten unterschiedlichen hochdimensionalen Quadraturformeln separat betrachtet werden. Die Monte Carlo Methode als auch die auf der arithmetischen Folge basierende Quasi-Monte Carlo Formel zeigen in allen untersuchten Dimensionen gute Konvergenzeigenschaften und scheinen von dem Dimensionszuwachs unbetroffen zu sein, siehe Abbildungen 53 und 54. Die auf der Halton-Punktfolge basierende Quasi-Monte Carlo Formel b¨ußt aus den oben genannten Gr¨unden mit steigender Dimension zunehmend ein, siehe Abbildung 55.
71
Abbildung 51: R¨aumliche Materialverteilungen, generiert aus dem 153., 1211. und 1987. Halton-Punkt (von links nach rechts) der Dimension 150.
1 Verteilung der maximalen Verschiebung bei 10 Ellipsen Verteilung der maximalen Verschiebung bei 20 Ellipsen Verteilung der maximalen Verschiebung bei 30 Ellipsen
0.9 0.8 0.7
#Verteilung
0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.01
0.02
0.03 0.04 #maximale Verschiebung u
0.05
0.06
0.07
Abbildung 52: Verteilungen der maximalen Verschiebung bei 10, 20 und 30 Ellipsen.
72
0.1 Dimension 50 Dimension 100 Dimension 150 0.01
relativer Fehler
0.001
0.0001
1e-05
1e-06
1e-07 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 53: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte durch die Monte Carlo Methode f¨ur die Dimensionen 50, 100 und 150.
73
0.1 Dimension 50 Dimension 100 Dimension 150
relativer Fehler
0.01
0.001
0.0001
1e-05 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 54: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte durch die auf der arithmetischen Folge basierenden Quasi-Monte Carlo Methode f¨ur die Dimensionen 50, 100 und 150.
74
1 Dimension 50 Dimension 100 Dimension 150 0.1
relativer Fehler
0.01
0.001
0.0001
1e-05
1e-06 0
10000
20000
30000
40000 50000 60000 #Integrationspunkte
70000
80000
90000
100000
Abbildung 55: relativer Fehler des Erwartungswertes der maximalen Verschiebung u¨ ber der Anzahl der ausgewerteten Integrationspunkte durch die auf der Halton-Punktfolge basierenden Quasi-Monte Carlo Methode f¨ur die Dimensionen 50, 100 und 150.
75
7
Zusammenfassung
Das in dieser Diplomarbeit entwickelte System zur statistischen Auswertung numerischer Simulationen basiert auf der Approximation von Integralen u¨ ber dem Raum Ω, siehe Abschnitt 3. Das System besteht aus zwei Teilsystemen, den Quadraturalgorithmen und der objektiven Funktion. Die Teilsysteme sind voneinander entkoppelt, so dass eines dieser leicht ausgetauscht werden kann ohne das andere anpassen zu m¨ussen. Die Quadraturalgorithmen liefern das Fundament zur numerischen Quadratur und statistischen Auswertung. Die objektive Funktion f , verkn¨upft mit einer Statistikfunktion, bildet den Integranden. Als Quadraturverfahren sind die Monte-Carlo Methode, Quasi-Monte Carlo Formeln, der Smolyak Algorithmus sowie ein modifizierter Smolyak Algorithmus zur Ausnutzung von Symmetrien implementiert, siehe Abschnitte 2 und 5.1. Die konkreten Beispiele objektiver Funktionen aus Abschnitt 6.2 umfassen jeweils die Simulation eines 2D FE Modells mit anschießender Extraktion der maximalen Verschiebung. Aus einem Integrationspunkt ~ω generiert eine objektive Funktion f die r¨aumliche Materialverteilung des 2D FE Modells, die durch n Einschl¨usse beschrieben wird. Die Simulation wird durch das FE Programm ParaFEP durchgef¨uhrt. Aus dem Simulationsergebnis wird die maximale Verschiebung extrahiert und als Funktionswert f (~ω ) an das Quadraturverfahren zur¨uckgegeben. Auf diese Weise konnte f¨ur unterschiedliche Symmetriegruppen und unterschiedliche Anzahlen an Einschl¨ussen der Erwartungswert der maximalen Verschiebung durch das in dieser Diplomarbeit implementierte System statistisch ausgewertet werden. Desweiteren wurden die verschiedenen implementierten Quadraturverfahren in den Dimensionen 6, 50, 100 und 150 innerhalb der konkreten Beispiele aus Abschnitt 6.2 auf Konvergenz gepr¨uft. Der Smolyak Algorithmus zeigte trotz Symmetrieausnutzung in der Dimension 6 das mit Abstand schlechteste Konvergenzverhalten. Ab der Dimension 50 war dieser sogar instabil, siehe Abschnitt 6.2.3. Die Quasi-Monte Carlo Formel mit Halton-Punktfolge hatte in der Dimension 6 gute Konvergenzeigenschaften, ließ aber mit Dimensionszuwachs zunehmend nach und war ab einer Dimension von 150 g¨anzlich unbrauchbar, siehe Abschnitte 6.2.3 und 6.2.4. Hingegen konvergierten die Monte-Carlo Methode als auch die auf der arithmetischen Folge basierende Quasi-Monte Carlo Formel gut und zeigten sich von dem Dimensionszuwachs unbeeindruckt, siehe Abschnitte 6.2.3 und 6.2.4. In dieser Diplomarbeit wurden mehr als 15000000 2D FE Modelle mit einer Vernetzung von 870 und 1024 Elementen simuliert und statistisch ausgewertet.
76
8
Ausblick
Derzeit werden noch keine beliebigen Statistikfunktionen, wie beispielsweise Kovarianzen, von dem in dieser Arbeit implementierten System (siehe Abschnitt 3) unterst¨utzt. Zwar hat der Experte generell die M¨oglichkeit dies zu umgehen und eine solche Funktionalit¨at umzusetzen, indem er die objektive Funktion entsprechend definiert. Dann w¨are allerdings die objektive Funktion von der statistischen Funktion nicht mehr entkoppelt. F¨ur ein flexibles System ist daher eine Erweiterung um die Unterst¨utzung beliebiger Statistikfunktionen w¨unschenswert. Desweiteren k¨onnte die Effizienz des Smolyak Algorithmus durch eine Berechnung der Gewichte mittels der in [15] vorgestellten Baumstruktur gesteigert werden. Zu guter letzt sei erw¨ahnt, dass geplant ist, die in dieser Arbeit implementierten Quadraturverfahren in das am Institut f¨ur Wissenschaftliches Rechnen der Technischen Universit¨at Braunschweig entwickelte Software-Framework PLATON [8] zu integrieren.
9
Danksagung
Im Rahmen dieser Diplomarbeit wurde ich von den Herren Professor Matthies und Dr. Niekamp unterst¨utzt, bei denen ich mich hiermit herzlichst bedanken m¨ochte. Die wertvollen Ratschl¨age des Herrn Professor Matthies waren sehr hilfreich. Ich habe davon profitiert. Das Engagement des Herrn Dr. Niekamp wirkte zus¨atzlich motivierend. Er war stets geduldig, opferte viel Zeit und war jederzeit zu einem Gespr¨ach bereit. Sein umfangreiches Wissen und seine hohe Kompetenz habe ich zu sch¨atzen gelernt.
Literatur [1] Helmut Brass. Quadraturverfahren. Vandenhoeck & Ruprecht in G¨ottingen und Z¨urich, 1977. [2] Barry W. Brown and James Lovato. RANDLIB - Library of C Routines for Random Number Generation. [3] H.-J. Bungartz and S. Dirnstorfer. Multivariate Quadrature on Adaptive Sparse Grids. Computing, 71:89–114, 2003.
77
[4] Chin-Yun Chen. Adaptive numerische Quadratur und Kubatur mit automatischer Ergebnisverfikation. PhD thesis, Institut f¨ur Angewandte Mathematik, Universit¨at Karlsruhe, 1998. [5] Thomas Gerstner and Michael Griebel. Numerical Algorithms, volume 18, pages 209–232. Kluwer Academic Publishers, 1998. Numerical integration using sparse grids. [6] Thomas Gerstner and Michael Griebel. Dimension-Adaptive TensorProduct Quadrature. Computing, 71, September 2003. [7] Stefan Heinrich and Alexander Keller. Quasi-Monte Carlo Methods in Computer Graphics, Part I: The QMC-Buffer. Technical Report 242, Fachbereich Informatik, AG Numerische Algorithmen, Universit¨at Kaiserslautern, Postfach 3049, D-67653 Kaiserslautern, 1994. [8] Krosche, Markus and Niekamp, Rainer and Matthies, Hermann G. PLATON: A Problem Solving Environment for Computational Steering of Evolutionary Optimisation on the Grid. In EUROGEN 2003 ISBN: 84-9599933-1 (published on CD). [9] Hermann G. Matthies. Computational Aspects of Probability in Non-Linear Mechanics. Technical report, Institute of Scientific Computing, Technische Universit¨at Braunschweig, Brunswick, Germany, 2004. [10] Giovanni Monegato. Numerical Algorithms, volume 26, pages 173–196. Kluwer Academic Publishers, 2001. An overview of the computational aspects of Kronrod quadrature rules. [11] R. Niekamp. Verteilte Algorithmen f¨ur h-, p- und d-adaptive Berechnungen in der nichtlinearen Strukturmechanik. Prof. Dr.-Ing P. Wriggers, Institut f¨ur Baumechanik und Numerischer Mechanik, Universit¨at Hannover, 2000. [12] Erich Novak. Numerische Verfahren f¨ur Hochdimensionale Probleme und der Fluch der Dimension. Jahresbericht der DMV 101 (1999), pages 151– 177, Juli 1999. [13] Knut Petras. Software - smolpack and finpack. bs.de:8080/ petras/software.html.
http://www-public.tu-
[14] Knut Petras. Asymptotically minimal Smolyak Cubature. Technical report, Institut f¨ur Angewandte Mathematik, Technische Universit¨at Braunschweig, Pockelsstraße 14, D-38106 Braunschweig, Germany, 2000.
78
[15] Knut Petras. Fast Calculation of Coefficients in the Smolyak Algorithm. Technical report, Institut f¨ur Angewandte Mathematik, Technische Universit¨at Braunschweig, Pockelsstraße 14, D-38106 Braunschweig, Germany, 2000. [16] Knut Petras. Numerische Methoden in der Finanzmathematik. Skript zur vorlesung, Institut f¨ur Angewandte Mathematik, Abteilung f¨ur Numerische Mathematik, Technische Universit¨at Braunschweig, Pockelsstrasse 14, 38106 Braunschweig, 2001/2002. [17] Leszek Plaskota and Grzegorz W. Wasilkowski. Smolyak’s Algorithm for Integration and L1 -Approximation of Multivariate Functions with Bounded Mixed Derivatives of Second Order. Technical report, Department of Mathematics, Informatics and Mechanics, University of Warsaw, Banacha 2, 02-097 Warsaw, Poland, Department of Computer Science, University of Kentucky, 773 Anderson Hall, Lexington, KY 40506-0046, USA, January 2004. [18] Prof. Dr. Thomas Sauerbier. Statistik f¨ur Wirtschaftswissenschaftler. R. Oldenbourg Verlag M¨unchen Wien, 2000. [19] Edwin Miles Stoudenmire. Monte Carlo Integration: an Overview with Examples. Technical Report MATH 4255, Georgia Tech, April 2004. [20] Andreas Keese und Hermann G. Matthies. Fragen der numerischen Integration bei stochastischen finiten Elementen f¨ur nichtlineare Probleme. Informatikbericht 4, Institut f¨ur Wissenschaftliches Rechnen, Technische Universit¨at Braunschweig, Braunschweig, Germany, Mai 2003. [21] Andreas Keese und Hermann G. Matthies. Numerical Methods and Smolyak Quadrature for Nonlinear Stochastic Partial Differential Equations. Informatikbericht 5, Institute of Scientific Computing, Technical University Braunschweig, Brunswick, Germany, Mai 2003. [22] Eric Weisstein. mathworld. http://mathworld.wolfram.com. [23] Saul A. Teukolsky William H. Press, Brian P. Flannery and William T. Vetterling. Numerical Recipes in C: the Art of Scientific Computing. Cambridge University Press, 1988-1992.
79
