Source: https://porlasbarbasdekutta.wordpress.com/2012/11/05/resolucion-de-ecuaciones-diferenciales-partei-resolviendo-sistemas-gauss-seidel/
Timestamp: 2017-08-19 03:33:18+00:00

Document:
Resolución de ecuaciones diferenciales. Parte(I) Resolviendo sistemas, Gauss Seidel. | Por las barbas de Kutta
← Matrices de coeficientes de los sistemas de OpenFOAM. Salida y representación visual con Python.
Resolución de ED. Parte(II) Gauss Seidel, cosas a tener en cuenta y matrices dinámicas. →
Resolución de ecuaciones diferenciales. Parte(I) Resolviendo sistemas, Gauss Seidel.
Publicado el 5 de noviembre de 2012	por Samuel
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.
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) y a partir de ella se construye una nueva aproximación a la solución
.Siguiendo este esquema se construye una sucesión de vectores con el objetivo de que la solución converja, dicho de otro modo:
Es decir, al cabo de un número finito de iteraciones en ausencia de errores de redondeo obtendremos que es la solución del sistema inicial .
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 tendremos:
Se modifican sucesivamente los valores de las variables hasta llegar a , 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
*      Contiene las funciones basicas para la primera version del programa GaussSeidel
#ifndef FUNCIONES_H_
#define FUNCIONES_H_
#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;
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;
#endif /* FUNCIONES_H_ */
Ahora el fichero GaussSeidelV1.cpp donde desarrollamos la versión básica del algoritmo:
// Name        : GaussSeidel.cpp
// Description : Gauss-Seidel solver V1, Ansi-style
//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
for (unsigned int iter = 0; iter < iteraciones; ++iter) {
for (unsigned int fila = 0; fila < rowSize; ++fila) {
for (unsigned int j = 0; j < rowSize; ++j) {
double elemento = coefMatrix[ind(fila,j)];
if(fila != j)
suma = suma + elemento*C0[j];
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]);
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.

References: Resolución 

Resolución 
 Resolución 
 resolución 
 resolución 
 resolución