Resolución de sistemas de ecuaciones lineales de grado "n"

[fzp]

Dejo aquí un método -que creo que funciona- para la resolución de un sistema de "n" ecuaciones lineales con "n" incógnitas.
Está basado en el algoritmo de Gauss (no sé si otros lo llaman de Gauss-Jordan u otros de Gauss-Jacobi). Sea como sea que lo llamen consiste en dejar la matriz de coeficientes en una matriz triangular con 0 (ceros) por debajo de la diagonal principal de a matriz de coeficientes. Y luego ir resolviendo despejando in cógnitas de atrás hacia delante.

Sé que es un método poco eficiente. Para empezar los errores de cálculos sucesivos (especialmente en sistema muy grandes con muchas ecuaciones) pueden producir desbordamientos (overflow) o que, si hay términos relativamente pequeños respecto a otros, se produzcan reducciones a "1" que desvirtúen el resultado final, para sistemas de ecuaciones muy particulares. Como mi intención es simplemente una primera aproximación para reslover sistemas del tipo de los que se dan en Ingeniería (cálculo de estructuras, mallas de tuberías, etc...) de momento éso no me importa demasiado.

También la forma de introducir los datos es totalmente ineficiente -a mano-. Es sólo para pruebas; lo suyo es que sea a partir de datos de matrices que previamente se hayan almacenado en disco. Es más, esas matrices no se introducirán a mano -ni siquiera en disco- sino que provendrán de otros programas que las habrán guardado en disco a partir de otros datos (matrices de rigidez de una estructura -de construcción- a partir de sus datos de longitud de barras, inercia, material, etc; o longitud de tubería, diámetro, material, persión, etc...).

Por lo tanto, el hecho de la forma de introducción de datos no debe de ser valorado; solo el funcionamiento del algoritmo de resolución del sistema.

Por lo tanto, respecto al programa en general, me gustaría tener opiniones sobre:
- Sobre todo: ¡fallos que pueda tener! - He probado varios sistemas, algunos con bastantes incógnitas, y los ha resuelto. He probado varios sistemas con fallos "evidentes" (filas = 0; columnas = 0; filas/columnas = combinación lineal de otras filas/columnas;...) y me los ha detectado. En ese sentido, mi miedo es más bien -no que no me detecte sistemas sin solución-, sino que sistemas que SÍ tienen solución, me los detecte como que no la tienen. ¿Veis algún fallo? ¿Alguna forma de mejorarlo?

- ¿Se podría mejorar, mediante un algoritmo iterativo, la solución? He visto cosas, pero ninguna lo suficientemente clara -para mi- sobre cómo mejorar una solución a partir de una primera solución. En concreto he leído sobre los métodos de Richardson, de Jacobi, de Gauss-Seidel, que incluso empiezan con vectores/matrices de soluciones "alejadas" de la solución verdadera; cuando se podría mplementar desde una primera solución "muy aproximada" - la que ofrece el algoritmo que aquí se da-; pero no sé como implementarla.

- Por último: sé que hay un método muy potente (y con reducción de operaciones y de tiempo de cálculo y de disminución por errores de operaciones (¡justo lo que quiero disminuir!) que es la "condensación de ecuaciones"..., pero no encuentro información viable. Si alguien me puede ofrecer literatura sobre ese tema, pues estaría genial.

Bueno, aquí dejo el código que he hecho, para las críticas/comentarios/ayudas/mejoras/etc... que creáis convenientes (cualquier cosa será bien recibida):

Código (c) [Seleccionar] Expandir

#include <stdio.h>
#include <stdlib.h>

void intro_sistema (int n, double ** a);
void combina_fila (int j, int n, double ** a); // Para dejar coeficientes de la columna a cero
void intercambia_fila (int i, int j, int n, double ** a); // Para que no haya coefic. == 0 en la diagonal
void calcula_incognitas (int n, double ** a, double * solucion); // A partir de una matriz ya triangular

int main ()
{
int n; // numero ecuaciones lineales
int i, j;
double ** a = NULL; // matriz ampliada: coeficientes | terminos independientes
double * solucion = NULL;

printf ("Introducir no. de ecuaciones: ");
scanf ("%d", &n);

a = (double **) malloc ( n * sizeof (double *) ); // filas
for (j = 0; j < n; ++j)
a[j] = (double *) malloc ( (n+1) * sizeof (double) ); // columnas
solucion = (double *) malloc ( n * sizeof (double) );

intro_sistema (n, a);

for (j = 0; j < n-1; j++)  // Procesa elementos diagonal matriz coeficientes
{
if (a[j][j] == 0)
{
for (i = j+1; i < n; i++) // Busca un coeficiente de la misma columna != 0
{
if (a[i][j] != 0)
{
intercambia_fila (i, j, n, a);
break;
}
printf ("\n\nEl sistema no tiene solucion unica"); // No hay coficiente != 0 sobre el que pivotar
return 0;
}
}
combina_fila (j, n, a);
}

if (a[n-1][n-1] == 0)
{
printf ("\n\nEl sistema no tiene solucion unica\n"); // Toda la fila inferior es de ceros
return 0;
}
else
{
calcula_incognitas (n, a, solucion);
printf ("\n\n");
for (j = 0; j < n; j++)
printf ("soluc[%d] = %lf\n", j, solucion[j]);
}
return 0;
}

void intro_sistema (int n, double ** a)
{
int i, j;
double var;

printf ("\nIntroducir matriz de coeficientes del sistema de ecuaciones lineales\n");
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
{
printf ("a[%d][%d] = ", i, j);
scanf ("%lf", &var);
a[i][j] = var;
}
printf ("\nIntroducir terminos independientes del sistema\n");
for (i = 0; i < n; i++)
{
printf ("b[%d] = ", i);
scanf ("%lf", &var);
a[i][n] = var;
}
}

void combina_fila (int j, int n, double ** a)
{
double factor;
int i, k;

for (i = j+1; i < n; i++)
{
if (a[i][j] != 0)
{
factor = -a[j][j] / a[i][j];
for (k = j; k < n+1; k++)
a[i][k] = (factor * a[i][k]) + a[j][k];
}
}
}

void intercambia_fila (int i, int j, int n, double ** a)
{
double var_intercamb;
int k;

for (k = j; k < n+1; k++)
{
var_intercamb = a[i][k];
a[i][k] = a[j][k];
a[j][k] = var_intercamb;
}
}

void calcula_incognitas (int n, double ** a, double * solucion)
{
double acumulador;
int i, j;

solucion[n-1] = a[n-1][n] / a[n-1][n-1];
for (i = n-2; i >= 0; i--)
{
acumulador = a[i][n];
for (j = i+1; j < n; j++)
acumulador = acumulador - (a[i][j] * solucion[j]);
solucion[i] = acumulador / a[i][i];
}
}


Pues nada, lo dicho, que se agradece cualquier comentario.