Simulador de reactores químicos basado en Excel (página 3)

Tabla 1. Valores recomendados para .

La ecuación NRTL es, en realidad, un conjunto de 3 ecuaciones, como se muestra en seguida entre (2.32) y (2.34).

(2.32)

(2.33)

(2.34)

A partir de la derivación de la ecuación (2.32) pueden obtenerse las expresiones (2.35) y (2.36), mediante las cuales es posible calcular los coeficientes de actividad de las sustancias presentes en la mezcla.

(2.35)

(2.36)

Los parámetros están definidos como se aprecia en (2.37) y (2.38).

(2.37)

(2.38)

A su vez, los parámetros presentan las equivalencias observadas en (2.39) y (2.40).

(2.39)

(2.40)

Se deduce de lo anterior que el modelo de actividad NRTL posee 2 parámetros ajustables ( y ), aunque contiene también una tercera constante () que depende de la naturaleza de los componentes.

Si se continúa aceptando como cierta la hipótesis de que las interacciones entre 2 moléculas generan los mayores aportes a las propiedades en exceso, la ecuación NRTL puede generalizarse a sistemas multicomponentes, obteniéndose las ecuaciones (2.41) a (2.43).

(2.41)

(2.42)

(2.43)

Las pruebas realizadas por Renon y Prausnitz en mezclas binarias revelaron que, en sistemas simétricos (), la energía de Gibbs en exceso disminuye conforme crece, lo cual se debe a la reducción del número de interacciones entre las moléculas 1 y las 2, causadas por el aumento de la aleatoriedad; por otra parte, el parámetro mostró ser un buen indicativo del grado de alejamiento de la idealidad de la mezcla, puesto que, entre mayor sea su valor, mayor es tal distanciamiento. Además, los autores mencionados notaron que, en los sistemas asimétricos, el coeficiente de actividad de la sustancia 1 a dilución infinita depende principalmente de , mientras que depende fundamentalmente de . De otro lado, se determinó que el alto grado de asimetría de un sistema (lo cual se traduce en la obtención de un prominente sesgo en la curva de energía de Gibbs en exceso contra la fracción molar de un componente) se relaciona con la diferencia entre y , mientras que en sistemas moderadamente asimétricos, la no idealidad se relaciona con su suma, al tiempo que la curva es más o menos achatada dependiendo del valor de . Por último, Renon y Prausnitz concluyeron que y son funciones débiles de la temperatura, sin embargo, en la mayoría de los casos es posible emplear datos isobáricos para generar parámetros promedio que sean útiles en diferentes rangos de temperatura.

2.3.4. Ventajas comparativas. En su artículo, Renon y Prausnitz comparan el desempeño de la ecuación NRTL con los de las ecuaciones de Wilson, Heil y van Laar. Este análisis permitió concluir que la ecuación de Wilson es especialmente precisa para representar el equilibrio líquido-vapor de sistemas alcohol-hidrocarburo, a pesar de que no es capaz de predecir sistemas en los que haya separación de fases, lo cual sí realiza la ecuación de Heil, aunque sus predicciones no mejoran significativamente las del modelo de van Laar (siendo este último de menor complejidad que el primero). Por su parte, NRTL generó los mejores ajustes para todos los tipos de mezclas analizadas cuando fue fijado adecuadamente, de acuerdo con las recomendaciones de la Tabla 1. En cuanto al equilibrio líquido-líquido, los autores mencionados observaron que la ecuación de Wilson no es aplicable en este caso, notaron también que NRTL es aplicable toda vez que sea menor a 0.426 y concluyeron que el modelo de Heil es aplicable siempre que las sustancias empleadas no formen una mezcla altamente no ideal con bajas solubilidades. Los experimentos mostraron también que los parámetros del NRTL varían en este caso linealmente con la temperatura de modo que es posible predecir el equilibrio líquido-vapor extrapolando linealmente con la temperatura los valores obtenidos del equilibrio líquido-líquido a una temperatura menor.

2.4. MODELO TERMODINÁMICO PARA UN SIMULADOR DE REACTORES DE TANQUE AGITADO

Como se verá en el capítulo 5, el volumen de control del reactor simulado en el presente trabajo de grado contiene únicamente y en todo momento, líquido comprimido o incluso en su punto de burbuja, sin tener en cuenta la fase gaseosa que pueda haber asociada a dicho líquido. De acuerdo con esta idea, el paquete seleccionado debe ser robusto para el cálculo de propiedades en la fase líquida, puesto que ella será siempre la fase de operación. En este caso, se ha seleccionado como modelo termodinámico la ecuación de estado de Patel y Teja (PT), como regla de mezclado se escogió la de Wong y Sandler (WS), dando como resultado el modelo PTWS, y para determinar las capacidades caloríficas se optó por el modelo de Bureš. Las razones que justifican las anteriores elecciones se exponen en las secciones siguientes.

2.4.1. Ecuación de estado de Patel y Teja, (PT, 1982). Las ecuaciones de estado pueden distinguirse principalmente, como se vio antes, según el número de parámetros que contengan. Las ecuaciones con más de 3 parámetros representan de forma precisa datos volumétricos de las sustancias, sin embargo, rara vez se emplean en cálculos de equilibrios de fase debido no solo a la dificultad que presentan para la determinación de expresiones generalizadas que sean apropiadas para cálculos en mezclas, sino también al elevado tiempo computacional requerido para resolverlas. En esta área la ecuación PT busca aminorar tales dificultades requiriendo únicamente la temperatura y la presión críticas, así como 2 parámetros adicionales para caracterizar cada sustancia.

Por otra parte, La ecuación PT predice adecuadamente las propiedades volumétricas de los líquidos y adicionalmente pretende que, en el caso de sustancias no polares, el resultado sea comparable con los de las ecuaciones PR y SRK, con la ventaja complementaria de predecir valores ajustados para compuestos polares. La ecuación propuesta corresponde a la ecuación (2.44), análoga a la ecuación (2.21).

(2.44)

Cabe anotar que si se hace que el parámetro se hace igual a en (2.44) se obtiene la ecuación PR, mientras que si se hace que sea igual a cero en (2.44), se obtiene la ecuación SRK. Las ecuaciones de SRK y PR asumen que el valor del factor de compresibilidad crítico permanece constante y por tal motivo las densidades en fase líquida predichas por estas ecuaciones difieren considerablemente de sus valores experimentales. Para lograr una predicción adecuada del comportamiento de las sustancias a bajas y altas presiones se requiere que el factor de compresibilidad crítico () sea tratado como un parámetro empírico diferente, por lo general, del valor de ZC.

Luego de un proceso derivativo sujeto a unas cuantas constricciones, se deduce que las constantes , y , poseen las estructuras señaladas por (2.45), (2.9) y (2.46).

(2.45)

(2.9)

(2.46)

Las definiciones de y se presentan a continuación en (2.47) y (2.48).

