Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Sistemas de ecuaciones lineales

Consideraciones computacionales

Uno de los problemas más importantes del cálculo técnico es la solución de sistemas de ecuaciones lineales simultáneas.

En notación matricial, el problema general se plasma de la siguiente manera: Dadas dos matrices A y b, ¿existe una matriz única x de manera tal que Axb o xAb?

Resulta didáctico considerar un ejemplo de 1 por 1. Por ejemplo, la ecuación

7x = 21

¿tiene una solución única?

Naturalmente, la respuesta es sí. La solución única de la ecuación es x = 3. La solución se obtiene fácilmente a través de la división:

x = 21/7 = 3.

La solución no se obtiene comúnmente mediante el cálculo del inverso de 7, es decir 7–1= 0.142857..., y luego multiplicando 7–1 por 21. Este proceso sería más arduo, y si 7–1 se representa como un número finito de dígitos, también sería menos preciso. Pueden aplicarse consideraciones muy semejantes a los conjuntos de ecuaciones lineales con más de una incógnita. MATLAB® resuelve dichas ecuaciones sin calcular inversas de matrices.

Aunque no es una notación matemática estándar, MATLAB utiliza la terminología de la división, común en el caso de escalares, para describir la solución de un sistema general de ecuaciones simultáneas. Los dos símbolos de división, barra diagonal, /, y barra invertida, \, corresponden a las dos funciones mrdivide y mldivide de MATLAB. Estos operadores se utilizan para las dos situaciones en las que la matriz desconocida aparece a la izquierda o a la derecha de la matriz de coeficientes:

x = b/A

Indica la solución a la ecuación de matrices xA = b, que se obtiene usando mrdivide.

x = A\b

Indica la solución a la ecuación de matrices Ax = b, que se obtiene usando mldivide.

Considere “dividir” ambos lados de la ecuación, Ax = b o xA = b por A. La matriz de coeficientes A siempre está en el “denominador”.

Las condiciones de compatibilidad de dimensiones para x = A\b requieren que las matrices A y b tengan la misma cantidad de filas. Entonces, la solución x tiene la misma cantidad de columnas que b y su dimensión de filas es igual a la dimensión de columnas de A. Para x = b/A, los roles de las filas y las columnas se intercambian.

En la práctica, las ecuaciones lineales del tipo Ax = b se dan con mayor frecuencia que aquellas del tipo xA = b. En consecuencia, la barra invertida se usa con mucha mayor frecuencia que la barra diagonal. El resto de esta sección se concentra en el operador barra invertida; las propiedades correspondientes del operador barra diagonal se pueden inferir de la identidad:

(b/A)' = (A'\b').

No es necesario que la matriz de coeficientes A sea cuadrada. Si A tiene un tamaño de m por n, entonces, existen tres casos:

m = n

Sistema cuadrado. Busque una solución exacta.

m > n

Sistema sobredeterminado, con más ecuaciones que incógnitas. Encuentre una solución de mínimos cuadrados.

m < n

Sistema subdeterminado, con menos ecuaciones que incógnitas. Encuentre una solución básica con un máximo de m componentes distintos de cero.

El algoritmo de mldivide

El operador mldivide emplea diferentes mecanismos de solución para tratar distintos tipos de matrices de coeficientes. Los diversos casos se diagnostican de forma automática mediante el análisis de la matriz de coeficientes. Para obtener más información, consulte la sección “Algoritmos” de la página de referencia de mldivide.

Solución general

La solución general para un sistema de ecuaciones lineales Axb describe todas las soluciones posibles. Se puede encontrar la solución general mediante:

  1. La solución del sistema homogéneo correspondiente Ax = 0. Para hacerlo, use el comando null, escribiendo null(A). Esto devuelve una base para el espacio de solución de Ax = 0. Cualquier solución es una combinación lineal de vectores de base.

  2. La búsqueda de una solución particular para el sistema no homogéneo Ax =b.

