Práctica 5
Resolución de sistemas lineales

Sistemas triangulares

Sistemas triangulares superiores

Resolvamos el sistema Ax=b, donde la matriz, A, se introduce como una tabla rectangular y el vector de términos independientes, b, como una lista, escribiendo los elementos entre llaves y separándolos por comas, tal  y como se hace a continuación:

Clear["Global`*"]

A={{1,2,-1,3,-4},{0,6,0,0,1},{0,0,2,2,0},{0,0,0,1,-4},{0,0,0,0,5}};
b={25,24,-15,40,61};
MatrixForm[A]

Para ello hallamos el orden del sistema y definimos un vector, x, que contendrá finalmente a la solución, y que se inicializa, por ejemplo, con ceros.

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

Hallamos la solución mediante sustitución regresiva y la mostramos:

For[i=n,i>=1,i--,
       x[[i]]=(b[[i]]-Sum[A[[i,k]] x[[k]], {k,i+1,n}])/A[[i,i]]
    ]

La solución es

x

Mathematica y la resolución de sistemas de ecuaciones lineales

Mathematica incorpora el comando LinearSolve para resolver sistemas lineales. Comprobemos la solución.

A={{1,2,-1,3,-4},{0,6,0,0,1},{0,0,2,2,0},{0,0,0,1,-4},{0,0,0,0,5}};
b={25,24,-15,40,61};
LinearSolve[A,b]

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}}; \)

L.U == A

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

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:

LinearSolve[A, b]

Métodos iterativos de resolución de sistemas lineales

En esta práctica nos ocuparemos de tres métodos iterativos de resolución de sistemas lineales: Jacobi y Gauss-Seidel.

Jacobi

Definimos la matriz de coeficientes, que no debe tener ningún elemento diagonal nulo, y el vector de términos independientes.

A={{2,-1,0},{1,6,-2},{4,3,-8}};
b={2,-4,5};

Especificamos el vector inicial., el número de iteraciones y el orden de la matriz del sistema, y producimos las aproximaciones deseadas

x={0.,0.,0.};
iteraciones=10;
n=Dimensions[A][[1]];
For[m=1,m<=iteraciones,m++,
    y=Table[0,{i,n}];
    For[i=1,i<=n,i++,
        y[[i]]=(b[[i]]-Sum[A[[i,j]]x[[j]],{j,i-1}]-
                       Sum[A[[i,j]]x[[j]],{j,i+1,n}])/A[[i,i]]
       ];
    x=y;Print[x];
   ]

La solución exacta es

LinearSolve[A,b]

Gauss-Seidel

Las siguientes órdenes proporcionan la aproximación de la única solución del sistema inicialmente considerado tras definir la matriz y el vector de términos independientes del método:

x={0.,0.,0.};
iteraciones=10;
For[m=1,m<=iteraciones,m++,
    y=Table[0,{i,n}];
    For[i=1,i<=n,i++,
        y[[i]]=(b[[i]]-Sum[A[[i,j]]y[[j]],{j,i-1}]-
                       Sum[A[[i,j]]x[[j]],{j,i+1,n}])/A[[i,i]]
       ];
    x=y;Print[x];
   ]

Ejercicios

1.- Resuelva mediante una orden directa, un sistema lineal, cuya matriz de coeficientes y vector de términos independientes son, respectivamente,  a1={{10,1,2,3},{1,-5,4,2,4},{2,3,10,4},{3,3,3,25,4}} y  b1={1,12.1,13.2,573}.

2.- 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, estudiando si admite una descomposición de Cholesky.

3.- Apoyándose en el ejercicio 2, resuelva el sistema Ax = b, donde  b = (((2, 6, 1, 9))^T).

4.- Obtenga aproximaciones por Jacobi y por Gauss-Seidel el sistema cuya matriz de coeficientes es 
                     A={{2, 2, 3, 1, 3, dm},{2, 1, 2, 2, 5, dm},{1, 1, 1, 1, 4, dm},
                            {0, 0, 1, 1, 5, dm},{-1, 0, -1, 1, -1, -3},{1, -1, 3, -4, 5, ds}}
y cuyos términos independientes son  b={5.1,-5.1,4.1,-4.1,111.1,3}, partiendo de la aproximación inicial  x={d2, d2, d2, d2, d2, d2}  y efectuando 100 iteraciones.

5.- Obtenga las soluciones de los ejercicios 3 y 4 usando la sentencia LinearSolve y también invirtiendo la matriz de coeficientes.