(2.47)

(2.48)

Por otro lado, el valor de corresponde a la menor raíz positiva de la ecuación cúbica mostrada en (2.49).

(2.49)

A su vez, puede determinarse a partir de la expresión (2.50).

(2.50)

Además, para la función en la ecuación (2.45), se define la expresión observada en la ecuación (2.51), similar a la ecuación (2.13), aunque es claro que en (2.51) se emplea el factor empírico en lugar de la función del factor acéntrico.

(2.51)

– Cálculo de los parámetros y . El cálculo de los parámetros de la ecuación PT se realiza por medio de un procedimiento de ensayo y error que contiene los pasos que se exponen a continuación.

a) Se fija el valor inicial para como el valor más cercano a ZC entre 0.307 y 1.1 ZC.

b) Se calculan ,y con las expresiones (2.48) a (2.50).

c) Se halla un valor que satisfaga la condición de equilibrio a lo largo de la curva de saturación para cada temperatura* .

d) Se calcula por el método de los mínimos cuadrados el valor de empleando asimismo los valores obtenidos en c) aplicados en la ecuación (2.52).

(2.52)

En la ecuación (2.52), es el número del punto de la curva de equilibrio. Cabe mencionar que la minimización de (2.52) es un proceso derivativo que conduce a un polinomio cúbico en , del cual se emplea, en los pasos siguientes del algoritmo, la menor raíz positiva encontrada.

e) Con todos los datos recopilados hasta este momento se calculan las densidades líquidas y se comparan con las reportadas por la literatura mediante la desviación absoluta promedio. Si esta desviación supera la tolerancia, se repite el algoritmo con un nuevo modificado en 0.001. Si no la supera, los valores hallados corresponden a los buscados.

La ecuación PT resulta ser, por la forma en que se ajustan los parámetros empíricos, de gran precisión para estimar densidades de fases líquidas y equilibrios líquido-vapor. Por otro lado, la ecuación PT genera mejores resultados para hidrocarburos pesados y compuestos polares que las ecuaciones SRK y PR; además, genera resultados comparables con las de estas ecuaciones para mezclas de hidrocarburos livianos empleando reglas de mezclado clásicas. Valderrama reporta, asimismo, que la ecuación PT es acertada para los cálculos de volúmenes saturados de líquidos polares o no polares y de volúmenes de vapores saturados.

La expresión desarrollada por Patel y Teja puede generalizarse con el factor acéntrico para sustancias no polares utilizando las correlaciones (2.53) y (2.54), las cuales evitan recurrir al algoritmo antes expuesto.

(2.53)

(2.54)

Sin embargo, es oportuno mencionar que la ecuación falla en regiones cercanas al punto crítico al estimar la densidad de la fase líquida, aunque esta debilidad puede contrarrestarse asumiendo que , en la región en la que , es una función lineal de la temperatura, o alternativamente, cambiando la regla de mezclado empleada por una de mayor desempeño, como la de Wong y Sandler. Por su parte, Patel modificó el término de atracción en la ecuación PT, el cual fue inicialmente definido como se aprecia en la ecuación (2.13). En tal modificación eliminó la constante aunque introdujo 4 nuevas constantes, mejorando con ello las predicciones de la presión de vapor, de la densidad en fase líquida y de la capacidad calorífica molar en fase líquida de un elevado número de sustancias en el rango práctico entre 273K y 523K. La nueva función implantada por Patel se muestra en la ecuación (2.55).

(2.55)

Es preciso en este punto mencionar que, a pesar de que esta modificación mejora el desempeño de la ecuación PT, el hecho de requerir 4 constantes dificulta enormemente su utilización y por tal motivo, el paquete termodinámico empleado en el presente trabajo de grado opera con la ecuación (2.13) en lugar de la (2.55).

2.4.2. Regla de mezclado de Wong y Sandler, (WS, 1992). Como se mencionó anteriormente, las ecuaciones de estado se utilizan ampliamente en cálculos de equilibrios de fase, aunque en la práctica se ha encontrado que para modelar comportamientos complejos de mezclas altamente no ideales se requieren reglas de mezclado diferentes a las de vdW de un solo fluido. Numerosos autores (como Panagiotopoulos y Reid) han propuesto reglas de este estilo con parámetros de interacción binarios dependientes de la composición que han sido empleadas con éxito en algunos casos particulares, aunque tales reglas no son generalizables a cualquier sistema puesto que son inconsistentes en el límite de baja densidad con la mecánica estadística, la cual prevé que el segundo coeficiente virial debe ser una función cuadrática de la composición. Las correcciones realizadas permiten el ajuste con la teoría, sin embargo, dado que son intencionales y no inherentes al modelo, no mantienen la naturaleza cúbica de la ecuación de estado cuando se aplica en mezclas.

Por otro lado, las reglas en cuestión poseen el llamado síndrome de Michelsen-Kistenmacher, lo cual dificulta en gran medida la representación de sistemas conformados por compuestos de una misma familia y con pesos moleculares similares.

Una alternativa a estas reglas de mezclado fue propuesta inicialmente por Hurón y Vidal. El método consiste en igualar la energía libre de Gibbs en exceso a presión infinita calculada con una ecuación de estado con la determinada mediante un modelo de actividad (como NRTL), suponiendo que el volumen en exceso es cero cuando la presión es infinita, lo cual conduce a una regla de mezclado lineal para el parámetro . Sin embargo, Sheng et al mostraron en sus trabajos que al igualar las energías libres en exceso de Helmholtz en lugar de las de Gibbs, el término puede calcularse con una gran variedad de reglas. A partir de los trabajos antes mencionados y teniendo en cuenta que la regla de Hurón y Vidal no es acorde con la mecánica estadística, Wong y Sandler desarrollaron un nuevo tipo de reglas de mezclado basado en la energía libre en exceso de Helmholtz, el cual, a pesar de ser independiente de la densidad, cumple con los límites de alta y baja densidad y no contradice las premisas de la mecánica estadística.

Tomando como base la ecuación de estado vdW y considerando los postulados de la mecánica estadística, puede deducirse que los parámetros de esta ecuación de estado son los que se muestran a continuación entre (2.56) y (2.58).

(2.56)

(2.57)

(2.58)

Cabe anotar que se calcula por medio de un modelo de actividad que puede ser tanto aleatorio como de composición local, basándose en la suposición de que la expresión (2.59) es cierta, lo cual es justificable si se toma en consideración el hecho de que el volumen en exceso suele ser despreciable y que la energía libre de Helmholtz en exceso es una función débil de la presión.

(2.59)

Análisis posteriores demuestran que la regla de Hurón y Vidal es un caso particular de la regla de WS cuando la función () se aproxima mediante una serie truncada en el término de orden cero y cuando todos los equivalen a cero.

