sábado, 11 de marzo de 2017

métodos de Gauss y Gauss-Jordan

      • Introducción y antecedentes (para qué sirve y qué resuelven) los métodos
      Los metodos de eliminacion de Gauss y Gauss-Jordan, son métodos que se utilizan para resolver sistemas de ecuaciones de cualquier número.

      • Requisitos de los métodos
      Para aplicar alguno de estos métodos, se necesita, primero que nada, un sistema de ecuaciones y de ella generar una matriz de coeficientes y una matriz de constantes. Posteriormente, se crea una matriz aumentada, esto mediante la unión de la matriz (o vector) de constantes a la de coeficientes. Para asegurar el funcionamiento, el sistema de ecuaciones debe tener solución, por lo que no deben ser paralelas o que los coeficientes de una ecuación sean múltiplos de los otros coeficientes de las demás ecuaciones.

      También es importante considerar que para ambos métodos, las ecuaciones del sistema de ecuaciones inicial necesitan estar acomodadas de manera que los coeficientes más grandes estén en la diagonal.

      Esto significa que si comenzaramos con el siguiente sistema de ecuaciones:
      Lo primero que tendríamos que hacer es reacomodarlo de la siguiente manera:
      • Diferencias entre uno y otro método
      La diferencia entre los métodos de Gauss y Gauss-Jordan es que en Gauss necesitas hacer backwards substitution y en Gauss-Jordan no. Esto, ya que el método de Gauss-Jordan va un poco más allá que Gauss y no solo se concentra en hacer los coeficientes debajo de la diagonal cero, sino también los coeficientes arriba de la diagonal.

      En el método de Gauss hace eliminación hacia adelante N-1 veces, siendo N la cantidad de ecuaciones.

      • ¿Es un método iterativo? (¿si?, ¿no?, ¿por qué?)
      A pesar de que haces iteraciones para la eliminación hacia adelante y sustitución hacia atrás en el caso de la eliminación de Gauss; no es un método iterativo. Esto porque las iteraciones que hace están definidas por el tamaño de la matriz y se realizan como parte de un proceso numérico. Por otra parte, un método iterativo se detiene a partir de la evaluación de un tipo de error cada vez que se repite un proceso que empieza con un valor semilla... esta característica no está presente en ninguno de estos dos métodos.

      • Pivoteo, y escalamiento. (¿Por qué son útiles?)

      El pivoteo es el proceso de reacomodo de la matriz para que el coeficiente más grande sea el pivote antes de usar uno de estos métodos. Se le llama pivoteo parcial si sólo se intercambian renglones, y pivoteo completo si se cambian tanto renglones como columnas. Esto con el motivo de disminuir la probabilidad de un error de redondeo cuando el elemento pivote es cercano a 0 y provoca una división indefinida. El pitoe completo es raramente usado porque agrega más complejidad a la computadora al momento de intercambiar filas y columnas.

      El escalamiento se puede definir como un proceso para simplificar los coeficientes de un sistema de ecuacion; por ejemplo, dado un sistema de ecuaciones, realizar un escalamiento implica, para cada renglón o ecuación, dividir por alguna constante tal que el coeficiente más pequeño de cada renglón sea uno. Esto, con el propósito de disminuir errores de redondeo que se presentan usualmente cuando los coeficientes de alguna ecuación es relativamente más grande que el de otra. Un ejemplo es el siguiente:

      2x1 + 100,000x2 = 100,000
      x1 + x2 = 2
      Usando escalamiento el sistema queda de la siguiente forma
      0.00002x1 + x2 = 1
      x1 + x2 = 2

      • Diagrama de flujo Gauss Elimination

      La mitad de este diagrama de flujo es exactamente igual a la primera mitad del diagrama de GaussJordan, para la segunda mitad es el primero proceso pero en sentido contrario, donde se busca crear una matriz identidad.
      • Código fuente para métode Gauss en Scilab

      //Declara la matriz cuadrada de coeficientes como A
//Declara el vector columna de constantes como b

//Generar matriz aumentada y desplegar
AUM = [A b]
X = AUM
[rX cX] = size(X)
//Eliminación hacia adelante
for p = 1 : 1 : rX-1   //Fila del pivote
 pivote = X(p,p);
 for i = p + 1 : 1 : rX   //Fila a procesar
     prim_fila = X(i,p)
     X(i,:) = X(i,:) - (X(p, :)/pivote)*prim_fila
 end
end

x(rX) = X(rX, cX)/X(rX, cX - 1)
for i = rX : -1 : 1
 s = 0
 for j = i + 1 : rX
     s = s + X(i, j) * x(j)
 end
 x(i) =( X(i,cX) - s ) / X(i,i)
