init

P2/CMakeLists.txt

@@ -0,0 +1,12 @@
+# Set the minimum required version of CMake
+cmake_minimum_required(VERSION 3.10)
+
+# Set the project name and specify the C++ as the project language
+project(HM_P2)
+
+# Specify the C++ standard
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+# Add an executable with the above sources
+add_executable(HM_P2 main.cpp CMyVektor.cpp CMyMatrix.cpp)

P2/CMyMatrix.cpp

@@ -0,0 +1,145 @@
+#include "CMyMatrix.h"
+#include "CMyVektor.h"
+#include <iostream>
+#include <iomanip>
+#include <math.h>
+
+// Konstruktor
+CMyMatrix::CMyMatrix(int z, int s) {
+    zeilen = z;
+    spalten = s;
+    matrix.resize(zeilen);
+    for (int i = 0; i < zeilen; ++i)
+        matrix[i].resize(spalten);
+}
+
+double CMyMatrix::getZeilen() {
+    return zeilen;
+}
+
+double CMyMatrix::getSpalten() {
+    return spalten;
+}
+
+void CMyMatrix::setWert(double wert, int zeile, int spalte) {
+    matrix[zeile][spalte] = wert;
+}
+
+double CMyMatrix::getWert(int zeile, int spalte) {
+    return matrix[zeile][spalte];
+}
+
+double& CMyMatrix::operator()(int zeile, int spalte) {
+    return matrix[zeile][spalte];
+}
+
+CMyMatrix CMyMatrix::invers() {
+    if (this->getZeilen() != 2 || this->getSpalten() != 2)
+        throw std::runtime_error("Keine 2x2 Matrix");
+
+    double a = this->getWert(0, 0);
+    double b = this->getWert(1, 0);
+    double c = this->getWert(0, 1);
+    double d = this->getWert(1, 1);
+
+    double determinate = a * d - b * c;
+
+    if (determinate == 0)
+        throw std::runtime_error("detA == 0");
+
+    double kehrwertDeterminante = 1 / determinate;
+    CMyMatrix inverse(2,2);
+
+    inverse.setWert(kehrwertDeterminante * d, 0, 0);
+    inverse.setWert(kehrwertDeterminante * -c, 0, 1);
+    inverse.setWert(kehrwertDeterminante * -b, 1, 0);
+    inverse.setWert(kehrwertDeterminante * a, 1, 1);
+
+    return inverse;
+}
+
+// Matrix-Vektor Multiplikation
+CMyVektor operator*(CMyMatrix A, CMyVektor x) {
+    // M Spalten == V Zeilen
+    if (A.getSpalten() != x.getDimension())
+        throw std::runtime_error("Matrix(m x n) * Vektor(i) n != i");
+
+    CMyVektor ergebnis(A.getZeilen());
+
+    for (int i = 0; i < A.getZeilen(); i++) {                                            // Zeile
+        for (int j = 0; j < A.getSpalten(); j++) // Spalte
+            ergebnis[i] += A(i, j) * x[j];       // Aufsummierung
+    }
+    return ergebnis;
+}
+
+std::ostream& operator<<(std::ostream& os, CMyMatrix x) {
+
+    for (int i = 0; i < x.getZeilen(); ++i) {
+        os << "\t\t\t(";
+        for (int j = 0; j < x.getSpalten(); ++j)
+            j == x.getSpalten() - 1 ? os << std::fixed << std::setprecision(10) << x(i, j) : os << std::fixed << std::setprecision(10) << x(i, j) << "|";
+
+        os << ")\n";
+    }
+    return os;
+}
+
+CMyMatrix jacobi(CMyVektor x, CMyVektor(*f)(CMyVektor)) {
+    double const h = pow(10, -4);
+    CMyVektor f_x = f(x);
+    CMyVektor tmp{ x };
+
+    CMyMatrix jacobiMatrix(f_x.getDimension(), x.getDimension());
+
+    for (int i = 0; i < f_x.getDimension(); ++i) { // Zeile bzw. Teilfunktion f_i
+        for (int j = 0; j < x.getDimension(); ++j) { // Spalte bzw. Variable x_j
+            tmp[j] += h;
+            jacobiMatrix(i, j) = (f(tmp)[i] - f_x[i]) / h;
+            tmp[j] -= h;
+        }
+    }
+    return jacobiMatrix;
+};
+
+void newtonverfahren(CMyVektor x, CMyVektor(*f)(CMyVektor)) {
+    int schritt = 0;
+    CMyVektor f_x = f(x);
+    CMyMatrix jacobi_f_x = jacobi(x, f);
+    CMyMatrix jacobiInvers = jacobi_f_x.invers();
+
+    double funktionLaenge = f_x.length();
+
+    CMyVektor dx = -1 * (jacobiInvers * f_x);
+
+    while (schritt != 50 && funktionLaenge >= 1e-5) {
+        std::cout << "Schritt " << schritt << ":"
+        << "\n\t            𝐱 = " << x
+        << "\n\t         𝑓(𝐱) = " << f_x << "\n"
+        << "\n\t        ∂𝑓(𝐱) = \n " << jacobi_f_x
+        << "\n\t      ∂𝑓(𝐱)⁻¹ = \n" << jacobiInvers
+        << "\n\t           Δ𝐱 = " << dx
+        << "\n\t     ||𝑓(𝐱)|| = " << funktionLaenge << std::endl << std::endl;
+
+        schritt++;
+        x = x + dx;
+        f_x = f(x);
+        jacobi_f_x = jacobi(x, f);
+        jacobiInvers = jacobi_f_x.invers();
+        dx = -1 * (jacobiInvers * f_x);
+        funktionLaenge = f_x.length();
+
+        if (schritt == 50) {
+            std::cout << "Ende wegen Schrittanzahl = 50 bei";
+            std::cout << "\n\t             x = " << x
+                << "\n\t          𝑓(𝐱) = " << f_x
+                << "\n\t      ||𝑓(𝐱)|| = " << funktionLaenge << std::endl;
+        }
+        else if (funktionLaenge < 1e-5) {
+            std::cout << "Ende wegen ||𝑓(𝐱)||<1e-5 bei";
+            std::cout << "\n\t             x = " << x
+                << "\n\t          𝑓(𝐱) = " << f_x
+                << "\n\t      ||𝑓(𝐱)|| = " << funktionLaenge << std::endl;
+        }
+    }
+};