Cuando Wong y Sandler compararon los resultados de su regla de mezclado con los obtenidos utilizando la regla de Panagiotopoulos-Reid (ambas reglas aplicadas a la ecuación PRSV, una modificación de PR hecha por Stryjek y Vera) y con NRTL, concluyeron que su regla es mejor que la de Panagiotopoulos-Reid puesto que no sufre del síndrome de Michelsen-Kistenmacher y notaron que sus resultados son comparables o levemente mejores que los de NRTL para cálculos en fase líquida. La regla fue probada no solo en equilibrios líquido-vapor y líquido-líquido de sistemas altamente no ideales binarios y ternarios reconocidos en la literatura por ser vagamente descritos por los modelos existentes, sino que también fue probada en gases cerca del punto crítico, alcanzándose resultados satisfactorios en todos los casos.

2.4.3. Extensión de la regla de mezclado de Wong-Sandler a la ecuación de estado Patel-Teja (PTWS, 1997). Desde que Hurón y Vidal introdujeron el modelo de energía libre en exceso en las ecuaciones cúbicas de estado, la creación de reglas de mezclado del tipo de composición local para sistemas complejos se ha convertido en un constante objeto de investigación por parte de numerosos autores. Entre todas las reglas desarrolladas, la de WS, antes expuesta, ha recibido una especial atención por ser independiente de la densidad y por permitir la extensión a sistemas sometidos a altas presiones, de los valores determinados de los parámetros en sistemas a bajas presiones por medio de la extrapolación. Como se mostró anteriormente, la regla de Wong y Sandler fue aplicada con éxito a ecuaciones de 2 parámetros como vdW ó PRSV; sin embargo, fueron Yang et al quienes extendieron esta regla de mezclado a la ecuación de estado de 3 parámetros de Patel y Teja.

Sabiendo que la ecuación de partida es (2.44), las ecuaciones que definen la regla de mezclado WS para la ecuación PT en particular, se muestran en (2.60) y (2.61).

(2.60)

(2.61)

Aparte de las igualdades (2.60) y (2.61), la ecuación PTWS requiere las definiciones que se observan a continuación entre (2.62) y (2.66)* .

(2.62) (2.63)

(2.64) (2.65)

(2.66)

La energía en exceso de Helmholtz necesaria para calcular algunos términos de la regla de mezclado se calcula usualmente con NRTL (de hecho, el simulador desarrollado emplea este modelo) y Yang et al recomiendan que los parámetros para su aplicación se consulten en DECHEMA data series.

Por otro lado, los resultados de PTWS en los sistemas binarios no polares analizados por Yang et al, no presentan desviaciones relativas en la presión superiores al 3.29%, al tiempo que las desviaciones absolutas en la fracción molar no superan la magnitud de 0.028. Asimismo, las mezclas binarias con un componente polar no presentan desviaciones en la presión superiores al 4.63% ni superiores a 0.0545 en la fracción molar. De la misma manera, los sistemas binarios conformados por 2 sustancias polares no presentan desviaciones por encima de 7.25% en la presión ni de 0.042 en la fracción molar. Además de lo anterior, se analizaron algunos sistemas binarios con sulfuro de hidrógeno y las divergencias en la presión fueron en todos los casos inferiores a 2.06% y a 0.091 en la fracción molar, registrándose una alta precisión en las predicciones en la región cercana a la crítica (región problemática cuando PT se emplea con otras reglas de mezclado). Para el caso de los sistemas ternarios estudiados, Yang et al observaron buenas aproximaciones con los datos de la composición de las fases líquido-vapor coexistentes que habían sido reportados anteriormente, no superando desviaciones de 0.032 en las fracciones molares. En cuanto al cálculo de las densidades de fase en equilibrio líquido-vapor, las discrepancias en el líquido no superaron el 5.50% y los investigadores concluyeron que el modelo opera correctamente en sistemas altamente asimétricos incluso en la región cercana a la crítica. Adicional a esto, el cálculo de densidades de fluidos comprimidos no reportó divergencias por encima del 2%. Por último, los factores de compresibilidad determinados no mostraron desviaciones superiores al 3.01%.

Un análisis global de los datos anteriormente expuestos permite afirmar que PTWS es aplicable con un buen nivel de confiabilidad a la predicción del equilibrio líquido-vapor de numerosos sistemas, incluyendo sistemas altamente asimétricos; igualmente, puede aplicarse satisfactoriamente no solo a la determinación de algunas propiedades volumétricas de líquidos y gases sino también de algunas variables de estado como la presión, sobre amplios rangos operativos que incluyen las cercanías de la región crítica.

2.4.4. Modelo de Bureš para las capacidades caloríficas molares de los gases. Las capacidades caloríficas de las sustancias puras en fase gaseosa ideal se requieren en las industrias de procesos químicos para cálculos de equilibrio químico y de fase y en balances de entalpía, como en el presente trabajo de grado.

Como consecuencia del mejoramiento en las técnicas de cómputo y el desarrollo de las bases de datos, las ecuaciones se emplean cada vez con mayor frecuencia que las tablas. Las ecuaciones parabólicas e hiperbólicas comúnmente utilizadas para evaluar las capacidades caloríficas molares como funciones de la temperatura no permiten la extrapolación confiable por fuera de su intervalo de construcción, lo cual se constituye en una desventaja sustancial cuando es necesario calcular el equilibrio químico a altas temperaturas bajo condiciones adiabáticas o cuando se desea diseñar un reactor adiabático (bien sea que la reacción sea endotérmica o exotérmica). Cuando se extrapola, debido a la naturaleza misma de las ecuaciones, las expresiones parabólicas e hiperbólicas convencionales tienden a algún infinito, el negativo o el positivo, brindando resultados alejados de la realidad. Esta dificultad puede superarse con una ecuación como la observada en (2.67).

(2.67)

La igualdad (2.67) está fundamentada en postulados de la termodinámica estadística, por lo cual representa la frecuencia de vibración característica para un grupo correspondiente de átomos; además, si se sustituye un valor promedio por una serie de vibraciones, (2.67) describe la dependencia de la capacidad calorífica de toda una molécula con la temperatura.

Los comportamientos límites de la expresión (2.67) son de notoria importancia puesto que la acotan tanto superior como inferiormente y evitan el problema de los infinitos al extrapolar que se mencionó antes. Por tales razones, se desarrollan a continuación las demostraciones de los resultados a los que convergen los valores límites, las cuales no aparecen en el artículo original (solo se registran los valores finales). Cuando la temperatura es cercana a cero, mediante la definición matemática de límite, se deduce la expresión (2.68).

(2.68)

Si se realiza el cambio de variable , se asume que cuando tiende a cero por derecha, tiende a infinito positivo y se obtiene la expresión (2.69).

(2.69)

Cuando se invierte el exponencial en el numerador y se expande el término cuadrático del denominador, resulta la expresión (2.70).

(2.70)

Realizando el producto indicado en el denominador se estructura la expresión (2.71).

(2.71)

