Héctor Jorquera González - Métodos numéricos aplicados a Ingeniería

El macro siguiente en Matlab ®permite resolver el sistema de ecuaciones 1.7 y 1.8 usando la rutina tridisolve de Moler [1], y los resultados se presentan en la Figura 1.2.

% macro que resuelve el ejemplo 1.3

F1=100; F2=10; xentra=0; yentra=0.05; K=4;

% factor común que se emplea en los balances de masa

alpha=F2*K/F1;

% construcción del lado derecho del sistema de ecuaciones

b=zeros(10,1);

b(1)=-yentra;

b(end)=-F2*xentra/F1;

% construcción de las tres diagonales del sistema de ecuaciones

e=ones(10,1);

f=-(1+alpha)*e;

g=alpha*e;

% solución para las concentraciones {y_i}

y=tridisolve(e,f,g,b);

% solución para las concentraciones {x_i}

x=K*y;

% gráfico de las concentraciones de equilibrio en cada fase

eje=1:10;

plot(eje,x,’k-o’,eje,y,’k:^’);

title(‘Extractor líquido-líquido’);

xlabel(‘Número de etapa’);

ylabel({‘concentraciones de soluto en equilibrio’;’kg soluto/kg fase’});

ylim([0 0.25]);

legend(‘fase orgánica’,’fase acuosa’,’Location’,’Best’);

r=A*x-b;

Nótese que basta con poner cinco etapas de equilibrio para que ya no haya más cambios en las concentraciones del soluto que se puede extraer en la fase orgánica; la salida de la fase acuosa es con una concentración y 10= 0.03 kg/kg. Se propone al lector construir una función de Matlab ®que resuelva el caso general de N etapas de equilibrio para un extractor líquido-líquido, considerando que puede ingresar algo de soluto en el solvente (x ENTRA> 0).

FIGURA 1.2. Solución del sistema de ecuaciones 1.17 y 1.18

1.1.3 Número de operaciones requeridas

El trabajo computacional en términos de multiplicaciones y adiciones es importante como una medida de la eficiencia del algoritmo. Para el caso de la eliminación de Gauss-Jordan en un sistema de necuaciones, tenemos que el número de operaciones está dado por:

• Fase de eliminación: n 3/3 + O(n 2)

• Fase de sustitución hacia atrás: n 2/2 + O(n)

Es decir, el número de operaciones crece con el cubo de la dimensión del sistema. Si pensamos que hoy es típico en una aplicación compleja que n sea del orden de un millón de elementos, podemos apreciar que la mayor complejidad nos lleva rápidamente a costos computacionales muy elevados para este tipo de algoritmos. Retomaremos esta problemática cuando veamos los métodos iterativos.

1.1.4 Métodos directos implementados en Matlab ®

En Matlab ®, hay varias opciones para resolver los sistemas de ecuaciones lineales:

a) La alternativa más obvia es usar la función que calcula la inversa de una matriz y entonces la solución está dada por: x = inv(A)· b. Sin embargo, esta alternativa está limitada solamente a matrices cuadradas que posean inversa y estén bien condicionadas (ver sección 1.3); además, no es el método más eficiente cuando la matriz posee cierta estructura como ser triangular, simétrica, pocas bandas diagonales y el resto elementos cero, etcétera.

b) Matlab ®dispone del operador ‘\’ que corresponde a dividir por la izquierda la ecuación básica A· x = b, la que queda expresada como: x = A\b. Este operador posee mayor funcionalidad que la inversa, ya que se puede aplicar a sistemas sobre- y subdeterminados que poseen más (menos) ecuaciones que incógnitas, por lo que en tal caso se obtiene una solución en el sentido de mínimos cuadrados. El algoritmo intenta primero explorar alguna característica especial de la matriz A; por ejemplo, si es simétrica, triangular, etcétera, antes de realizar el cálculo numérico; verbigracia, si Aes rala, entonces se llama a rutinas de Matlab ®especializadas en matrices ralas.

c) La factorización LU permite resolver de manera eficiente el sistema de ecuaciones 1.1 cuando se necesita evaluar la solución xpara múltiples lados derechos b.

Sin embargo, todos estas funciones trabajan basadas en un proceso de eliminación de Gauss-Jordan, el que sabemos que se pone costoso cuando el número de ecuaciones aumenta. Luego, es necesario considerar esquemas alternativos de solución de 1.1 que no pasen por formar una inversa de la matriz A.

Ejemplo 1.3. Comparación de métodos de solución

Resuelva el sistema de ecuaciones del ejemplo 1.1 usando las funciones inv , ‘\’ o lu de Matlab ®.

El siguiente macro aplica las tres metodologías para resolver el sistema de ecuaciones; los tres métodos encuentran la misma solución, y solo el error de redondeo producido en el cálculo de la inversa de la matriz Aproduce un residuo pequeño de norma 8.9·10 -16, siendo los restantes residuos exactamente cero. En el caso de la factorización LU, la rutina entrega una matriz de permutación P, de modo que P· A= L· U; esta matriz Prepresenta los intercambios de filas realizados en la estrategia de pivoteo parcial, para prevenir ceros (o casi ceros) en la diagonal.

% este macro resuelve el ejemplo 1.2

%% ingreso de la matriz y vector lado derecho

A=[2 1 -3; -1 3 2; 3 1 -3];

b=[-1 12 0]’;

%% solución aplicando la inversa de la matriz A

x1=inv(A)*b;

r1=A*x1-b;

%% solución aplicando el operador ‘\’

x2=A\b;

r2=A*x2-b;

%% solución aplicando factorización lu

[L,U,P]=lu(A);

b3=P*b;

y3=L\b3;

x3=U\y3;

r3=A*x3-b;

1.2 Métodos iterativos

Existe una amplia variedad de métodos iterativos aplicados a los sistemas de ecuaciones lineales. A menudo estos métodos se usan en la solución de sistemas lineales con un gran número de variables (n >> 1), los que son costosos de resolver vía eliminación gaussiana. La mayoría de estos métodos son iteraciones estacionarias del tipo (Axelson, [4]):

Comúnmente la matriz My el vector gson independientes del índice de iteración k . Una condición suficiente para la convergencia de la iteración 1.9 es que la norma de la matriz Msea menor que 1. En particular una condición necesaria y suficiente para la convergencia es que

Donde ρ( M) es el radio espectral de la matriz M, y corresponde al mayor valor absoluto de los valores propios de M. En Matlab ®, lo calcularíamos como ρ= max ( abs ( eig (M))).

La forma de iteración dada por la ecuación 1.9 se obtiene de la ecuación original A· x= bmediante la descomposición de la matriz original A. En particular, podemos descomponer la matriz Ade la siguiente forma:

Hay que hacer notar que en Matlab ®podemos generar la descomposición de la siguiente forma: D= diag ( diag ( A)); L= tril ( A,-1) y U = triu ( A,1). A partir de esta descomposición es posible obtener varios esquemas que se han utilizado con frecuencia como métodos iterativos (Nakamura, [5]; Gerald, [6]), y que se describen a continuación.