Práctica 5. Sistemas de ecuaciones lineales I (Método de Gauss)

En la práctica 2 vimos cómo expresar un sistema de ecuaciones lineales en forma matricial, cómo comprobar que un vector es la solución de un sistema Ax=b, y cómo usar la 'división izquierda' para hallar la solución de un sistema Ax=b con matriz A regular.

También vimos ejemplos de sistemas mal condicionados, y sistemas donde la 'división izquierda' no da una solución porque el sistema es incompatible.

En esta práctica veremos cómo resolver sistemas triangulares y el método de Gauss para resolver sistemas de ecuaciones lineales.

Resolución de un sistema triangular superior.
Resolución de un sistema triangular inferior.
Método de Gauss (Ejemplo del proceso de hacer ceros).
Método de Gauss (Programación del método).

Resolución de un sistema triangular superior.

Consideremos un sistema triangular de tres ecuaciones y tres incógnitas como el siguiente:

A11*x1 + A12*x2 + A13*x3 = b1,
         A22*x2 + A23*x3 = b2,
                  A33*x3 = b3.

Para hallar su solución podemos despejar x3 en la última ecuación, x3=b3/A33. Con este valor de x3, podemos despejar x2 en la segunda ecuación x2=(b2-A23*x3)/A22. Y con los valores hallados de x3 y x2 podemos despejar x1 en la primera ecuación x1=(b1-A12*x2-A13*x3)/A11.

Emplearemos este método para resolver un sistema de n ecuaciones lineales y n incógnitas cuya matriz de coeficientes tiene todos los coeficientes por debajo de la diagonal iguales a cero, es decir, es triangular superior.

Consideramos un sistema como el siguiente:

A(1,1)*x(1) + A(1,2)*x(2) + ... +  A(1,n-1)*x(n-1)  +  A(1,n)*x(n)  = b(1),
              A(2,2)*x(2) + ... +  A(2,n-1)*x(n-1)  +  A(2,n)*x(n)  = b(2),

...

                                  A(n-1,n-1)*x(n-1) + A(n-1,n)*x(n) = b(n-1),
                                                        A(n,n)*x(n) = b(n).

Podemos despejar x(n) de la última ecuación con lo que x(n)=b(n)/A(n,n). Conociendo x(n) podemos despejar x(n-1) de la penúltima ecuación.

Siguiendo este proceso, supongamos que ya hemos hallado x(n), x(n-1), ... x(i+1) y queremos despejar x(i) de la ecuación i-ésima

       A(i,i)*x(i) + A(i,i+1)*x(i+1) + ... + A(i,n)*x(n) = b(i),

entonces tenemos

       x(i)=(b(i) - S )/A(i,i),

donde S = A(i,i+1)*x(i+1) + ... + A(i,n)*x(n) .

Nos proponemos expresar el valor de S como un producto escalar de un vector fila por un vector columna.

S = A(i,i+1)*x(i+1) + ... + A(i,n)*x(n) .

En Matlab, el vector fila A(i,i+1:n) contiene los elementos que están en la fila i, y en las columnas desde la i+1 hasta la n, es decir A(i,i+1), A(i,i+2), A (i,i+3), ... A(i,n-1), A(i,n).

El vector columna x(i+1:n) contiene los elementos de x que van desde el i+1 hasta el n, es decir x(i+1), x(i+2), ..., x(n-1), x(n). Por tanto el valor de S se puede calcular como S=A(i,i+1:n)*x(i+1:n), con lo que usamos las operaciones entre vectores de Matlab.

Podemos escribir el pseudocódigo que permite resolver el sistema.

Entrada A, b.
n = número de filas de A.
Inicializamos x.
Para i=n, n-1, n-2, ..., 1
   x(i)= (b(i)- A(i,i+1:n)*x(i+1:n))/A(i,i)
Siguiente i.

El siguiente código de Matlab permite resolver un sistema triangular superior.

% Matriz del sistema
A=[1 3 2 2 1; 0 1 2 4 1; 0 0 2 1 -1; 0 0 0 3 -2; 0 0 0 0 2];
% Término independiente
b=[-2;-1;-2;-9;6];
% Hallamos el número de ecuaciones del sistema.
n=length(A);
% Definimos un vector-columna para almacenar la solución.
x=zeros(n,1);
% El bucle va desde la última ecuación hasta la primera
for i=n:-1:1
   x(i)= (b(i)- A(i,i+1:n)*x(i+1:n))/A(i,i);
end
disp('La solución del sistema es ')
disp(x)

La solución del sistema es 
     1
    -2
     1
    -1
     3

Ej 1. Empleando el código que hemos visto escriba un script que permita resolver el sistema

x1 + 2*x2 + 3*x3 +   x4 +   x5 = 7,
       x2 + 4*x3 +   x4 -   x5 = 6,
          -   x3 + 2*x4 -   x5 = -5,
                   3*x4 +   x5 = -2,
                          2*x5 = 2.

Compruebe el resultado con la 'división izquierda' (Sol. x1=1, x2=0, x3=2, x4=-1, x5=1).

