Prof. Dr. Marc Steinbach M.Sc. Jan Thiedau, Dipl.-Math. Jens Hübner Hannover, 7. Mai 2015 1. Programmierblatt zu Numerische Mathematik II – Nichtlineare Ausgleichsprobleme – Das Ziel dieses Programmierblatts ist die Implementation einer Lösungsroutine für nichtlineare Ausgleichsprobleme. Dazu implementieren Sie zuächst eine QR-Faktorisierung und damit eine Lösungsroutine für lineare Ausgleichsprobleme, die Sie dann für den Löser für nichtlineare Ausgleichsprobleme verwenden können. (a) Implementieren Sie für eine Matrix A ∈ Rm×n , m ≥ n eine QR-Faktorisierung A = QR mit Q ∈ O(m) und R ∈ Rm×n obere Dreiecksmatrix (s. Vorlesung). Zur Berechnung der Faktorisierung können Sie entweder Householder-Spiegelungen oder Givens-Rotationen verwenden. Ihre QR-Faktorisierung sollte folgenden Funktionskopf haben: [Q,R] = factorize_qr(A) Hinweise: Sie sollen Q explizit berechnen. Sie können Ihre Implementation mit der MATLAB-Routine qr vergleichen. (b) Implementieren Sie eine Lösungsroutine für lineare Ausgleichsprobleme der Form min kAx − ak22 x∈Rn unter Cx − c = 0 (LAP) mit A ∈ Rna ×n , na ≥ n, a ∈ Rna , C ∈ Rnc ×n , nc ≤ n, rank(C) = nc , c ∈ Rnc . Berechnen Sie dazu in einem ersten Schritt die Faktorisierungen der benötigten Matrizen mit einer Funktion [Lc,Qc,A1,A2,Ra,Qa] = factorize_LAP(A,C) Es soll die LQ-Faktorisierung von C, die Aufteilung von A in A1 , A2 und die QR-Faktorisierung von A2 zurückgegeben werden. Verwenden Sie hierzu Ihre QR-Faktorisierung. Im zweiten Schritt soll dann mit Hilfe der Faktorisierungen das lineare Ausgleichsproblem gelöst werden: x = solve_LAP(a,c,Lc,Qc,A1,A2,Ra,Qa) Hinweise: Zum Lösen der Dreieckssysteme können Sie die MATLAB-Routine linsolve mit den geeigneten Optionen verwenden. Zur Berechnung der LQ-Faktorisierung der Beschränkungsmatrix C können Sie die QR-Faktorisierung auf die Transponierte C T anwenden – Sie benötigen dafür keine zusätzliche Faktorisierungsmethode. Zum Überprüfen der Richtigkeit Ihrer Implementation können Sie Ihre Methode mit der MATLAB-Routine lsqlin vergleichen. (c) Implementieren Sie das Gauß-Newton-Verfahren für nichtlineare Ausgleichsprobleme min kFa (x)k22 x∈Rn unter Fc (x) = 0, (NLAP) mit Fa ∈ C 1 (Rn , Rna ) und Fc ∈ C 1 (Rn , Rnc ). Verwenden Sie dabei zum Lösen der linearisierten Probleme (LAP(x)) Ihren Löser solve_LAP. Der Funktionskopf Ihrer Routine sollte folgende Form haben: [x, dx, iter] = GGN_method(Fa, Fc, x0, tol) Dabei sind Fa und Fc die Funktionen zum Auswerten der Zielfunktion beziehungsweise der Beschränkungsfunktion, x0 der Startwert und tol die Fehlertoleranz für das Abbruchkriterium. Die Rückgabewerte sind zum einen die Lösung x des nichtlinearen Ausgleichsproblems (NLAP), der letzte Schritt sowie die Anzahl Iteration bis zum Erreichen des Abbruchkriteriums. Hinweise: • Benutzen Sie als Abbruchkriterium für das Gauß-Newton-Verfahren die Norm des aktuellen Schrittes ∆x. • Für die Implementation von nichtlinearen Funktionen können Sie MATLAB-Funktionen schreiben und mit Hilfe des function handle (@) an Ihr Gauß-Newton-Verfahren übergeben. Implementieren Sie die Funktionen so, dass jeweils Funktionswert und Ableitung gemeinsam berechnet werden. [f, jacobian] = func(x) Mögliche Beispiele • Probleme aus Übungsaufgabe 3.1 (Beschränkte lineare Probleme) • Beispiel aus dem Skript Numerik 1 Kapitel 4.5 in der nichtlinearen und linearen Form. Die Daten der Messpunkte sind im Stud.IP verfügbar. • In der euklidischen Ebene seien die Messpunkte (xi , yi ) ∈ R2 mit i xi yi 1 0 1.9900 2 0.125 2.4840 3 0.25 2.6614 4 1.5 2.8660 gegeben. Es soll nun ein Kreis bestimmt werden, so dass die Messpunkte möglichst nahe an der Kreislinie liegen. a) Formulieren Sie diese Aufgabe als nichtlineares Ausgleichsproblem und lösen Sie es mit Ihrer Implementierung des Gauß-Newton-Verfahrens Als Startwert für die Berechnung können Sie einen Kreis mit Mittelpunkt 12 , 32 und Radius R = 43 verwenden. b) Nehmen Sie als zusätzliche Beschränkung an, dass der Mittelpunkt des Kreises auf einem Kreis um (0, 0) mit Radius 2 liegt. Hinweise zur Abgabe • Die Abgabe bis spätestens 01.06.2015 erfolgt per Email an [email protected] mit dem Betreff „Abgabe Programmieraufgabe 1“. In die Email gehören Name, Matrikelnummer, eine vollständige Liste der Dateinamen, die Angabe der verwendeten Programmiersprache (MATLAB, OCTAVE, ...) und gegebenenfalls Kommentare zur Abgabe. Sämtliche Dateien, die Sie abgeben möchten gehören in den Anhang der Email. • Kommentieren Sie Ihre Lösung ausführlich im erstellten Quellcode!