Lösung nichtlinearer Gleichungssysteme

1. Lösung nichtlinearer Gleichungssysteme Präsentation in Rahmen des Seminars „Parallele Programmierung“ SS 2003 Referent: Roman Achziger Datum: 11.06.2003

2. Einleitung: Motivation • Anwendungsgebiete • Gebiete der Technik, Wirtschafts- und Naturwissenschaften • Lösung nichtlinearen Abhängigkeiten • Approximation durch lineare Modelle • Einsatz numerischen Algorithmen • Die Standardform nichtlinearer Gleichungssysteme: = oder kurz F(x)

3. Grundlagen der Iterationsverfahren • Fixpunktiteration: Betrachtung des skalaren Falls • Transformation der vorliegende Gleichung f(x)=0 in eine äquivalente Fixpunktform (x) = x Definition: Ein Element x* R heißt Fixpunkt der Abbildung : R R falls (x*)=x* gilt. • Wähle die Abbildung so, dass x* die Abbildung (x*) = x* und f(x*) = 0 erfüllt • Ausgehend vom Startwert x(0) gilt die Rekursion x(k+1) = (x(k)), • Abbildung verkleinert sukzessive den Abstand zwischen x(k+1) und x*

4. Grundlagen der Iterationsverfahren • Spinnweb-Diagramm für (x)= exp(-x2) • Die Folge oszilliert um den Fixpunkt

5. Grundlagen der Iterationsverfahren • Die Wahl einer gute Iterationsfunktion ist entscheidend Definition: Eine auf einer Teilmenge A R definierte Abbildung : A A heißt Kontraktion, wenn eine Lipschitz-Konstante 0<L<1 existiert, so dass für alle Paare x,y A gilt: | (x) - (y)|L |x-y|. L = |´(Z)|, mit einem Zwischenwert Z zwischen x und y  Die Abbildung konvergiert, wenn dieses Maximum der ersten Ableitung kleiner als 1

6. Grundlagen der Iterationsverfahren: Beispiel • Gesucht ist die Nullstelle von f(x) = x3+x-1,5 • Einsetzen: f(0,8) = -0,188 und f(0,9) = 0,129  Da die Funktion stetig ist, existiert eine Nullstelle x* mit 0,8<x*<0,9 1) 1(x)= -x3+1,5 ; | 1´(x)|=3x2 | 1´(0,8)| <| 1´(x*)|<| 1´(0,9)|<=> 1,92 <| 1´(x*)|< 2,43  abstoßender Fixpunkt: x(0)=0,9; x(1)=0,77; x(2)=1,04; x(3)=0,37 2) 2(x)=(-x3+x+1,5)/2 ; | 2´(x)| =| (-3x2+1)/2| | 2´(0,8) |<| 2´(x*)|<| 2´(0,9)| d.h. 0,46<| 2´(x*)|< 0,715  Die Folge konvergiert gegen x* für Startwerte aus [0,8; 0,9] x(0)=0,9; x(1)=0,8355; x(2)=0,8761; x(3)=0,8518; x(4)=0,8667; x(5)=0,8577; x(5)=0,8633

7. Grundlagen der Iterationsverfahren Banachsche Fixpunktsatz: • Es sei : A A eine Kontraktion über der abgeschlossenen Menge AR mit der Kontraktionskonstanten L< 1. Dann : 1) existiert ein eindeutiger Fixpunkt x* = (x*) in A; 2) konvergiert die durch x(k+1) = (x(k)) definierte Folge für jeden Startwert x(0) A gegen den Fixpunkt; 3) gelten die Abschätzungen:

8. Grundlagen der Iterationsverfahren: Einzugsbereich Definition: Die Menge E  R heißt Einzugsbereich der Lösung x*, wenn x(k)= x* für alle x(0) E gilt. • globaleKonvergenz: • Ist der gesamte Definitionsbereich von Einzugsbereich  Iterationsfolge konvergiert für alle Startwerte • lokaleKonvergenz: • Iterationsfolge konvergiert nur auf einer Umgebung um x*

9. Grundlagen der Iterationsverfahren: Konvergenzgeschwindigkeit • Der Fehler (k) := x(k) -x* stellt auf einer kleinen Umgebung von x* das Konvergenzverhalten der Folge dar • Lineare Konvergenz: Es gibt eine Konstante q (0; 1), Konvergenzordnung p=1, so dass für alle k=0, 1, 2, . . . gilt. der Fehler (k+1) nimmt linear mit dem Faktor q ab • Quadratische Konvergenz: Es gibt eine Konvergenzfaktor K > 0 und eine Konvergenzordnung p=2, so dass für alle k=0, 1, 2, . . . gilt. Fehler(k+1) ist proportional zum Quadrat von (k) mit dem Faktor K

10. Grundlagen der Iterationsverfahren: Abbruchkriterien • Kompromiss zwischen der Genauigkeit und der Effizienz • Abbruch der Verfahren beim Erreichen einer geforderten Genauigkeit • Fehler-Kriterium: Prozess stoppt, wenn || x(k+1)-x(k)|| < abs • Residuums-Kriterium: Prozess stoppt, wenn ||F(x)|| < f • Abbruch der Verfahren bei Nicht-Konvergenz • Kriterien: ||x(k)||>xmax oder || F(x(k))|| > fmax oder ||F(x(k))|| > ||F(x(k-1))||> ||F(x(k-2))|| • sicherste Methode: Beschränkung der maximalen Anzahl der Iterationsschritte  Verhinderung einer unendlichen Folge

11. Lösung nichtlinearer Gleichungen: Bisektionsverfahren • Nullstellensatz für stetige Funktionen • Sei eine stetige Funktion f : [a,b] R mit f(a)*f(b) < 0 gegeben. Dann existiert ein x* (a,b), so dass f(x*) = 0 gilt. • Algorithmus: Startintervall: I0 =[x(k-1),x(k)] mit x(k-1) < x(k) und f(x(k-1))f(x(k))< 0 while (Abbruchkriterium nicht erfüllt) { (x(k+1) = x(k-1) + x(k)) /2 if (f(x(k-1))f(x(k+1)) 0 ) x(k)= x(k+1)else x(k-1)= x(k+1) } return x(k) • Lineare Konvergenz mit dem Konvergenzfaktor ½ • Global konvergent  

12. Lösung nichtlinearer Gleichungen: Bisektionsverfahren • Parallelen Realisierung • Unterteilung des Intervalls in in p+1 Teilintervalle • parallelen Berechnung der p Funktionswerte • Dasjenige mit einen Vorzeichenwechsel wird aus den p+1 Teilintervallen ausgewählt • wiederhole bis Abbruchkriterium erfüllt • Nach k-Schritt wird das Anfangsintervall auf das ( p+1)-k –fache reduziert  speed-up von Sp= O(log2 p)

13. Lösung nichtlinearer Gleichungen: Linearisiertes Newton-Verfahren • Ersetzung der Funktion durch eine lineare Modellfunktion • Geometrisch: Der Schnittpunkt der Tangente mit der x-Achse gilt als Näherungswert für die Nullstelle x* • Iterationsvorschrift: x(k+1)= x(k) - , k =0,1,2,..., • Fixpunktiterationen: (x(k))= x(k)- • Newton-Iteration: x(k+1) = (x(k)), k = 0, 1, 2,...., • Eigenschaften: • quadratische Konvergenz • lokale Konvergenz • selbstkorrigierend

14. Lösung nichtlinearer Gleichungssysteme: Fixpunktiteration • Übertragung des skalaren Falls der Fixpunktiteration auf ein n-dimensionales Problem vorliegende Problem: sind : Rn Rn und F: Rn Rn gegeben, finde x* Rn, so dass (x*) = x*. Für einen Fixpunkt x* von  gilt F(x*) = 0 • Ausgehend von x(0) Rn konvergiert die Iterationsvorschrift x(k+1) = ( x(k)), für k = 0, 1, 2..... gegen x*

15. Lösung nichtlinearer Gleichungssysteme:Fixpunktiteration • parallele Realisierung • vergleichbar zur der parallelen Realisierung eines Iterations-verfahrens zur Lösung linearer Gleichungssysteme • Aufteilung der Komponenten i(x(k)) der Kontraktionsabbildung  auf n Prozessoren • blockweise Aufteilung Fixpunktabbildung  auf die Prozessoren jeder einzelne Prozessor berechnet n/p Spalten des nächsten Iterationsvektors • Verteilung des Iterationsvektor per Multi-Broadcastoperation