P2/CMyMatrix.h

@@ -0,0 +1,31 @@
+#include <vector>
+#include <iostream>
+#include "CMyVektor.h"
+
+class CMyMatrix {
+private:
+    std::vector<std::vector<double>> matrix;
+    int zeilen = 0;
+    int spalten = 0;
+
+public:
+    CMyMatrix(int zeilen, int spalten);
+
+    double getZeilen();
+    double getSpalten();
+
+    void setWert(double wert, int zeile, int spalte);
+
+    double getWert(int zeile, int spalte);
+
+    double& operator()(int zeile, int spalte);
+
+    CMyMatrix invers();
+
+    friend CMyVektor operator*(CMyMatrix A, CMyVektor x);
+
+    friend std::ostream& operator<<(std::ostream& os, CMyMatrix x);
+};
+
+CMyMatrix jacobi(CMyVektor x, CMyVektor(*funktion)(CMyVektor));   //como friend-> error de compiler "C3767"
+void newtonverfahren(CMyVektor x, CMyVektor(*funktion)(CMyVektor));  //libro pagina 221

P2/CMyVektor.cpp

@@ -0,0 +1,199 @@
+#include "CMyVektor.h"
+#include <math.h>
+#include <iostream>
+
+int CMyVektor::getDimension() {
+    return dimension;
+}
+
+// Vektor Element bei index i
+double& CMyVektor::operator[](int i) {
+    return werte.at(i);
+}
+
+// Betrag des Vektors
+double CMyVektor::length() {
+    double l = 0.0;
+    for (int i = 0; i < werte.size(); i++)
+        l += pow(werte.at(i), 2);
+
+    return(sqrt(l));
+}
+
+// Ausgabe im Vektorformat
+std::ostream& operator<< (std::ostream& os, CMyVektor& vektor) {
+    os << "(";
+    for (int i = 0; i < vektor.getDimension(); i++) {
+        os << vektor[i];
+        if (i < vektor.getDimension() - 1) os << "|";
+    }
+    os << ")";
+
+    return os;
+}
+
+/* Vektor Addition
+ a1	+	b1	=	a1+b1
+ a2	+	b2	=	a2+b2
+ :		:		  :
+ ai	+	bi	=	ai+bi
+*/
+CMyVektor operator+(CMyVektor a, CMyVektor b) {
+    // Vektoren addierbar?
+    if (a.getDimension() != b.getDimension())
+        return a;
+
+    CMyVektor sum(a.getDimension());
+
+    for (int i = 0; i < a.getDimension(); i++)
+        sum[i] = a[i] + b[i];
+
+    return sum;
+}
+
+/* Skalare Multiplikation
+ a1	*	s	=	v1
+ a2	*	s	=	v2
+ :		:		 :
+ ai	*	s	=	vi
+*/
+CMyVektor operator*(double s, CMyVektor a) {
+    CMyVektor vektor(a.getDimension());
+
+    for (int i = 0; i < a.getDimension(); i++)
+        vektor[i] = s * a[i];
+
+    return vektor;
+}
+
+CMyVektor gradient(double f(CMyVektor x), CMyVektor x) {
+    CMyVektor grad(x.getDimension());
+    CMyVektor xh(x.getDimension()); 	// x mit "Wackelei" an x_i
+    xh = x;
+
+    double fx = f(x); 
+    double h = 1e-8;
+
+    // Berechnung der partiellen Ableitung nach x_i und Zusammenfassung in grad
+    for (int i = 0; i < x.getDimension(); i++) {
+        //  xh[0]	            xh[1]              ...  xh[n]  
+        //  x_1+h,x_2,...,x_n   x_1,x_2+h,...,x_n       x_1,x_2,...,x_n+1
+        xh[i] += h;
+
+        // grad[i] = df/dx_i => numerische Ableitung nach dem jeweilige x_i
+        grad[i] = (f(xh) - fx) / h;
+
+        // Zurücksetzen der "Wackelei"
+        xh[i] -= h;
+    }
+    return grad;
+}
+
+CMyVektor gradientenverfahren(double f(CMyVektor x), CMyVektor x, double lambda) {
+    int schritt_zaehler = 0;
+    CMyVektor gradfx(x.getDimension());
+    CMyVektor x_neu(x.getDimension());
+    gradfx = gradient(f, x);
+
+    while (gradient(f, x).length() >= 1e-5 && schritt_zaehler < 25) {
+        x_neu = x + lambda * gradfx;
+
+        std::cout
+            << "Schritt " << schritt_zaehler << ": " << std::endl
+            << "               x = " << x << std::endl
+            << "               λ = " << lambda << std::endl
+            << "            f(x) = " << f(x) << std::endl
+            << "       grad f(x) = " << gradfx << std::endl
+            << "   ||grad f(x)|| = " << gradfx.length() << std::endl << std::endl
+            << "           x_neu = " << x_neu << std::endl
+            << "        f(x_neu) = " << f(x_neu) << std::endl << std::endl;
+
+        // Halbierung
+        if (f(x_neu) <= f(x)) {
+            double lamda_test = lambda * 0.5;
+            CMyVektor x_test = x + lamda_test * gradfx;
+
+            std::cout
+                << "     ↯ f(x) = " << f(x) << " ≥ f(x_neu) = " << f(x_neu) << std::endl << std::endl
+                << "     ? Test mit halbierter Schrittweite (λ = " << lamda_test << "): " << std::endl
+                << "          x_test = " << x_test << std::endl
+                << "       f(x_test) = " << f(x_test) << std::endl << std::endl;
+
+            if (f(x) < f(x_test))
+                std::cout
+                << "     ✓ f(x) = " << f(x) << " < f(x_test) = " << f(x_test) << std::endl
+                << "     ! Übernehme Schrittweite." << std::endl << std::endl;
+
+            // weiter Halbieren wenn f(x_test) <= f(x)
+            while (f(x_test) <= f(x)) {
+                std::cout
+                    << "     ↯ f(x) = " << f(x) << " ≥ f(x_test) = " << f(x_test) << std::endl << std::endl;
+
+                lamda_test *= 0.5;
+                x_test = x + lamda_test * gradfx;
+
+                std::cout
+                    << "     ? Test mit halbierter Schrittweite (λ = " << lamda_test << "): " << std::endl
+                    << "          x_test = " << x_test << std::endl
+                    << "       f(x_test) = " << f(x_test) << std::endl << std::endl;
+
+                // Ausgabe beim letzten Durchlauf der Schleife
+                if (!(f(x_test) <= f(x)))
+                    std::cout
+                    << "     ✓ f(x) = " << f(x) << " < f(x_test) = " << f(x_test) << std::endl << std::endl
+                    << "     ! Übernehme Schrittweite λ = " << lamda_test << std::endl << std::endl;
+            }
+
+            lambda = lamda_test;
+            x_neu = x_test;
+        }
+
+        // Verdopplung
+        else {
+            double lamda_test = lambda * 2;
+            CMyVektor x_test = x + lamda_test * gradfx;
+
+            std::cout
+                << "     ? Test mit doppelter Schrittweite (λ = " << lamda_test << "): " << std::endl
+                << "          x_test = " << x_test << std::endl
+                << "       f(x_test) = " << f(x_test) << std::endl << std::endl;
+
+            if (f(x_test) > f(x_neu)) {
+                x_neu = x_test;
+                lambda = lamda_test;
+                std::cout
+                    << "     ! Übernehme verdoppelte Schrittweite λ = " << lamda_test << std::endl
+                    << "     ✓ f(x_neu) = " << f(x_neu) << " < f(x_test) = " << f(x_test) << std::endl << std::endl;
+            }
+            else {
+                std::cout
+                    << "     ↯ f(x_neu) = " << f(x_neu) << " ≥ f(x_test) = " << f(x_test) << std::endl
+                    << "     ! behalte aktuelle Schrittweite λ = " << lambda << std::endl << std::endl;
+            }
+        }
+        std::cout << std::endl;
+
+        x = x_neu;
+        gradfx = gradient(f, x);
+        schritt_zaehler++;
+    }
+
+    // Zusammenfassung des Endes 
+    if (gradfx.length() < 1e-5)
+        std::cout
+        << "Ende wegen ||grad f(x)|| < 10^-5 bei" << std::endl;
+
+    else
+        std::cout
+        << "Ende wegen 25. Schritt" << std::endl;
+
+    std::cout
+        << "               x = " << x << std::endl
+        << "               λ = " << lambda << std::endl
+        << "            f(x) = " << f(x) << std::endl
+        << "       grad f(x) = " << gradfx << std::endl
+        << "   ||grad f(x)|| = " << gradfx.length() << std::endl
+        << "-----------------------------------------------------" << std::endl << std::endl;
+
+    return x_neu;
+}

