DESCRIPCIÓN MATEMÁTICA Y ANÁLISIS DE ESTABILIDAD DE PROCESOS FERMENTATIVOS

MATHEMATICAL DESCRIPTION AND STABILITY ANALYSIS OF FERMENTATIVE PROCESSES

 

VÍCTOR MANUEL TREJOS
Ingeniero Químico, Universidad Nacional de Colombia. Sede Manizales. vmtrejosm@gmail.com

JAVIER FONTALVO ALZATE
Profesor Asociado, Universidad Nacional de Colombia. Sede Manizales. jfontalvoa@unal.edu.co

MIGUEL ÁNGEL GÓMEZ GARCIA
Profesor Asociado, Universidad Nacional de Colombia. Sede Manizales. magomez@unal.edu.co

 

Recibido para revisar junio 23 de 2008, aceptado septiembre 18 de 2008, versión final octubre 20 de 2008

 

---

RESUMEN: El auge de los procesos fermentativos como fuente de etanol carburante conduce a plantear estudios teóricos para el entendimiento detallado del proceso. En el presente trabajo se evalúa la producción de etanol en un sistema continuo con células de levadura. Los modelos de crecimiento microbiano Monod y Haldane son analizados en un CSTR reactor (Continuous Stirred Tank Reactor). Se analiza la estabilidad del sistema, el comportamiento dinámico, la caracterización de puntos de equilibrio, planos de fase y puntos de bifurcación para cada uno de los dos modelos de crecimiento.

PALABRAS CLAVE: Puntos de equilibrio, comportamiento dinámico, puntos de bifurcación.

ABSTRACT: The advances in the fermentative processes as source of ethanol fuel leads to outline theoretical studies for the detailed understanding of that process. In this paper, the ethanol production in a continuous system with yeast cells is evaluated. The behavior of Monod and Haldane models for microbial growth is analyzed for a CSTR (Continuous Stirred Tank Reactor). The stability of the system, the dynamic behavior, the characterization of fixed points, phase planes and bifurcation points for each one of the two models of microbial growth are analyzed.

KEYWORDS: Fixed points, dynamic behavior, bifurcation points.

---

 

1. INTRODUCCIÓN

En la actualidad se fermentan mieles o jugos azucarados con el fin de obtener alcoholes cuyo uso va desde la industria farmacéutica y de cosméticos, hasta el alcohol combustible para vehículos de motor. El etanol es una de las alternativas más prometedoras de combustible utilizándose como aditivo de la gasolina [1]. En Colombia, actualmente se está implementado el uso de etanol para oxigenar la gasolina y existe un creciente aumento en su producción [2,3].

En términos generales la fermentación implica el empleo de microorganismos para llevar a cabo la transformación de materia orgánica [3]. La producción comercial de productos de fermentación ha empleado diversas especies de bacterias, levaduras y hongos. Estos procesos fermentativos son autocatalíticos ya que el microorganismo no solo cataliza la reacción sino que además crece y se reproduce en el medio de reacción. Los microorganismos, en un medio de cultivo apropiado, se dividen activamente empleando los nutrientes que le aporta el medio para "fabricar" nuevos microorganismos. Este tipo de reacciones generalmente presentan inhibición por producto.

Existen diferentes configuraciones de reactores para llevar a cabo la reacción de fermentación siendo el reactor en continuo el más conveniente [4]. Las fermentaciones en continuo pueden considerarse como un sistema abierto en el que se va añadiendo biomasa y sustrato continuamente al biorreactor y se va eliminando simultáneamente igual volumen de caldo de fermentación. El reactor continuo más utilizado para llevar a cabo este tipo de reacciones, es el CSTR (de sus siglas en inglés Continuous Stirred Tank Reactor) o reactor de mezcla [4].

En la literatura se encuentran variantes de este reactor en las que, por ejemplo, el producto es removido utilizando una membrana, previniendo la inhibición por producto, mejorando la actividad celular, aumentando la productividad volumétrica y requiriendo volúmenes más pequeños de fermentador [7].

En los modelos dinámicos propuestos para las fermentaciones la concentración de biomasa, sustrato y producto son variables dependientes del tiempo que permiten evaluar el rendimiento de la fermentación. Se propone entonces la utilización de modelos matemáticos que sirvan de herramientas para la simulación de procesos fermentativos continuos, con el fin de obtener información detallada del comportamiento del reactor en cualquier instante de tiempo.