Aplicando la regla de L’Hôpital al límite indicado en (2.71), se deduce la expresión (2.72).

(2.72)

Aplicando de nuevo la regla de L’Hôpital a la expresión mostrada en (2.72), se obtiene el límite (2.73).

(2.73)

Si se realiza el cambio de variable , se deriva de ello que cuando tiende a infinito positivo, tiende asimismo a infinito positivo y se obtiene la expresión (2.74).

(2.74)

Operando sobre la expresión (2.59) puede llegarse a la expresión (2.75).

(2.75)

El límite presentado en (2.75) corresponde al tipo , por lo cual, al dividir por la mayor potencia presente () y evaluando el límite se obtiene la expresión (2.76), que corresponde al límite inferior de la ecuación de Bureš.

(2.76)

Por otro lado, cuando la temperatura tiende a infinito positivo se tiene la expresión (2.77).

(2.77)

Si se realiza el cambio de variable , puede determinarse que cuando tiende a infinito positivo, tiende a cero por derecha y se obtiene la expresión (2.78).

(2.78)

Luego de algunas manipulaciones idénticas a las expuestas para el caso del límite inferior, puede llegarse a la expresión (2.79).

(2.79)

Dado que en la expresión (2.79) no existen indeterminaciones de ningún tipo, el límite puede evaluarse sin inconveniente de la manera que se muestra en (2.80).

(2.80)

Consecuentemente, el límite superior de la ecuación de Bureš corresponde a .

Una ventaja adicional de la ecuación (2.67) es que resulta fácilmente integrable analíticamente, tanto para cálculos de entalpías como de entropías, lo cual evita recurrir a soluciones numéricas que solo brindarían valores aproximados. Debido a la importancia de estas expresiones dentro del desarrollo del presente trabajo de grado, estas integraciones se muestran detalladamente a continuación, dado que en el artículo original solo se reportan sus resultados, obviando el procedimiento.

Las entalpías requieren la integración de la capacidad calorífica mostrada en (2.81), para su cálculo.

(2.81)

La integral en el miembro derecho de la igualdad (2.81) equivale, en términos de la igualdad (2.67) a la expresión (2.82).

(2.82)

Si se realiza el cambio de variable , los diferenciales poseen la relación , con lo cual se obtiene la expresión (2.83) (cabe anotar que los límites de la integral no se cambiarán porque al final de la integración se retornará a la variable ).

(2.83)

Realizando un nuevo cambio de variable , los diferenciales tienen la siguiente relación , y resulta la expresión (2.84).

(2.84)

La integral del miembro izquierdo de la igualdad (2.84) se evalúa fácilmente para obtener la expresión (2.85).

(2.85)

Revirtiendo los cambios de variable realizados, se obtiene el resultado en función de la temperatura que se observa en (2.86).

(2.86)

En cuanto a la entropía, la expresión para calcularla en términos de la capacidad calorífica corresponde a la igualdad (2.87).

(2.87)

La integral en el miembro derecho de la igualdad (2.87) equivale, en términos de la igualdad (2.67) a la expresión observada en (2.88).

(2.88)

Teniendo en cuenta que la integral es un operador lineal, puede separarse la integral del miembro derecho de (2.88) en dos partes, de acuerdo con la aditividad del operador lineal. Realizando nuevamente el cambio de variable en la segunda parte de la integral, los diferenciales se relacionan según , y con esto, se llega a la expresión (2.89) (nuevamente se menciona que los límites de la integral no se cambiarán puesto que al final de la integración se retornará a la variable ).

(2.89)

Evaluando la primera integral del miembro derecho de (2.89) y efectuando en la segunda el cambio de variable , en el que los diferenciales tienen la relación , se llega a la igualdad indicada en (2.90).

(2.90)

Para poder evaluar la integral restante en (2.90) se debe efectuar inicialmente una integración por partes que conduzca a la expresión (2.91) y, posteriormente, debe realizarse un proceso de descomposición en fracciones parciales de la integral resultante, para obtener la expresión (2.92).

(2.91)

(2.92)

Reemplazando lo encontrado en la expresión (2.92) en la igualdad (2.90) se deduce la expresión (2.93).

(2.93)

Restituyendo los cambios de variable efectuados anteriormente y aplicando algunas propiedades de los exponenciales y de los logaritmos, se llega a la expresión (2.94).

(2.94)

Además de discutir los comportamientos límites de la ecuación (2.67), es igualmente importante mencionar que dicha ecuación fue empleada por Bureš para ajustar el de 250 sustancias entre hidrocarburos y sus principales derivados con oxígeno, nitrógeno, azufre y cloro, así como algunas sustancias inorgánicas. Bureš evaluó el valor de las constantes , y , con base en el método iterativo de Gauss-Newton y estimó que la ecuación (2.67) puede emplearse en el amplio rango entre 200K y 5000K sin obtenerse errores mayores a 1.38%, en promedio, con respecto a los valores reportados en la literatura, mientras que las ecuaciones clásicas solo operan, generalmente, entre 273K y 1500k. Por otra parte, este autor notó en sus investigaciones que al extrapolar las ecuaciones polinómicas y racionales convencionalmente empleadas por encima de 1500K, se originan valores poco precisos que no pueden ser utilizados en cálculos subsiguientes. A su vez, Henao resalta la versatilidad del modelo de Bureš mencionando el hecho de que con tan solo 3 constantes genera buenos resultados mientras que los modelos más conocidos requieren 4, 5 ó más constantes para su operación, con el agravante de que tales modelos no permiten la extrapolación. Asimismo, Henao presenta en su libro un extenso compendio de constantes con cerca de 1700 sustancias para el modelo de Bureš, que incluye no solo hidrocarburos sino también derivados halogenados, alcoholes y fenoles, aldehídos, cetonas, ácidos carboxílicos, ésteres, éteres, aminas y nitrilos, entre otros.

3. MÉTODOS NUMÉRICOS PARA RESOLVER SISTEMAS DE ECUACIONES DIFERENCIALES

3.1. CONTEXTUALIZACIÓN,

La simulación digital es una poderosa herramienta para resolver las ecuaciones que describen sistemas de ingeniería química, aunque posee dos dificultades: la solución simultánea de ecuaciones algebraicas no lineales y la resolución numérica de ecuaciones diferenciales ordinarias. La primera se resuelve empleando algún método iterativo y la segunda utilizando ecuaciones en diferencias finitas. La precisión y la estabilidad de estas aproximaciones debe ser tenida en cuenta porque el método para obtenerlas o el algoritmo de solución empleado afectan notoriamente la convergencia. En la actualidad existen numerosos algoritmos cuya eficacia depende del tipo de problema y aunque infortunadamente no existe un algoritmo que opere de manera adecuada para todo tipo de situaciones, algunos autores recomiendan el algoritmo explícito simple de primer orden de Euler para un gran número de aplicaciones de ingeniería.

