Resolución de sistemas compatibles determinados con wxMaxima (y III): métodos iterativos

En las entradas anteriores de la serie "Resolución de sistemas compatibles determinados con wxMaxima" expliqué como aplicar en wxMaxima algunos métodos directos y algunos métodos basados en factorización para resolver sistemas de ecuaciones lineales (compatibles determinados).

En esta entrada explico cómo aplicar algunos de los métodos denominados iterativos: el método de Jacobi, el método de Gauss-Seidel y el método de relajación.

Podemos decir, simplificando la cuestión, que la principal diferencia entre los métodos iterativos y los métodos vistos en las entradas anteriores es que los métodos iterativos consisten en empezar con un valor inicial que se le da al vector solución x y se aplica un algoritmo a las sucesivas aproximaciones que se van obteniendo. Es decir, la entrada de una aplicación del algoritmo es la salida de la aplicación anterior del algoritmo. Estos algoritmos pueden acercarse más a la solución del sistema en cada iteración (en ese caso se dice que el algoritmo converge) o alejarse de ella (o entrar en un bucle) por lo que no obtendríamos la solución del sistema.

En el caso de que el algoritmo converja a la solución debemos decidir cuándo consideramos que la solución es lo "suficientemente buena". Ese parámetro lo llamamos tolerancia y, de nuevo simplificando la explicación, podemos decir que es el margen de error que vamos a permitir (o una medida de lo lejos que puede estar la aproximación de la solución).

Recordemos que estamos resolviendo el siguiente sistema:

10x+4y+4z+3t+5s=1
  4x+8y+2z+2t+5s=1
  4x+2y+6z+3t+4s=1
  3x+2y+3z+8t+5s=1
  5x+5y+4z+5t+6s=1

Y que en la primera entrada de la serie comprobamos que se trata de un sistema compatible determinado.

Al aplicar en wxMaxima cualquiera de los 3 métodos iterativos que se explican en esta entrada, se definen una serie de parámetros previos (que son comunes a los 3 métodos):

D:zerofor(A);
for i:1 thru matrix_size(A)[2] do D[i,i]:A[i,i];
D;

x0:zerofor(b);

n:matrix_size(A)[2];
tol: 1e-10;
Kmax: 10000;

El parámetro tol es la tolerancia comentada anteriormente y el parámetro Kmax es el número máximo de iteraciones que vamos a permitir al algoritmo (para evitar que entre en un bucle infinito). En nuestro ejemplo permitimos un máximo de 10000 iteraciones. Vamos entonces con los métodos:

Método de Jacobi:

MJ:D;
NJ:MJ-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MJ).(NJ.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

Al ejecutar la instrucción anterior wxMaxima devuelve un error de desbordamiento, lo que significa que, para el sistema que estamos resolviendo, el método de Jacobi no converge.

Método de Gauss-Seidel:

MGS:copymatrix(A);

for i:1 thru n do
 for j:1 thru n do
   if i

MGS;

NGS:MGS-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MGS).(NGS.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

x;

float(%);

Método de relajación:

En este ejemplo escogemos un parámetro de relajación de 1.5 (sobre-relajación).

w:1.5;

MRA:copymatrix(A);

for i:1 thru n do
 for j:1 thru n do
   if i<=j then MRA[i,j]:0;

MRA;

MR:D/w+MRA;

NR: MR-A;

r: transpose(x0);
x: r;

for k:1 thru Kmax do
   ( x: invert(MR).(NR.x+b), r: A.x-b,
      if (mat_norm(r,1) < tol*mat_norm(transpose(b),1)) then k:Kmax );

x;

Con esta entrada doy por finalizada la serie "Resolución de sistemas compatibles determinados con wxMaxima". Espero que pueda ser de utilidad.