add num8

2020-06-26 20:19:29 +02:00
parent 333cd04e58
commit 0c6d5e130e
8 changed files with 499 additions and 0 deletions
@@ -0,0 +1,91 @@
+#include <iostream>
+#include <vector>
+#include "hdnum.hh"
+
+
+template<class N>
+class PolynomialProblem
+{
+public:
+  // Exportiere Größentyp
+  typedef std::size_t size_type;
+
+  // Exportiere Zahlentyp
+  typedef N number_type;
+
+  // Dimension der untenstehenden Funktion
+  std::size_t size () const
+  {
+    return 1;
+  }
+
+  // Funktionsauswertung
+  void F (const hdnum::Vector<N>& x, hdnum::Vector<N>& result) const
+  {
+    result[0] = pow(x[0], 4) - 3.0*pow(x[0], 3) + 5.0*x[0] + 3.0;
+  }
+
+  // Jacobimatrix
+  void F_x (const hdnum::Vector<N>& x, hdnum::DenseMatrix<N>& result) const
+  {
+    result[0][0] = 4.0*pow(x[0], 3) - 9.0*pow(x[0], 2) + 5.0;
+  }
+};
+
+
+int main ()
+{
+  typedef std::complex<double> Number; // Der Zahlentyp soll hier komplex sein
+
+  typedef PolynomialProblem<Number> Problem;
+  Problem problem;
+
+
+  // Wir schauen Startpunkte auf einem Raster an, das später in ein Bild umgesetzt werden soll
+  for (double start_value_i = -5.0; start_value_i <= 5.0; start_value_i += 1e-2) {
+
+    for (double start_value = -5.0; start_value <= 5.0; start_value += 1e-2) {
+
+      hdnum::Vector<Number> u(problem.size()); // Vektor für Startpunkt bzw Lösung
+      u[0] = Number(start_value, start_value_i); // Startpunkt
+
+      hdnum::Newton newton;
+      newton.set_maxit(20);
+      newton.set_verbosity(0);
+      newton.set_reduction(1e-10);
+      newton.set_abslimit(1e-20);
+      newton.set_linesearchsteps(3);
+
+      newton.solve(problem,u); // Lösen
+
+      double id = 0;
+      if(newton.has_converged()) {
+
+        // Falls unser Newton-Löser konvergiert ist, stellen wir fest welche Nullstelle
+        // Nullstelle wir getroffen haben und legen einen Index für jede Nullstelle fest.
+
+        if (u[0].real() < 0 && u[0].imag() < 0)
+          id = 1;
+        else if (u[0].real() < 0 && u[0].imag() > 0)
+          id = 2;
+        else if (u[0].real() > 0 && u[0].imag() < 0)
+          id = 3;
+        else if (u[0].real() > 0 && u[0].imag() > 0)
+          id = 4;
+        else
+          std::cout << "Unbekannte Nullstelle!" << std::endl;
+
+        // Damit auch der "Abstand" zur Nullstelle (im Sinne von benötigten Iterationen
+        // um sie zu erreichen) sichtbar wird, addieren wir noch die Iterationszahl
+        // um einen Faktor abgeschwächt
+        id += .1 * (double)newton.iterations();
+      }
+
+      std::cout << id << ",";
+    }
+
+    std::cout << std::endl;
+  }
+
+  return 0;
+}