Luego, se puede escribir cualquier solución de Axb como la suma de la solución particular para Ax =b, del paso 2, más una combinación lineal de los vectores de base del paso 1.

En el resto de esta sección, se describe cómo usar MATLAB para buscar una solución particular para Ax =b, como en el paso 2.

Sistemas cuadrados

La situación más común incluye una matriz de coeficientes cuadrada A y un único vector columna b a la derecha.

Matriz de coeficientes no singular

Si la matriz A no es singular, entonces la solución, x = A\b, tiene el mismo tamaño que b. Por ejemplo:

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

Puede confirmarse que A*x es exactamente igual a u.

Si A y b son cuadradas y tienen el mismo tamaño, x= A\b también tiene ese tamaño:

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

Puede confirmarse que A*x es exactamente igual a b.

Estos dos ejemplos tienen soluciones exactas y de números enteros. Esto se debe a que se escogió como matriz de coeficientes pascal(3), que es una matriz de rango completo (no singular).

Matriz de coeficientes singular

Una matriz cuadrada A es singular si no tiene columnas linealmente independientes. Si A es singular, la solución para Ax = b o bien no existe o no es única. El operador barra invertida, A\b, emite una advertencia si A es casi singular o si detecta una singularidad exacta.

Si A es singular y Ax = b tiene solución, se puede encontrar una solución particular que no es única. Para ello, escriba

P = pinv(A)*b

pinv(A) es una pseudoinversa de A. Si Ax = b no tiene solución exacta, entonces pinv(A) devuelve una solución de mínimos cuadrados.

Por ejemplo:

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

es singular, como puede verificar si escribe

rank(A)

ans =

     2

Debido a que A no es de rango completo, tiene algunos valores singulares iguales a cero.

Soluciones exactas Para b =[5;2;12], la ecuación Ax = b tiene una solución exacta, dada por:

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Para comprobar que pinv(A)*b es una solución exacta, escriba

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Soluciones de mínimos cuadrados. Sin embargo, si b = [3;6;0], Ax = b no tiene solución exacta. En este caso, pinv(A)*b devuelve una solución de mínimos cuadrados. Si escribe

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

no recupera el vector original b.

Puede determinar si Ax =b tiene una solución exacta buscando la forma escalonada reducida por filas de la matriz aumentada [A b]. Para hacerlo en este ejemplo, introduzca

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Dado que la fila inferior contiene solo ceros, con excepción de la última entrada, la ecuación no tiene solución. En este caso, pinv(A) devuelve una solución de mínimos cuadrados.

Sistemas sobredeterminados

Este ejemplo muestra cómo con frecuencia se encuentran sistemas sobredeterminados en diversos tipos de ajustes de curvas para datos experimentales.

Una cantidad y se mide en varios valores de tiempo t diferentes para producir las siguientes observaciones. Para introducir los datos y verlos en una tabla, utilice las siguientes instrucciones.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)
B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Intente modelar los datos con una función exponencial decreciente

y(t)=c1+c2e-t.

La ecuación anterior indica que el vector y se debe aproximar mediante una combinación lineal de otros dos vectores. Uno es un vector constante que contiene solo unos, y el otro es el vector con componentes exp(-t). Los coeficientes desconocidos, c1 y c2, se pueden calcular haciendo un ajuste de mínimos cuadrados que minimice la suma de los cuadrados de las desviaciones de los datos respecto del modelo. Existen seis ecuaciones en dos incógnitas, representadas mediante una matriz de 6 por 2.

E = [ones(size(t)) exp(-t)]
E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

Use el operador barra invertida para obtener la solución de mínimos cuadrados.

c = E\y
c = 2×1

    0.4760
    0.3413

En otras palabras, el ajuste de mínimos cuadrados para los datos es el siguiente:

y(t)=0.4760+0.3413e-t.

