Читать книгу Métodos numéricos aplicados a Ingeniería - Héctor Jorquera González - Страница 9

	1.Sistemas de ecuaciones lineales

1. SISTEMAS DE ECUACIONES LINEALES

La resolución numérica de sistemas de ecuaciones lineales se puede dividir en los siguientes dos métodos:

a) Método directo: Si el número de etapas para resolver las ecuaciones es finito. Por ejemplo, la resolución del sistema de ecuaciones lineales A · x = b cuando existe la inversa de la matriz A.

b) Método iterativo: Si se requiere infinitas etapas para resolver las ecuaciones en forma exacta. Por ejemplo, resolver un sistema de ecuaciones lineales mediante el método de Jacobi.

1.1 Métodos de solución directa

Un método directo es un algoritmo con un número finito y predefinido de pasos, al final de los cuales se obtiene la solución.

1.1.1 Eliminación de Gauss-Jordan

Este algoritmo consiste en producir una serie de transformaciones del sistema lineal original, hasta obtener un sistema triangular superior. Supongamos que nuestro punto de partida consiste en el sistema de ecuaciones lineales

Se genera una secuencia de operaciones hasta que se transforma el sistema original de ecuaciones a un sistema con matriz triangular superior U:

La solución de esta ecuación es la misma que la de la ecuación 1.1.

Para conseguir la estructura triangular superior se procede a eliminar los elementos bajo la diagonal principal, haciéndolos cero a través de operaciones (sumas y restas) con las filas de la matriz A. Al final del procesamiento se obtiene la siguiente estructura de matriz:

Esto genera el sistema de ecuaciones U · x = y, el cual se puede resolver por sustitución hacia atrás para obtener la solución x.

Ejemplo 1.1. Eliminación de Gauss-Jordan

Utilice Matlab^® para realizar la eliminación de Gauss-Jordan del siguiente sistema de ecuaciones A · x = b,

Solución: Se hace el pivoteo de la matriz A en conjunto con el vector b, de manera de producir directamente el resultado mostrado en la ecuación 1.2; para esto se opera con una matriz aumentada M = [A b]. El siguiente macro de Matlab^® resuelve el problema propuesto:

% este macro resuelve el Ejemplo 1.1

%% ingreso de la matriz y vector lado derecho

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

b=[-1 12 0]’;

x=zeros(size(b));

%% definición matriz aumentada M y eliminación Gauss-Jordan

M=[A b];

M(1,:)=M(1,:)/M(1,1); % se normaliza fila 1 con primer pivote

M(2,:)=M(2,:)-M(1,:)*M(2,1); % se genera un ‘0’ en posición (2,1)

M(3,:)=M(3,:)-M(1,:)*M(3,1); % se genera un ‘0’ en posición (3,1)

M(2,:)=M(2,:)/M(2,2); % se normaliza fila 2 con segundo pivote

M(3,:)=M(3,:)-M(2,:)*M(3,2); % se genera un ‘0’ en posición (3,2)

%% etapa de solución de la ecuación por sustitución hacia atrás

x(3)=M(3,4)/M(3,3);

x(2)=M(2,4)-M(2,3)*x(3);

x(1)=M(1,4)-M(1,3)*x(3)-M(1,2)*x(2);

%% verificación de la solución

r=A*x-b;

El resultado de aplicar este macro es:

Notas

1) El método como aquí se ha descrito (eliminación de Gauss-Jordan) falla si cualquiera de los pivotes (elementos de la diagonal de la matriz A) se hace cero o si es muy pequeño por causa de errores de redondeo; por ejemplo, resta de cantidades de similar magnitud.

2) En la práctica es necesario usar una estrategia de pivoteo para escoger el “mejor” pivote, de manera de minimizar los problemas numéricos. El mejor pivote es usualmente el elemento con mayor valor absoluto en la columna que queda bajo la diagonal principal durante la eliminación hacia adelante.

3) En Matlab^®, el algoritmo está implementado como la factorización LU (rutina lu), y también en el operador ‘\’, que usa la eliminación gaussiana para resolver sistemas lineales. Moler [1] ha desarrollado una interfaz gráfica lugui que permite visualizar el proceso de factorización de una matriz cuadrada, para distintas opciones: sin pivote, pivoteo parcial o pivoteo total.