En este trabajo se evalúa el comportamiento dinámico y la estabilidad de un sistema fermentativo en continuo. Se caracteriza el sistema con sus puntos fijos, puntos de bifurcación y planos de fase para los modelos de crecimiento microbiano de Haldane y Monod. Este trabajo busca mostrar la importancia del análisis de estabilidad de un sistema en continuo y cómo éste puede definir las regiones de interés operativo, en este caso para la producción de etanol.

 

2. MODELOS CINÉTICOS EN ESTUDIO

Investigaciones en procesos fermentativos han resultado en un alto número de diferentes ecuaciones que describen el crecimiento microbiano. La más famosa de ellas es la expresión propuesta por Monod [8]. Aunque cada uno de estos modelos puede ser descrito por una ecuación flexible, en general de tres parámetros, la falta de consistencia con los datos experimentales ha conducido a desarrollar ecuaciones alternas como las propuestas por Teissier, Moser, Haldane y Hinshelwood, entre otras [9-11]. En la Tabla 1 se presentan las expresiones de los modelos de crecimiento microbiano mencionados anteriormente.

Tabla 1. Modelos de crecimiento celular [8-11]
Table 1. Models for microbial growth [8-11]

Se puede observar de la Tabla 1 la simplicidad que ofrece la ecuación de Monod a diferencia de las ecuaciones restantes Teissier, Moser, Haldane y Hinshelwood, que poseen mínimo dos o más parámetros ajustables. Aunque estas últimas ecuaciones son el resultado de modelos propuestos con el fin de mejorar el ajuste de los datos experimentales, no tienen un principio bioquímico. La ecuación de Monod está originalmente basada en el concepto de una única velocidad enzimática limitante con control de crecimiento de acuerdo con la cinética de Michaelis-Menten. Una discusión más extensa de estos modelos puede encontrarse fácilmente en la literatura [9].

Para este trabajo se escogieron los modelos de crecimiento microbiano de Haldane (inhibición por sustrato) y Monod debido a su amplia

utilización en la literatura para la descripción de procesos fermentativos [2,7,12].

 

3. MODELO MATEMÁTICO

El modelo matemático propuesto para describir el proceso tiene en cuenta una entrada y una salida en el biorreactor, comportándose como un sistema CSTR. El esquema del proceso, así como la nomenclatura de las corrientes, se presenta en la Figura 1.

Figura 1. Esquema de un CSTR para un proceso fermentativo
Figure 1. Scheme of a CSTR for a fermentative process

La reacción llevada a cabo en el reactor se puede expresar de la forma:

De donde () es sustrato, () biomasa y () concentración de producto. El balance de materia global para el sistema de la Figura 1 viene dado por [12]:

Entra – Sale + Genera = Acumula (1)

Considérese el sistema presentado en la Figura 1 donde solo hay una entrada y una salida. Asumiendo volumen de reacción constante el balance para la biomasa según la ecuación (1) toma la forma:

(2)

Además, los balances para el sustrato y el producto se expresan como:

(3)

(4)

Reagrupando las ecuaciones (2), (3) y (4) en función del factor de dilución (D=F/V) y de la productividad definida como gramos de producto/gramos de biomasa (Y_PX) y gramos de biomasa / gramos de sustrato (Y_XS) se tiene que:

(5)

(6)

(7)

El factor de dilución, en unidades de <h^-1>, corresponde al inverso del tiempo de residencia . En estado estable las expresiones en función del tiempo se vuelven cero.

(8)

Reemplazando (8) en las ecuaciones (5), (6) y (7) se obtiene que:

(9)

(10)

(11)

Las ecuaciones (9), (10) y (11) representan el comportamiento de la biomasa, sustrato y producto en estado estable.

 

4. PARÁMETROS DE LA SIMULACIÓN

Los parámetros utilizados en la simulación fueron tomados de un amplio intervalo de valores reportados en las referencias [13-19].

De acuerdo con los parámetros de simulación propuestos en la Tabla 2, se observa en la Figura 2 que cada modelo de crecimiento microbiano presenta un comportamiento diferente. Sin embargo, ambos modelos exhiben un valor máximo en la velocidad de crecimiento microbiano [4].

