Читать книгу Métodos numéricos aplicados a Ingeniería - Héctor Jorquera González - Страница 10
Оглавление| 2.Ecuaciones no lineales |
2. ECUACIONES NO LINEALES
Para resolver ecuaciones algebraicas en general se usan métodos iterativos de solución. Un método iterativo consta de las siguientes etapas:
i) Estimar un valor inicial para la solución buscada.
ii) Una fórmula para actualizar la solución aproximada que se obtiene.
iii) Un criterio para detener el proceso de actualización (chequeo de convergencia).
Notas:
a) Es importante distinguir entre el proceso (o algoritmo) iterativo completo y la fórmula de iteración (o de actualización).
b) La verificación del término “satisfactorio” del proceso es esencial y debe anticipar todas las posibles salidas del método iterativo. Debe ser capaz de notar cuando el algoritmo falla.
c) La etapa inicial requiere a menudo incluir pruebas para validar la consistencia de los datos iniciales (concentraciones positivas, etc.).
FIGURA 2.1. Diagrama de flujo típico de un proceso de iteración
A continuación vamos a presentar los métodos más comunes para resolver ecuaciones no lineales en una variable.
2.1 Método del punto fijo
Este método solo utiliza el último valor estimado para obtener la siguiente aproximación a la solución. Por esto se denomina a veces iteración funcional. Consideremos la solución de la ecuación escalar:
Ahora calculamos una secuencia {x(k)}; k = 1, 2,..., que si es exitosa, converge a la solución x*. El usuario debe proveer la primera estimación x(1), lo cual no es un paso trivial. Las siguientes aproximaciones están dadas por:
Donde x = φ(x) posee la misma solución que f(x) = 0, es decir, x* = φ(x*).
La fórmula de iteración (o de actualización) x = φ(x) se obtiene por reordenamiento de términos o adición de términos a ambos lados de la ecuación original f(x) = 0. Este proceso de pasar de f(x) = 0 a x = φ(x) no es único, como se ve en el siguiente ejemplo.
Ejemplo 2.1. Método del punto fijo
Considere la ecuación: f(x) = x3 - 7x - 6 = 0. Algunas funciones de iteración posibles son:
φ1(x) = (x3 - 6)/7
φ2(x) = (7x + 6)/x2
φ3(x) = (2x3 + 6) /(3x2 - 7)
Como sabemos a priori que f(x) = (x+1)·(x+2)·(x-3), lo que hacemos es probar el desempeño de las fórmulas arriba detalladas cuando escogemos como valor inicial x(1) = -1.1 y x(1) = -2.2, respectivamente.
Haciendo los cálculos con las distintas funciones de iteración, para un error relativo máximo de 10-5 en la solución se obtiene la siguiente tabla de resultados:
TABLA 2.1. Comparación de resultados x*(iteraciones) para distintas funciones de actualización
A continuación se presenta el código Matlab® utilizado para generar las iteraciones de punto fijo.
function [raíz,err,nit]=punto_fijo(fun,x1,tol,maxit)
% función que ejecuta la iteración de punto fijo (o iteración funcional)
% variables de entrada:
% fun: nombre de la función de actualización
% x1: estimador inicial de la solución (o raíz)
% tol: tolerancia en el error relativo de la raíz; defecto: tol = 1e-5
% maxit: máximo número de iteraciones permitidas; defecto: maxit = 100
% variables de salida:
% raíz: solución de la ecuación original
% err: error relativo estimado para la raíz
% nit: número de iteraciones funcionales realizadas
if nargin<4,
maxit=100;
end
if nargin<3,
tol=1e-5;
end
if nargin<2,
error(‘no se puede calcular nada’);
end
xr=x1;
iter=1;
while (1)
xrold=xr;
xr=fun(xrold);
if xr~=0,
ea=abs((xrold-xr)/xr);
end
iter=iter+1;
if ea<tol || iter>maxit,
break;
end
end
raíz=xr;
err=ea;
nit=iter;
Definición: Una función φ: R → R es contracción en un conjunto D ∈ R si existe una constante 0 < L < 1 tal que:
Gráficamente, esto significa que un intervalo dado del dominio D se transforma, mediante la aplicación de φ, en un subconjunto de D. Veamos esto en el siguiente ejemplo.
Ejemplo 2.2. Convergencia del método del punto fijo
Volvamos a considerar φ1(x) = (x3-6)/7 en torno al origen x = 0. Gráficamente, empleando la función fplot, podemos apreciar que el intervalo [0,2] se transforma en el intervalo [-6/7, 2/7], el que claramente no es subconjunto de [0,2], por lo cual φ1 no es una contracción allí. Sin embargo, también es posible apreciar gráficamente que
x ∈ [-2, 0] → φ1(x) ∈ [-2, -6/7], ∴ φ1(x) es contracción en D ≡ [-2, 0]
FIGURA 2.2. Función de actualización φ1(x) para el ejemplo 2.2
% figura 2.2
f1=@(x) x;
f2=@(x) (x^3-6)/7;
fplot(f1,[-2 2],’-k’);
ylim([-2 2]);
axis equal;
hold on
fplot(f2,[-2 2],’k.’);
xlabel(‘x’);
ylabel(‘y’);
legend(‘y=x’,’y=\phi_1(x)’);
Hay que hacer notar que el concepto de contracción está íntimamente asociado al dominio de la función φ1, y este dominio se supone que contiene la solución buscada de antemano, por lo que nuestra intuición física nos puede dar alguna estimación de cuál debería ser este intervalo, de manera de poder aplicar la definición.
A continuación estableceremos el resultado matemático fundamental en que se basa este método iterativo.
2.2 Teorema de la función contractante (o del punto fijo)
Sea φ una función continua φ: R → R y sea D ∈ R tal que: x∈ D → φ (x) ∈ D. Supongamos además que existe una constante L (0 < L < 1)
tal que
Luego, para todo x(1) ∈ D, la secuencia 2.2 converge a un único punto fijo x* ∈ D.
Gráficamente, la secuencia de iteraciones 2.2 hace que la distancia entre un estimador x(k) y la solución x* sea cada vez más pequeña, hasta que tiende a cero, debido a que la constante de Lipschitz L es menor que 1.0. No obstante, para que haya convergencia, se deben cumplir todas las condiciones del teorema anterior, las que se pueden verificar fácilmente empleando la derivada de φ(x) como se ve en el siguiente ejemplo.
Ejemplo 2.3. Verificación de convergencia del método del punto fijo
Consideremos f(x) = x3 - 7x - 6 = 0 y φ1(x) = (x3 - 6)/7 para hallar la raíz x* = -1
a) Por el ejemplo previo, sabemos que φ1 es una contracción en D = [-2,0] y x* ∈ D.
b) A través del teorema del valor medio se tiene:
Por lo que tenemos que L ≡ máx.{ξ ∈ Δ} | φ1’(ξ) |
La condición de contracción se cumple si consideramos la intersección de ambos intervalos, luego la región de convergencia es: [-√7/3, 0]; como esta región contiene x* = -1, la iteración de punto fijo con φ1(x) convergerá a x* = -1, para cualquier valor inicial x(1) ∈ [-√7/3, 0].
Verificar esto utilizando la función punto_fijo presentada en el ejemplo 2.1. Adviértase que φ1(x) no es capaz de calcular la raíz x* = -2. ¿Por qué?
2.2.1 Representación gráfica de la iteración de punto fijo
Para el caso de una ecuación escalar, el comportamiento de la iteración x(k+1) = φ(x(k)) se puede representar como una trayectoria que une los gráficos de y = φ(x) e y = x, como se aprecia en la siguiente figura. Es posible notar que la magnitud y signo de la derivada φ’(x) determinan el éxito o fracaso de la iteración de punto fijo.
FIGURA 2.3. Resultados posibles para una iteración de punto fijo
2.3 Métodos de interpolación
La idea consiste en que la función f(x) = 0 se puede expandir en torno a x = x1 mediante serie de Taylor; se trunca la serie y a esta aproximación a f(x) se le busca(n) el(los) cero(s).
2.3.1 Interpolación lineal (método de Newton)
Expandimos f(x) en torno a x = x(k); llamando P(x) a la aproximación lineal:
Y la nueva aproximación x(k+1) es la raíz de P(x) = 0:
Ejemplo 2.4. Visualización de convergencia del método de Newton
El método de Newton corresponde a una secuencia de interpolaciones lineales, tal como se ilustra en la siguiente figura para el caso de f(x) = x-x2+0.5x3, comenzando con x(1) = 3.5.
FIGURA 2.4. Visualización del método de Newton
A continuación se presenta el código en Matlab® empleado para generar la figura anterior, excepto por los títulos x(k) que se añadieron con un editor de imágenes.
% macro para construir la Figura 2.4
f1=@(x) x-x^2+0.5*x^3; % función no lineal f1(x)
f1p=@(x) 1-2*x+1.5*x^2; % derivada de la función f1’(x)
fplot(f1,[-2 4],’k-’);
axis square;
xlabel(‘x’);
ylabel(‘y’);
hold on;
plot([-2 4],[0 0],’k-’)
% secuencia iterativa del método de Newton
x(1)=3.5;
f(1)=f1(x(1));
for i=2:4,
% iteración del método de Newton
x(i)=x(i-1)-f(i-1)/f1p(x(i-1));
f(i)=f1(x(i));
% secuencia de gráficos del método de Newton
xpl=[x(i-1) x(i-1)];
fpl=[0 f(i-1)];
plot(xpl,fpl,’k:’);
xpl=[x(i-1) x(i)];
fpl=[f(i-1) 0];
plot(xpl,fpl,’k:’);
end
Notas:
i) La función de actualización es: φ(x) = x - f(x)/f’(x)
ii) El método es de segundo orden: cuadráticamente convergente. Esto significa que el número de cifras decimales correctas se duplica en cada iteración cuando x(k) → x*
Podemos apreciar en la ecuación 2.7 que define la actualización del método de Newton, que se pueden presentar problemas si cerca de la solución buscada la derivada f’ se hace cero o muy pequeña. El siguiente teorema nos aclara qué se puede esperar de la iteración 2.7.
Teorema: Si la solución x* de f(x) = 0 está en un intervalo I en donde f’(x) ≠ 0 y además f’ es continua en I y se cumple la condición de convexidad:
Entonces para todo x(0) ∈ I, la iteración de Newton produce una secuencia {x(k)} que satisface exactamente una de las siguientes condiciones:
i) x(k) ∈ I y la secuencia {x(k)} converge monótonamente a x*
ii) f(x(0)) f(x(1)) < 0 y la secuencia {x(k)} converge monótonamente a x*
iii) x(1) I (si I = R esta opción no es válida)
Corolario 1: La condición de convexidad se puede cambiar a una de concavidad:
Y el teorema sigue siendo válido.
Corolario 2: En vez de condiciones de convexidad o concavidad, se puede utilizar la condición: f” (x) ≠ 0 y es continua ∀ x ∈ I.
2.3.2 Interpolación cuadrática
Dados los problemas del método de Newton clásico cuando la derivada cambia de signo y la iteración se torna inestable, se han propuesto esquemas para evitar dividir por la derivada f’(x).
En la metodología de interpolación cuadrática se aproxima f(x) en torno a x = x(k) con un polinomio de segundo orden:
y x(k+1) se obtiene haciendo P(x(k+1)) = 0:
Puesto que se supone que x(k+1) ≈ x(k), se escoge la raíz con menor valor absoluto en la fórmula de iteración.
Comentarios:
i) La cancelación de términos de similar valor absoluto, pero distinto signo, puede conducir a pérdida de cifras significativas en la aritmética.
ii) Dado que no se puede garantizar en general que el discriminante γ(k) en la ecuación 2.11 sea positivo, lo que se hace es remplazar uno de los factores (x-x(k)) por la expresión del método de Newton clásico 2.7, con lo que se obtiene el método de Newton extendido:
El siguiente ejemplo muestra que este método permite reducir las oscilaciones que ocurren con el método de Newton clásico, aunque no se puede concluir que el algoritmo 2.12 sea siempre robusto con respecto al estimador inicial x(1).
Ejemplo 2.5. Comparación del método de Newton clásico y su versión extendida
Compare el desempeño del método de Newton clásico con su variante extendida para el caso de f(x) = x3-x-1 partiendo con x(1) = -1.
Se obtiene la siguiente figura, donde se aprecia que el método clásico oscila en las primeras 10 iteraciones sin acercarse a la raíz correcta, mientras que al método extendido le toma solamente 4 o 5 iteraciones llegar a la raíz correcta, aun cuando la derivada cerca de x(1) = -1 va cambiando de signo.
FIGURA 2.5. Desempeño de la iteración de Newton clásica y extendida
A continuación se presenta el código de Matlab® que permite hacer los cálculos y la figura anterior.
% macro que resuelve el ejemplo 2.4
f=@(x) x^3-x-1;
fp=@(x) 3*x^2-1;
fpp=@(x) 6*x;
% secuencia iterativa del método de Newton y de su variante cuadrática
xN=zeros(10,1);
fN=zeros(10,1);
xN(1)=-1;
fN(1)=f(xN(1));
xC=zeros(10,1);
fC=zeros(10,1);
xC(1)=-1;
fC(1)=f(xC(1));
for i=2:10,
% iteración del método de Newton
xN(i)=xN(i-1)-fN(i-1)/fp(xN(i-1));
fN(i)=f(xN(i));
% iteración por método de Newton extendido
der=fp(xC(i-1));
xC(i)=xC(i-1)-fC(i-1)/(der-0.5*fC(i-1)*fpp(xC(i-1))/der);
fC(i)=f(xC(i));
end
% gráfico de comparación de resultados
plot(1:10,xN,’k-o’,1:10,xC,’k--^’);
axis square;
title(‘f(x) = x^3 - x - 1’);
xlabel(‘Número de iteraciones’);
ylabel(‘Raíz estimada’);
legend(‘Newton clásico’,’Newton extendido’);
dc(3) = k2*c(2);
2.3.3 Rutinas implementadas en Matlab® para ecuaciones escalares
En Matlab®, tenemos disponible como rutina de propósito general la rutina fzero. Esta rutina es híbrida, ya que combina varias técnicas en forma heurística para lograr convergencia. La rutina decide si realiza una bisección, el método de la secante o una interpolación cuadrática inversa, dependiendo de las condiciones locales. Se recomienda practicar su uso empleando la rutina fzerogui desarrollada por Moler [1]. Esta función está limitada solamente al caso de funciones continuas y no es capaz de hallar raíces donde la función no cruza el eje cambiando de signo, como es el caso de x2 en x = 0.
Para el caso particular de los polinomios, Matlab® tiene la función roots que calcula todas las raíces. En caso de que se quiera hacer cero una ecuación cuando hay datos reales, con errores de medición, se puede plantear el problema original 2.1 como la minimización de (f)2 y aplicar entonces la función de minimización escalar fminbnd.
2.4 Sistemas de ecuaciones: el método de Newton y sus variantes
La naturaleza vectorial de un sistema de ecuaciones, en contraste con el carácter escalar de las ecuaciones en una incógnita vistas anteriormente, no modifica mucho la manera en que se resuelve iterativamente el problema. No obstante, veremos algunos detalles que se deben considerar en el momento de implementar los distintos algoritmos.
Ya hemos visto que el método de Newton corresponde a una interpolación lineal local (porque la serie de Taylor def(x) era local). De modo análogo procedemos para los sistemas de ecuaciones: aproximamos
Donde el gradiente, o jacobiano J, de f viene dado por
Por lo tanto, para obtener el nuevo valor x(k+1) tenemos que resolver el sistema lineal
Para la corrección d(k) = x(k+1) -x(k). Este es el método de Newton para n variables. Nótese que es necesario que la inversa del jacobiano exista (no debe haber singularidad). La siguiente tabla muestra el algoritmo estándar del método en n variables.
TABLA 2.2. Algoritmo típico del método de Newton
2.4.1 Variaciones del método de Newton
En la ecuación 2.15, J(k) tiene que ser evaluado en cada iteración a partir de expresiones para las derivadas parciales. Para sistemas con muchas ecuaciones, esto no es práctico. Existen varias alternativas para reducir la cantidad de trabajo o simplificar la técnica de solución.
a) Podemos remplazar J(k) por una matriz A constante, con lo que se obtiene el método de la cuerda paralela:
b) Si escogemos A = I, tenemos una iteración de punto fijo en cada variable por separado, con pérdida de rapidez de convergencia. Una elección más apropiada es A = J(0), con lo cual la iteración queda:
Esta se utiliza frecuentemente en resolución de ecuaciones diferenciales ordinarias por métodos implícitos (capítulo 3). Generalmente, el jacobiano J se revalúa cada m iteraciones, para mejorar la rapidez de convergencia.
c) Cuando las derivadas parciales no se entregan explícitamente, el jacobiano se puede determinar por diferencias finitas. Esto conduce al método de Newton discreto:
Comentarios:
De las aproximaciones recién comentadas se desprenden dos desventajas:
i) Hay que evaluar el jacobiano (a veces por diferencias finitas). Esto implica n2 operaciones por paso (se puede reducir si se actualiza el jacobiano cada m iteraciones).
ii) La solución para d(k) implica n3/3 operaciones por iteración, aunque si el jacobiano se fija, esto significa n2 operaciones/iteración.
En general, las funciones f son caras de evaluar en problemas de ingeniería de procesos, por lo que los métodos de Newton son costosos en términos de tiempo computacional. Para disminuir el número de operaciones/iteración, se recurre a los métodos de cuasi-Newton, que usan una aproximación al jacobiano, los que son actualizados en cada paso a fin de reducir el número de operaciones. Para más detalles acerca de estos métodos, se recomienda consultar el artículo de Sargent [2] y el de Dennis y Mori [3].
El siguiente ejemplo muestra un código Matlab® para resolver un sistema de n ecuaciones usando el método de Newton discreto 2.18. Se calcula cada columna ‘j’ del jacobiano introduciendo una perturbación a la variable ‘j’ mediante diferencias finitas.
Ejemplo 2.6. Aplicación del método de Newton discreto en N variables
Se desea resolver el sistema de ecuaciones no lineales:
Se usa el método de Newton discreto. El método programado debe controlar que el número de iteraciones no exceda un valor máximo, y como criterio de error utilizar la norma del vector f: err = ||f(x(k))||= norm(f) que debe ser menor a 1·10-5. ¿Cuántas raíces distintas encuentra usted?
La siguiente tabla muestra los resultados de aplicar el método de Newton discreto al sistema de ecuaciones anterior, con distintos vectores iniciales x(1). En este caso encontramos tres raíces diferentes. Esto se debe a que si eliminamos x2 y x3 en las ecuaciones del sistema, hallamos una ecuación cúbica para x1, la cual posee tres raíces reales distintas. La función que implementa el método de Newton y el sistema de ecuaciones se entregan a continuación de la tabla; hay que hacer notar que el vector f(x) debe ser un vector columna.
TABLA 2.3. Solución del ejemplo 2.6
A continuación el listado del macro de Matlab® que resuelve el ejemplo 2.6.
function [raiz,f,err,nit]=newton_dis_N(fun,x1,tol,maxit)
% Esta rutina resuelve un sistema de ecuaciones no lineales usando el método clásico de
% Newton-Raphson con jacobiano discreto. Variables de entrada:
% fun: función externa que calcula el vector f(x)
% x0: vector inicial estimado para la solución
% tol: tolerancia en error relativo de x
% maxit: máximo número de iteraciones aceptable
% variables de salida:
% raíz: vector solución estimada
% f: valor de la función en la solución estimada
% err: estimación del error alcanzado
% nit: número de iteraciones realizadas por el método
if nargin<4,
maxit=100;
end
if nargin<3,
tol=1e-5;
end
if nargin<2,
error(‘no se puede calcular nada’);
end
%% inicialización del cálculo
xr=x1; iter=0;
ep=1e-4; % incremento para calcular diferencias finitas en el jacobiano
N=length(x1); I=eye(N); J=zeros(N,N);
%% ciclo principal del cálculo
while (1)
xrold=xr;
f=fun(xrold);
for j=1:N,
J(:,j)=(fun(xrold+ep*I(:,j))-f)/ep;
end
xr=xrold-J\f;
ef=norm(f); % error en la norma de f(x)
iter=iter+1; % incremento del número de iteraciones
if ef<tol || iter>maxit,
break;
end
end
%% asignación de variables de salida
raiz=xr; err=ef; nit=iter;
function y = f1(x)
%sistema de ecuaciones no lineales del ejemplo 2.5
y(1) = x(1)^2+x(2)^2-x(3)-2;
y(2) = x(1)+5*x(2)+1;
y(3) = x(1)*x(3)-2*x(1)+1;
y = y(:);
2.4.2. Rutinas implementadas en Matlab® para sistemas de ecuaciones
En el caso de sistemas de ecuaciones lineales, Matlab® tiene como herramienta básica de resolución la rutina fsolve. Esta rutina recibe como argumentos mínimos una función de Matlab® donde se expresa el sistema de ecuaciones como función del vector x, que es el segundo argumento mínimo requerido por fsolve. La función puede aplicar tres algoritmos diferentes: Gauss-Newton, Levenberg-Marquardt y Región de Confianza, y por defecto calcula el jacobiano del sistema de ecuaciones por diferencias finitas.
Si esta función tiene dificultades para hallar la solución, sobre todo cuando la dimensión de x aumenta, entonces se puede intentar minimizar la función escalar f(x)·f(x), que es mínima cuando f(x) = 0. Así se puede probar con fminunc que minimiza una función escalar de un vector x sin restricciones, o bien probar fmincon que opera con restricciones y por tanto reduce el espacio de búsqueda de la solución.
2.5 Problemas propuestos
2.5.1 Método del punto fijo para ecuaciones escalares
1) Considere la resolución de la ecuación cúbica:
Se propone encontrar una raíz positiva (>1) de la ecuación anterior, utilizando las siguientes funciones de actualización:
a) Despejando la raíz cúbica, llámela φ1(x) = (x2+2)1/3
b) Despejando la raíz cuadrada, llámela φ2(x) = (x3-2)1/2
Realice iteraciones del punto fijo usando las dos funciones de actualización ya descritas. Indique si las iteraciones están convergiendo o no, y por qué se produce tal comportamiento.
2) Para el flujo turbulento en una cañería lisa, el factor de fricción c viene dado por la solución de la ecuación algebraica:
Donde NRe es el número de Reynolds. Utilice un método de punto fijo para hallar el valor de c para los siguientes valores del número de Reynolds: NRe = 104, 105, 106. Un punto de partida para el coeficiente de fricción puede ser la fórmula de Blasius: c = 0.316(NRe)–0.25
3) Se desea calcular la densidad del CO2 supercrítico a la presión de 104 kPa y una temperatura de 340 K. En estas condiciones una ecuación de estado apropiada para caracterizar las propiedades p-v-T del fluido es la ecuación de estado de Peng-Robinson, dada por:
En donde a = 350 m6kPa/kmol2 y b = 0.07 m3/kmol.
Proponga una fórmula de iteración de punto fijo y resuelva para el volumen molar del sistema.
4) R. DeSantis ha deducido la siguiente ecuación para el factor de compresibilidad Z:
Donde b es la constante de Van der Waals y v el volumen molar. Si b = 0.08 L/mol y Z = 0.8, proponga una iteración de punto fijo y encuentre el volumen molar del sistema.
2.5.2 Métodos de interpolación para ecuaciones escalares
1) Una partícula en caída libre alcanza una velocidad terminal v en m/s que está dada por la siguiente ecuación algebraica:
a) Resuelva esta ecuación en forma numérica mediante el método de Newton.
b) Establezca un criterio de convergencia para detener la iteración e incorpórelo dentro de una macro de Matlab®.
c) Verifique si obtuvo o no convergencia, si el método funciona o no, etcétera.
2) La ecuación no lineal
posee tres soluciones distintas en el intervalo [-1, 4]. Programe una función en Matlab® que resuelva la ecuación anterior mediante el método de Newton para distintos valores iniciales x(1). Luego aplíquela para hallar las tres raíces de f(x).
3) Resuelva el problema anterior empleando el método de interpolación cuadrática para hallar las tres raíces de la ecuación algebraica.
4) Para un problema de equilibrio químico, hay que resolver la siguiente ecuación no lineal:
Donde 0 < x < 1. Escriba un macro que utilice el método de interpolación cuadrática para resolver la ecuación (1). ¿Cómo podría Ud. decidir si su iteración está convergiendo o no? Explique.
5) Para el flujo turbulento en una cañería con diámetro D y espesor e, el factor de fricción c viene dado por la solución de la ecuación de Colebrook:
Donde NRe es el número de Reynolds.
a) Programe una función en Matlab® para hallar c que tome como argumentos los valores de NRe y de e/D, y que encuentre c mediante el método de Newton.
b) Ahora haga un macro en Matlab® que calcule el valor de c para distintas combinaciones NRe = 104 a 108, e/D = 0.0001 a 0.01 y que luego haga un gráfico que represente la función c(NRe, e/D) en el dominio donde se han hecho los cálculos.
6) Considere la sedimentación de una partícula esférica sólida en un medio líquido en reposo. La partícula está sometida a la fuerza de gravedad (Fg) y la fuerza de fricción con el líquido (Fd). La ecuación para la velocidad terminal Vt de una partícula de diámetro DP viene dada por la expresión:
Donde ρ y ρP son las densidades del líquido y de la partícula, respectivamente; g es la aceleración de gravedad (9.8 m/s2); y el factor de fricción CD está determinado por el número de Reynolds de la partícula: Re = Vt·ρ·DP/μ, (μ es la viscosidad del líquido) por la siguiente expresión:
Construya una función de Matlab® que reciba como argumentos de entrada ρ, μ, ρP, DP, y entregue como argumentos de salida: a) Número de Reynolds Re; b) Coeficiente de fricción CD; y c) La velocidad de sedimentación Vt (en m/s). Pruebe su función con los siguientes datos de entrada:
i) ρ = 1000 kg/m3, μ = 9·10-4 kg/m·s, ρP = 2600 kg/m3 , DP = 0.02 m.
ii) ρ = 1000 kg/m3, μ = 9·10-4 kg/m·s, ρP = 2000 kg/m3 , DP = 2·10–4 m.
iii) ρ = 1000 kg/m3, μ = 9·10-4 kg/m·s, ρP = 1500 kg/m3 , DP = 2.5·10–6 m.
2.5.3 Sistemas de ecuaciones
1) Los problemas de equilibrio químico en mezclas de compuestos producen sistemas de ecuaciones no lineales con múltiples soluciones. Considere por ejemplo, el siguiente sistema de reacciones en equilibrio:
Si inicialmente se combinan 3 moles de H2 con 1 mol de CO hasta que se alcanza el equilibrio, las concentraciones de equilibrio se obtienen al resolver el siguiente sistema de ecuaciones algebraicas:
