uni/sose2020/num/uebungen/rohrleitungsnetzwerk.cc

// Kompilieren mit:
//
// g++ -I../hdnum/ -o rohrleitungsnetzwerk rohrleitungsnetzwerk.cc
//
// Das setzt voraus, dass ihr Programm in einem Ordner parallel zur
// entpackten HDNum Bibliothek liegt.

#include <iostream>
#include <cstdlib>
#include "hdnum.hh"


// Funktion zum Aufstellen der Matrix
template<class NumberType>
void flussMatrix( hdnum::DenseMatrix<NumberType> &A ) {
    int M( A.rowsize() );
    int N( A.colsize() );
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    // Numerierung wie auf Zettel 3 nur mit 0 beginnend
    // also v_0, ... ,v_(N^2-1)
    // der Referenzknoten v_r hat Druck 0

    // berechnung der kantenlaenge
    int n = floor(sqrt(N+1));

    for(int i = 0; i < N; i++) {
        int edges = 0;
        if ((i+1)%n != 0) { // nicht linker rand
            edges++;
            if (i-1 >= 0) A(i, i-1) = -1; // falls nicht der referenzknoten
        }
        if ((i+2)%n != 0) { // nicht rechter rand
            edges++;
            A(i, i+1) = -1;
        }
        if (i+n < N) { // nicht unterer rand
            edges++;
            A(i, i+n) = -1;
        }
        if (i-n >= -1) { // nicht oberer rand
            edges++;
            if (i-n >= 0) A(i, i-n) = -1; // falls nicht der referenzknoten
        }
        A(i, i) = edges;
    }
}


// Funktion zur Berechnung der Frobenius-Norm einer Matrix
template<class NumberType>
NumberType frobeniusNorm(const hdnum::DenseMatrix<NumberType> &A) {
    // Error checking
    int M(A.rowsize());
    int N(A.colsize());
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    NumberType result=0.0;

    // iteriere ueber alle zeilen und spalten, quadriere die elemente und summiere
    for (int i=0; i < N; i++) {
        for (int j=0; j < N; j++) {
            result += pow(A(i,j), 2);
        }
    }

    // ziehe wurzel aus summe
    return sqrt(result);
}


// Funktion zur Berechnung des betragsgrößten Eigenwertes mit Potenzmethode
template<class NumberType>
NumberType maxEigenwert(const hdnum::DenseMatrix<NumberType> &A) {
    // Error checking
    int M(A.rowsize());
    int N(A.colsize());
    if(M!=N)
        HDNUM_ERROR("Matrix muss quadratisch sein!");

    // start vektor
    hdnum::Vector<NumberType> r(N);
    r[0] = 1;
    // work copy
    hdnum::Vector<NumberType> r_tmp(N);
    hdnum::Vector<NumberType> diff(N);
    // fuehre iterationsschritt 10000 mal aus
    for (int k=0; k<10000; k++) {
        A.mv(r_tmp, r); // r_tmp = Ar
        r_tmp /= r_tmp.two_norm(); // normiere r_tmp
        r = r_tmp;
    }
    A.mv(r_tmp, r);
    // berechne eigenwert mit rayleigh quotient
    return (r * r_tmp)/r.two_norm_2();
}

// Hauptprogramm
int main(int argc, char ** argv)
{

  // Anzahl der Knoten
  const int N(10);
  std::cout << "Knotenanzahl N: " << N << std::endl;

  // Größe der Matrix
  const int n(N*N-1);

  // Datentyp für die Matrix
  typedef double REAL;

  // Matrix initialisieren
  hdnum::DenseMatrix<REAL> A(n,n);

  // Pretty-printing einmal setzen für alle Matrizen
  A.scientific(false);
  A.width(15);

  flussMatrix(A);
  if (N<=4)
    std::cout << A << std::endl;

  // Bei Schwierigkeiten mit Teilaufgabe a) können Sie Teilaufgaben b)
  // und c) mit folgender Matrix testen
  int size_b = 3;
  hdnum::DenseMatrix<REAL> B(size_b,size_b);
  for (std::size_t i=0; i<size_b; ++i)
  {
    for (std::size_t j=0; j<size_b; ++j)
    {
      B[i][j] = i+j;
    }
  }
  std::cout << "Zeilen-Summen-Norm: " << A.norm_infty() << std::endl;
  std::cout << "Spalten-Summen-Norm: " << A.norm_1() << std::endl;
  std::cout << "Frobenius-Norm: " << frobeniusNorm(A) << std::endl;
  std::cout << "Maximaler Eigenwert: " << maxEigenwert(A) << std::endl;

  return 0;
}