Tabla 2. Parámetros de la simulación
Table 2. Simulation Parameters

Figura 2. Modelos de crecimiento microbiano.
Figure 2. Models for microbial growth

De acuerdo con la Figura 2 el modelo de Haldane presenta velocidades de crecimiento altas a bajas concentraciones de sustrato y velocidades de crecimiento bajas a altas concentraciones de sustrato, mostrando el efecto negativo que tiene sobre las células el sustrato a altas concentraciones. Este fenómeno es conocido como inhibición por sustrato.

4.1 Condición De Lavado
La primera solución de las ecuaciones (9) y (10) se conoce como solución trivial o condición de lavado (washout), ya que no hay reacción química. Esto equivale a decir que no hay biomasa cargada dentro del biorreactor y por tanto la concentración inicial de sustrato es igual a la concentración de sustrato en el alimento. Teniendo en cuenta que para todos los casos no hay biomasa en la corriente de entrada (), la ecuación (9) toma la forma:

(12)

Resolviendo simultáneamente las ecuaciones (9-11) para X₁=P₁=0 se obtiene la condición de lavado válida para ambos modelos, Monod y Haldane.

4.2 Solución No Trivial
Para este caso se tiene que la concentración de biomasa es diferente de cero , por lo que las ecuaciones (9) y (10) tienen solución de la forma:

(13)

(14)

(15)

La ecuación (13) es implícita para la concentración de sustrato S₁ y por lo tanto su solución variará de acuerdo al modelo de crecimiento microbiano. A continuación se muestra el análisis de cada modelo con sus restricciones y limitaciones.

4.2.1 Modelo Monod
De acuerdo al modelo Monod y teniendo en cuenta que la ecuación (13) tiene solución algebraica, para la concentración de sustrato se tiene que:

(16)

En la ecuación (16) se observa que es una condición para que la concentración de sustrato sea numéricamente positiva. Según la ecuación (16) el valor máximo que puede alcanzar la concentración de sustrato S₁ es cuando sea igual a la concentración de alimento. Esta condición garantiza que el factor de dilución D también alcanza un valor máximo. Para el caso de estudio:

(17)

Teniendo en cuenta la ecuación (16) y la condición se puede concluir que:

(18)

La ecuación (18) es una condición necesaria para obtener las soluciones en estado estable con sentido físico. En la Figura 3 se muestra el punto de bifurcación de la función, resultado de la solución dada por la ecuación (17), además de la restricción impuesta por la ecuación (18).

Figura 3. Diagrama de bifurcación del modelo Monod
Figure 3. Bifurcation diagram for the

Los puntos donde el factor de dilución D hace que el número de soluciones y las características de estabilidad del sistema cambien se conocen como puntos de bifurcación. De acuerdo con la Figura 3 se observan las soluciones con sentido físico para el modelo Monod y se muestran las regiones de estabilidad e inestabilidad explicadas más adelante en la sección 6.

Monod model

4.2.2 Modelo de Haldane
Para el modelo de Haldane la solución de la ecuación (13) tiene forma cuadrática. Planteando la ecuación en función de la concentración de sustrato se tiene:

(19)

La ecuación (19) tiene solución cuadrática de la forma:

(20)

La función que está dentro del radical se define como el discriminante dado por:

(21)

La ecuación (21) para el modelo de Haldane, tiene las siguientes posibilidades:

:Se tienen dos raíces positivas que tienen significado físico.

:Las raíces no tienen significado físico

:Existe una única solución que corresponde a la velocidad máxima de dilución.

Este es el punto donde la función sufre inflexión para constituir tres posibles soluciones. Si el discriminante es igual a cero la ecuación (20) constituye un resultado directo del punto de inflexión de la función dado por:

(22)

De acuerdo con la ecuación (22) el punto de inflexión (turning point) se da para un factor de dilución de correspondiente a una concentración de sustrato . Las ecuaciones (14-15) no se modifican para este modelo.

El factor de dilución se comporta como el parámetro de bifurcación del sistema, separando las regiones que presentan estabilidad e inestabilidad. La Figura 4 presenta la variación de la concentración de sustrato en función del factor de dilución D.

Figura 4. Diagrama de bifurcación del modelo de Haldane
Figure 4. Bifurcation diagram for the Haldane model

