Eliminación de Gauss-Jordan en C#

Eliminación de Gauss-Jordan en C#

Código de la técnica Eliminación de Gauss-Jordan para resolver sistemas de ecuaciones lineales de reales y complejos

Por Harvey Triana, 22 de noviembre de 2006

Introducción

La “Eliminación de Gauss-Jordan”, es una de mis soluciones preferidas, el código lo escribí inicialmente en FORTRAN en 1992, Visual Basic clásico en 1996, y ahora en C# en 2006. Si bien la solución aplicada a números reales es aceptable en Visual Basic, hablo de Visual Basic clásico, la solución para números complejos produce un código Visual Basic necesariamente basado en métodos aritméticos con nombre, lo que limita las expresiones a una forma difícil de construir, y difícil de depurar en términos intuitivamente lógicos. El tratamiento de números complejos con C, o C#, es la aplicación clásica de sobrecarga de operadores, en donde al final, tenemos un código natural y tenazmente equivalente a la solución con números reales.

El articulo de Fahad Gilani (ver referencia 4), es una interesante exposición acerca del rendimiento de C# sobre la cuestión, y una inspiración para mi.

Eliminación de Gauss-Jordan Aplicada a Números Reales

En la matemática, la eliminación Gaussiana o eliminación de Gauss-Jordan, llamada así debido a Carl Friedrich Gauss y Wilhelm Jordan, es un algoritmo del álgebra lineal para determinar las soluciones de un sistema de ecuaciones lineales, y encontrar matrices e inversas. El algoritmo lo puede consultar de las referencias 1 y 2. Aquí les presento la loable versión en C#.

Creé un programa de consola, Article_GJE, y agregué una clase de nombre LinearEquationsSolver.cs. El código de la solución por eliminación de Gauss-Jordan para números reales es como sigue:

Eliminación de Gauss-Jordan para números reales

using System;

namespace Article_GJE