P2/CMyVektor.h

@@ -0,0 +1,34 @@
+#pragma once
+#include <vector>
+#include <iostream>
+
+class CMyVektor {
+private:
+	int dimension;
+	std::vector<double> werte;
+public:
+	// Vektor mit Dimension anlegen
+	CMyVektor(int n) : dimension{ n } {
+		werte.resize(n);
+	}
+
+	CMyVektor(std::vector<double> x) : dimension{ (int)x.size() } {
+		werte.resize(x.size());
+		for (int i = 0; i < x.size(); i++)
+			werte[i] = x[i];
+	}
+
+	int getDimension();
+	double length();
+
+	// Komponente von Vektor übergeben
+	double& operator[](int index);
+};
+
+CMyVektor operator+(CMyVektor a, CMyVektor b);
+CMyVektor operator*(double lambda, CMyVektor a);
+
+std::ostream& operator<< (std::ostream& os, CMyVektor& vektor);
+
+CMyVektor gradient(double f(CMyVektor x), CMyVektor x);
+CMyVektor gradientenverfahren(double f(CMyVektor x), CMyVektor x, double lambda = 1.0);

P2/main.cpp

@@ -0,0 +1,41 @@
+#include <math.h>
+#include "CMyVektor.h"
+#include "CMyMatrix.h"
+
+using namespace std;
+
+// f( [x,y] ) = sin(x * y) + sin(x) + cos(y)
+double ff(CMyVektor X) {
+	double x = X[0];
+	double y = X[1];
+
+	return (sin(x * y) + sin(x) + cos(y));
+}
+
+// g( [x,y,z] ) = -(2x^2 - 2xy + y^2 + z^2 - 2x - 4z)
+double gg(CMyVektor X) {
+	double x = X[0];
+	double y = X[1];
+	double z = X[2];
+
+	return -(2 * pow(x, 2) - 2 * x * y + pow(y, 2) + pow(z, 2) - 2 * x - 4 * z);
+}
+
+CMyVektor f(CMyVektor x) {
+	CMyVektor result(3); // da f: R^4 -> R^3
+	result[0] = x[0] * x[1] * exp(x[2]);
+	result[1] = x[1] * x[2] * x[3];
+	result[2] = x[3];
+	return result;
+}
+
+int main() {
+
+	cout << "Aufgabe_2 :\n" << jacobi(CMyVektor{ {1, 2, 0, 3} },
+		[](CMyVektor x) -> CMyVektor { return { {x[0] * x[1] * exp(x[2]), x[1] * x[2] * x[3], x[3]} }; })
+		<< endl;
+
+	newtonverfahren({ {1, 1} },
+        [](CMyVektor x) -> CMyVektor { return { {pow(x[0], 3) * pow(x[1], 3) - 2 * x[1], x[0] - 2} }; });
+    return 0;
+}