En la Figura 4 se observa que a bajos factores de dilución el modelo de Haldane se comporta de igual manera que el modelo de Monod. Se muestra que para altos factores de dilución, ambos modelos tienen solamente un estado estable factible y se presentan las regiones de estabilidad e inestabilidad explicadas más adelante en la sección 6.

4.3 Solución En Estado Estable
Para ejemplificar los resultados obtenidos en las Figuras 3 y 4, considérese un factor de dilución de . La solución de la ecuación (13) para los dos modelos de crecimiento microbiano puede hacerse algebraicamente como se explicó anteriormente o por medio de un método numérico multivariable como el Newton-Raphson [20]. Los resultados son presentados en las Tabla 3.

Tabla 3. Solución para el modelo Monod y Haldane
Table 3. Answer for the Monod and Haldane model

Los valores reportados en la Tabla 3 corresponden a la solución en estado estable para cada uno de los modelos de crecimiento microbiano para un factor de dilución dado.

El modelo de Haldane presenta tres estados estables a diferencia del modelo Monod que exhibe solo dos estados estables. Se observa que las concentraciones de biomasa y sustrato varían de acuerdo al tipo de expresión que describe el crecimiento microbiano.

 

5. DINÁMICA Y ESTABILIDAD

Para comprobar los estados estables presentados en la Tabla 3 se modeló el comportamiento dinámico del sistema resolviendo las ecuaciones diferenciales (5) y (6). Para la solución de las ecuaciones diferenciales del modelo dinámico se utilizó un Runge-Kutta-Felberg con tamaño de paso variable [20].

5.1 Modelo Monod
La solución del modelo dinámico para diferentes alimentaciones de biomasa y sustrato arrojó los resultados mostrados en la Figura 5.

Figura 5. Estabilidad del sistema utilizando el modelo Monod
Figure 5. Stability of the system using the model Monod

En la Figura 5 se observa que para diferentes condiciones iniciales de concentración de biomasa, sustrato y biomasa se obtienen diferentes curvas utilizando el modelo Monod. Se puede apreciar como el sistema siempre alcanza uno de los dos posibles estados estables del sistema.

5.2 Modelo de Haldane
La solución del modelo dinámico arrojó los resultados mostrados en la Figura 6.

Figura 6. Estabilidad del sistema utilizando el modelo Haldane
Figure 6. Stability of the system using the Haldane model

En la Figura 6 se observa que para alimentaciones cercanas de biomasa, sustrato y producto cada una de las curvas tiende a estados estables diferentes. Sin embargo se puede apreciar con claridad que cuando se alimenta a diferentes concentraciones de biomasa y sustrato solo se llega a dos de los tres posibles estados estables del sistema.

 

6. ANÁLISIS DE ESTABILIDAD

Con el fin de mejorar la comprensión de la estabilidad del sistema utilizando modelos de crecimiento microbiano, se hizo una caracterización de los puntos fijos para los dos modelos utilizados. El método consiste en encontrar los eigenvalores de la función que permitan determinar si es un nodo estable, inestable o punto silla. Una revisión más detallada se presenta en las referencias [5, 6, 21].

La estabilidad de cada una de las soluciones de los estados estables es determinada a partir de los eigenvalores de la matriz jacobiana. Se desea entonces determinar las raíces de la ecuación:

(23)

La matriz se define como la matriz jacobiana de las funciones variantes con el tiempo, definidas por las ecuaciones (5), (6) y (7).

₍₂₄₎

Resolviendo la ecuación (26) se tiene que:

(25)

La ecuación (25) es una expresión cúbica, con tres eigenvalores resultantes de su solución. De acuerdo con la ecuación (25) se obtienen dos soluciones completamente diferentes para la condición de lavado y para la solución no trivial. La derivada de la velocidad específica de crecimiento presentada para cada modelo de crecimiento microbiano viene dada como:

Modelo Monod

(26)

Modelo Haldane (27)

De acuerdo con las ecuaciones (26) y (27) los resultados de estabilidad serán diferentes para cada modelo de crecimiento microbiano.

6.1 Condición De Lavado
La ecuación característica para la condición de lavado del sistema viene dada como:

(28)

La ecuación (28) tiene tres raíces dadas como: y .

6.2 Solución No Trivial
El análisis de estabilidad de la condición no trivial es dada como:

(29)

