Implementación de Métodos Numéricos para Autovalores y Sistemas Lineales

Algoritmo de Jacobi para Autovalores y Autovectores

Este segmento de código implementa el algoritmo de Jacobi para el cálculo de autovalores y autovectores de una matriz simétrica. El método es iterativo y busca diagonalizar la matriz mediante una serie de rotaciones.

Inicialización de la Matriz de Autovectores

Se inicializa la matriz de autovectores como una matriz identidad, ya que inicialmente representa la base canónica.

Autovectores = 0.;
for(int i = 0; i < dim; i++) Autovectores[i][i] = 1.;

Bucle Principal del Algoritmo de Jacobi

El algoritmo itera hasta que la matriz se considera diagonalizada (elementos fuera de la diagonal menores que una tolerancia) o se alcanza el número máximo de iteraciones.

for (Niter = 0; Niter < NMaxIter; Niter++){
    int p, q;
    real max = 0;

Cálculo del Elemento Mayor Fuera de la Diagonal

En cada iteración, se busca el elemento con el mayor valor absoluto fuera de la diagonal principal. Este elemento determinará el plano de rotación.

    for (int i = 0; i < dim; i++){
        for (int j = i + 1; j < dim; j++){
            if (fabs(A[i][j]) > max){
                max = fabs(A[i][j]);
                p = i;
                q = j;
            }
        }
    }

Criterio de Convergencia

Si el elemento más grande fuera de la diagonal es menor que una tolerancia predefinida (TOL), se considera que el algoritmo ha convergido y se devuelven los autovalores (elementos de la diagonal).

    if (max < TOL) {
        Array1D a(dim);
        for(int l = 0; l < dim; l++) a[l] = A[l][l];
        return (a);
    }

Cálculo del Ángulo de Rotación

Se calcula el ángulo alpha necesario para la rotación de Jacobi, que anulará el elemento A[p][q].

    real alpha = 0.5 * atan((2 * A[p][q]) / (A[q][q] - A[p][p]));

Modificación de la Matriz A

Se aplican las rotaciones de Jacobi a la matriz A para acercarla a su forma diagonal. Esto implica actualizar los elementos de las filas y columnas p y q.

    for (int k = 0; k < dim; k++){
        if(k != p && k != q){
            real D = A[p][k];
            A[k][p] = cos(alpha) * D - sin(alpha) * A[q][k];
            A[p][k] = A[k][p];
            A[k][q] = sin(alpha) * D + cos(alpha) * A[q][k];
            A[q][k] = A[k][q];
        }
    }
    real app = A[p][p];
    A[p][p] = cos(alpha) * (cos(alpha) * A[p][p] - sin(alpha) * A[q][p]) - sin(alpha) * (cos(alpha) * A[p][q] - sin(alpha) * A[q][q]);
    A[q][q] = cos(alpha) * (cos(alpha) * A[q][q] + sin(alpha) * A[p][q]) + sin(alpha) * (cos(alpha) * A[q][p] + sin(alpha) * app);
    A[p][q] = 0;
    A[q][p] = 0;

Modificación de la Matriz de Autovectores

Simultáneamente, se actualiza la matriz de autovectores aplicando las mismas rotaciones para obtener los autovectores correspondientes.

    for(int k = 0; k < dim; k++){
        real r_kp = Autovectores[k][p];
        Autovectores[k][p] = cos(alpha) * Autovectores[k][p] - sin(alpha) * Autovectores[k][q];
        Autovectores[k][q] = sin(alpha) * r_kp + cos(alpha) * Autovectores[k][q];
    }
}

Manejo de Error por Exceso de Iteraciones

Si el bucle de iteraciones finaliza sin alcanzar la convergencia, se imprime un mensaje de error.

printf ("Error, número de iteraciones excedido\n");
return(Array1D());

Nota: El fragmento de código finaliza con un cierre de llave } que sugiere el fin de una función.

Método de Relajación (SOR/Gauss-Seidel)

Este fragmento de código implementa un método iterativo de relajación, comúnmente utilizado para resolver sistemas de ecuaciones lineales. La estructura sugiere una implementación del método de Gauss-Seidel con posible sobre-relajación (SOR) si w es el factor de relajación.