A través de los años han sido desarrollados numerosos paquetes para la simulación que han eximido al ingeniero de conocer los métodos numéricos empleados, porque los programas detectan automáticamente los errores y la estabilidad y ajustan parámetros de solución como el intervalo y el tamaño de paso de manera que la solución cumpla con la tolerancia especificada. La solución numérica de las ecuaciones diferenciales ordinarias se puede realizar con métodos explícitos como el algoritmo de Euler o el de Runge-Kutta, o con métodos implícitos como el de Euler.

El modelo dinámico de un reactor de tanque agitado genera un sistema de ecuaciones algebraico-diferenciales de alta complejidad que no permite una resolución analítica, por lo que es necesario recurrir a los métodos numéricos para obtener una solución.

Los problemas de valores iniciales para una sola ecuación diferencial ordinaria poseen la estructura de la ecuación (3.1).

(3.1)

Con y siendo números reales positivos o cero.

La condición inicial del problema es y es la aproximación de la solución en el tiempo i.

Teniendo en cuenta que N es el número de particiones realizadas en el intervalo comprendido entre y , el tamaño de paso se calcula, en este caso, como se observa a continuación en (3.2).

(3.2)

3.2. MÉTODO DE EULER

El método de Euler surge de expandir la solución de la ecuación diferencial en términos de la serie de Taylor truncada en el primer término, de acuerdo con esto la aproximación de la solución a la ecuación diferencial es la ecuación (3.3).

(3.3)

Esto quiere decir que para poder hallar el valor de la aproximación en ti+1 se requiere conocer el valor de la aproximación () y la función () en ti.

Este método genera una aproximación a la solución en varios valores en el intervalo de hasta , estos valores se denominan puntos de red o nodos. Los puntos de red tienen una distribución uniforme en el intervalo de tiempo y están separados entre sí por el tamaño de paso, cuando el tamaño de paso disminuye debe haber una mayor precisión en las aproximaciones, sin embargo reducir excesivamente el tamaño de paso requiere un número mayor de cálculos y esto aumenta el error debido al redondeo, por lo cual se busca un equilibrio al seleccionar el tamaño de paso adecuado teniendo en cuenta dicho error.

El algoritmo del método de Euler consiste en definir el intervalo de la solución (fijar los valores de y de ), la condición inicial de la solución () y el tamaño de paso () y luego evaluar la aproximación de la solución desde hasta cada unidades de tiempo.

3.3. MÉTODOS DE TAYLOR DE ORDEN SUPERIOR

Estos métodos son una ampliación del método de Euler en los cuales se trunca la serie de Taylor utilizada al aproximar la solución de la ecuación diferencial en un término , siendo el orden del método de Taylor (el método de Euler corresponde al método de Taylor de orden 1). De acuerdo con esto, los métodos de Taylor tienen la solución general mostrada en (3.4).

(3.4)

Donde : Orden del método de Taylor que se va a emplear.

El procedimiento de solución de estos métodos numéricos es análogo al del método de Euler, no obstante en estos métodos es necesario evaluar las derivadas de la función en cada uno de los pasos. En este último punto es donde se encuentra la principal dificultad de aplicación de los métodos de orden superior, dado que el cálculo y la evaluación de las derivadas es un procedimiento lento y complicado, aunque a medida que aumenta el orden del método de Taylor más precisa es la aproximación solución.

3.4. MÉTODOS DE RUNGE-KUTTA

La aplicación de los métodos de Taylor de orden superior para obtener la solución de una ecuación diferencial ordinaria con un valor inicial suele complicarse por la necesidad de determinar las derivadas de orden superior y evaluarlas a través del tiempo. Los métodos de Runge-Kutta resultan de modificar los métodos de Taylor para que el orden de las cotas del error se conserve, pero evitan la necesidad de emplear derivadas de órdenes elevados en la solución. Las técnicas mencionadas utilizan el desarrollo de Taylor de ƒ, la función que aparece en el miembro derecho de la ecuación diferencial que se va a solucionar. Dado que ƒ es una función de dos variables (t, y) debe emplearse el teorema generalizado de Taylor a funciones de este tipo que resulta más complicado que el de una sola variable debido a que aparecen en él todas las derivadas parciales de ƒ. Este método requiere que ƒ y todas sus derivadas hasta la n+1 sean continuas dentro de los valores empleados para t y y, de modo que el teorema de Taylor sea aplicable.

La aproximación de las técnicas de Runge-Kutta podría incrementar el error, pero se hace de manera que el error introducido sea del mismo orden del que ya presenta el método de Taylor y, en consecuencia, los nuevos errores no afectan significativamente los cálculos.

El método de Runge-Kutta más básico es el de segundo orden, también conocido como método del punto medio, el cual se presenta a continuación en (3.5) y (3.6).

(3.5)

(3.6)

Con

El método de Runge-Kutta de segundo orden puede ser depurado si se mejora la aproximación que lleva incorporada y de esta manera se obtienen el método de Euler modificado y el método de Heun. El método de Euler modificado se presenta a continuación en (3.5) y (3.7).

(3.5)

(3.7)

El método de Heun está representado por las ecuaciones (3.5) y (3.8).

(3.5)

(3.8)

Las fórmulas de Taylor de orden superior pueden convertirse en técnicas de Runge-Kutta de una forma análoga aunque las manipulaciones algebraicas son complicadas. De estos métodos de orden superior, el más común es el de orden 4 que contiene cuatro evaluaciones de la función y que resulta de resolver un sistema de ecuaciones con 12 incógnitas. El método consiste en el sistema mostrado entre (3.9) y (3.14).

  (3.9)

(3.10)

(3.11)

(3.12)

(3.13)

(3.14)

El principal esfuerzo computacional en la aplicación de los métodos de Runge-Kutta consiste en la evaluación de f en cada paso. En los métodos de segundo orden el error en cada paso es función de h al cubo a costa de dos evaluaciones de la función, mientras que en los de orden 4 se requieren cuatro evaluaciones por paso y el error local es función de h a la quinta potencia. De ahí en adelante el relativo decrecimiento del orden del error hace que sean preferibles los métodos de orden menor que 5 con tamaños de paso menores con respecto a métodos de orden superior con mayor tamaño de paso.

La comparación entre los métodos de Runge-Kutta de orden bajo se hace con base en el número de evaluaciones por paso. Así, si el de cuarto orden requiere de 4 evaluaciones por paso y el de Euler requiere una sola evaluación, se considera que el primero debe dar respuestas más precisas que el segundo cuando el segundo emplea un tamaño de paso equivalente a la cuarta parte del primero. En todas las comparaciones posibles el método de cuarto orden ha probado ser el más preciso y eficiente y por ello es el de mayor aplicación.

3.5. MÉTODOS MULTIPASO

La base de estos métodos consiste en emplear un predictor-corrector. Estas técnicas de resolución numérica emplean la información proveniente de más de uno de los puntos de red precedentes para determinar la aproximación del siguiente punto.

