podemos apreciar que la mayor complejidad nos lleva rápidamente a costos computacionales muy elevados para este tipo de algoritmos. Retomaremos esta problemática cuando veamos los métodos iterativos.
1.1.4 Métodos directos implementados en Matlab®
En Matlab®, hay varias opciones para resolver los sistemas de ecuaciones lineales:
a) La alternativa más obvia es usar la función que calcula la inversa de una matriz y entonces la solución está dada por: x = inv(A)·b. Sin embargo, esta alternativa está limitada solamente a matrices cuadradas que posean inversa y estén bien condicionadas (ver sección 1.3); además, no es el método más eficiente cuando la matriz posee cierta estructura como ser triangular, simétrica, pocas bandas diagonales y el resto elementos cero, etcétera.
b) Matlab® dispone del operador ‘\’ que corresponde a dividir por la izquierda la ecuación básica A·x = b, la que queda expresada como:x = A\b. Este operador posee mayor funcionalidad que la inversa, ya que se puede aplicar a sistemas sobre- y subdeterminados que poseen más (menos) ecuaciones que incógnitas, por lo que en tal caso se obtiene una solución en el sentido de mínimos cuadrados. El algoritmo intenta primero explorar alguna característica especial de la matriz A; por ejemplo, si es simétrica, triangular, etcétera, antes de realizar el cálculo numérico; verbigracia, si A es rala, entonces se llama a rutinas de Matlab® especializadas en matrices ralas.
c) La factorización LU permite resolver de manera eficiente el sistema de ecuaciones 1.1 cuando se necesita evaluar la solución x para múltiples lados derechos b.
Sin embargo, todos estas funciones trabajan basadas en un proceso de eliminación de Gauss-Jordan, el que sabemos que se pone costoso cuando el número de ecuaciones aumenta. Luego, es necesario considerar esquemas alternativos de solución de 1.1 que no pasen por formar una inversa de la matriz A.
Ejemplo 1.3. Comparación de métodos de solución
Resuelva el sistema de ecuaciones del ejemplo 1.1 usando las funciones inv, ‘\’ o lu de Matlab®.
El siguiente macro aplica las tres metodologías para resolver el sistema de ecuaciones; los tres métodos encuentran la misma solución, y solo el error de redondeo producido en el cálculo de la inversa de la matriz A produce un residuo pequeño de norma 8.9·10-16, siendo los restantes residuos exactamente cero. En el caso de la factorización LU, la rutina entrega una matriz de permutación P, de modo que P·A = L·U; esta matriz P representa los intercambios de filas realizados en la estrategia de pivoteo parcial, para prevenir ceros (o casi ceros) en la diagonal.
% este macro resuelve el ejemplo 1.2
%% ingreso de la matriz y vector lado derecho
A=[2 1 -3; -1 3 2; 3 1 -3];
b=[-1 12 0]’;
%% solución aplicando la inversa de la matriz A
x1=inv(A)*b;
r1=A*x1-b;
%% solución aplicando el operador ‘\’
x2=A\b;
r2=A*x2-b;
%% solución aplicando factorización lu
[L,U,P]=lu(A);
b3=P*b;
y3=L\b3;
x3=U\y3;
r3=A*x3-b;
1.2 Métodos iterativos
Existe una amplia variedad de métodos iterativos aplicados a los sistemas de ecuaciones lineales. A menudo estos métodos se usan en la solución de sistemas lineales con un gran número de variables (n >> 1), los que son costosos de resolver vía eliminación gaussiana. La mayoría de estos métodos son iteraciones estacionarias del tipo (Axelson, [4]):
Comúnmente la matriz M y el vector g son independientes del índice de iteración k. Una condición suficiente para la convergencia de la iteración 1.9 es que la norma de la matriz M sea menor que 1. En particular una condición necesaria y suficiente para la convergencia es que
Donde ρ(M) es el radio espectral de la matriz M, y corresponde al mayor valor absoluto de los valores propios de M. En Matlab®, lo calcularíamos como ρ=max(abs(eig(M))).
La forma de iteración dada por la ecuación 1.9 se obtiene de la ecuación original A·x = b mediante la descomposición de la matriz original A. En particular, podemos descomponer la matriz A de la siguiente forma:
Hay que hacer notar que en Matlab® podemos generar la descomposición de la siguiente forma: D = diag(diag(A)); L = tril(A,-1) y U = triu(A,1). A partir de esta descomposición es posible obtener varios esquemas que se han utilizado con frecuencia como métodos iterativos (Nakamura, [5]; Gerald, [6]), y que se describen a continuación.
1.2.1 Método de Jacobi (Desplazamientos simultáneos)
En este método se aplica la descomposición 1.11 y se deja la diagonal de la matriz original en el lado izquierdo de la ecuación original, de modo que
Por lo tanto, en cada iteración se calculan los nuevos componentes del vector de solución, usando la información de los componentes de la etapa anterior, por lo que en la ecuación 1.12 el componente ‘i’ del vector solución actualizado queda como:
1.2.2 Método de Gauss-Seidel (Desplazamientos sucesivos)
En esta variante se pasa solo la diagonal superior de la matriz A al lado derecho de la ecuación original, obteniéndose:
En este método, la componente i-ésima del vector solución en la etapa (k+1) se actualiza usando la información de los componentes del vector x(k) y los nuevos componentes del vector x(k+1) hasta el componente i-1:
Es decir, se emplea información más actualizada que en el caso del método de Jacobi.
1.2.3 Método