4) En el caso particular de matrices que sean definidas positivas o que posean diagonal dominante¹, el pivoteo siempre funciona y no hay problemas de encontrar pivotes iguales a cero.

1.1.2 Caso de matrices tridiagonales

Un ejemplo simple de aplicación de la eliminación de Gauss-Jordan es el caso de una matriz tridiagonal, de la forma

Si se cumple la condición de diagonal dominante para la matriz: |α_i| > |γ_i| + |β_i-1|, (para todo i), entonces el proceso de eliminación de Gauss-Jordan queda dado por el algoritmo de Thomas:

Y la etapa de sustitución hacia atrás queda como

Notas

1) Este algoritmo no solo evita operaciones aritméticas innecesarias con elementos nulos de la matriz, sino que asimismo reduce la cantidad de almacenamiento requerido de O(n²) a O(5·n) aproximadamente, ya que almacena la matriz A como tres vectores α, β y γ. Este algoritmo está disponible en la rutina tridisolve desarrollada por Moler [1].

2) Este enfoque se puede implementar de modo eficiente cuando la matriz original posee una estructura regular. Esta matriz tridiagonal aparece frecuentemente cuando se modelan balances de masa en procesos de separación, tales como destilación (Wang y Henke, [2]), absorción, extracción por solventes (Hanna y Sandall, [3]); igualmente aparece en problemas de contorno (ecuaciones diferenciales ordinarias o a derivadas parciales) cuando se usan métodos de diferencias finitas (capítulos 4 y 5).

Ejemplo 1.2. Solución de un sistema tridiagonal

Se remueve anilina desde una corriente acuosa usando tolueno como solvente extractor; los flujos de agua (F₁) y tolueno (F₂) son constantes (en base libre de soluto). La unidad de separación es una torre en contracorriente con 10 etapas de extracción, como se muestra en la Figura 1.1, con la siguiente notación:

x_i= kg de anilina en fase orgánica/kg de tolueno en fase orgánica.

y_i = kg de anilina en fase acuosa/kg de agua en fase acuosa.

FIGURA 1.1. Esquema de la extracción líquido-líquido del ejemplo 1.2

Los subíndices ‘i’ corresponden a las composiciones de las corrientes que salen de la etapa ‘i’ de extracción. Las ecuaciones del balance de masa en cada etapa ‘i’ se pueden poner como sigue:

Donde K es el coeficiente de reparto de soluto entre ambas fases, supuesto constante en todas las etapas. Para las etapas 1 y 10 hay que considerar:

Construya un macro que resuelva este sistema de ecuaciones y que grafique las concentraciones de anilina en ambas fases, en función del número de etapas, si se considera que K = 4 = constante.

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₁₀ = 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 n ecuaciones, tenemos que el número de operaciones está dado por:

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

• Fase de sustitución hacia atrás: n²/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 A es 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 x para 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 A produce 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 P representa 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 M y el vector g son 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 M sea 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 = b mediante la descomposición de la matriz original A. En particular, podemos descomponer la matriz A de 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.

1.2.1 Método de Jacobi (Desplazamientos simultáneos)

En este método se aplica la descomposición 1.11 y se deja la diagonal de la matriz original en el lado izquierdo de la ecuación original, de modo que

Por lo tanto, en cada iteración se calculan los nuevos componentes del vector de solución, usando la información de los componentes de la etapa anterior, por lo que en la ecuación 1.12 el componente ‘i’ del vector solución actualizado queda como:

1.2.2 Método de Gauss-Seidel (Desplazamientos sucesivos)

En esta variante se pasa solo la diagonal superior de la matriz A al lado derecho de la ecuación original, obteniéndose:

En este método, la componente i-ésima del vector solución en la etapa (k+1) se actualiza usando la información de los componentes del vector x^(k) y los nuevos componentes del vector x^(k+1) hasta el componente i-1:

Es decir, se emplea información más actualizada que en el caso del método de Jacobi.

1.2.3 Método de relajaciones sucesivas

La mayor velocidad de convergencia del método de Gauss-Seidel (con respecto al método de Jacobi) se debe a que actualiza toda la información disponible del vector x^(k+1) para calcular los restantes componentes de dicho vector, y esto se efectúa al incluir el componente L de la matriz original A en el lado izquierdo de la ecuación 1.14.