La solución general de los métodos multipaso es la ecuación (3.15).

(3.15)

Donde: y : Constantes que dependen del método específico que se va a emplear.

Los métodos multipaso requieren de m valores iniciales y como únicamente se dispone de un valor inicial los restantes m-1 valores se obtienen usando un método de un paso, luego se procede a evaluar las funciones en cada punto sucesivo empleando un tamaño de paso adecuado.

Los métodos multipaso se dividen en dos categorías: la primera se conoce como método explícito y ocurre cuando el coeficiente βm es igual a cero. En este caso la solución en el punto que se va a evaluar se puede obtener directamente. La segunda, en caso de que el coeficiente βm sea diferente de cero, se denomina método implícito. Para aplicar directamente este método se debe resolver la ecuación implícita para ti+1, esta solución es particular a cada problema y no se asegura la consecución de una solución única para ti+1. En la práctica se utiliza una combinación de ambos métodos, una aproximación dada por un método explícito que es corregida evaluando este valor en el miembro derecho de la ecuación implícita; la combinación de ambas técnicas se denomina método predictor-corrector.

La ecuación explícita que debe utilizarse para evaluar un método-predictor corrector es la ecuación (3.16).

(3.16)

Esta expresión se evalúa en y se utiliza para obtener el valor final de ti+1 con el método implícito, es decir, se corrige el valor inicial derivado del método explícito empleando la ecuación (3.17).

(3.17)

Las principales ventajas de los métodos que emplean un predictor-corrector frente a los métodos de un solo paso radican en que son más estables y el error acumulado crece más lentamente, por lo que son preferibles cuando el número de puntos que es necesario evaluar es elevado, dado que se aprovechan los cálculos anteriores y por ello solo hay que evaluar en cada iteración, sin necesidad de calcular varias constantes en cada paso.

En los algoritmos de los métodos predictor-corrector se pueden incluir técnicas de extrapolación para mejorar la precisión de las aproximaciones a las soluciones de problemas de valor inicial. Esto se logra incluyendo en el algoritmo un paso en el que se realiza la extrapolación, el cual puede incluirse tanto en el paso del predictor como en el del corrector.

El principal inconveniente que tiene la utilización de estos métodos consiste en su dependencia de los valores iniciales que se les suministren y que en la práctica deben obtenerse a partir de un método de un solo paso, con lo cual, a pesar de tener un buen control del error, la precisión depende del método empleado para generar los primeros puntos.

Existen distintos tipos de métodos multipaso, la diferencia entre ellos radica en los valores fijados para las constantes; la forma más común de hallar dichos valores se basa en el método de coeficientes indeterminados, y las constantes halladas dependen de los coeficientes que cada método específico fija en cero. A continuación se muestra la ecuación general para determinar los coeficientes yde los métodos multipaso.

(3.18)

Donde: = 0,1…n

En esta ecuación se asume que 00 equivale a 1; además, en algunos de los métodos se fijan previamente los valores de algunas de las constantes y.

Los métodos más comúnmente empleados son los de orden 4 porque equilibran la precisión del método con la simplicidad de las fórmulas lo cual disminuye los errores derivados del redondeo.

3.5.1. Método predictor-corrector de Adams-Bashforth-Moulton. Esta técnica se basa en el método de Adams-Bashforth como predictor, en el cual los valores de y son iguales a cero con i desde 2 hasta n, y el método de Adams-Moulton como corrector, en el que los valores de equivalen a cero con i desde 2 hasta n.

El algoritmo del método predictor-corrector de Adams-Bashforth-Moulton de cuarto orden es el conformado por (3.19) y (3.20).

(3.19)

(3.20)

Donde: : Resultado del método predictor. : Resultado del método corrector

3.5.2. Método predictor-corrector de Milne-Simpson. En el método de Milne-Simpson a diferencia del método de Adams-Bashforth-Moulton se tiene el valor de igual a cero y diferente de cero. Con estas condiciones se reduce el error local de truncamiento. Este método presenta problemas de estabilidad y por lo tanto no es recomendable su aplicación.

El algoritmo del método predictor-corrector de Milne-Simpson de cuarto orden es el formado por las ecuaciones (3.21) y (3.22).

(3.21)

(3.22)

3.5.3. Método predictor-corrector de Hamming. Este método se basa en el predictor-corrector de Milne-Simpson, pues emplea el mismo predictor y un corrector en el cual se corrigen los problemas de estabilidad del método anterior.

El algoritmo del método predictor-corrector de Hamming de cuarto orden consta de las ecuaciones (3.23) y (3.24).

(3.23)

(3.24)

3.6. TÉCNICAS ADAPTABLES

Las técnicas adaptables emplean un tamaño de paso variable que permite hacer aproximaciones más eficientes en cuanto al número de cálculos que se deben realizar en cada paso, pero tienen la dificultad de que su aplicación es complicada de implementar en el proceso de solución. Estos métodos adoptan el número y la posición de los nodos que se utilizan en la aproximación para mantener el error local dentro de una cota especificada inicialmente. Adicionalmente, tienen una notable ventaja que consiste en que evitan calcular las derivadas de orden superior de la función debido a que el procedimiento de selección del tamaño de paso proporciona una estimación del error local que es netamente algebraico. Las técnicas adaptables parten de suponer que un método ideal sería aquel en el que dada una tolerancia () mayor que cero, utilizaría el menor número de nodos con el que fuera posible garantizar que el error global no sería nunca mayor que la tolerancia en ningún punto de la partición, lo cual sería inconsistente con un tamaño de paso constante. Para acercarse a este caso ideal se emplean métodos de diferentes órdenes consecutivos que permitan estimar el error local (dado que el global generalmente no se puede determinar) y, empleando esa estimación, elegir un tamaño de paso que controle el error global. Si q es un múltiplo del tamaño de paso, es la aproximación con el método de orden superior (n+1) y es la aproximación de orden inferior (n), entonces se tiene la desigualdad (3.25) en la que se muestra el intervalo de valores de q que aseguran que el nuevo tamaño de paso (qh) mantenga el error global dentro de la cota.

  (3.25)

Cuando q<1 se rechaza la elección inicial de h para el paso i-ésimo y se repiten los cálculos empleando qh como tamaño de paso. Si se acepta el valor calculado para el paso i-ésimo con tamaño de paso h y se modifica el tamaño a qh para el paso (i+1)-ésimo.

Una de las técnicas adaptables que se emplea con mayor frecuencia es la de Runge-Kutta-Fehlberg (conocida como RKF45), que utiliza los métodos de Runge-Kutta de quinto y cuarto orden como se muestra en (3.26) y (3.27).

(3.26)

(3.27)

El valor de q y de las constantes se determina como se muestra entre (3.28) y (3.34).

(3.28)

(3.29)

(3.30)

(3.31)

(3.32)

(3.33)

(3.34)

