Sistemas de ecuaciones lineales con MATLAB

Determinante

Un determinante es un número asociado a una matriz cuadrada. Para una matriz cuadrada de dimensión 2×2

A=[ a 11 a 12 a 21 a 22 ]| A |= a 11 a 22 a 12 a 21

El determinante de una matriz m×m se desarrolla en términos de una combinación de determinates de matrices de dimensión m-m-1 y así, sucesivamente, hasta llegar a los determinantes de matrices 2×2.

Por ejemplo, el determinante de una matriz 3×3 es

| a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 |= a 11 | a 22 a 23 a 32 a 33 | a 12 | a 21 a 23 a 31 a 33 |+ a 13 | a 22 a 3 a 32 a 33 |

El determinante se calcula mediante la siguiente fórmula

det(A)= j=1 n (1) 1+j a 1j det( M 1j )

Donde M1j es una submatriz obtenida eliminando la fila 1 y la columna j de la matriz A

En la página titulada Matrices hemos visto como se accede a los elementos de una matriz

Tomemos la matriz cuadrada A de dimensión 4×4, la submatriz M12 que se obtiene eliminando la primera fila y la segunda columna es

>> A=[1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]
A =
     1     2     3     4
     5     6     7     8
     9    10    11    12
    13    14    15    16

>> M12=[A(2:4,1),A(2:4,3:4)]
M12 =
     5     7     8
     9    11    12
    13    15    16

Definimos la función determinante que toma la matriz A de dimensión n×n y produce matrices M1j de dimensión n-1×n-1 y calcula el determinante de acuerdo con la fórmula anterior, de forma recursiva

function d = determinante(A)
    n=length(A);
    if n==1;
        d=A(1,1);
    else
        d=0; 
        sgn=1;
        for j=1:n
            M1j=[A(2:n,1:j-1) A(2:n,j+1:n)];
            d=d+sgn*A(1,j)*determinante(M1j);
            sgn=-sgn;
        end
    end
end

Calculamos el determinante de esta matriz utilizando la función determinante

| 1 0 3 2 1 0 3 1 2 |=13

>> A=[-1 0 3; 2 -1 0; 3 1 -2]
A =
    -1     0     3
     2    -1     0
     3     1    -2
>> determinante(A)
ans =    13

En MATLAB la función det(A) calcula el determinante de la matriz cuadrada A.

>> det(A)
ans = 13

Matriz inversa

Se denomina matriz identidad I a aquella matriz cuadrada de dimensiones m×m en la cual los elementos de la diagonal valen 1 y el resto de los elementos vale cero. MATLAB dispone de la función eye(m) para crear una matriz cuadrada de dimensión m con los elementos de la diagonal principal unos y el resto ceros.

El producto de la matriz identidad I por otra matriz A nos da la matriz A

A=[ a 11 a 12 a 13 a 21 a 21 a 23 ]I=[ 1 0 0 0 1 0 0 0 1 ] A*I=[ a 11 a 12 a 13 a 21 a 21 a 23 ]

Si la matriz A es cuadrada, se obtiene el mismo resultado efectuando el producto A*I o I*A

Si A es una matriz cuadrada, B es su matriz inversa si el producto A*B=B*A=I

En MATLAB se puede obtener la matriz inversa de A elevando A a la potencia -1, A-1 o bien, mediante la función inv(A)

>> A=[-1 0 3; 2 -1 0; 3 1 -2]
A =
    -1     0     3
     2    -1     0
     3     1    -2
>> B=inv(A)
B =
    0.1538    0.2308    0.2308
    0.3077   -0.5385    0.4615
    0.3846    0.0769    0.0769
>> A*B
ans =
    1.0000         0         0
   -0.0000    1.0000   -0.0000
         0   -0.0000    1.0000

Rango de una matriz

El rango de una matriz es el máximo número de filas linealmente independientes. La función rank calcula el rango de una matriz. Sea la matriz A

A=[ 0 1 2 1 1 1 0 1 3 1 2 0 2 3 2 1 ]

>> A=[0 -1 2 1; 1 -1 0 -1; 3 1 2 0; 2 -3 2 -1];
>> rank(A)
ans =     3

Como podemos apreciar la fila cuatro es combinación lineal de la fila 1 y la fila 2:. a4j=a1j+2·a2j, j=1...4

División por la izquierda y por la derecha

La división por la izquierda se utiliza para resolver la ecuación AX=B. En esta ecuación X es el vector columna de las incógnitas, B es el vector columna de los términos independientes y A es una matriz cuadrada.

a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 [ a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 23 ][ x 1 x 2 x 3 ]=[ b 1 b 2 b 3 ] AX=B

A-1AX=IX=X

De modo que

X=A-1B