Ej 2. Modifique el programa para que imprima (usando disp o fprintf) un mensaje como 'El elemento 3 de la diagonal es cero' si algún elemento de la diagonal vale cero, pero que no interrumpa el método. Aplique el programa a la resolución del sistema A*x=b donde A=[1,2,3,1,1;0,2,4,-1,-1;0,0,0,2,-1;0,0,0,3,1;0,0,0,0,-2], b=[2;3;-5;-5;-2]. (Sol. x5=1, x4=-2.) (El sistema planteado tiene infinitas soluciones, es compatible indeterminado. Resolviendo a mano se podrían dar las soluciones dependiendo de un parámetro : x5=1, x4=-2, x3=x1-1, x2=-2*x1+3 )

Resolución de un sistema triangular inferior.

Consideremos ahora un sistema de n ecuaciones lineales y n incógnitas cuya matriz de coeficientes tiene todos los coeficientes por encima de la diagonal iguales a cero, es decir, es triangular inferior.

A(1,1) *x(1)                                                        =b(1),
A(2,1) *x(1) +  A(2,2) *x(2)                                        =b(2),
...
A(n-1,1)*x(1) + A(n-1,2)*x(2) + ... + A(n-1,n-1)*x(n-1)              =b(n-1),
A(n,1) *x(1) +   A(n,2)*x(2) + ... +   A(n,n-1)*x(n-1) + A(n,n)*x(n)=b(n).

En este sistema se puede despejar x(1) de la primera ecuación. Una vez hallado el valor de x(1), se puede hallar x(2) despejando en la segunda ecuación. Supongamos que ya hemos hallado x(1), x(2), ..., x(i-1). Considerando la ecuación i-ésima

( A(i,1)*x(1) + A(i,2)*x(2) + ... + A(i,i-1)*x(i-1) ) + A(i,i)*x(i) =b(i)

podemos despejar x(i).

Ej 3. a) Despeje (en papel) x(i) en la última ecuación. b) Exprese (en papel) el valor de x(i) usando un producto escalar de un vector fila por un vector columna. c) Programe la resolución de un sistema triangular inferior. d) Compruebe el programa resolviendo el sistema A*x=b donde A=[1,0,0,0,0; 0,2,0,0,0; 1,1,1,0,0; 0,-1,0,3,0; -2,0,1,0,-2], b=[1;4;6;10;-9]. (Sol. x1=1, x2=2, x3=3, x4=4, x5=5 )

Método de Gauss (Ejemplo del proceso de hacer ceros).

El método de Gauss permite resolver sistemas de ecuaciones lineales. Este método consiste en realizar, en primer lugar, determinadas operaciones elementales por filas sobre la matriz ampliada del sistema hasta que ésta se transforme en otra en la que todas sus columnas salvo la última formen una matriz triangular superior (proceso de hacer ceros). Y, a continuación, resolver el sistema triangular superior asociado a la matriz resultante (que es equivalente al original).

Veamos un ejemplo de la primera parte del método de Gauss (proceso de hacer ceros):

Consideramos el siguiente sistema de tres ecuaciones lineales y tres incógnitas

x + 2*y + 3*z = 5,
2*x +   y       = 1,
x +   y + 2*z = 4.

Para hacer las operaciones elementales sobre las filas de la matriz ampliada [A b] lo que haremos será realizar las mismas operaciones elementales sobre la matriz A y sobre el vector b de su formulación matricial.

Definimos la matriz de coeficientes del sistema y el vector de términos independientes.

A=[1 2 3; 2 1 0; 1 1 2];
b=[5;1;4];

Hacemos ceros en la primera columna debajo del elemento pivote A(1,1) .

A la segunda fila le restamos dos veces la primera

A(2,:)=A(2,:)-2*A(1,:) ;
b(2)=b(2)-2*b(1) ;
% A la tercera fila le restamos la primera
A(3,:)=A(3,:)-A(1,:) ;
b(3)=b(3)-b(1) ;
% Vemos que hemos hecho ceros en la primera columna
disp([A b])

     1     2     3     5
     0    -3    -6    -9
     0    -1    -1    -1

Para hacer cero en A(3,2) debemos restar a la tercera fila un múltiplo de la segunda.

A(3,:)=A(3,:)-(-1/-3)*A(2,:) ;
b(3)=b(3)-(-1/-3)*b(2) ;
% Vemos que hemos hecho ceros debajo de la diagonal, con lo que el sistema
% que resulta es triangular superior.
disp([A b])

     1     2     3     5
     0    -3    -6    -9
     0     0     1     2

Método de Gauss (Programación del método).

En lo que sigue vamos a implementar el método de Gauss para sistemas en los que en el proceso de hacer ceros no es preciso permutar filas. La programación del método de Gauss cuando se precisa permutar filas en el proceso de hacer ceros es algo más complicada y no la veremos aquí.

Siguiendo la notación dada en clase al explicar el método de Gauss, la variable i indicará la columna en la que vamos a hacer ceros. La variable j indica la fila en la que estamos haciendo ceros.

Para simplificar la notación y ahorrar memoria, no usaremos superíndices, llamaremos A0 a la matriz original del sistema y b0 al término independiente original. La matriz A y el vector b contendrán los sucesivos sistemas equivalentes, si se acaban de hacer ceros en la columna i, A indicará la matriz con superíndice i+1.

Dado i, suponemos que ya se han hecho ceros en las columnas 1, ..., i-1 y hacemos ceros en la columna i en las filas j>=i+1 mediante la operación por filas 'sumar a la fila j la fila i multiplicada por lj ' donde lj es el valor lj=-A(j,i)/A(i,i) .

El programa ProgGauss.m implementa el proceso de hacer ceros del método de Gauss para resolver sistemas en los que en el proceso de hacer ceros no es preciso permutar filas.

Dicho programa incluye las instrucciones tic y toc que sirven para mostrar el tiempo que emplea Matlab en ejecutar un conjunto de instrucciones. Basta con poner la instrucción tic antes del conjunto de instrucciones y toc al final. Al ejecutarse la instrucción toc se muestra el tiempo que emplea Matlab en ejecutar el conjunto de instrucciones.

Definimos en forma matricial el sistema que queremos resolver.

A0=[1,2,3,2 ; 3,-3,-6,1 ; 2,1,7,0 ; 2,3,2,2];
b0=[12;-10;22;9];
% Inicializamos la matriz A y el vector b en los que efectuaremos las
% operaciones elementales por filas.
A=A0;
b=b0;
n=length(A);

Haremos ceros en la columna i, donde i va desde la primera columna a la penúltima.

for i=1:n-1
   % j va desde la fila i+1 hasta la última.
   for j=i+1:n
   lj=- A(j,i)/A(i,i);
   % Sumamos a la fila j un múltiplo de la fila i.
   A(j,:)= A(j,:)+lj * A(i,:);
   % Efectuamos la misma operación en el vector b.
   b(j)=b(j)+lj * b(i);
   end
   % fprintf('Tras hacer ceros en la columna %d\n',i);
   % disp('queda el sistema equivalente:')
   % disp([A b])
end

Una vez transformado el sistema original en otro equivalente pero triangular superior, hay que resolver éste último y para ello se puede utilizar el método que se vió al principio de esta práctica para resolver sistemas triangulares superiores o también, por comodidad, la división izquierda.

x=A\b

Podemos comprobar que la solución del sistema triangular A*x=b obtenida es solución del sistema original A0*x=b0.

A0*x-b0

Ej 4. Comprueba que el sistema A0*x=b0 dado por A0=[1,0,3,2,-1 ; 3,0,-2,1,1 ; 2,7,1,0,3 ; 1,2,3,-1,2; 1,2,3,-1,0] b0=[5;4;3;2;1] es compatible determinado y que A0\b0 calcula la solución del sistema, Comprueba que, sin embargo, el método de Gauss implementado no obtiene la solución del sistema. ¿A qué se debe? (Indicación : observe las operaciones que se realizan con los elementos A(i,i). )

Ej 5. Escriba un programa que aplique el método de Gauss al sistema A0*x=b0 donde A0=[1,0,3,2,-1 ; 3,1,-2,1,1 ; 2,7,1,0,3 ; 1,2,3,-1,2; 1,2,3,-1,0]; b0=[17;16;44;22;20] de forma que transforme el sistema a uno triangular superior equivalente, y luego aplique el método visto para resolver sistemas triangulares superiores. (Sol. x1=5, x2=4, x3=3, x4=2, x5=1).

Las operaciones que se aplican en esta implementación del método de Gauss actuan sobre toda la fila a pesar de que los primeros elementos de cada fila serán cero, con lo cual se están efectuando operaciones innecesarias.

Así al hacer ceros en la columna i, en las filas j=i+1,...,n los elementos A(j,i) son cero. Se proponen dos formas de evitar operaciones innecesarias modificando la línea A(j,:)= A(j,:)+lj * A(i,:);

i) Emplear un bucle for con variable k que vaya operando en las columnas de la i+1 hasta la n. ii) Hacer la operación por filas indicando que deben considerarse las columnas de la i+1 hasta la n, es decir, debe asignarse sólamente A(j,i+1:n).

Ej 6. Compruebe (con papel) que, en el método de Gauss implementado, los elementos A(j,i) para j=i+1,...,n son cero.

Ej 7. Implemente el método de Gauss sustituyendo la línea A(j,:)= A(j,:)+lj * A(i,:) por un bucle for como se indica en i) y asignando A(j,i)=0.

Puede comprobar el método con el mismo ejemplo A0=[1,2,3,2 ; 3,-3,-6,1 ; 2,1,7,0 ; 2,3,2,2]; b0=[12;-10;22;9];

Ej 8. Implemente el método de Gauss sustituyendo la línea A(j,:)= A(j,:)+lj * A(i,:) por una línea que comienze por A(j,i+1:n)= ... y asignando A(j,i)=0.