Implementación de Algoritmos de Jacobi y Método de la Potencia en C++

Implementación de Algoritmos de Álgebra Lineal

1. Método de Jacobi

El siguiente código implementa el método de Jacobi para la diagonalización de matrices simétricas:

Array1D<real> jacobi(Array2D<real>& A, Array2D<real>& Autovectores, const int NMaxIter, const real TOL, int &Niter) {
  /** COMPROBAMOS QUE LAS DIMENSIONES SON COHERENTES */
  // CODIGO HEDIONDO

  /** HACER ALUMNO */
  Array1D<real> autovalores(N);
  Autovectores = Array2D<real>(N, N, 0.);
  for(int i = 0; i < N; i++) Autovectores[i][i] = 1.;

  for (Niter = 0; Niter < NMaxIter; Niter++) {
    int p = 0;
    int q = 1;
    real Max = 0;
    for (int i = 0; i < N; i++) {
      for (int j = i + 1; j < N; j++) {
        if (fabs(A[i][j]) > Max) {
          Max = fabs(A[i][j]);
          p = i;
          q = j;
        }
      }
    }
    if (Max < TOL) {
      for(int l = 0; l < N; l++) autovalores[l] = A[l][l];
      return (autovalores);
    }
  }

  real alfa = 0.5 * atan((2. * A[p][q]) / (A[q][q] - A[p][p]));
  real co = cos(alfa);
  real si = sin(alfa);

  for(int k = 0; k < N; k++) {
    real temp = Autovectores[k][p];
    Autovectores[k][p] = co * temp - si * Autovectores[k][q];
    Autovectores[k][q] = co * Autovectores[k][q] + si * temp;
    if(k == p || k == q) continue;
    A[p][p] = co * (co * App - si * A[p][q]) - si * (co * A[p][q] - si * A[q][q]);
    A[q][q] = co * (co * A[q][q] + si * A[p][q]) + si * (co * A[p][q] + si * App);
    A[p][q] = A[q][p] = 0.;
  }
  return(Array1D<real>());
}

2. Método de la Potencia

Función para calcular el autovalor dominante y su autovector asociado:

int mn_metodo_potencia(Array2D<real> &A, Array1D<real> &u_max, real &l_max, int Nmax, real TOL) {
  /// HACER ALUMNO
  l_max = mn_norma_euclidea(u_max);
  u_max = u_max / l_max;
  for(int i = 0; i < Nmax; i++) {
    Array1D<real> u_new = A * u_max;
    real l_new = mn_norma_euclidea(u_new);
    if(l_new == 0) return(-2);
    if(mn_producto_escalar(u_new, u_max) < 0.) l_max = -l_new;
    else l_max = l_new;
    u_new = u_new / l_max;
    if(mn_norma_euclidea(u_new - u_max) < TOL) {
      u_new = u_max;
      return(i);
    }
    u_max = u_new;
  }
  return(-3);
}