Una discusión detallada sobre la estabilidad de la ecuación (29) es dada en las referencias [5, 6, 21]. En las Tablas 5 y 6 se muestran los resultados obtenidos de la caracterización de puntos fijos para los dos modelos de crecimiento microbiano en estudio para un factor de dilución de . Se define el tipo de nodo de cada punto de estabilidad con lo que se establece que tan estable o inestable es cada uno de ellos.

Tabla 4. Caracterización de los puntos fijos para el modelo Monod
Table 4. Characterization of the fixed points for the Monod model

Tabla 5. Caracterización de los puntos fijos para el modelo de Haldane
Table 5. Characterization of the fixed points for the Haldane model

 

7. DIAGRAMAS DE FASE

El análisis de plano de fase proporciona una mejor comprensión de la estabilidad o inestabilidad de cada uno de los nodos de los dos modelos de crecimiento microbiano.

7.1 Plano De Fases Del Modelo Monod
De acuerdo con la Tabla 5, al modelarse el sistema con la ecuación del modelo Monod existen dos estados estables posibles, el primero inestable y el segundo estable, lo que indica la presencia de un nodo tipo silla y otro nodo al que llegan las líneas de tendencia. En la Figura 7 se presenta el plano de fase para el modelo Monod que confirma las observaciones hechas en la caracterización de los puntos fijos presentadas en la Tabla 5.

Figura 7. Plano de fase modelo Monod
Figure 7. Phase plane for the Monod model

7.2 Plano De Fases Del Modelo Haldane
El modelo de Inhibición por sustrato, según los resultados mostrados en la Tabla 6, presenta tres posibles estados estables, dos de ellos estables y uno inestable tipo silla.

En la Figura 8 se presenta el plano de fase para el modelo de inhibición por sustrato donde se observa con claridad las líneas de tendencia de cada uno de los nodos caracterizados en la Tabla 6.

Figura 8. Plano de fase modelo Haldane
Figure 8. Phase plane for the Haldane model

 

8. CONCLUSIONES

En este trabajo se analizó el comportamiento de un sistema fermentativo en continuo con células de levadura. Se desarrollaron modelos matemáticos con dos expresiones cinéticas de crecimiento microbiano Monod y Haldane. Se mostró la forma de utilizar modelos cinéticos con el fin de aplicarlos al análisis de procesos fermentativos: determinación de los estados estables de un biorreactor tipo CSTR, análisis del comportamiento dinámico y caracterización de los puntos fijos del sistema presentados en planos de fase. A partir del caso de estudio, la producción de etanol, se demostró que existen diferentes regiones de estabilidad. Cada una de ellas conduce a condiciones de operación con características específicas diferentes. El no considerar las particularidades de dichas regiones, en la selección de las condiciones de operación, puede conducir a alcanzar condiciones de operación de bajo rendimiento y productividad, que no son de interés. Para el caso particular de estudio para la producción de etanol, la Región I es la región que presenta el mejor intervalo de condiciones de operabilidad para ambas cinéticas de reacción.

 

NOMENCLATURA

: Velocidad de crecimiento microbiano < h^-1 >.
: Velocidad de crecimiento máxima < h^-1 >.
: Concentración de Sustrato a la salida <g L^-1>.
: Concentración de Sustrato a la entrada < g L^-1 >.
: Concentración de biomasa a la salida < g L^-1 >.
: Concentración de biomasa a la entrada < g L^-1 >.
: Volumen del biorreactor < L >.
: Velocidad de formación de biomasa < g L^-1s^-1 >.
: Velocidad de consumo de sustrato < g L^-1s^-1 >.
: Constante cinética de velocidad correspondiente a cada modelo <g L^-1>.
: Constante cinética del modelo de Haldane < L g^-1 >.
: Velocidad de flujo volumétrico < L s^-1 >.
: Gramos de producto / gramos de biomasa
: Gramos de biomasa / gramos de sustrato

 

REFERENCIAS