En el método de relajaciones sucesivas se considera la posibilidad de acelerar la convergencia de la iteración de Gauss-Seidel usando un parámetro de relajación ω, de manera que se hace una actualización de la iteración de la forma

El parámetro ω es mayor que 0 y menor que 2. Si ω es mayor que 1.0, hablamos de sobrerrelajación, mientras que si 0 < ω < 1, el método se denomina subrelajación. Existen técnicas para estimar el valor óptimo del parámetro ω para una situación dada (Axelson, [4]), incluso modificando su valor a lo largo de la iteración (esquemas adaptivos), que básicamente consisten en minimizar el radio espectral (ρ) de la matriz M = (D/ω+L)^-1·((1/ω−1)·D-U) que aparece en la ecuación anterior.

Notas

1) Los métodos de Jacobi y de Gauss-Seidel son útiles solamente cuando la matriz es diagonal dominante, lo cual garantiza la convergencia de ambos métodos.

2) Usualmente se tiene, en orden creciente de rapidez de convergencia que: Jacobi < Gauss-Seidel < Relajaciones, ya que en este último caso se puede escoger ω para minimizar ρ(M).

3) Se puede demostrar que si 0 < ω < 2 y A es definida positiva, el método de relajaciones sucesivas converge (Sewell, [7]).

4) Se puede demostrar que si 0 < ω < 1 y A es diagonal dominante, el método de relajaciones sucesivas converge (Sewell, [7]).

Ejemplo 1.4. Comparación de métodos iterativos

Para los siguientes sistemas de ecuaciones lineales, calcule el número de iteraciones necesarias para conseguir un valor menor a 10^-5 en la norma del residuo r^(k) = A·x^(k)-b usando los siguientes métodos: i) Jacobi; ii) Gauss-Seidel; iii) relajaciones aplicadas a Gauss-Seidel, escogiendo ω de forma de minimizar ρ(M) en el método de Gauss-Seidel.

Aplicando los distintos esquemas iterativos se obtienen los siguientes resultados:

TABLA 1.1. Resultados del ejemplo 1.4

* Calculado usando el valor óptimo de ω que minimiza ρ (M) en 1.16.

Nótese que el caso b) no converge porque la matriz no es diagonal dominante ni tampoco se puede reordenar las filas para que lo sea. Para los casos a) y c), ambas matrices poseen diagonal dominante, pero solo c) es definida positiva (todos sus valores propios son estrictamente positivos), de ahí que el método de relajaciones alcance su mayor eficacia de implementación.

Por lo tanto, una conclusión relevante es que todos los métodos iterativos de la forma general 1.9 deben pasar por un chequeo previo de que las condiciones de convergencia se cumplen, antes de intentar resolverlos por cualquier método iterativo.

A continuación presentamos un listado con el código Matlab^® de resolución iterativa mediante el método de Jacobi. El lector puede modificar fácilmente el código para implementar el método de Gauss-Seidel o el de relajaciones sucesivas.

% función que resuelve A*x=b según método de Jacobi

function [X,it,err]=Jacobi(A,b,tol,maxit)

if nargin<4,

maxit=100;

end

if nargin <3,

tol=1e-5;

end

D=diag(diag(A));

L=tril(A,-1);

U=triu(A,1);

M=-inv(D)*(L+U);

% chequeo que los valores propios cumplen la condición de convergencia

l=max(abs(eig(M)));

if l>= 1.0,

error(‘Matriz de iteración no garantiza convergencia’);

end

g=D\b;

X=zeros(size(b));

it=0;

err=norm(b);

while err > tol && it < maxit,

X=M*X+g;

err=norm(A*X-b);

it=it+1;

end

1.2.5 Estimación del error en métodos iterativos

Sea x* la solución de la ecuación genérica 1.9; se puede demostrar que una cota del error después de la k-ésima etapa de iteración viene dado por:

Por lo tanto, para que el método presente convergencia, es necesario que la norma de la matriz M sea menor que 1.0. Esto es equivalente a exigir que el radio espectral de la matriz M, definido en la ecuación 1.10, sea menor que 1.0.

