Resolución de ED. Parte(III) Gauss Seidel, sistema y solver mediante clases y representación gráfica.

Creo que ya tenemos lo necesario para poder seguir en la siguiente entrada pasando a la discretización de las ecuaciones y al planteamiento del sistema.

Esta es la tercera de las entradas sobre el método de Gauss-Seidel pero no es la última, volveré sobre ella para las optimizaciones. Y más adelante también para tratar las matrices dispersas y los algoritmos específicos para resolver matrices con 3 o 5 bandas. Hubiese querido hacerlas seguidas pero el compilador da un error de arquitectura en el Mac -lo haremos pues en los Xeon..😛 – las haré esta semana.

Para ahorrar código en el programa principal y poder centrarnos en entender los desarrollos matemáticos he planteado el sistema de ecuaciones en forma de clase -con templates para después reciclarlo en GPU en simple precisión si llegamos a ello.

Se ha separado dentro de la clase la matriz del coeficiente del vector B y del vector de soluciones, he preferido hacerlo así en lugar de hacerlo mediante clases heredadas para que el código sea más claro y más ligero aunque a mí me guste más mediante herencia para sobrecargar los operadores más fácilmente y tener bastante más juego [seguro que también lo vemos].

Al final os he puesto un ejemplo donde se construye una matriz tribanda que nos va a ser bastante útil en la siguiente entrada.

Definición de la clase:

/*
 * LinearSystems.h
 *
 *  Created on: 06/11/2012
 *      Author: Samuel Rodriguez Bernabeu
 */

#ifndef LINEARSYSTEMS_H_
#define LINEARSYSTEMS_H_

template
class LinSystem {
public:
	LinSystem(int rowSize);
	LinSystem(const char* TypeofMatrix, int rowSize);
	~LinSystem() {};

	//------Funciones de la clase----//
	void showAMatrix();
	void showBVector();
	void showSolution();
	void saveSolution();
	void plotSolution();
	void inputA(Type valor, int posicion);
	void inputB(Type valor, int posicion);
	int getSize();
	int dominancia();
	void InfoDominancia();
	void solveSeidel(double Tolerance, double x0);

protected:
	const char* MatrixType;
	int RS;//rowSize
	int N;//numero de elementos de la matriz de coeff

	//------punteros-------//
	Type *coefMatrix;
	Type *bVector;
	Type *SVector;//vector de soluciones
};

//--------------DEFINICION DE LAS FUNCIONES DE LA CLASE---------//
template
LinSystem::LinSystem(int rowSize)
{
	RS = rowSize;
	N = RS*RS;
	MatrixType = "UserDefined";

	//--------Creamos los arrays-------------//
	coefMatrix = new Type[N];
	bVector = new Type[RS];

	//--------Inicializamos el vector B------//
	srand(time(0));

	for (int i = 0; i < RS; ++i) {
		bVector[i] = 0.0;
	}

	//-----Inicializamos la matriz de coef---//
	for (int i = 0; i < N; ++i)
	{
			coefMatrix[i] = 0.0;
	}
}

template
LinSystem::LinSystem(const char* TypeofMatrix, int rowSize)
{
	RS = rowSize;
	N = RS*RS;
	MatrixType = TypeofMatrix;

	//--------Creamos los arrays-------------//
	coefMatrix = new Type[N];
	bVector = new Type[RS];

	//--------Inicializamos el vector B------//
	srand(time(0));

	for (int i = 0; i < RS; ++i) {
		bVector[i] = rand()%10-5;
	}

	//-----Inicializamos la matriz de coef---//

	if (strcmp (MatrixType, "Zeros") == 0){
		for (int i = 0; i < N; ++i) {
			coefMatrix[i] = 0.0;
		}
	}//Este es un ejemplo, se pueden implementar muchas...

	else
	{
		for (int i = 0; i < RS; ++i) {
			for (int j = 0; j < RS; ++j) {
				if(i!=j){
				coefMatrix[i*RS+j] = rand()%2-1 ;
				}
				else{
					coefMatrix[i*RS+j] = rand()%80-40;
					//nos aseguramos de que los coeficientes de la diagonal
					//tengan mas "suerte" con el azar.
				}
			}
		}

		//dentro de esto quiero comprobar las filas dominantes, eso lo hare despues
	}
}

