Resolución de ecuaciones diferenciales. Parte(I) Resolviendo sistemas, Gauss Seidel.

Hoy me apetece hacer algo diferente. Voy a escribir una serie de post, unos 3 o 4, en los que voy a tratar con ejemplos los pasos básicos para resolver ecuaciones diferenciales sobre un dominio con condiciones de contorno. Espero que os guste.

Creo que las explicaciones van a seguir este orden:

1. Resolución de sistemas lineales de ecuaciones. Método de Gauss-Seidel.

2. Teoría de discretización de una ecuación diferencial. Esquemas y métodos. Práctica; ecuación del calor unidimensional, planteamiento de la matriz de coeficientes, resolución con el esquema del punto 1 y visualización de resultados. Mismo desarrollo que el anterior.

3. Ecuación del calor bidimensional.

4. Algo más complejo, espero.

En este primer post vamos a escribir un programa para resolver sistemas de ecuaciones lineales de las que nos interesarán más adelante. Digo de las que nos van a interesar porque vamos a resolver un tipo concreto de sistemas en las que no nos vamos a tener que preocupar por problemas de divisiones por cero y otras dificultades que se pueden plantear si lo que queremos es hacer un “solver” más general.

Soy consciente -y espero que vosotros también, sino asimiladlo- de cualquier programa que podamos hacer en el campo del álgebra no será ni la mitad de eficiente que cualquiera de los que hay en librerías especializadas, por lo que si os veis en la situación de resolver un producto de matrices, un escalar por matriz, vector… etc, tiréis de librerías como las que ofrece Numpy o equivalentes en su respectivo lenguaje.

Dicho esto me justifico, a pesar de que sea lento creo que es bueno comprender el esquema de resolución y plantearlo al menos una vez. Hace ver las cosas desde otro punto de vista.

Método de Gauss-Seidel.

El método de Gauss-Seidel forma parte de los métodos indirectos o iterativos de resolución de sistemas de ecuaciones. En estos métodos se comienza con un valor inicial para cada una de las variables (semilla) x^{0}=\left(x_{1}^{0},x_{2}^{0},x_{3}^{0},...,x_{n}^{0}\right) y a partir de ella se construye una nueva aproximación a la solución
x^{1}=\left(x_{1}^{1},x_{2}^{1},x_{3}^{1},...,x_{n}^{1}\right).Siguiendo este esquema se construye una sucesión de vectores \left\{ x^{k}\right\} con el objetivo de que la solución converja, dicho de otro modo:

\underset{x\rightarrow\infty}{lim}x^{k}=x^{*}

Es decir, al cabo de un número finito de iteraciones en ausencia de errores de redondeo obtendremos x^{*} que es la solución del sistema inicial Ax=b.

Cada una de las iteraciones del método cuenta con n sub-iteraciones en las que se modifica el valor de una de las variables de la semilla de la iteración, así, para x_{1} tendremos:

x_{1}^{1}=\frac{b_{1}-(a_{12}x_{2}^{0}+a_{13}x_{3}^{0}+...+a_{1n}x_{n}^{0}}{a_{11}}\,,\, x_{i}^{1}=x_{i}^{0}\,,\, i=1,2,..,n

Se modifican sucesivamente los valores de las variables hasta llegar a x_{n}, punto en el cual tendremos completo nuestro vector de valores para iniciar la siguiente iteración.

En la práctica los errores de redondeo van a ser inevitables, por lo que hemos de imponer una condición de parada para detener el algoritmo, lo cual haremos imponiendo un valor de tolerancia. Cuando la variación de el valor de una variable en la iteración t+1 respecto a su valor en t sea menor que la tolerancia tendremos satisfecho nuestro criterio de parada. Tened en cuenta que el criterio de parada se puede aplicar a todas las variables del sistema o sólo a unas determinadas […].

Implementando el algoritmo.

Una vez tenemos las cosas algo más claras vamos con algo de código. Lo he escrito en C++ sin muchas florituras, vamos a partir de este código y vamos a ver qué podemos hacer con él. No he querido entrar en clases y lo he hecho lo más simple posible para que todo el mundo lo siga.

Fichero “Funciones.h”

/*
 * Funciones.h
 *
 *  Created on: 05/11/2012
 *      Author: samuelrodriguezbernabeu
 *
 *      Contiene las funciones basicas para la primera version del programa GaussSeidel
 */

#ifndef FUNCIONES_H_
#define FUNCIONES_H_

#include
#include <math.h>
#include <time.h>
#include
using namespace std;

#define rowSize 3 //Longitud de las filas y las columas (matriz cuadrada m=n)

int ind(int i, int j)
{
	//La funcion nos da el indice traducido a la notacion de vector lineal a partir de la notacion i,j
	int posicion = i*rowSize+j;
	return posicion;
};

int diagPosition(int i)
{
	//Nos da la posicion del elemento diagonal dentro de la fila que le pasemos como argumento (la fila es i)
	int posicion = i*rowSize + i;
	return posicion;
};

double diffclock(clock_t clock1,clock_t clock2)
{
	double diffticks=clock1-clock2;
	double diffms=(diffticks*1000)/CLOCKS_PER_SEC;
	return diffms;
};

#endif /* FUNCIONES_H_ */

Ahora el fichero GaussSeidelV1.cpp donde desarrollamos la versión básica del algoritmo:

//============================================================================
// Name        : GaussSeidel.cpp
// Author      : Samuel Rodriguez Bernabeu
// Version     : 1.0
// Description : Gauss-Seidel solver V1, Ansi-style
//============================================================================

#include "Funciones.h"

int main() {

	//Parametros necesarios para el sistema
	double coefMatrix[rowSize*rowSize] = {3, -1, 1, 1, -5, 1, 1, -1, 4};
	double xVector[rowSize] = {1,8,11};
	double C0[rowSize] = {0}; 			//vector para la semilla
	double C1[rowSize] = {0}; 			//vector para la copia de la semilla
	unsigned int iteraciones = 20; 		//numero de iteraciones

	double errorVector[rowSize]; //este vector va a contener los datos de error

	cout << "\tResultados de las iteraciones:\n";

	//Variables auxiliares necesarias para el programa

	double suma = 0;		//almacena el sumatorio de los productos

        for (unsigned int iter = 0; iter < iteraciones; ++iter) {
		for (unsigned int fila = 0; fila < rowSize; ++fila) {

			suma = 0;

			for (unsigned int j = 0; j < rowSize; ++j) {
				double elemento = coefMatrix[ind(fila,j)];
				if(fila != j)
				{

					suma = 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[fila] = (xVector[fila]-suma)/coefMatrix[diagPosition(fila)];

		}
		//Hacemos el calculo del error en cada iteracion. Despues implementaremos una condicion de parada
		for (unsigned int var = 0; var < rowSize; ++var) {

			errorVector[var] = fabs(C0[var]-C1[var]);
                        //C1 es el anterior, C0 el actual
			C1[var] = C0[var];
                        //Hacemos una copia para la siguiente iteracion
		}
		printf("Iter: %d x:%.6f y:%.6f z:%.6f\n", iter,C0[0],C0[1],C0[2]);
	}
        return 0;
}

Los coeficientes de la matriz y de b los he copiado de la entrada de la Wiki para no pensar mucho =). Podéis ver el resultado de las iteraciones:

Iter: 0 x:0.333333 y:-1.533333 z:2.283333
Iter: 1 x:-0.938889 y:-1.331111 z:2.651944
                     ...
Iter: 9 x:-0.980000 y:-1.260000 z:2.680000

Sencillo ¿verdad?. Un paso interesante ahora sería plantear el programa de forma que crease matrices y vectores de una dimensión dada para que podamos probar nuestro solver en unas condiciones más exigentes. De nuevo lo hacemos de forma simple (lo suyo es usar templates y así tener la biblioteca para otras condiciones, pero por ahora está bien así).

Vamos entonces a crear matrices pero, adelantándonos a los tiempos vamos a crear matrices bandeadas. Cuando planteemos la discretización de las ecuaciones diferenciales veremos que son éste tipo de matrices las que nos van a aparecer.

Por ahora no tengo más tiempo, seguiré ampliando esta entrada esta semana, aún hay muchas cosas por hacer, ¿sugerencias? =)

Esta entrada fue publicada en básico, C/C++, matemáticas. 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