Las siguientes instrucciones evalúan el modelo en incrementos uniformemente espaciados en t, y luego representan gráficamente el resultado junto con los datos originales:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*c no es exactamente igual a y, pero la diferencia bien podría ser menor que los errores de medición en los datos originales.

Una matriz rectangular A es de rango deficiente si no tiene columnas linealmente independientes. Si A es de rango deficiente, entonces, no hay una solución única de mínimos cuadrados para AX = B. A\B emite una advertencia si A es de rango deficiente y produce una solución de mínimos cuadrados. Puede utilizar lsqminnorm para encontrar la solución X que tiene la norma mínima entre todas las soluciones.

Sistemas subdeterminados

En este ejemplo se demuestra cómo la solución para los sistemas subdeterminados no es única. Los sistemas lineales subdeterminados incluyen más incógnitas que ecuaciones. En MATLAB, la operación de división matricial izquierda encuentra una solución de mínimos cuadrados básica que, como máximo, tiene m componentes distintos de cero para una matriz de coeficientes de m por n.

El siguiente es un breve ejemplo al azar:

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)
R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

El sistema lineal Rp = b incluye dos ecuaciones en cuatro incógnitas. Como la matriz de coeficientes contiene enteros pequeños, es apropiado usar el comando format para visualizar la solución en formato racional. La solución particular se obtiene de la siguiente manera:

format rat
p = R\b
p =

       0       
      17/7     
       0       
     -29/7    

Uno de los componentes distintos de cero es p(2), debido a que R(:,2) es la columna de R con la norma de mayor tamaño. El otro componente que no es cero es p(4), ya que R(:,4) domina tras la eliminación de R(:,2).

La solución general completa al sistema subdeterminado se puede caracterizar mediante la adición de p a una combinación lineal arbitraria de los vectores del espacio nulo, que se pueden encontrar usando la función null con una opción que solicita una base racional.

Z = null(R,'r')
Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

Se puede confirmar que R*Z es igual a cero y que el residuo R*x - b resulta pequeño para cualquier vector x, donde

x = p + Z*q

Dado que las columnas de Z son los vectores del espacio nulo, el producto Z*q es una combinación lineal de dichos vectores:

Zq=(x1x2)(uw)=ux1+wx2.

Para ilustrarlo, escoja un q arbitrario y construya x.

q = [-2; 1];
x = p + Z*q;

Calcule la norma del residuo.

format short
norm(R*x - b)
ans =

   2.6645e-15

Cuando hay infinitas soluciones disponibles, la solución con la norma mínima resulta de especial interés. Puede utilizar lsqminnorm para calcular la solución de mínimos cuadrados de norma mínima. Esta solución tiene el menor valor posible para norm(p).

p = lsqminnorm(R,b)
p =

    -207/137   
     365/137   
      79/137   
    -424/137  

Resolución de varios lados derechos

Algunos problemas se relacionan con la resolución de sistemas lineales que tienen la misma matriz de coeficientes A, pero diferentes lados derechos b. Cuando los diferentes valores de b están disponibles al mismo tiempo, puede construir b como matriz con varias columnas y resolver a la vez todos los sistemas de ecuaciones mediante un único comando de barra invertida: X = A\[b1 b2 b3 …].

No obstante, a veces, todos los valores de b no están disponibles al mismo tiempo, lo que significa que tendrá que resolver varios sistemas de ecuaciones de forma consecutiva. Cuando resuelve uno de estos sistemas de ecuaciones mediante la barra diagonal (/) o la barra invertida (\), el operador factoriza la matriz de coeficientes A y utiliza esta descomposición de matrices para calcular la solución. Sin embargo, a partir de entonces, cada vez que resuelva un sistema de ecuaciones similar con una b diferente, el operador calculará la misma descomposición de A, lo cual es un cálculo redundante.