Bucle Principal de Iteraciones

El algoritmo itera hasta un número máximo de iteraciones (nMax) o hasta que el error entre iteraciones consecutivas caiga por debajo de una tolerancia.

Array1D v(u.dim());
for(int n = 1; n <= nMax; n++){
    real error = 0.; // Variable para almacenar el error en la iteración

Recorrido de Incógnitas y Actualización

Se recorren las incógnitas del sistema de forma ascendente, actualizando cada componente del vector solución u.

    for (int i = 0; i < u.dim(); i++){
        real temp = b[i]; // Variable para guardar el nuevo valor de u[i]
        for(int j = 0; j < i; j++){
            temp -= A[i][j] * u[j];
        }
        for(int j = i + 1; j < u.dim(); j++){
            temp -= A[i][j] * u[j];
        }
        temp = w * temp / A[i][i] + (1 - w) * u[i];
        v[i] = u[i];
        u[i] = temp;
    }

Criterio de Convergencia

Después de cada iteración completa, se calcula el error entre el vector solución actual (u) y el anterior (v). Si este error es menor que la tolerancia (TOL), el algoritmo ha convergido y se devuelve el número de iteraciones.

    error = mn_error_vectores(u, v);
    if (error < TOL){
        return n;
    }
}

Manejo de Error por Exceso de Iteraciones

Si el algoritmo no converge dentro del número máximo de iteraciones, se devuelve -1.

return -1;

Nota: El fragmento de código finaliza con un cierre de llave } que sugiere el fin de una función.

Método de la Potencia para el Autovalor Dominante

Este código implementa el método de la potencia, un algoritmo iterativo para encontrar el autovalor de mayor magnitud (autovalor dominante) y su autovector asociado de una matriz.

Inicialización del Autovector y Autovalor Máximo

Se inicializa el autovalor máximo con la norma del autovector inicial y se normaliza este último.

l_max = u_max.norma();
u_max = u_max / l_max;

Bucle Principal de Iteraciones

El algoritmo itera hasta un número máximo de iteraciones (Nmax) o hasta que el autovector converge.

for(int i = 0; i < Nmax; i++){

Generación de un Nuevo Vector

Se crea un nuevo vector multiplicando la matriz A por el autovector actual (u_max).

    Array1D u_new = A * u_max;

Cálculo del Nuevo Autovalor y Normalización

Se aplica la norma euclídea al nuevo vector para obtener una estimación del autovalor dominante. Se incluye una lógica para manejar el signo del autovalor.

    real l_new = u_new.norma();
    Array1D temp = u_max * u_new; // Producto cartesiano (¿producto punto?)
    for(int k = 0; k < temp.dim(); k++){
        if(temp[k] < 0){
            l_max = -l_new;
        } else {
            l_max = l_new;
        }
    }
    u_new = u_new / l_max; // Normalización del nuevo autovector

Criterio de Convergencia

Si la norma de la diferencia entre el autovector actual y el nuevo es menor que una tolerancia (TOL), se considera que el algoritmo ha finalizado y se devuelve el número de iteraciones.

    if((u_max - u_new).norma() < TOL){
        u_max = u_new;
        return i;
    }
    u_max = u_new; // Actualización del autovector para la siguiente iteración
}

Manejo de Error por Exceso de Iteraciones

Si el algoritmo no converge dentro del número máximo de iteraciones, se devuelve -1.

return -1;

Nota: El fragmento de código finaliza con un cierre de llave } que sugiere el fin de una función.

Método de Eliminación Gaussiana con Pivoteo

Este segmento de código implementa el método de eliminación Gaussiana con pivoteo parcial para resolver sistemas de ecuaciones lineales de la forma Ax = b. Es una parte fundamental para el método de la potencia inversa.

Inicialización del Vector de Pivoteo

Se inicializa un vector piv para llevar un registro de los intercambios de filas durante el pivoteo.

for(int k = 0; k < b.dim(); k++) piv[k] = k;

Aplicación del Método de Gauss (Fase de Eliminación)