end


disp (X)
disp(x)

      • Código fuente para métode Gauss-Jordan en Scilab

       
       1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
      // para resolver esta ecuación simultánea en SCILAB:
 // 2x -y +1 = 7
 // x +2y -z = 6
 // x + y + z = 12

a=[2,-1,1;1,2,-1;1,1,1]
 // matriz de coeficientes, dimensiones 3x3:
  // 2 -1  1
  // 1  2 -1
  // 1  1  1

b=[7;6;12]
 // matriz de constantes, dimensiones 3x1:
  // 7
  // 6
  // 12

aumAB=[a b]
 // matriz aumentada, que es así:
  // 2 -1  1  7
  // 1  2 -1  6
  // 1  1  1  12

// y para hacer el proceso de Gauss-Jordan elimination...
rref(aumAB)
 // y nos da:
  // 1  0  0  3
  // 0  1  0  4
  // 0  0  1  5

 // entonces, la respuesta para nuestra ecuación simúltanea es
 // x=3, y=4, z=5

      • Casos de falla y lo que tu programa hace en estos casos
      En una matriz cuadrada, el método de Gauss fallará si la determinante de la matriz es cero. Otra manera en que se puede ver que esto ocurre es cuando alguna de las ecuaciones del sistema de ecuaciones es múltiplo escalar de otra.

      • Pruebas y resultados con sistemas de más de 3 incógnitas
      Podemos incluso probar estos métodos con ecuaciones de más de 3 incógnitas. Con 4 incógnitas:
      2w+2x-5y+z=-16
      -w+x+6y-z=15
      2w-x+y+6z=3
      w+x+2y-z=7

      Ambos métodos dan como respuesta:

      w=2, x=-2, y=3, z=-1

      Lo cual es correcto.

      Y con 6 incógnitas podemos observar lo siguiente:
      a +  b - 2c +  d + 3e -  f  =    4
         2a -  b +  c + 2d +  e - 3f  =   20
          a + 3b - 3c -  d + 2e +  f  =  -15
         5a + 2b -  c -  d + 2e +  f  =  - 3
        -3a -  b + 2c + 3d +  e + 3f  =   16
         4a + 3b +  c - 6d - 3e - 2f  =  -27

      Ambos métodos dan como respuesta:

      a=1, b=-2, c=3, d=4, e=2, f=-1

      Lo cual es correcto.

      • Conclusiones

      Estos dos métodos de eliminación, a diferencia del método de Cramer, son menos complejos para la computadora cuando el tamaño de la matriz crece; esto es porque no realizan operaciones grandes, como lo hace el cálculo de una determinante, y solamente recorren cada elemento de la matriz para irla modificando Funcionan mejor para matrices cuadradas mayores a 4, por lo dicho anteriormente, y su implementación es muy fácil en scilab (véase código de GaussJordan).
      Dado que no son iterativos, la respuesta no tiende a un porcentaje de error, simplemente encuentra la solución si es que existe. Sin embargo, al ser implementados en la computadora, siguen existiendo los errores de redondeo. A pesar de esto, antes de aplicar alguno de estos métodos, existe el procedimiento de pivoteo y escalamiento para disminuir la probabilidad de sufrir dicho error.
Publicado por en No hay comentarios:

martes, 7 de marzo de 2017

análisis comparativo de métodos

Tabla Comparativa para
Bisección, Newton-Raphson, Secante y Bairstow
Previamente hemos hablado a detalle de ya varios métodos numéricos para encontrar las raíces de un polinomio, por lo que surge la necesidad de comparar uno con otro para tomar la mejor decisión cuando se requiera. Es así, que a continuación presentaremos una tabla comparativa tomando en cuenta diferentes aspectos de cada uno de los métodos, que cómo podrás ver, varían considerablemente. El porqué de cada uno de estos datos puede ser entendido más claramente leyendo nuestros posts anteriores para cada uno de los métodos.
Mas información sobre los métodos...
Bisección:
Newton-Raphson:
Secante:
Bairstow


Publicado por en No hay comentarios:

viernes, 3 de marzo de 2017

método Bairstow

Otro método numérico: Bairstow.
Llamado en honor a su creador, Leonard Bairstow, tras publicarlo por primera vez en el libro Applied Aerodynamics de 1920.