Este método tiene la ventaja de solo requerir seis evaluaciones de ƒ en cada paso mientras que si se emplearan independientemente serían cuatro evaluaciones para el de cuarto orden y seis para el de quinto orden, siendo el RKF45, además de esto, de mayor precisión que los métodos antes expuestos. El método RKF45 puede ser mejorado si al implementarlo se evitan grandes o pequeñas modificaciones del paso y en algunas ocasiones si el paso nunca se aumenta de tamaño y solo se reduce cuando sea necesario controlar el error.

Otra técnica adaptable es el predictor–corrector de Adams con tamaño de paso variable que emplea el método explícito de Adams-Bashforth de cuatro pasos para hacer la predicción y el método implícito de Adams-Moulton de tres pasos como corrector de la aproximación. Debido a que este método es multipaso y adicionalmente tiene paso variable, es necesario calcular los valores de partida en nuevos nodos igualmente espaciados y cualquier cambio en el tamaño de paso requiere que se calculen nuevos valores de partida en ese punto, lo cual hace que el algoritmo sea más complicado que los métodos de un solo paso.

3.7. ECUACIONES DIFERENCIALES RÍGIDAS

Todos los métodos para la aproximación numérica de la solución de un problema de valores iniciales, bien sea una sola ecuación o un sistema de ellas, poseen fórmulas para el error que incluyen una derivada de orden superior al de la solución de la ecuación. Uno de los supuestos más importantes de estas técnicas es que el error puede mantenerse bajo control, sin embargo, surgen con frecuencia problemas en los que el error crece tanto que llega a dominar los cálculos; tales problemas se componen de ecuaciones diferenciales rígidas, que son aquellas en cuya solución analítica aparece un término de la forma , en donde c es una constante positiva de gran magnitud. Generalmente, este término exponencial es tan solo una parte de la solución, llamada solución transitoria, que se suma a la solución permanente, la parte más importante de la solución.

La parte transitoria de la solución decae con rapidez a medida que t crece debido a su carácter exponencial, pero su enésima derivada , no decae tan rápidamente e incluso puede que crezca conforme lo haga t, desviando el error por encima de la cota. Esta dificultad se supera comúnmente porque las condiciones físicas del problema del cual se deriva una ecuación de este estilo permiten predecir su rigidez y tomar las medidas apropiadas para controlar el error, lo cual se hace con una restricción sobre el tamaño del paso.

Para establecer el comportamiento de un método numérico concreto cuando se aplica a un sistema rígido, se realiza una prueba preliminar con la ecuación (3.35), la cual es una ecuación rígida con parte permanente nula, en la que λ es un número real negativo, como se muestra a continuación.

(3.35)

La solución analítica de la ecuación (3.35) es la ecuación (3.36), una vez que se ha evaluado la condición inicial.

(3.36)

En general, debe existir una función Q tal que al aplicarse a la ecuación de prueba (3.35), se obtenga la aproximación indicada en (3.37).

(3.37)

La precisión del método, y por ello su aplicabilidad a sistemas rígidos dependerá del grado de acercamiento de a ; además, el tamaño del paso debe ser escogido de modo que se satisfaga la condición (3.38).

(3.38)

La desigualdad mostrada en (3.38) debe ser satisfecha puesto que, de lo contrario, el error crecerá sin límite y no se alcanzará la aproximación de la solución buscada.

3.8. SOLUCIÓN DE SISTEMAS DE ECUACIONES DIFERENCIALES EN UN SIMULADOR DE REACTORES DE TANQUE AGITADO

Para determinar la solución numérica de un sistema de ecuaciones con valores iniciales se pueden emplear los métodos revisados anteriormente para una sola ecuación. El método seleccionado para la solución del sistema de ecuaciones diferenciales resultante del modelo del reactor, el cual se expondrá en el capítulo 5, es la técnica de Runge-Kutta de cuarto orden (RK4) puesto que puede extenderse con facilidad a sistemas de ecuaciones, es un método estable, no requiere el cálculo de la derivada de la función ni de un alto número de evaluaciones de la misma en cada paso. Todo esto, sumado al hecho de que no necesita un tamaño de paso excesivamente pequeño para obtener aproximaciones precisas, permite afirmar que el método es adecuado para realizar los procedimientos requeridos. Sin embargo, es claro que los demás métodos de aproximación de un solo paso también pueden extenderse a los sistemas de ecuaciones en cuestión. Asimismo, los métodos multipaso y las técnicas de predicción-corrección también pueden ampliarse a dichos sistemas, aunque el nivel de cálculos es elevado. De nuevo, si se emplea el control del error, cada componente del conjunto de aproximaciones debe ser lo suficientemente precisa, ya que de otra forma, debe calcularse de nuevo la solución numérica. La extensión de la técnica de extrapolación también puede realizarse, aunque en la práctica no es utilizada con frecuencia debido a que su notación es sumamente complicada.

Como se verá en el capítulo 5, el modelo de un reactor de tanque agitado en estado transitorio genera un sistema de balances de masa (donde es el número de sustancias que atraviesan el volumen de control) y un balance de energía, lo cual genera un sistema de ecuaciones, con variables dependientes, las moles de cada componente en el interior del reactor y la energía interna del volumen de control, más una variable independiente: el tiempo (). Un esquema general del conjunto de ecuaciones en cuestión puede observarse en (3.39).

(3.39)

Este sistema de ecuaciones diferenciales de primer orden se encuentra sujeto al conjunto de condiciones iniciales (3.40), cuando .

(3.40)

Para el método de Runge-Kutta de cuarto orden h y N tienen idénticos significados que anteriormente. Se utilizará la notación para denotar una aproximación de y se definirá el tiempo en cada nodo como se muestra en (3.41).

(3.41)

Es necesario calcular cuatro constantes para cada función ui en cada tj, transformando las condiciones iniciales, como se muestra entre (3.42) y (3.46).

(3.42)

(3.43)

(3.44)

(3.45)

(3.46)

Las aproximaciones sucesivas de las funciones que integran el sistema se calculan empleando la ecuación (3.47).

(3.47)

El método RK4 es una técnica que ofrece una excelente relación entre la precisión de la aproximación y el número de cálculos que es necesario realizar para determinarla, sin embargo, al extender técnicas adaptables como RKF45 a sistemas de ecuaciones pueden obtenerse valores más precisos de la aproximación de la solución, aunque el número de cálculos y, por tanto, el tiempo de cómputo aumentan significativamente y, en ocasiones, como se vio antes, la notación del método puede resultar sumamente compleja.

4. MÉTODOS NUMÉRICOS PARA RESOLVER ECUACIONES DE UNA SOLA VARIABLE

Como se observa en el capítulo 6, en el cual se explica el algoritmo de solución del problema planteado, el presente trabajo de grado únicamente precisa resolver una ecuación en lugar de un sistema de ellas, empleando un método numérico. Por tal razón, se exponen a continuación los principales métodos numéricos utilizados para aproximar la solución de una ecuación algebraica de una sola variable y se describe también el algoritmo computacional que emplea el simulador desarrollado.