template
void LinSystem::showAMatrix()
{
	cout << "\n";
	for (int i = 0; i < RS; ++i) {
		for (int j = 0; j < RS; ++j) {
			cout << coefMatrix[i*RS + j] << "\t";
		}
		cout << "\n";
	}
	cout << "\n";
}

template
void LinSystem::showBVector()
{
	cout << "\n";
	for (int i = 0; i < RS; ++i) {
		cout << bVector[i] << "\t";
	}
	cout << "\n";
}

template
void LinSystem::showSolution()
{
	cout << "\n";
	for (int i = 0; i < RS; ++i) {
		cout << SVector[i] << "\t";
	}
	cout << "\n";
}

template
void LinSystem::plotSolution()
{
	system("gnuplot -persist ./plotData.gp");
	//de esta manera se ejecuta el script permitiendo que continue el programa
}

template
void LinSystem::saveSolution()
{
	FILE *out;
	out = fopen("solveData.dat", "w");
	fprintf(out,"# solucion al sistema de ecuaciones\n");

	for (int i = 0; i < RS; ++i) {
		fprintf(out,"%d\t%f\n", i+1,SVector[i]);
	}

	fclose(out);

	cout << "... ... ... Guardado de datos completo" << endl;
}

template
void LinSystem::inputA(Type valor, int posicion)
{
	//modificamos el tipo de la matriz
	MatrixType = "UserDefined";

	coefMatrix[posicion] = valor;
}

template
void LinSystem::inputB(Type valor, int posicion)
{
	bVector[posicion] = valor;
}

template
int LinSystem::getSize()
{
	return RS;
};

template
int LinSystem::dominancia()
{
	double sumaElementos;
	unsigned int counter = 0;
	for (int i = 0; i < RS; ++i) {
		sumaElementos = 0;
		for (int j = 0; j < RS; ++j) { 			if(i!=j) 			{ 				sumaElementos += fabs( coefMatrix[i*RS+j] ); 			} 		} 		if (fabs(coefMatrix[i*RS+i]) > sumaElementos)
		{
			counter = counter+1;
		}
	}
	//cout << "Filas dominantes" << counter << endl;
	return counter;
}

template
void LinSystem::InfoDominancia()
{
	cout << "\nINFO; Filas diag. dominante:  "<< dominancia() << " entre " << RS << " filas totales. Porcentaje: " <<
			(dominancia()*100)/RS << endl;
}

template
void LinSystem::solveSeidel(double Tolerance, double x0 = 0)
{
	double suma = 0;		//almacena el sumatorio de los productos
	double elemento = 0;

	//------Declaramos vectores para semillas-----//
	Type *C0 = new Type[RS]; 				//vector para la semilla
	Type *C1 = new Type[RS]; 				//vector para la copia de la semilla
	SVector = new Type [RS];

	double error; //variable para calcular el error
	int convergenceCounter; //variable auxiliar
	unsigned iter = 0;

	//iniciamos con el valor de la semilla si lo hay
	if(x0 != 0)
	{
		for (int i = 0; i < RS; ++i) {
			C0[i] = x0;
		}
	}

	while (true)
	{
		for (int i = 0; i < RS; ++i) {

			suma = 0;

			for (int j = 0; j < RS; ++j) {
				elemento = coefMatrix[i*RS+j];
				if(i != j)
				{

					suma += elemento*C0[j];
				/*
				 * Como las matrices que vamos a tratar siempre van a tener una diagonal no nula
				 * no vamos a tener problemas con esto...
				 */
				}

			}
			//Actualizamos el valor de la semilla
			C0[i] = (bVector[i]-suma)/coefMatrix[i*RS+i];

		}
		//Hacemos el calculo del error en cada iteracion.
		//Vamos con la condicion de parada
		convergenceCounter = 0; //variable auxiliar
		for (int var = 0; var < RS; ++var) {

			error = fabs(C0[var]-C1[var]); 	//C1 es el anterior, C0 el actual
			C1[var] = C0[var]; 							//Hacemos una copia para la siguiente iteracion

			if(error <= Tolerance){
				++convergenceCounter;
			}
		}
		printf("Iter: %d x:%.6f y:%.6f z:%.6f\n", iter+1,C0[0],C0[1],C0[2]);

		//Parada por tolerancia y salida del bucle
		if (convergenceCounter == RS)
		{
			cout << "INFO; Se ha alcanzado el criterio de parada tras ."<< iter << " iteraciones." << endl;
			break;
		}

		//-----Contador de iteraciones y criterio de parada por iteraciones maximas
		++iter;
		if (iter == MAX_ITER) {
			cout << "INFO; Se ha alcanzado el numero maximo de iteraciones: " << MAX_ITER << endl;
		}
	}
	//copiamos las soluciones a los vectores de la clase
	for (int i = 0; i < RS; ++i) {
		SVector[i] = C0[i];
	}
}

