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).