4.1. CONTEXTUALIZACIÓN

El problema de solucionar una ecuación algebraica en la variable , consiste básicamente en encontrar un valor de dicha variable que satisfaga la condición (4.1).

(4.1)

Cuando se trata de la solución exacta de la ecuación, se obtiene un valor que satisface por completo la igualdad presente en (4.1), mientras que al tratarse de una solución numérica, se habla de un valor que satisface la condición (4.1) con una determinada precisión; de tal forma que el valor es una aproximación de la solución , que hace que el miembro derecho de la expresión (4.1) sea un valor cercano a cero.

El nivel de aproximación a cero está directamente relacionado con el grado de aproximación entre y , razón por la cual es necesario aproximarse lo suficiente, aunque es igualmente adecuado no excederse en la precisión para evitar problemas de convergencia y elevados tiempos de cálculo. Para esto, los métodos tienen incorporados diversos criterios que detienen su ejecución. Tales criterios se basan generalmente en el concepto de tolerancia (), el cual es un valor especificado al principio del método, que se compara con la diferencia absoluta de la aproximación de la solución hallada en un paso y la aproximación del paso anterior: cuando tal diferencia es menor que el valor fijado inicialmente, el criterio se satisface y el método se detiene, de otra forma, continúan las iteraciones. Otro tipo de tolerancia se compara con una diferencia absoluta entre y cero, con la consecuente detención del método al satisfacerse que la diferencia absoluta sea menor que el valor fijado. Por otra parte, cuando no se presentan variaciones significativas entre la aproximación de un paso y el del anterior, un criterio de parada alternativo consiste en definir un número máximo de iteraciones que al cumplirse interrumpa el método, con lo cual se evitan ciclos infinitos.

4.2. MÉTODO DE LA BISECCIÓN

El método de la bisección está basado en el teorema del valor intermedio, el cual establece que para una función continua , definida en el intervalo , si la función definida en el punto tiene el signo contrario a la función definida en , esto es , entonces deberá existir al menos una solución para dentro de dicho intervalo. Para hallar esta solución se van reduciendo sistemáticamente a la mitad los subintervalos de y en cada paso se usa la mitad en la cual exista un cambio de signos al evaluar la función. El algoritmo es el siguiente:

a) Se determina , como el punto medio del intervalo .

b) Se evalúa en , y , para obtener , y , respectivamente.

c) Si , entonces la raíz buscada está en el intervalo , si , entonces la raíz buscada se encuentra en el intervalo . Deben renombrarse las variables para el nuevo intervalo y continuar subdividiendo dicho intervalo hasta cumplir el criterio de parada.

d) Criterio de parada: se continúa hasta satisfacer alguna de las condiciones mostradas entre (4.2) y (4.4).

(4.2)

(4.3)

(4.4)

Cabe anotar que, en cada paso, la aproximación se calcula de acuerdo con la ecuación (4.5).

(4.5)

Se recomienda usar está expresión en lugar de su simplificación para evitar problemas asociados al redondeo cuando esté cerca de la precisión de la máquina en la que se realicen los cálculos.

Existe un teorema que permite predecir el número de iteraciones necesarias para alcanzar un valor dado de tolerancia, el cual se muestra en la expresión (4.6). Aunque este valor en la mayoría de los casos es mayor de lo que realmente se requiere, el resultado opera como un buen estimativo o también como un límite superior que puede servir para restringir el número máximo de iteraciones.

(4.6)

Al despejar de la ecuación (4.6) la variable y teniendo en cuenta que es la tolerancia, se deduce la desigualdad mostrada en (4.7).

(4.7)

Las principales desventajas del método radican en que debe conocerse el intervalo en el que se encuentra la solución, solo sirve para hallar raíces reales y el proceso de cálculo es enormemente laborioso.

4.3. MÉTODO DE NEWTON-RAPHSON

Este método es una de las técnicas numéricas con mayor velocidad de convergencia. La técnica parte del polinomio de Taylor truncado en el primer término alrededor del punto , como se muestra en (4.8).

(4.8)

Mediante manipulaciones matemáticas sencillas es posible obtener una expresión para a partir de la aproximación mostrada en (4.8), teniendo en cuenta que es cero puesto que es la solución de la ecuación (4.1). El proceso mencionado concluye con la deducción la aproximación mostrada en (4.9).

(4.9)

Esta aproximación es la base del método de Newton-Raphson, en la cual se comienza por una aproximación inicial , suministrada por quien utiliza el método, hasta encontrar una solución , dentro de la tolerancia especificada

El algoritmo se basa en la expresión (4.10) para el cálculo sucesivo de las aproximaciones.

(4.10)

En la ecuación (4.10) se observa que, si la derivada de la función en es igual a cero para algún valor de , el método debe interrumpirse y no se alcanza una solución apropiada.

La obtención de este método a partir de las series de Taylor resalta la importancia de una aproximación inicial () cercana a la solución para garantizar la convergencia del método.

Por otro lado, la necesidad de calcular la derivada de la función en cada paso es una condición con un enorme inconveniente: tiene asociada la generación de errores significativos, puesto que los métodos numéricos para el cálculo de derivadas son, en general, poco precisos. Además, cuando la derivada de la función tiende a cero en un punto o en una región, la convergencia puede llegar a ser inalcanzable (este problema persiste al extender el método de Newton-Raphson a sistemas de ecuaciones, traduciéndose en la condición de que el determinante del jacobiano tiende a cero).

4.4. MÉTODO DE LA SECANTE

El método de Newton-Raphson es una técnica robusta de rápida convergencia, sin embargo, requiere conocer el valor de la derivada de la función en cada paso, lo cual, como se mencionó antes, es un proceso con considerables errores incorporados. Para evitar este problema, se utiliza el método de la secante, el cual utiliza dos puntos (2 aproximaciones iniciales) para aproximar la derivada sin necesidad de emplear diferenciales. Esta aproximación parte del límite observado en (4.11).

(4.11)

Si el valor de se hace igual a y se omite el límite, la igualdad mostrada (4.11) se transforma en la expresión (4.12).

(4.12)

Cuando la aproximación presente en (4.12) se reemplaza en la expresión (4.10), base del algoritmo de Newton-Raphson, se obtiene la expresión (4.13), la cual permite aproximar la solución por el método de la secante.

(4.13)

La convergencia de este método es un poco menor que la del método de Newton-Raphson. El método de la secante y el de Newton-Raphson por lo general se emplean para refinar las respuestas conseguidas con otra técnica más robusta, como la bisección. Un inconveniente del método de la secante consiste en que, a pesar de tener dos valores iniciales, la raíz no necesariamente se encuentra entre ellos, lo cual no permite asegurar la convergencia del método (como sí sucede en la bisección).