No obstante lo anterior, si la norma de la matriz M es mayor que 1.0, es posible resolver la ecuación 1.7 de manera indirecta usando los residuos r^(j) de dicha iteración, a través de métodos del gradiente conjugado (Axelsson, [4]; Sewell, [7]), los cuales son más complejos de implementar y cuyo desempeño depende fuertemente de la estructura del espectro de valores propios de la matriz M, por lo que no es fácil predecir su desempeño.

Cuando tengamos que resolver en forma iterativa un caso como el recién descrito, habría que probar con los métodos iterativos que trae Matlab^® y que se describen en la siguiente sección.

1.2.6 Métodos iterativos implementados en Matlab^®

Los siguientes métodos iterativos se encuentran implementados en el ambiente Matlab^®:

• gmres: rutina que utiliza el método de minimización generalizada de residuos.

• cgs: usa el método de gradientes conjugados cuadrados.

• bicg: emplea el método de gradiente biconjugado.

• pcg: utiliza el método de gradientes conjugados precondicionados.

• qmr: maneja el método de cuasi minimización de residuos.

• bicgstab: aplica el método de gradiente biconjugado estabilizado.

Todas estas rutinas permiten reiniciar las iteraciones, incorporar criterios de tolerancia, número máximo de iteraciones, etcétera. Para mayores detalles de cada método, consultar por ejemplo Axelsson [4] y la ayuda en Matlab^®.

1.3 Análisis del error

Consideremos la posibilidad de que haya errores en los datos de entrada del problema y, por lo tanto, necesitamos estimar la magnitud del error en la solución de las ecuaciones. Para estos fines, la matriz inversa A^-1 es fundamental en la determinación del error de la solución.

Cuando los pequeños errores en la matriz de coeficientes A producen grandes cambios en la solución, entonces se dice que el sistema de ecuaciones está mal condicionado. Esto normalmente se debe a que las filas de la matriz A son casi linealmente dependientes entre sí.

Consideremos el caso cuando haya error en el lado derecho de las ecuaciones, es decir:

y cuya solución es x+δx. Luego se tiene que

y además se cumple

Por lo tanto, de las ecuaciones anteriores se deduce que una cota superior para el error relativo viene dada por (Moler, [1])

En donde κ(A) es el número de condición² de la matriz A. Si κ >> 1, se dice que la matriz A está mal condicionada. Por ejemplo, si κ(A) es 10⁸, un error relativo de 10^-8 en el lado derecho, ¡podría producir un error tan grande como 100% en el error relativo de la solución!

Un punto importante que conviene destacar aquí es que la característica de mal condicionamiento de un matriz a menudo proviene del mismo usuario que la genera en forma inadvertida. Una manera fácil de producir matrices mal condicionadas consiste en tener ecuaciones con coeficientes con órdenes de magnitud distintos entre una ecuación y otra. No existen procedimientos numéricos que nos garanticen que una matriz no presente este problema.

Una solución que se puede plantear cuando el sistema es mal condicionado, es aplicar alguna de las rutinas iterativas que se han mencionado en la sección 1.2, ya que así no se utiliza el procedimiento de inversión de la matriz (que es el que genera normalmente los problemas de mal condicionamiento y ambigüedad de los resultados), sino un método iterativo que no trabaja con la inversa de la matriz original.

1.4 Problemas propuestos

1) Considere las matrices:

Ingrese las matrices al ambiente Matlab^® e indique los resultados al aplicar las siguientes operaciones:

i) x = A(3,:)

ii) y = B(:,1)

iii) C = A(2:4;1:2)

iv) D = diag(A)

v) E = B’

vi) C = diag(diag(A))

vii) A(3,:)*B(:,2)

viii) A(2,3)*B(4,1)

ix) eig(A)

x) [U, R] = eig(A)

xi) A+4*ones(4,4)-eye(4)

2) Considere los objetos:

Encuentre lo que se obtendría al efectuar las siguientes operaciones en Matlab^®:

i) A*B

ii) C = diag(diag(A))

iii) A(1:2,3:4)*B(3:4)

iv) A(3,:)*B

v) A(:,2)’*B

vi) A(2,3)*B(1)

vii) B’* diag(A)

viii) [U, R] = eig(A*A’)

ix) A-5*ones(4,4)+3*eye(4)