Se realiza la fase de eliminación hacia adelante, transformando la matriz en una forma triangular superior. Se incluye pivoteo parcial para mejorar la estabilidad numérica.

for(int k = 0; k < b.dim(); k++){
    // Detección del máximo en la diagonal hacia abajo (pivoteo)
    real max = fabs(A1[piv[k]][k]);
    int kmax = k;
    for(int m = k + 1; m < b.dim(); m++){
        if(fabs(A1[piv[m]][k]) > max){
            max = fabs(A1[piv[m]][k]);
            kmax = m;
        }
    }
    if(kmax != k){
        int paso = piv[kmax];
        piv[kmax] = piv[k];
        piv[k] = paso;
    }
    // Conversión a cero de los elementos bajo la diagonal
    for(int j = k + 1; j < b.dim(); j++){
        real mul = -A1[piv[j]][k] / A1[piv[k]][k];
        A1[piv[j]][k] = 0.;
        for(int n = k + 1; n < b.dim(); n++) A1[piv[j]][n] += mul * A1[piv[k]][n];
        b1[piv[j]] += mul * b1[piv[k]];
    }
}

Fase de Remonte (Sustitución Hacia Atrás)

Una vez que la matriz está en forma triangular superior, se resuelve el sistema mediante sustitución hacia atrás.

if(A1[piv[b.dim() - 1]][b.dim() - 1] == 0.) return Array1D(); // Manejo de matriz singular
u[b.dim() - 1] = b1[piv[b.dim() - 1]] / A1[piv[b.dim() - 1]][b.dim() - 1];
for(int k = b.dim() - 2; k >= 0; k--){
    real sum = 0;
    for(int l = k + 1; l < b.dim(); l++) sum += A1[piv[k]][l] * u[l];
    u[k] = (b1[piv[k]] - sum) / A1[piv[k]][k];
}
return(u);

Método de la Potencia Inversa para el Autovalor Mínimo

Este código implementa el método de la potencia inversa, un algoritmo iterativo para encontrar el autovalor de menor magnitud (autovalor mínimo) y su autovector asociado de una matriz. Este método requiere la resolución de un sistema lineal en cada iteración, a menudo utilizando un método como la eliminación Gaussiana.

Inicialización del Autovector y Autovalor Mínimo

Se inicializa el autovalor mínimo con la norma del autovector inicial y se normaliza este último.

l_min = u_min.norma();
u_min = u_min / l_min;

Bucle Principal de Iteraciones

El algoritmo itera hasta un número máximo de iteraciones (Nmax) o hasta que el autovector converge.

for(int i = 0; i < Nmax; i++){

Generación de un Nuevo Vector mediante Eliminación Gaussiana

Se crea un nuevo vector resolviendo el sistema lineal A * u_new = u_min, lo que se logra implícitamente a través de la función mn_gauss.

    Array1D u_new = mn_gauss(A, u_min);

Cálculo del Nuevo Autovalor y Normalización

Se aplica la norma euclídea al nuevo vector para obtener una estimación del autovalor mínimo. Se incluye una lógica para manejar el signo del autovalor.

    real l_new = u_new.norma();
    Array1D temp = u_min * u_new; // Producto cartesiano (¿producto punto?)
    for(int k = 0; k < temp.dim(); k++){
        if(temp[k] < 0){
            l_min = -l_new;
        } else {
            l_min = l_new;
        }
    }
    u_new = u_new / l_min; // Normalización del nuevo autovector

Criterio de Convergencia

Si la norma de la diferencia entre el autovector actual y el nuevo es menor que una tolerancia (TOL), se considera que el algoritmo ha finalizado. El autovalor mínimo se calcula como el inverso del valor obtenido, y se devuelve el número de iteraciones.

    if((u_min - u_new).norma() < TOL){
        l_min = 1. / l_min; // El autovalor es el inverso del obtenido
        u_min = u_new;
        return i;
    }
    u_min = u_new; // Actualización del autovector para la siguiente iteración
}

Manejo de Error por Exceso de Iteraciones

Si el algoritmo no converge dentro del número máximo de iteraciones, se devuelve -1.

return -1;

Nota: El fragmento de código finaliza con un cierre de llave } que sugiere el fin de una función.