#endif /* LINEARSYSTEMS_H_ */

Fichero principal:

//============================================================================
// Name        : GaussSeidel_linux.cpp
// Author      : Samuel Rodriguez Bernabeu
// Version     : 4
// Copyright   : GNU
// Description : Hello World in C++, Ansi-style
//============================================================================

#include "Includes.h"

int main() {
	cout << "Gauss Seidel V2" << endl;

	LinSystem  LS("Diag",20);
	if (LS.getSize() <= 10) {LS.showAMatrix();}
	LS.InfoDominancia();
	LS.solveSeidel(10e-6);

	LinSystem  UserMatrix(6);
	if (UserMatrix.getSize() <= 10) {UserMatrix.showAMatrix();}

	for (int i = 0; i < 6; ++i) {
		for (int j = 0; j < 6; ++j) {
			UserMatrix.inputA(2, 6*i+i+1); 		//banda superior
			UserMatrix.inputA(7, 6*i+i);			//diagonal
			UserMatrix.inputA(3, 6*(1+i)+i); 	//banda inferior
		}
	}

	if (UserMatrix.getSize() <= 10) {UserMatrix.showAMatrix();}

	return 0;
}

Fichero con bibliotecas

/*
 * Includes.h
 *
 *  Created on: 06/11/2012
 *      Author: samuelrodriguezbernabeu
 */

#ifndef INCLUDES_H_
#define INCLUDES_H_

#define MAX_ITER 500

#include
#include <time.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include 	//depuracion
#include 	//comparar cadenas
#include 	//fabs
using namespace std;

//------------Ficheros de declaracion de clases-------//
#include "LinearSystems.h"

#endif /* INCLUDES_H_ */

Notas adicionales:

Dentro de la clase hay implementado dos métodos para mostrar la gráfica de las soluciones. El primero guarda los datos en un fichero de texto plano y el segundo las representa mediante una llamada a Gnuplot vía shell. No he probado estas cosas en Windows nunca, así que no sé si funcionarán, Linux es más flexible en este sentido. Os agradecería información sobre el tema si alguien lo ejecuta en un entorno Windows, así aprendemos todos. =)

Hasta que no resolvamos sistemas coherentes no tiene mucho sentido representar las soluciones de sistemas aleatorios, pero podéis jugar, no obstante nos adelanto…

Os dejo el script para Gnuplot tambien. plotData.gp

set grid;
set ylabel 'Valor numerico';
set xlabel 'Nodo';
set title 'Solucion del sistema de ecuaciones';
plot 'solveData.dat' w lines t 'x'

También veremos cuando resolvamos la ecuación unidimensional algo así:

Distribución de temperaturas en los nodos. Ec. del calor unidimensional

Esta entrada fue publicada en básico, C/C++, Gnuplot, matemáticas y etiquetada , , , , . Guarda el enlace permanente.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s