Leonard Bairstow
Video explicativo:
¿Cómo funciona? Antes que nada, Bairstow es un método abierto: necesita 2 estimaciones iniciales, pero no es necesario que éstas encierren a la raíz en lo absoluto.
Todo comienza con un polinomio. Lo que hace este método es ir reduciendo el grado del polinomio 2 grados a la vez a través de divisiones sintéticas dobles.
Gracias a estas divisiones sintéticas, podremos formular un sistema de ecuaciones de 2 incógnitas donde a través del método Cramer (del cual hablaremos más a detalle en otro post), y obtener factores del polinomio original.
Con todos estos factores, se obtendrán las raíces del polinomio original.
Cuando encontremos un factor por el cual dividamos el polinomio original y nos resulte en un polinomio de grado 1 ó 0, significa que ya tenemos suficientes factores para encontrar todas las raíces.



Ahora bien, ¿bajo qué criterio determinamos que la división sintética resulte en un polinomio dos grados menor? Aquí es donde se encuentra nuestro margen de error: qué tan cercanos están los 2 valores de menor grado del resultado de la división sintética a cero. Cuando se acercan lo suficiente, los consideramos como 0 y entonces el resultado de la división sintética tiene es 2 grados menor. Por último, si el resultado es de grado 1 ó 0, significa que hemos encontrado todos los factores del polinomio original. Las raíces de cada uno de estos factores conforman nuestra respuesta final.
Código fuente en C++, con la función (x^5)-(3.5x^4)+(2.75x^3)+(2.125x^2)-(3.875x)+1.25:

 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <iostream>
#include <cmath>

using namespace std;
double EPSILON = 0.01;

void sinteticWeirdFormula( double *arr, double *arrOr, int size, double r, double s ) {
    arr[size-1] = arrOr[size-1];
    arr[size-2] = arr[size-1] * r + arrOr[size-2];
    
    for ( int i = size - 3; i >= 0; i-- ) {
        arr[i] = ( arr[i + 1] * r ) + ( arr[i + 2] * s ) + arrOr[i];
    }
    return;
}

int main() {
    double arrVal[] = {1.25, -3.875, 2.125, 2.75, -3.5, 1}; // Coeficientes del polinomio
                    //  0       x1     x2    x3    x4   x5
    int size = sizeof(arrVal) / sizeof(arrVal[0]);
    int count = 1;
    
    double arrB[size]; //Initial array B
    double arrC[size]; //Initial array C
    
    double r = -1;
    double s = -1;
    
    double dr; //delta r
    double ds; //delta s
    
    do {
        sinteticWeirdFormula( arrB, arrVal, size, r, s );
        sinteticWeirdFormula( arrC, arrB, size, r, s );
        /** Cramer
         * A = arrC[2]
         * B = arrC[3]
         * C =  - arrB[1]
         * P = arrC[1] 
         * Q = arrC[2]
         * R = - arrB[0]
         * dr = (CQ - BR) / (AQ - PB)
         * ds = (AR - CP) / (AQ - PB)
        */
        dr = ( -arrB[1] * arrC[2] - arrC[3] * -arrB[0] ) / ( arrC[2] * arrC[2] - arrC[1] * arrC[3] );
        ds = ( -arrB[0] * arrC[2] - arrC[1] * -arrB[1] ) / ( arrC[2] * arrC[2] - arrC[1] * arrC[3] );
        
        double errorR = ( dr / r ) * 100;
        double errorS = ( ds / s ) * 100;
        
        if ( abs(errorR) <= EPSILON && abs(errorS) <= EPSILON ) {
            break;
        }
        r += dr;
        s += ds;
        count++;
    } while (true);
    
    double x1 = ( r + sqrt( r * r + 4 * s ) ) / 2;
    double x2 = ( r - sqrt( r * r + 4 * s ) ) / 2;
    
    cout << "Root 1: " << x1 << endl
    << "Root 2: " << x2 << endl;
    
    cout << "r: " << r << endl
    << "s: " << s << endl;
    
    for ( int i = 0; i < size; i++ ) {
        cout << "Coeficientes de la funcion con menor grado " << arrC[i] << endl;
    }
    /*for ( int i = 0; i < size; i++ ) {
        cout << "b: " << arrB[i] << endl;
    }*/   
    return 0;
}


Ahora bien, debemos también de considerar que este método no funciona para funciones trascendentales, ni polinomios con coeficientes complejos. También, si se le da una precisión demasiado pequeña, Bairstow puede tener divergencia.

 Ahora que hemos analizado a detalle este método, podemos ver cómo aunque este método cuenta con ventajas muy grandes como no necesitar la derivada de la función como en Newton-Raphson y permitirnos encontrar más de una raíz (incluso imaginarias), sigue sin ser un método perfecto para todos los casos. Lo importante es tomar la mejor decisión sobre qué método usar para cada caso específico.
Fuentes:

Publicado por en No hay comentarios:
Suscribirse a: Entradas (Atom)