Métodos numéricos aplicados a Ingeniería. Héctor Jorquera González
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.