16. Lösung nichtlinearer Gleichungssysteme: Newton-Verfahren • Linearisierung der Funktion F: Rn Rn • Iterationsvorschrift: x(k+1)= x(k)-[DF(x(k))]-1 F(x(k)) für k=0, 1, 2, . . .. • Fixpunktproblem: (x) = x(k)-[DF(x(k))]-1 F(x(k)) • Fixpunktiteration x(k+1) = (x(k)) DF(x)= • Eigenschaften: • quadratische Konvergenz • lokale Konvergenz

17. Lösung nichtlinearer Gleichungssysteme: Newton-Verfahren • Einführung eines Korrekturterms (k): 1) x(k+1) = x(k) + (k) 2) DF(x(k))*(k) = - F(x(k)) Überführung des Problem des Lösens eines Systems nichtlinearer Gleichungen in das Problem des Lösens eines linearen Gleichungs- systems • Rechenaufwand pro Iterationsschritt: • n2 +n Komponenten-Funktionsauswertung • 2n3/3+O(n2) arithmetischen Operationen • n2 + O(n) Speicheroperationen

18. Lösung nichtlinearer Gleichungssysteme: Modifiziertes Newton-Verfahren • Zyklisches Aktualisieren der Jacobimatrix • LU-Faktorisierung der Jakobimatrix nur alle n Iterationsschritte • Differenzenapproximation der Jacobimatrix • Bestimmung einer Näherung der Jacobimatrix durch n-dimensionale Differenzen • Berechne für jeden Eintrag aus DF(x(k)): (x(k))  Möglichkeit der Wahl für rj und ej: ||F(x(k) + rjej)- F(x(k))|| F(x(k))

19. Lösung nichtlinearer Gleichungssysteme: Parallele Implementierung des Newton-Verfahrens • Verwendung des Ergebnisvektors des k-ten Iterationsschrittes als Datum bei der Berechnung des nächsten Iterationsschrittes  die Iterationen müssen nacheinender durchgeführt werden • Teilaufgaben pro Schritt: • Auswertung der Funktion F • Berechnung der Jacobimatrix • Durchführung der Gauß-Elimination • Berechnung des nächsten Approximationsvektors • Kontrolle der Konvergenz des Verfahrens mit Bestimmung der Norm des Korrekturterms

20. Lösung nichtlinearer Gleichungssysteme: Parallele Implementierung des Newton-Verfahrens • Parallelisierung der Teilaufgaben des Iterationsschrittes unter Beachtung der Datenabhängigkeiten • Datenabhängigkeit in der Newton-Iteration:

21. Lösung nichtlinearer Gleichungssysteme: Parallele Implementierung des Newton-Verfahrens • Bei der parallelen Ausführung auf eine gemeinsame Datenverteilung achten • Gauß-Elimination verursacht den meisten Rechenaufwand  Datenverteilung bzgl. Gauß-Elimination optimieren • Durchführung der Gauß-Elimination auf mehreren Prozessoren: • zeilenzyklische Berechnung • gesamtzyklische Berechnung

22. Lösung nichtlinearer Gleichungssysteme:Parallele zeilenzyklische Implementierung • Verteilung der Jacobimatrix zyklisch auf die vorhandenen Prozessoren  Berechnung der Jacobimatrix in der entsprechenden Verteilung • Approximation der Jacobimatrix durch: (x(k))  • Berechnung der Jacobimatrix: • Ein Prozessor Pq speichert und berechnet alle Zeilen i mit q = ((i-1) mod p) +1  maximal Zeilen • Vektor x(k) wird repliziert gehalten  keine Kommunikation • Gesamtaufwand: T1(n,p)= (n+1) Teval(F)

23. Lösung nichtlinearer Gleichungssysteme:Parallele zeilenzyklische Implementierung • Durchführung der Gauß-Elimination • Zeilenzyklische Ablage der Jacobimatrix • Zeilenzyklische Ablage der rechten Seite des zu lösenden Gleichungssystems  zeilenzyklischen Variante der Gauß-Elimination ohne Umverteilung der Daten • Ein Prozessor speichert die benötigten Zeilen der Jacobimatrix und den entsprechenden Funktionswert ab

24. Lösung nichtlinearer Gleichungssysteme:Parallele zeilenzyklische Implementierung • Berechnung des neuen Iterationsvektors: • Replizierte Ablage des Iterationsvektors x(k) der Newton-Iteration  die Prozessoren führen die Berechnung ohne Kommunikation durch • Jeder Prozessor pberechnet: • Werte der gespeicherten Funktion • Zugehörige Zeilen der Jacobimatrix • Zugehörige Werte des Vektors (k) • Zugehörige Werte des Vektors x(k+1) • Zugehörige Komponenten der Norm ||(k)||

25. Lösung nichtlinearer Gleichungssysteme:Parallele zeilenzyklische Implementierung double *newton_zyklisch (double *F( )) { double *x_neu, *x_alt, *temp, *y, f[MAX_SIZE], fg[MIN_SIZE] ; int i, j ; x_neu = x_buf1; x_old= x_buf2; initialize (x_neu); do { temp=x_neu ; x_neu= x_alt ; x_alt=temp ; for (i=me; i<n; i+=p) { fg[ i ] = -F(i, x_alt); for (j=0; i<n; i++) { x_alt[ j ] += r[ j ]; Berechnung der Jacobimatrix f [ j ] = F(i, x_alt); x_alt[ j ] - = r[ j ]; DF[ i ][ j ] = ( fg[ i ] + f [ j ] ) /r [ j ]; } } y = gauss_zyklisch (DF, fg); zeilenzyklische Gauß-Elimination for ( j=0 ; j<n; j++ ) Berechnung des neuen Iterationsvektors x_neu[ j ] = x_alt[ j ] + y[ j ]; } while (norm(x_alt, x_neu ) > (1-L) /L ); Überprüfung des Abbruchkriteriums return x_neu ; }

26. Lösung nichtlinearer Gleichungssysteme:Parallele gesamtzyklische Implementierung • Berechnung der Jacobimatrix: • Blockzyklische schachbrettartige Datenverteilung derMatrix Bei gleichen Anzahl der Prozessoren in jeder Dimension berechnet ein Prozessor etwa Einträge • Der benötigte Iterationsvektor wird repliziert abgelegt  lokales Arbeiten mit der Funktion F

27. Lösung nichtlinearer Gleichungssysteme:Parallele gesamtzyklische Implementierung • Vermeidung von Kommunikation bei der Berechnung der Jacobimatrix zwischen den Prozessoren  Fi(x(k)) mindestens Mal errechnen, da sich an der Berechnung einer Zeile der Jacobimatrix Prozessoren beteiligen • Es ergibt sich der Gesamtaufwand von: T2(n,p) = n/p( )Teval(F) • Durchführung der Gauß-Elimination: • Allokierung der Jacobimatrix DF vergleichbar mit der zeilenzyklischen Variante • Eine allokierte Zeile enthält nur noch +1 Einträge

28. Lösung nichtlinearer Gleichungssysteme :Vergleich der beiden Varianten • Zeilenzyklische Implementierung: • geringere Anzahl der Funktionsauswertungen • Gesamtzyklische Implementierung: • besseres Kommunikationsverhalten der Gauß-Elimination • Wahl abhängig davon, welcher der beiden Aspekte überwiegt • Auswertungsaufwand der Funktion F  Zeilenzyklische Implementierung • Kommunikationsverhalten der Zielmaschine  Gesamtzyklische Implementierung

29. Zusammenfassung • Vielzahl von unterschiedlichen Iterationsverfahren • Beurteilung anhand des Effizienzvergleich • Anzahl der Iterationsschritte • Rechenaufwand pro Iterationsschritt • Bisektionsverfahren: • für skalaren Fall • einfaches Verfahren mit globaler, lineare Konvergenz • Newton-Verfahren: • für skalaren und auch mehrdimensionalen Fall • quadratische, lokale Konvergenz • Parallele Realisierung: • Verteilung des Berechnungsaufwands auf mehrere Prozessoren

30. Herzlichen Dank für Ihre Aufmerksamkeit! Noch Fragen?