Resolución de sistemas lineales por métodos directos

Descomposición LU

Si la matriz de coeficientes, A, del sistema Ax = b se expresa como producto, A = L U, de una matriz triangular inferior, L, y de una superior, U, la resolución del mismo se reduce a hallar las soluciones de los sistemas Ly = b y Ux = y. Veamos un ejemplo. La matriz

A = {{1, 2, 4}, {2, 9, 7}, {-1, 8, -2}};

es el producto de las matrices

\(L = {{1, 0, 0}, {2, 1, 0}, {\(-1\), 2, 4}}; 

    U = {{1, 2, 4}, {0, 5, \(-1\)}, {0, 0, 1}}; \)

Entonces, la solución del sistema  que tiene a A como matriz de coeficientes y

b = {10, 9, -2};

como vector de términos independientes se resuelve en dos etapas: primero un sistema triangular inferior y luego uno superior.

\(n = \(Dimensions[A]\)[[1]]; 

    y = Table[0, {i, n}]; 

    For[i = 1, i ≤ n, \(i++\), 
      y[[i]] = 
        \(b[[i]] - 
            ∑\+\(k = 1\)\%\(i - 1\)L[[i, k
                  ]]\ 
                y[[k]]\)\/L
            [[i, i]]]; 

    x = Table[0, {i, n}]; 

    For[i = n, i ≥ 1, \(i--\), 
      x[[i]] = 
        \(y[[i]] - 
            ∑\+\(k = i + 1\)\%n
                 U[[i, k]]\ 
                x[[k]]\)\/U
            [[i, i]]]; 
x\)

Es fácil, y rápido, verificar que la solución es correcta:

A . x == b

Casos usuales

El problema es, evidentemente, hallar la descomposición LU de la matriz.

Factorización de Doolittle

La diagonal de la matriz L es unitaria. Veamos un ejemplo.

A = {{60, 30, 20}, {30, 20, 15}, {20, 15, 12}};

\(n = \(Dimensions[A]\)[[1]]; 

    U = Table[0, {i, n}, {j, n}]; 
L = Table[0, {i, n}, {j, n}]; 

    For[i = 1, i ≤ n, \(i++\), 
      L[[i, i]] = 1]; 

    For[j = 1, j ≤ n, \(j++\), 
      U[[1, j]] = 
        A[[1, j]]]; 

    For[i = 2, i ≤ n, \(i++\), 
      L[[i, 1]] = 
        A[[i, 1]]\/U
            [[1, 1]]]; 

    For[k = 2, k ≤ n, \(k++\), 
\ \ \ \ 
      U[[k, k]] = 
        A[[k, k]] - 
          ∑\+\(r = 1\)\%\(k - 1\)L[[k, r
                ]]\ 
              U[[r, k]]; 
      
\ \ \ \ 
      For[j = k + 1, j ≤ n, \(j++\), 
        U[[k, j]] = 
          A[[k, j]] - 
            ∑\+\(r = 1\)\%\(k - 1\)L[[k, r
                  ]]\ 
                U[[r, j]]]; 
      
\ \ \ \ \ 
      For[i = k + 1, i ≤ n, \(i++\), 
        L[[i, k]] = 
          \(A[[i, k]] - 
              ∑\+\(r = 1\)\%\(k - 1\)L[[i, r
                    ]]\ 
                  U[[r, k]]\)\/U
              [[k, k]]]
        
\ \ ]; \)

MatrixForm[L]MatrixForm[U]

Es inmediato comprobar que la factorización es correcta, pues

L . U == A

Factorización de Crout

La diagonal de la matriz U es unitaria. Se trata de forma análoga a la factorización de Crout.

Factorización de Cholesky

La matriz A que se pretende factorizar debe ser simétrica y definida positiva. La factorización de A debe ser de la forma A = L (L^T). Veamos un ejemplo.

\(A = {{2., \(-1.\), 0., 0.}, {\(-1.\), 2., \(-1.\), 0.}, 
        
\ \ {0., \(-1.\), 2., \(-1.\)}, {0., 0., \(-1.\), 
          2.}}; 
MatrixForm[A]\)

Es simétrica. Además, las submatrices diagonales principales tienen determinantes positivos:

\(n = \(Dimensions[A]\)[[1]]; 

    detprincipales = 
      Table[Det[
          Table[Table[
              A[[i, j]], {j, k}], {i, 
              k}]], 
{k, n}]\)

Por tanto, es definida positiva. Es aplicable el método de Cholesky.

\(L = Table[0, {i, n}, {j, n}]; 

    L[[1, 1]] = 
      \@A[[1, 1]]; 

    For[i = 2, i ≤ n, \(i++\), 
      L[[i, 1]] = 
        A[[i, 1]]\/L
            [[1, 1]]]; 

    For[k = 2, k ≤ n, \(k++\), 
\ \ \ \ \ 
      L[[k, k]] = 
        \@\(A[[k, k]] - 
            ∑\+\(r = 1\)\%\(k - 1\)L[[k, r
                  ]]\^2\); 
\ \ \ \ \ 
      For[i = k + 1, i ≤ n, \(i++\), 
        L[[i, k]] = 
          \(A[[i, k]] - 
              ∑\+\(r = 1\)\%\(k - 1\)L[[i, r
                    ]]\ 
                  L[[k, r]]\)\/L
              [[k, k]]]]; \)

MatrixForm[L]

Es inmediato ver que la factorización es correcta, pues

L . Transpose[L] == A

Se puede modificar el procedimiento anterior para, como es lógico, no tener necesidad de estudiar previamente si la matriz es definida positiva o no.

Aprovechamos la factorización para resolver el sistema Ax = b, donde

b = {1, 2, 3, 4};

Su solución se calcula de la siguiente manera:

\(U = Transpose[L]; 
y = Table[0, {i, n}]; 

    For[i = 1, i ≤ n, \(i++\), 
      y[[i]] = 
        \(b[[i]] - 
            ∑\+\(k = 1\)\%\(i - 1\)L[[i, k
                  ]]\ 
                y[[k]]\)\/L
            [[i, i]]]; 

    x = Table[0., {i, n}]; 

    For[i = n, i ≥ 1, \(i--\), 
      x[[i]] = 
        \(y[[i]] - 
            ∑\+\(k = i + 1\)\%n
                 U[[i, k]]\ 
                x[[k]]\)\/U
            [[i, i]]]; 
N[x, 10]\)

Podemos comprobarlo:

A . N[x] == b

Finalmente, indiquemos que Mathematica dispone de una orden para obterner la solución del sistema lineal Ax = b, que es LinearSolve[A,b].

Ejercicios

1.- Halle la descomposición LU correspondiente a Crout de la matriz

		(Subscript[A, 1])=(({{(39/8), (39/4), (117/8), -(117/4)}, {1, 3, 1, 0}, {2, 6, 1, -1}, {0, 1, -1, 10}})) .

2.- i) Apoyándose en la descomposición anterior, resuelva el sistema (Subscript[A, 1]) x = (Subscript[b, 1]),  donde

		(Subscript[b, 1]) = (((2, -6, 1, 9))^T).

     ii) Utilice, asimismo, los comandos LinearSolve e Inverse para resolver el mismo sistema.

3.- Obtenga la descomposición de Doolittle de (Subscript[A, 1]) y, apoyándose en ella, vuelva a obtener la solución del sistema (Subscript[A, 1]) x = (Subscript[b, 1]).

4.- Estudie si la matriz

		(Subscript[A, 2])=(({{(3025/64), (165/4), (55/8), 0}, {(165/4), 37, 7, 0}, {(55/8), 7, 6, 18}, {0, 0, 18, 90}}))

es definida positiva.

5.- Apoyándose en el ejercicio 4, resuelva el sistema (Subscript[A, 2]) x = (Subscript[b, 2]), donde

		(Subscript[b, 2]) = (((2, 6, 1, 9))^T).

6.- Aproveche la estructura tridiagonal y diagonalmente dominante de la matriz

		(Subscript[A, 3])=(({{40, 2, 0, 0, 0, 0}, {6, 79, 1, 0, 0, 0}, {0, 9, 117, 2, 0, 0}, {0, 0, -2, 117, 6, 0}, {0, 0, 0, -6, 117, 1}, {0, 0, 0, 0, -1, 38}}))

para hallar la correspondiente descomposición de Crout con el programa del ejercicio 1, y aplíquela para resolver el sistema (Subscript[A, 3]) x = (Subscript[b, 3]), siendo

		\(TraditionalForm\`b\_3 = \((2, \(-6\), 2.45, 39\/8, 1.365  .0)\)\^T\).