La solución para este problema es realizar un cálculo previo con la descomposición de A y, a continuación, reutilizar los factores para resolver los diferentes valores de b. Sin embargo, en la práctica, puede que sea difícil realizar el cálculo previo de la descomposición de esta manera, ya que se debe saber qué descomposición calcular (LU, LDL, Cholesky, etc.) y cómo multiplicar los factores para resolver el problema. Por ejemplo, con la descomposición LU, se deben resolver dos sistemas lineales para resolver el sistema original Ax = b:

[L,U] = lu(A);
x = U \ (L \ b);

En cambio, para resolver sistemas lineales con varios lados derechos consecutivos, lo más recomendado es usar objetos decomposition. Estos objetos le permiten aprovechar las ventajas de rendimiento que ofrece el cálculo previo de la descomposición de matrices, pero no es necesario saber cómo usar los factores de matrices. Puede reemplazar la descomposición de LU anterior de la siguiente manera:

dA = decomposition(A,'lu');
x = dA\b;

Si no sabe con certeza qué descomposición utilizar, decomposition(A) escoge el tipo correcto en función de las propiedades de A, algo similar a lo que hace la barra invertida.

Aquí le mostramos una prueba simple sobre las posibles ventajas de rendimiento de este enfoque. La prueba resuelve el mismo sistema lineal disperso 100 veces usando la barra invertida (\) y la función decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc
Elapsed time is 9.006156 seconds.
% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc
Elapsed time is 0.374347 seconds.

Para este problema, la solución decomposition es mucho más rápida que la ofrecida por el uso exclusivo de la barra invertida. Sin embargo, la sintaxis sigue siendo simple.

Métodos iterativos

Si la matriz de coeficientes A es grande y dispersa, generalmente, los métodos de factorización no son eficientes. Los métodos iterativos generan una serie de soluciones aproximadas. MATLAB ofrece varios métodos iterativos para manipular matrices de entrada grandes y dispersas.

FunciónDescripción
pcg

Método del gradiente conjugado precondicionado. Este método es adecuado para la matriz de coeficientes definida positiva hermítica A.

bicg

Método del gradiente biconjugado

bicgstab

Método del gradiente biconjugado estabilizado

bicgstabl

Método BiCGStab(l), variante del método del gradiente biconjugado estabilizado

cgs

Método del gradiente conjugado cuadrado

gmres

Método del residuo mínimo generalizado

lsqr

Método LSQR

minres

Método del residuo mínimo. Este método es adecuado para la matriz de coeficientes hermítica A.

qmr

Método del residuo cuasi-mínimo

symmlq

Método LQ simétrico

tfqmr

Método QMR sin trasposición

Multiproceso computacional

MATLAB admite multiprocesos computacionales para varias funciones de álgebra lineal y funciones numéricas que actúan elemento por elemento. Estas funciones se ejecutan automáticamente en varios procesos. Para que una función o expresión se ejecute más rápido en varias CPU, se deben cumplir una serie de condiciones:

  1. La función debe realizar operaciones que se dividan fácilmente en secciones ejecutadas de manera simultánea. Estas secciones se deben poder ejecutar con poca comunicación entre los procesos. Además, estas deben requerir pocas operaciones secuenciales.

  2. El tamaño de los datos debe ser lo suficientemente grande para que cualquier ventaja propia de la ejecución simultánea justifique el tiempo que se requiere para dividir los datos y administrar procesos de ejecución separados. Por ejemplo, la mayoría de las funciones se acelera solo cuando el arreglo contiene varios miles de elementos o más.

  3. La operación no debe estar limitada por la memoria; el tiempo de procesamiento no debe estar dominado por el tiempo de acceso de la memoria. Como regla general, las funciones complicadas se aceleran más que las funciones simples.

Las funciones inv, lscov, linsolve y mldivide muestran un significativo aumento de la velocidad en arreglos de doble precisión de gran tamaño (del orden de diez mil elementos o más) cuando está activado el multiproceso.

Consulte también

| | | |

Temas relacionados