[1] S.S.E.H. ELNASHAIE, A.H. FAKEEHA, E. HELAL AND E. ABASHAR. A mathematical model achieving the twin objectives of simplicity and accuracy for the simulation of immobilized packed bed fermentors, Math. Comput. Modelling 19 (5), 105-114, (1994).         [ Links ]
[2] VARGAS G. JAIME, M. PEÑUELA, M. ECHEVERRI, M. ORTIZ, M. C. ESCOBAR Y J. C. QUINTERO. Producción de alcohol por fermentación con levaduras libres e inmovilizadas. Revista Facultad de Ingeniería. Nº 22, (2001).         [ Links ]
[3] RIBAS M., HURTADO R., GARRIDO N., DÍAZ M. DE LOS RÍOS, F. DOMÉNECH, R. SABADÍ, T. GARCÍA Y D. RODRÍGUEZ. Modelación matemática y simulación de procesos fermentativos, Revista de Ingeniería Química 438, (2006).         [ Links ]
[4] WARD OWEN P. Biotecnología de la fermentación Editorial Acribia, S. A. Zaragoza España (1989).         [ Links ]
[5] KUZNETZOV YURI A. Elements of applied bifurcation theory. Springer, Verlag New York , (1995).
[6] PARKER T. S. AND CHUA L. O. Practical numerical algorithms for chaotic systems. Springer Verlag, 1989.         [ Links ]
[7] BOTERO A MAHECHA, GARHYAN P. AND S.S. E. H. ELNASHAIE. Bifurcation, stabilization and ethanol productivity enhancement for a membrane fermentor. Mathematical and computer modelling. Vol. 41 pages 391-406, (2005).         [ Links ]
[8] DORAN PAULINE M. Principios de ingeniería de los bioprocesos, Edición Acribia S.A. Saragoza España (1998).         [ Links ]
[9] HEIJNEN J.J. AND B. ROMEIN. Derivation of kinetic equations for growth on single substrates based on general properties of a simple metabolic network. Biotechnology. Prog. Vol. 11, pages 712-716, (1995).         [ Links ]
[10] J.B.S. HALDANE, Enzymes, Longmans Green, London , 1930.         [ Links ]
[11] C.N. HINSHELWOOD, Proc. Roy. Soc. London Ser. A 113 (1926), p. 230.         [ Links ]
[12] E. JAMES BAILEY AND DAVID OLLIS. Biochemical Engineering Fundamentals, second Editions, Mc Graw Hill Inc 1986.         [ Links ]
[13] ELNASHAIE S. E., PARAG GARHYAN MARCEL DEKKER. Conservation equation and modelling of chemical and biochemical processes. Said, Inc (2003).         [ Links ]
[14] DAUGULIS, J. ANDREW. MCLELLAN JAMES, LI JINGHONG. Experimental Investigation and modeling of oscillatory Behavior in the Continuous Culture of Zymomonas mobilis. Biothechnology and Bioengineering. 56, p. 99-105 (1997).         [ Links ]
[15] GARHYAN PRAG, ELNASHAIE S.S.E.H. Exploration and exploitation of bifurcation/chaotic behavior of a continuous fermentor for the production of ethanol. Chem. Eng. Sci. 58, 1479-1496, (2003).         [ Links ]
[16] GHOMMIDH C, VAIJA J, BOLARINWA S., NAVARRO J.M. Oscillatory Behaviour of Zymomonas in Continuous Cultures: A Simple Stochastic Model. Biotechnology Letters. Volume 2 No.9. pp. 659-664, (1989).         [ Links ]
[17] JARZEBSKI, ANDRZEJ B. Modeling Of Oscillatory Behavior in Continuous Ethanol Fermentation. Biotechnology Letters. Volume 14 No. 2. pp. 137-142, (1992).         [ Links ]
[18] JÖBSES I. , EGBERTS G., LUYBEN K. AND ROELS J. Fermentation Kinetics of Zymomonas Mobilis at High Ethanol Concentrations: Oscilations in Continuos Cultures. Biothechnology and Bioengineering. 28, p. 868-877, (1986).         [ Links ]
[19] GÜLNUR BIROL, PEMRA DORUKER, BETÜL KIRDAR, Z. ILSEN ÖNSAN AND KUTLU ÜLGEN. Mathematical description of ethanol fermentation by immobilised Saccharomyces Cerevisiae. Process Biochemistry. Vol 33 No. 7, pages 763-771, (1998).         [ Links ]
[20] MATHEWS JHON H. Numerical Methods for mathematics science and engineering., Second edition, Prentice Hall, Inc (1992).         [ Links ]
[21] B. W. BEQUETTE. Process Dynamics, Modeling, Analysis and Simulation. Prentice Hall international Series. (1998).         [ Links ]