{

class LinearEquationsSolver

{

/// <summary>

/// GaussianElimination()

/// Gaussian elimination is a method for solving matrix equations

/// By Harvey Triana

/// </summary>

/// <param name="a"> The matrix</param>

/// <param name="r"> The solution array</param>

/// <returns>Success function</returns>

public static bool GaussianElimination(double[,] a, double[] r)

{

double t, s;

int i, l, j, k, m, n;

try

{

n = r.Length - 1;

m = n + 1;

for (l = 0; l <= n - 1; l++)

{

j = l;

for (k = l + 1; k <= n; k++)

{

if (!(Math.Abs(a[j, l]) >= Math.Abs(a[k, l]))) j = k;

}

if (!(j == l))

{

for (i = 0; i <= m; i++)

{

t = a[l, i];

a[l, i] = a[j, i];

a[j, i] = t;

}

for (j = l + 1; j <= n; j++)

{

t = (a[j, l] / a[l, l]);

for (i = 0; i <= m; i++) a[j, i] -= t * a[l, i];

}

r[n] = a[n, m] / a[n, n];

for (i = 0; i <= n - 1; i++)

{

j = n - i - 1;

s = 0;

for (l = 0; l <= i; l++)

{

k = j + l + 1;

s += a[j, k] * r[k];

}

r[j] = ((a[j, m] - s) / a[j, j]);

}

return true;

}

catch

{

return false;

}

La función GaussianElimination toma como parámetros la matriz, que es el sistema de ecuaciones lineales, y un arreglo, en el cual se almacenará la posible solución. La función retorna el bolean “true” si el sistema de ecuaciones tiene solución y fue calculado. La función GaussianElimination es dinámica en cuanto al tamaño de la matriz, lo que la hace muy potente. La dimensión del arreglo que contiene la solución debe ser igual al número de variables del sistema de ecuaciones.

Para probar la función, dispuse un ejemplo, en la clase Program (predeterminada de la aplicación de consola), como sigue:

Un ejemplo de uso de la funcion GaussianElimination para numeros reales

using System;

namespace Article_GJE

{

class Program

{

static void Main(string[] args)

{

Sample1();

}

static void Sample1()

{

/* Solucionar el siguiente sistema de ecuaciones lineales

x + y + z = 6 | 1 1 1 6 |

x + z = 4 | 1 0 1 4 |

x + y = 3 | 1 1 0 3 | */

double[,] a = { { 1, 1, 1, 6 },

{ 1, 0, 1, 4 },

{ 1, 1, 0, 3 } };

double [] r = new double[3];

ShowMatrix(a, "Ejemplo 1");

if (LinearEquationsSolver.GaussianElimination(a, r))

ShowSolution(r);

else

Console.WriteLine("No es un sistema de ecuaciones lineales");

}

#region formated output

static void ShowMatrix(double[,] a, string Title)

{

Console.WriteLine(Title + '\n');

for (int i = 0; i <= a.GetUpperBound(0); i++)

{

Console.Write('|');

for (int j = 0; j <= a.GetUpperBound(1); j++)

{

Console.Write(ToStringSign(a[i, j]));

}

Console.Write(" | \n");

}

Console.WriteLine('\n');

}

static void ShowSolution(double[] r)

{

Console.WriteLine("Solución por Eliminación Gaussiana");

for (int i = 0; i <= r.GetUpperBound(0); i++)

{

Console.WriteLine(ToStringSign(r[i]));

}

Console.WriteLine("\n");

}

static private string ToStringSign(double v)

{

if (v < 0) return ' ' + v.ToString(); else return " " + v.ToString();

}

#endregion

}

Note que la asignación de la matriz con sintaxis de C# parece muy natural.

La salida del programa es la siguiente, x = 1, y = 2, z = 3:

Números Complejos

La solución de un sistema de ecuaciones lineales de números complejos es compleja, valga la redundancia. El trabajo algebraico es bastante laborioso por el gran número de operaciones matemáticas que requiere. El sistema de ecuaciones más pequeño, que es de dos por dos, con calculadora en mano, puede tomarnos mucho tiempo, y estar sujeto a errores humanos, imagínese uno grande.

La programación de la matemática de números complejos hace natural con la sobrecarga de operadores, lo cual es un mecanismo clásico del lenguaje C. Como alternativa a usar sobrecarga de operadores podemos usar métodos para cada operación, pero las líneas de código de las operaciones resultan muy confusas.

Bien, lo primero que se debe hacer es escribir una clase ComplexNumber. Hay muchas versiones en la red de tal clase, y todas a puntan a lo mismo, sobrecarga de operadores. La clase Articles_CGE.ComplexNumber es mi versión. Noten que use métodos de propiedad en vez de variables públicas. A mi parecer da un código más limpio.

La clase ComplexNumber

* ComplexNumber Class

* By Harvey Triana

namespace Article_GJE

{

using System;

public class ComplexNumber

{

// members

private double m_Real;

private double m_Imaginary;

// constructor

public ComplexNumber()

{

m_Real = 0;

m_Imaginary = 0;

}

public ComplexNumber(double Real, double Imaginary)

{

m_Real = Real;

m_Imaginary = Imaginary;

}

// properties

public double Real

{

get { return m_Real; }

set { m_Real = value; }

}

public double Imaginary

{

get { return m_Imaginary; }

set { m_Imaginary = value; }

}

// Equal method

public bool Equals(ComplexNumber a)

{

return (m_Real == a.Real && m_Imaginary == a.Imaginary);

}

// Let method

public void Let(ComplexNumber a)

{

m_Real = a.Real;

m_Imaginary = a.Imaginary;

}

// Absolute value of a complex number

public double Abs()

{

return (Math.Sqrt(m_Real * m_Real + m_Imaginary * m_Imaginary));

}

// Argument of a complex number in radians

public double Arg()

{

double r = 0;

if (m_Real != 0) r = Math.Atan(m_Imaginary / m_Real);

return (r);

}

// Argument of a complex number in degree

public double ArgDeg()

{

return (180 / Math.PI) * this.Arg();

}

// overridden ToString to return format: a + bi

public override string ToString()

{

string r;

if (m_Real >= 0)

r = ' ' + m_Real.ToString();

else

r = m_Real.ToString();

if (m_Imaginary >= 0)

r += " + " + ImaginaryToString(m_Imaginary);

else

r += " - " + ImaginaryToString(-m_Imaginary);

return r + "i";

}

string ImaginaryToString(double v)

{

if (v == 1) return ""; else return v.ToString();

}

// ToString Gaussian to return format: (a, b)

public string ToStringGaussian()

{

return string.Format("({0}, {1})",

m_Real.ToString(), m_Imaginary.ToString());

}

#region OVERLOAD OPERATORS

// overloaded binary + operator

public static ComplexNumber operator +(ComplexNumber a, ComplexNumber b)

{

return C(a.Real + b.Real, a.Imaginary + b.Imaginary);

}

// overloaded unary - operator

public static ComplexNumber operator -(ComplexNumber a)

{

return C(-a.Real, -a.Imaginary);

}

// overloaded binary - operator

public static ComplexNumber operator -(ComplexNumber a, ComplexNumber b)

{

return C(a.Real - b.Real, a.Imaginary - b.Imaginary);

}

// overloaded binary * operator

public static ComplexNumber operator *(ComplexNumber a, ComplexNumber b)

{

return C(a.Real * b.Real - a.Imaginary * b.Imaginary,

a.Real * b.Imaginary + b.Real * a.Imaginary);

}

// overloaded binary / operator

public static ComplexNumber operator /(ComplexNumber a, ComplexNumber b)

{

double c1, c2, d;

d = b.Real * b.Real + b.Imaginary * b.Imaginary;

if (d == 0)

{

return C(0, 0);

}

else

{

c1 = a.Real * b.Real + a.Imaginary * b.Imaginary;

c2 = a.Imaginary * b.Real - a.Real * b.Imaginary;

return C(c1 / d, c2 / d);

}

#endregion

//shortcut to return new ComplexNumber ...

static ComplexNumber C(double r, double i) {return new ComplexNumber(r, i);}

}

La matemática de complejos es bastante más amplia, otras operaciones y funciones transcendentales puedes ser agregadas a la clase.

Note que escribí el método Let, el cual se usa para asignar un número complejo a otro, ya que el operador “=” no se puede sobrecargar.

Eliminación de Gauss-Jordan Aplicada a Matrices de Complejos

Sobrecargando, y ajustando en algo el código de GaussianElimination, tengo la portentosa solución de ecuaciones lineales para números complejos. He aquí el resultado, el cual se agrega a la clase LinearEquationsSolver.cs.

Eliminación Gauss-Jordan para números complejos

/// <summary>

/// GaussianElimination()

/// Gaussian elimination is a method for solving matrix equations

/// By Harvey Triana

/// </summary>

/// <param name="a"> The matrix</param>

/// <param name="r"> The solution array</param>

/// <returns>Success function</returns>

public static bool GaussianElimination(ComplexNumber[,] a, ComplexNumber[] r)

{

ComplexNumber t = new ComplexNumber();

ComplexNumber s = new ComplexNumber();

int i, l, j, k, m, n;

try

{

n = r.Length - 1;

m = n + 1;

for (l = 0; l <= n - 1; l++)

{

j = l;

for (k = l + 1; k <= n; k++)

{

if (!(a[j, l].Abs() >= a[k, l].Abs())) j = k;

}

if (!(j == l))

{

for (i = 0; i <= m; i++)

{

t.Let(a[l, i]);

a[l, i].Let(a[j, i]);

a[j, i].Let(t);

}

for (j = l + 1; j <= n; j++)

{

t = (a[j, l] / a[l, l]);

for (i = 0; i <= m; i++) a[j, i] -= t * a[l, i];

}

r[n] = a[n, m] / a[n, n];

for (i = 0; i <= n - 1; i++)

{

j = n - i - 1;

s.Real = 0;

s.Imaginary = 0;

for (l = 0; l <= i; l++)

{

k = j + l + 1;

s += a[j, k] * r[k];

}

r[j] = ((a[j, m] - s) / a[j, j]);

}

return true;

}

catch

{

return false;

}

¿En que se diferencian? Note que la magnitud de un número complejo se mide con la función Absoluto para número complejo. Las asignaciones se hacen por el método Let, y de resto lo hacemos naturalmente gracias a la sobrecarga de operadores; quedando oculta la complejidad de las operaciones matemáticas.

El ejemplo de uso es similar al anterior. La salida formateada ahora cambia solo en los parámetros.

Ejemplo de la Eliminación Gauss-Jordan para números complejos

using System;

namespace Article_GJE

{

class Program

{

static void Main(string[] args)

{

Sample2();

}

static void Sample2()

{

/* Solucionar el siguiente sistema de ecuaciones lineales

(2+i)x - (i+1)y = 7i | (2,1) (-1,-1) (0, 7) |

(1+i)x + (2+3i)y = -2i | (1,1) ( 2, 3) (0,-2) | */

ComplexNumber[,] a = {{C(2, 1), C(-1, -1), C(0, 7)},

{C(1, 1), C( 2, 3), C(0, -2)}};

ComplexNumber[] r = new ComplexNumber[2];

ShowMatrix(a, "Ejemplo 2");

if (LinearEquationsSolver.GaussianElimination(a, r))

ShowSolution(r);

else

Console.WriteLine("No es un sistema de ecuaciones lineales");

}

static ComplexNumber C(double Real,double Imaginary)

{

return new ComplexNumber(Real, Imaginary);

}

#region formated output

static void ShowMatrix(ComplexNumber[,] a, string Title)

{

Console.WriteLine(Title + '\n');

for (int i = 0; i <= a.GetUpperBound(0); i++)

{

Console.Write('|');

for (int j = 0; j <= a.GetUpperBound(1); j++)

{

Console.Write(a[i, j].ToStringGaussian() + ' ');

}

Console.Write("| \n");

}

Console.WriteLine("\n");

}

static void ShowSolution(ComplexNumber[] r)

{

Console.WriteLine("Solución por Eliminación Gaussiana");

for (int i = 0; i <= r.GetUpperBound(0); i++)

{

Console.WriteLine(r[i].ToStringGaussian());

}

Console.WriteLine("\n");

}

#endregion

}

Un detalle para resaltar en este ejemplo es la asignación de la matriz, la cual se hizo con la función C, la cual retorna una variable de la clase ComplexNumber. Esto me permite escribir una asignación de variable más leíble, y simple. Es decir:

Modo largo de definir la matriz de complejos

ComplexNumber[,] a =

{{new ComplexNumber(2, 1), new ComplexNumber(-1, -1), new ComplexNumber(0, 7)},

{new ComplexNumber(1, 1), new ComplexNumber(2, 3), new ComplexNumber(0, -2)}};

Modo corto de definir la matriz de complejos

ComplexNumber[,] a = {{C(2, 1), C(-1, -1), C(0, 7)},

{C(1, 1), C( 2, 3), C(0, -2)}};

Pasar por Valor la Matriz

Bien, si estudio el código sabrá que las matrices, tanto en el caso de números reales como complejos, se pasan por referencia. El método de Gauss-Jordan normalmente hace una transformación de la matriz durante la solución. Si quisiéramos hacer alguna operación posterior con la matriz, por ejemplo comprobar la solución, no seria valido. Por ejemplo, la siguiente rutina la use para comprobar la validez de la función GaussianElimination aplicada a números complejos

Rutina que pueba la solucion Gauss-Jordan para números complejos

static void HardTest(ComplexNumber[,] a, ComplexNumber[] r)

{

ComplexNumber s = new ComplexNumber();

Console.WriteLine("Prueba de la solución");

for (int i = 0; i <= r.GetUpperBound(0); i++)

{

s.Real = 0;

s.Imaginary = 0;

for (int j = 0; j <= r.GetUpperBound(0); j++) s += a[i, j] * r[j];

Console.Write("Ecuación {0}: {1} \n", i.ToString(), s.ToString());

}

Por supuesto la función HardTest debe llamarse después de la solución. Para que HardTest retorne correctamente los valores, la matriz tiene que usar sus valores originales, o alternativamente ser pasada por valor. Para pasar una matriz por valor, C# entrega una elegante solución, el método Clone, el cual aplicado al caso es así:

Pasar una matriz por valor usando el método Clone

if (LinearEquationsSolver.GaussianElimination((ComplexNumber[,]) a.Clone(), r))

{

ShowSolution(r);

HardTest(a, r);

}

Conclusiones

Este articulo pone de manifiesto las potentes capacidades de C# para escribir aplicaciones que conciernen a problemas generales de matemática aplicada e ingeniería.

Si usted tiene una sugerencia para mejorar el código, por favor envíeme un e-mail.

Referencias

1. Gaussian Elimination, http://mathworld.wolfram.com/GaussianElimination.html

2. Eliminación de Gauss-Jordan, http://es.wikipedia.org/wiki/Eliminaci%C3%B3n_de_Gauss-Jordan

3. Solución de Ecuaciones Lineales Simultaneas, http://vexpert.mvps.org/articles/GJE_VB6.htm

4. Harness the Features of C# to Power Your Scientific Computing Projects, http://msdn.microsoft.com/msdnmag/issues/04/03/ScientificC/

-----