5. Raíces de ecuaciones 5.1 Métodos cerrados Parte II: Análisis Numérico 1 5.1.1 Métodos Gráficos Un método simple para obtener una aproximación a la raíz de la ecuación f(x)=0, consiste en graficar la función y observar donde cruza el eje x. Ejemplo: Utilizar gráficas por computadora para localizar las raíces de f(x) = x3 + x2 -3·x+5 Solución. Utilizando MATLAB, 0 THEN xl = xr fl = f(xl) il = 0 iu = iu+1 IF iu ≥ 2 THEN fu = fu/2 ELSE ea = 0 END IF IF ea < es OR iter ≥ imax THEN EXIT END DO ModFalsePos = xr END ModFalsePos Parte II: Análisis Numérico 19 Ejercicios Ejercicio 5.1 Determine las raíces reales de f(x) = -0.4x2 + 2.2x + 4.7: a. Gráficamente b. Usando el método de bisección para determinar la raíz más grande. Emplee como valores iniciales xl=5 y xu=10. Calcule el error estimado εa y el error verdadero εt para cada iteración. Ejercicio 5.2 Calcule la raíz real positiva de f(x)=x4-8x3-36x2+462x 1010 utilizando el método de la falsa posición. Use una gráfica para escoger el valor inicial y realice el cálculo con εs = 1.0 % Parte II: Análisis Numérico 20 Ejercicio 5.3 La concentración de saturación de oxígeno disuelto en agua se calcula con la ecuación 1.575701 × 10 5 6.642308 × 10 7 1.243800 × 1010 8.621949 × 1011 ln Osf = −139.34411 + − + − 2 3 4 Ta Ta Ta Ta donde Osf = concentración de saturación de oxígeno disuelto en agua a 1 atm (mg/L) y Ta = Temperatura absoluta (K). Recuerde que Ta = T + 273.15, donde T = temperatura (ºC). De acuerdo con ésta ecuación, la saturación disminuye con el incremento de la temperatura. Para aguas naturales típicas en climas templados, la ecuación sirve para determinar rangos de concentración de oxígeno desde 14.621 mg/L a 0ºC hasta 6.949 mg/L a 35ºC. Dado un valor de concentración de oxígeno, ésta fórmula y el método de bisección son útiles para resolver la temperatura en ºC. Si los valores iniciales se fijan en 0 y 35ºC, desarrolle y pruebe un programa de bisección para determinar T como una función de una concentración de oxígeno dada. Pruebe el programa para Osf=8, 10 y 14 mg/L. Compruebe sus resultados Parte II: Análisis Numérico 21 5.2 Métodos abiertos Parte II: Análisis Numérico 22 5.2.1 Iteración simple de punto fijo Los métodos abiertos utilizan una fórmula para predecir la raíz. Esta fórmula puede desarrollarse como una iteración simple de punto fijo (También llamada iteración de un punto o sustitución sucesiva o método de punto fijo), al reordenar la ecuación f(x)=0 de tal modo que x esté del lado izquierdo de la ecuación: x=g(x) x2 + 3 Por ejemplo, x2-2x+3 = 0, se reordena para obtener x = 2 Mientras que sen(x)=0, puede transformarse sumando x a ambos lados para obtener x=sen(x)+x Parte II: Análisis Numérico 23 De ésta manera, dado un valor inicial para la raíz xi , la ecuación anterior puede usarse para obtener una nueva aproximación xi+1, expresada por la fórmula iterativa xi+1=g(xi) El error aproximado se calcula usando el error normalizado: xi +1 − xi εa = 100% xi +1 Parte II: Análisis Numérico 24 Ejemplo Iteración simple de punto fijo Planteamiento del problema. Use una iteración simple de punto fijo para localizar la raíz de f(x) = e-x - x Solución. xi+1=e-xi i 1 2 3 4 5 6 7 8 9 10 xi 1 0.367879 0.692201 0.500473 0.606244 0.545396 0.579612 0.560115 0.571143 0.564479 Parte II: Análisis Numérico εa % 100.0 171.8 46.9 38.3 17.4 11.2 5.90 3.48 1.93 1.11 ετ % 76.3 35.1 22.1 11.8 6.89 3.83 2.20 1.24 0.705 0.399 25 Convergencia El error relativo porcentual verdadero en cada iteración del ejemplo anterior, es proporcional (por un factor de 0.5 a 0.6) al error de la iteración anterior. Esta propiedad se conoce como convergencia lineal. Parte II: Análisis Numérico 26 f(x) Un método gráfico alternativo consiste en separar la ecuación en dos partes, de esta manera f1(x)=f2(x) Entonces las dos ecuaciones y1 = f1(x) y y2 = f2(x) se grafican por separado. Así, los valores de x correspondientes a Las intersecciones de estas dos funciones representan las raíces de f(x)=0 Parte II: Análisis Numérico f(x) = e-x - x Raíz f(x) f2(x) = e-x f1(x) = x Raíz 27 y y1 = x y2= g(x) y y1 = x y2= g(x) x2 x1 x0 x y2= g(x) x0 x y y2= g(x) y1 = x y y1 = x x0 x x0 x Parte II: Análisis Numérico 28 FUNCTION Fixpt(x0, es, imax) xr = x0 iter = 0 DO xrold = xr xr = g(xrold) iter = iter+1 IF xr ≠ 0 THEN ea = xr − xrold ⋅ 100 xr END IF IF ea < es OR iter ≥ imax EXIT END DO Fixpt = xr END fixpt Parte II: Análisis Numérico 29 5.2.2 Método de Newton-Raphson A partir de la expansión en series f(x) de Taylor, se tiene: Pendiente = f ’(xi) f ( xi ) − 0 f ' ( xi ) = xi − xi +1 que se reordena para obtener f(xi) xi +1 f ( xi ) = xi − f ' ( xi ) 0 xi+1 xi x la cual se conoce como fórmula De Newton Raphson Parte II: Análisis Numérico 30 Ejemplo Método de Newton-Raphson Planteamiento del problema. Utilice el método de Newton Raphson para calcular la raíz de f(x)=e-x – x empleando como valor inicial x0 = 0 Solución. La primer derivada de la función es f ’(x)=-e-x -1 e − xi − x i que se sustituye para obtener xi +1 = xi − − xi −e −1 i 0 1 2 3 4 xi 0 0.500000000 0.566311003 0.567143165 0.567143290 Parte II: Análisis Numérico εt (%) 100 11.8 0.147 0.0000220 < 10-8 31 Algoritmo 1. 2. Se debe incluir una rutina de graficación en el programa Al final de los cálculos, se necesitará sustituir siempre la raíz final calculada en la función original, para determinar si el resultado se acerca a cero. Esta prueba protege el desarrollo del programa contra aquellos casos en los que se presenta convergencia lenta u oscilatoria, la cual puede llevar a valores pequeños de εa, mientras que la solución aún está muy lejos de una raíz. El programa deberá incluir siempre un límite máximo permitido del número de iteraciones para estar prevenidos contra soluciones oscilantes, de lenta convergencia o divergentes que podrían persistir en forma interminable. El programa deberá alertar al usuario para que tome en cuenta la posibilidad de que f ‘(x) sea igual a cero en cualquier momento durante el cálculo. Parte II: Análisis Numérico 32 3. 4. 5.2.3 El método de la secante Un problema potencial en la implementación del método de Newton Raphson es la evaluación de la derivada. En casos complejos, la derivada se puede aproximar mediante una diferencia finita dividida hacia atrás f ( xi −1 ) − f ( xi ) f ' ( xi ) ≡ xi −1 − xi Sustituyendo en la ecuación de Newton - Raphson f(x i) f(x i-1) x i-1 xi xi +1 f ( xi )( xi −1 − xi ) = xi − f ( xi −1 ) − f ( xi ) Parte II: Análisis Numérico 33 Ejemplo El método de la secante Planteamiento del problema. Con el método de la secante, calcule la raíz de f(x)=e-x –x. Comience los cálculos iniciales con los valores x-1=0 y x0 = 1.0. Solución. Primera iteración: x-1=0 f(x-1)=1 x0 =1 f(x0)=-0.63212 x1=1-((-0.63212)(0-1)/(1-(-0.63212)))=0.61270 Segunda iteración x0=1 f(x0)=-0.63212 x1 =0.61270 f(x1)=-0.07081 x2=0.61270-((-0.0708)(1-0.61270)/(-0.63212- … … (0.07081))) = 0.56384 Parte II: Análisis Numérico 34 Método de la secante modificada En lugar de considerar dos valores arbitrarios para aproximar la derivada, un método alternativo considera un cambio fraccionario de la variable independiente para estimar f’(x), f ( xi + δxi ) − f ( xi ) f ' ( x) ≅ δxi donde d es un pequeño cambio fraccionario. Ésta aproximación se sustituya en la ecuación de la secante para obtener la siguiente expresión iterativa: xi +1 δxi f ( xi ) = xi − f ( xi + δxi ) − f ( xi ) Parte II: Análisis Numérico 35 Ejercicios Ejercicio 5.4 Evaluar las raíces de las siguientes ecuaciones trascendentes a. sin x - 2exp(-x2) = 0 b. ax - ax = 0 para a = 2, e, or 3 c. ln(1 + x2) – x1/2= 0 d. e-x/(1 + cos x) - 1 = 0 Parte II: Análisis Numérico 36 Ejercicio 5.5 Para el flujo turbulento de un fluido a través de un tubo liso, es posible establecer la siguiente relación entre el factor de fricción cf y el número de Reynolds Re: Calcular cf para Re = 104, 105 y 106. Parte II: Análisis Numérico 37 Ejercicio 5.6 Desarrolle una función para calcular el volumen específico de un gas puro, dada la temperatura y la presión usando la ecuación de estado de Soave-Redlich-Kwong Las constantes a y b son obtenidas de las ecuaciones Parte II: Análisis Numérico 38 donde Pc y Tc son la presión crítica y temperatura crítica respectivamente. La variable α es una función empírica de la Temperatura El valor de S es una función del factor acéntrico ω Las propiedades físicas del n-butano son Parte II: Análisis Numérico 39 y la constante de los gases R es Calcule el volumen específico del vapor de n-butano a 500 K y en un rango de presiones de 1 a 40 atm. Compare los resultados gráficamente con aquellos que se obtienen de la ley de los gases ideales. ¿Qué conclusión obtiene de ésta comparación gráfica? Parte II: Análisis Numérico 40 Ejercicio 5.7 Repita el ejercicio 5.6 usando las ecuaciones de estado de Benedict-Webb-Rubin (BWR) y Patel-Teja (PT). Compare gráficamente los resultados con los obtenidos en el ejercicio 3. La ecuación de estado de Benedict-Webb-Rubin (BWR) es donde A0, B0, C0, a, b, c, α, y γ son constante. Donde P está en atmósferas, V está en litros por mol, y T está en kelvin, Los valores de las constantes para el n-butano son: Parte II: Análisis Numérico 41 La ecuación de estado de Patel-Teja es Donde a es función de la temperatura, y, b y c son constantes donde Parte II: Análisis Numérico 42 y Ωb es la más pequeña de las raíces positivas del polinomio cúbico F y ζc son funciones del factor acéntrico Parte II: Análisis Numérico 43 Ejercicio 5.8 La ecuación de Underwood para destilación multicomponente está dada por donde F = tasa de flujo molar de la alimentación n = numero de componentes en la alimentación zjF = fracción molar de cada componente en la alimentación q = calidad de la alimentación αj = volatilidad relativa de cada componente en condiciones promedio de la columna φ = raíz de la ecuación Parte II: Análisis Numérico 44 Underwood ha demostrado que (n-1) de la raíces de la ecuación se encuentran entre los valores de las volatilidades relativas como se Muestra Evalúe las n-1 raíces de ésta ecuación para el caso mostrado en la Tabla F=100 mol/h q=1 (líquido saturado) Parte II: Análisis Numérico 45