En MATLAB esta última expresión se escribe utilizando el operador división por la izquierda

X=A\B

La división por la derecha se utiliza para resolver la ecuación XC=D. En esta ecuación X y D son vectores fila y C es una matriz cuadrada.

XCC-1=DC-1

X=DC-1 o bien, X=D/C (división por la derecha)

Sistema de ecuaciones lineales

Sistema de n ecuaciones con n incógnitas. Regla de Cramer

Consideremos un sistema de n ecuaciones con n incógnitas cuya solución es única.

{ a 11 x 1 + a 12 x 2 +.... a 1n x n = b 1 a 21 x 1 + a 22 x 2 +.... a 2n x n = b 2 .... a n1 x 1 + a n2 x 2 +.... a nn x n = b n

Llamemos A a la matriz de los coeficientes, cuyo determinate d es distinto de cero

d=| a 11 a 12 ... a 1n a 21 a 22 ... a 2n ... ... ... ... a n1 a n2 ... a nn |0

La regla de Cramer nos dice que cada una de las incógnitas xj se puede obtener a partir del determinante de la matriz en la que se ha sustituido la columna j de la matriz A por el vector columna de los términos independientes b.

x j = 1 d | a 11 a 12 ... b 1 ... a 1n a 21 a 22 ... b 2 ... a 2n ... ... ... ... ... ... a j1 a j2 ... b j ... a jn ... ... ... ... ... ... a n1 a n2 ... b n ... a nn |

Sea el sistema

3 x 1 x 2 = 5 2 x 1 + x 2 + x 3 = 0 2 x 1 x 2 + 4 x 3 = 15 [ 3 1 0 2 1 1 2 1 4 ] [ x 1 x 2 x 3 ] = [ 5 0 15 ]

A=[3 -1 0; -2 1 1; 2 -1 4];
b=[5;0;15];
n=length(b);
d=det(A);
x=zeros(n,1);
for i=1:n
    Ab=[A(:,1:i-1),b,A(:,i+1:n)];
    x(i)=det(Ab)/d;
end
disp('Solución')
disp(x);
Solución
    2.0000
    1.0000
    3.0000

Sistema de m ecuaciones con n incógnitas

Sea el sistema

{ a 11 x 1 + a 12 x 2 +.....+ a 1n x n = b 1 a 21 x 1 + a 22 x 2 +.....+ a 2n x n = b 2 ... a m1 x 1 + a m2 x 2 +.....+ a mn x n = b m [ a 11 a 12 ... a 1n a 21 a 22 ... a 2n ... ... ... ... a m1 a m2 ... a mn ][ x 1 x 2 ... x n ]=[ b 1 b 2 ... b m ]

que podemos escribir Ax=b, donde A es una matriz de dimensión m×n, y x y b son dos vectores columna de longitudes n y m respectivamente. Tenemos un sistema de m ecuaciones con n incógnitas.

  1. El sistema tiene solución si y solo si el rango de la matriz A es igual al rango de la matriz ampliada A|b. Sistema compatible.
  2. Si el rango es igual al número n de incógnitas el sistema tiene una solución única. Sistema compatible determinado
  3. Si el rango es menor que el número n de incógnitas entonces hay un número infinito de soluciones. Sistema compatible indeterminado.

Vamos a ver algunos ejemplos:

1.-Sea el sistema de tres ecuaciones con dos incógnitas

{ 2 x 1 + x 2 =3 2 x 1 x 2 =0 x 1 2 x 2 =4

En forma matricial se escribe

[ 2 1 2 1 1 2 ][ x 1 x 2 ]=[ 3 0 4 ]

La matriz A y la matriz ampliada A|b son respectivamente

A=[ 2 1 2 1 1 2 ]A|b=[ 2 1 3 2 1 0 1 2 4 ]

>> A=[2 1; 2 -1; 1 -2];
>> b=[3;0;4];
>> Ab=[A b]
Ab =
     2     1     3
     2    -1     0
     1    -2     4
>> rank(A)
ans =     2
>> rank(Ab)
ans =     3

El sistema no tiene solución (primer caso). Cada una de las ecuaciones del sistema representa una recta, identificamos x1 con x y x2 con y. En la figura, vemos la representación gráfica de las rectas: y=3-2x, y=2x, y=x/2-2; que como vemos no se cortan en un punto.

line([-5,5],[13,-7])
line([-5,5],[-10,10])
line([-5,5],[-9/2,1/2])
xlabel('x')
ylabel('y')
title('No tiene solución')

2.-Sea el sistema de tres ecuaciones con dos incógnitas

{ 2 x 1 + x 2 =3 4 x 1 +2 x 2 =6 6 x 1 +3 x 2 =9 [ 2 1 4 2 6 3 ][ x 1 x 2 ]=[ 3 6 9 ]

>> A=[2 1; 4 2; 6 3];
>> b=[3;6;9];
>> Ab=[A b];
>> rank(A)
ans =     1
>> rank(Ab)
ans =     1

El sistema tiene solución, pero como el rango es menor que el número de incógnitas hay un número infinito de soluciones (tercer caso).

Tenemos tres ecuaciones iguales, la segunda es igual a la primera multiplicada por dos y la tercera es igual a la primera multiplicada por tres. En la figura, se representa la recta y=3-2x, que es la solución del sistema de ecuaciones.

line([-5,5],[13,-7])
xlabel('x')
ylabel('y')
title('Infinitas soluciones')

3.-Sea el sistema de tres ecuaciones con dos incógnitas y su representación matricial

{ 2 x 1 + x 2 =3 2 x 1 x 2 =5 x 1 2 x 2 =4 [ 2 1 2 1 1 2 ][ x 1 x 2 ]=[ 3 5 4 ]

>> A=[2 1; 2 -1; 1 -2];
>> b=[3;5;4];
>> Ab=[A b];
>> rank(A)
ans =     2
>> rank(Ab)
ans =     2

El sistema tiene solución única ya que el rango es igual al número de incógnitas (segundo caso).

>> X=A\b
X =    2.0000
   -1.0000

En la figura vemos la representación gráfica de las rectas: y=3-2x, y=2x-5, y=x/2-2; que como vemos se cortan en un el punto x=2, y=-1.

line([-5,5],[13,-7])
line([-5,5],[-15,5])
line([-5,5],[-9/2,1/2])
xlabel('x')
ylabel('y')
title('Solución única')

División por la izquierda \

Vamos a practicar la división por la izquierda \ resolviendo el sistema de tres ecuaciones con tres incógnitas:

3 x 1 x 2 = 5 2 x 1 + x 2 + x 3 = 0 2 x 1 x 2 + 4 x 3 = 15 [ 3 1 0 2 1 1 2 1 4 ] [ x 1 x 2 x 3 ] = [ 5 0 15 ] A X = B X = A \ B

>> A=[3 -1 0; -2 1 1; 2 -1 4];
>> b=[5;0;15];
>> Ab=[A b];
>> rank(A)
ans =     3
>> rank(Ab)
ans =     3
>> A\b
ans =
    2.0000
    1.0000
    3.0000
>> inv(A)*b
ans =
    2.0000
    1.0000
    3.0000

La función rref

Esta función nos puede ser útil para resolver sistemas de ecuaciones.

Ab=( 3 1 0 2 1 1 2 1 4 5 0 15 )

>> A=[3 -1 0; -2 1 1; 2 -1 4];
>> b=[5;0;15];
>> Ab=[A b];
>> R=rref(Ab)
R =
     1     0     0     2
     0     1     0     1
     0     0     1     3

Pasamos la matriz Ab a la función rref de MATLAB y obtenemos la matriz ampliada R del sistema equivalente de ecuaciones

R = ( 1 0 0 0 1 0 0 0 1 2 1 3 ) { x 1 = 2 x 2 = 1 x 3 = 3

{ x 1 + x 2 + x 3 =6 x 1 2 x 3 =4 2 x 1 + x 2 x 3 =18

Creamos la matriz A de los coeficientes y la matriz ampliada Ab

>> A=[1 1 1;1 0 -2;2 1 -1];
>> b=[-6;4;18];
>> Ab=[A b];
>> R=rref(Ab)
R =
     1     0    -2     0
     0     1     3     0
     0     0     0     1

El sistema equivalente es

{ x 1 2 x 3 =0 x 2 +3 x 3 =0 0=1

Se trata de un sistema incompatible, tal como podemos comprobar calculando el rango de la matriz A y de su ampliada Ab

>> rank(A)
ans =     2
>> rank(Ab)
ans =     3

{ x 1 + x 2 + x 3 =6 x 1 2 x 3 =4 2 x 1 + x 2 x 3 =10

Creamos la matriz A de los coeficientes y la matriz ampliada Ab

>> A=[1 1 1; 1 0 -2;2 1 -1];
>> b=[6;4;10];
>> Ab=[A b];
>> rref(Ab)
ans =
     1     0    -2     4
     0     1     3     2
     0     0     0     0

La última fila 0x1+0x2+0x3=0 se verifica para cualquier valor de x1, x2 y x3

El sistema equivalente es

{ x 1 2 x 3 =4 x 2 +3 x 3 =2 { x 1 =4+2 x 3 x 2 =23 x 3

Una solución sería por ejemplo, x1=4, x2=2 y x3=0.

>> rank(A)
ans =     2
>> rank(Ab)
ans =     2

El sistema es compatible e indeterminado.