3) Para las matrices A de los dos problemas anteriores, descompóngalas de manera que se puedan expresar de la siguiente forma: A =D + L + U, donde D es la diagonal principal de A, L una matriz triangular inferior y U una matriz triangular superior. Utilice los comandos diag, tril y triu de Matlab^®.

4) ¿Cuál de las siguientes matrices es singular?

Si la matriz no es singular, encuentre su inversa.

5) Encuentre la pendiente e intercepto de la recta de mínimos cuadrados que pasa por los puntos: (1, 1), (2, 1.5) y (3, 2.2). Grafique sus resultados usando plot.

6) ¿Tienen solución los siguientes sistemas de ecuaciones? Encuéntrela si ella existe.

7) Resuelva las siguientes ecuaciones lineales usando eliminación de Gauss-Jordan con diferentes opciones de pivoteo.

8) Emplee el método de Jacobi y el de Gauss-Seidel para resolver el sistema de ecuaciones del problema anterior, caso b); en ambos casos use como punto de partida x⁽⁰⁾ = [1,1,0.4]^T.

9) Considere el sistema lineal de ecuaciones A·x = b con

a) Resuelva el sistema anterior en Matlab^® usando: x = A\b, operador que realiza la eliminación gaussiana. ¿Se obtiene el mismo resultado que al aplicar: x=inv(A)*b?

b) Defina un nuevo vector del lado derecho como: b₁ = [30.00002; 42.00006; 22.00004] y calcule la nueva solución x₁ utilizando ‘\’ y la función inv. ¿Hay diferencias en los resultados? ¿Esperaba usted que hubiera diferencias en la solución numérica x?

c) Compare los errores relativos al modificar el lado derecho³: norm(b₁-b)/norm(b) y en la solución: norm(x₁-x)/norm(x). ¿Qué pasa?

d) Ahora aplique las rutinas iterativas pcg, bicg, gmres, qmr, bicgstab en su forma más simple, es decir, utilizando como argumentos solamente la matriz A y el vector b originales. Calcule x y x₁ para los casos en que se ingresa b y b₁ como vectores del lado derecho. ¿Qué resultados se obtienen? ¿Cómo se comparan estos resultados con los de las partes a) y b) anteriores? Nota: calcule los residuos de la solución x (o x₁) en cada caso.

e) Finalmente, concluya qué métodos de solución son más adecuados para abordar la resolución de este tipo de casos con matrices mal condicionadas.

10) Considere ahora el sistema lineal de ecuaciones A·x = b con

Repita el procedimiento del problema anterior, ahora usando b₁ = [ 5.00004 ; 10.0003], etc.

11) El siguiente sistema tridiagonal de ecuaciones lineales representa las ecuaciones de un modelo simple de un proceso de extracción líquido-líquido, donde la composición x(i) de un soluto de interés cambia a lo largo de las n etapas de extracción (por ejemplo, el sulfato de cobre se extrae de una solución acuosa pobre en cobre y se concentra en una corriente de solvente orgánico).

Tome n = 15 etapas y calcule el vector de concentraciones x. Grafíquelo usando plot(x).

Nota:

Para construir la matriz y el lado derecho del sistema de ecuaciones, se recomienda usar los comandos: zeros(n,1), ones(n,1) y spdiags. Emplee la función spy(M) para verificar que la matriz M posee una estructura tridiagonal.

1.5 Referencias

1. Moler, C. (2008). Numerical Computing with MATLAB. Disponible en la página web: http://www.mathworks.com/moler/; incluye una biblioteca de programas escritos en Matlab^®.

2. Wang, J.C. & Henke, G.E. (1966). Tridiagonal Matrix for Distillation. Hydrocarbon Process. Vol 45, N° 8, 155-163.

3. Hanna, O.W. & Sandall, O.C. (1995). Computational Methods in Chemical Engineering. Nueva Jersey: Prentice Hall.

4. Axelson, O. (1994). Iterative Solution Methods. Nueva York: Cambridge University Press.

5. Nakamura, S. (1992). Métodos numéricos aplicados con software México DF: Prentice Hall Hispanoamericana.

6. Gerald, C.F. Análisis numérico (2ª ed.). México DF: Alfaomega.

7. Sewell, G. (1988). The Numerical Solution of Ordinary and Partial Diferential Equations. Boston: Academic Press.