Resolución de ED. Parte(II) Gauss Seidel, cosas a tener en cuenta y matrices dinámicas.

Seguimos con lo que dejamos a medio en la entrada sobre resolución de sistemas de ecuaciones. Y vamos a ir adelantando algunas cosas más; como ya mencionamos nos interesan las matrices diagonales o bandeadas. Aún no vamos a entrar en ello pero diremos que hay un teorema que nos va a ayudar a encontrar la solución para matrices de este tipo:

Teorema; Si A (matriz de coeficientes) es estrictamente dominante por filas entonces el método de Gauss-Seidel converge para cualquier valor inicial de x.

Interesante, vamos a ver qué es eso. Una matriz cuadrada es de diagonal estrictamente dominante por filas si en cada fila el valor absoluto del elemento diagonal es mayor que la suma de los valores absolutos de los otros elementos de la fila.

Implementamos rápidamente esto en el código para tener información sobre este factor que, si hacemos algunas pruebas con números aleatorios para la matriz de coeficientes veremos que tiene efectivamente una gran trascendencia:

//loop para comprobar si es estrictamente dominante:
double sumaElementos = 0; 		//contador para la suma de elementos de la fila
unsigned int diagDomCounter = 0;	//contador para filas dominantes

for (int i = 0; i < rowSize; ++i) {
	sumaElementos = 0;
	for (int j = 0; j  sumaElementos)
	{
		++diagDomCounter;
	}
}
cout << "Tenemos " << diagDomCounter << " filas dominantes de " << rowSize << ".";
cout << "Porcentaje: " << (diagDomCounter*100)/rowSize << "%" << endl;

Seguramente si obtenemos un porcentaje como el siguiente tengamos bastantes papeletas para que el método llegue a una solución:

Tenemos 3 filas dominantes de 3. Porcentaje: 100%
7	-3	0
-1	6	0
-2	-1	9

	Resultados de las iteraciones:

Iter: 1 x:0.142857 y:1.357143 z:1.404762
Iter: 2 x:0.724490 y:1.454082 z:1.544785
Iter: 3 x:0.766035 y:1.461006 z:1.554786
Iter: 4 x:0.769002 y:1.461500 z:1.555501
Iter: 5 x:0.769214 y:1.461536 z:1.555552
Iter: 6 x:0.769230 y:1.461538 z:1.555555
Iter: 7 x:0.769231 y:1.461538 z:1.555556
Iter: 8 x:0.769231 y:1.461538 z:1.555556
Iter: 9 x:0.769231 y:1.461538 z:1.555556
Iter: 10 x:0.769231 y:1.461538 z:1.555556

El solver ha empleado 0.059 ms para N= 9 elementos.

Ahora que somos conscientes de la importancia de que nuestra matriz sea estrictamente dominante por filas vamos a pasar a jugar con matrices algo más grandes (si lo hubiese implementado antes posiblemente podríais haber obtenido errores de divergencia). Bien, como decía vamos a generar matrices de filas/columnas NxN «aleatoriamente» para asegurarnos un buen porcentaje de filas dominantes, así que vamos a alterar un poco el azar computacional.

double coefMatrix[nElementos] = {0.0}; 	//inicializada a cero
double BVector[rowSize] = {0};		   	//vector de soluciones
srand(time(0));							//Tiempo de referencia

for (int i = 0; i < rowSize; ++i) {
	for (int j = 0; j < rowSize; ++j) {
		if(i!=j){
		coefMatrix[ind(i,j)] = rand()%6-3 ;
		}
		else{
			coefMatrix[ind(i,j)] = rand()%20+1;
			//nos aseguramos de que los coeficientes de la diagonal
			//tengan mas "suerte" con el azar.
		}
	}
}

//generamos ahora los numeros aleatorios para el vector B
for (int i = 0; i < rowSize; ++i) {
	BVector[i] = rand()%10-5;
}

Haciéndolo de esta manera tendremos más posibilidades de que el sistema esté bien condicionado. No es hacer trampas, si vamos a jugar con algo es mejor que funcione!!

Atención; puesto que el valor absoluto del elemento diagonal ha de ser mayor que la suma de los valores absolutos de los elementos de la fila, cuanto mayor sea el número de filas más grande ha de ser el elemento diagonal […] tenedlo en cuenta al aumentar indiscriminadamente el número de elementos para tostar vuestro procesador =).

Tenemos 4 filas dominantes de 5. Porcentaje: 80%
17	-2	1	-1	0
-2	5	0	-3	-2
-2	-2	19	2	2
0	1	-3	20	-1
-1	-3	0	0	19

	Resultados de las iteraciones:

Iter: 1 x:-0.117647 y:0.752941 z:0.277399
Iter: 2 x:-0.048091 y:0.714008 z:0.295767
Iter: 3 x:-0.053763 y:0.710631 z:0.295096
Iter: 4 x:-0.054125 y:0.710120 z:0.295098
Iter: 5 x:-0.054186 y:0.710046 z:0.295096
Iter: 6 x:-0.054194 y:0.710036 z:0.295096
Iter: 7 x:-0.054196 y:0.710034 z:0.295095
Iter: 8 x:-0.054196 y:0.710034 z:0.295095
Iter: 9 x:-0.054196 y:0.710034 z:0.295095
Iter: 10 x:-0.054196 y:0.710034 z:0.295095
Iter: 11 x:-0.054196 y:0.710034 z:0.295095

El solver ha empleado 0.089 ms para N= 25 elementos.

Una vez hecho esto vamos a emplear matrices más grandes (n>100 elementos) y tengo curiosidad por saber si emplear clases o punteros varía los tiempos del programa. Por ahora dejo aquí una gráfica. Luego la subo en Gnuplot, porque esta es bastante cutre…