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 no viembre de 2006
Introducción
La “Eliminación de GaussGauss -Jordan”, es una de una de mis soluciones preferidas, el código lo escribí inicialmente en FORTRAN en 1992, Visual Basic clásico clásico en 1996, y ahora en C# 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 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 { ///
/// GaussianElimination() /// Gaussian elimination is a method for solving matrix equations /// By Harvey Triana /// ///
The matrix ///
The solution array ///
Success function 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 ///
/// GaussianElimination() /// Gaussian elimination is a method for solving matrix equations /// By Harvey Triana /// ///
The matrix ///
The solution array ///
Success function
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, {C(1, 1), C( 2,
7)},
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), 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,
new
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); }