METODOS NUMERICOS PARA INGENIEROS QUIMICOS CON MATLAB (I)

Francisco Muñoz Paba M.Sc.E-mail: f31paba @ yahoo.com. Departamento de Ingeniería Química, Grupo de

Simulación y Control de Procesos. Universidad del Atlántico, Barranquilla, Colombia

INTRODUCCION

Muchos planteamientos matemáticos sobre situaciones problémicas, en procesos químicos , son de difícil solución analítica y hacen que el ingeniero químico tenga que recurrir a los métodos numéricos para encontrar una respuesta a sus casos de estudio. Una necesidad muy frecuente es la de representar un conjunto de datos experimentales tomados en forma discreta ajustados a una expresión analítica que permita de forma mas fácil la estimación de, por ejemplo, valores intermedios, sumatorias o integrales y variaciones o razones de cambio entre ellos. El desarrollo de los métodos numéricos , la certidumbre de sus resultados y la posibilidad de ejecutarlos con la ayuda de códigos por computador hacen de ellos un recurso que ofrece ventajas con respecto a los métodos analíticos. En ésta revisión se presentan algunos métodos de ajuste de datos a ecuaciones con ejemplos a la ingeniería química que se resuelven con los pro cedimientos explicados y con la ayuda de un computador mediante la construcción de instrucciones cortas codificadas con MATLAB.

AJUSTE DE CURVAS PARA FUNCIONES POLINOMICAS.

Muchas funciones matemáticas incluyen términos como logarítmicos, exponenciales

o trigonométricos que las hacen de un manejo complejo. Una alternativa para afrontar tal dificultad la ofrecen los métodos numéricos permitiendo que una función se pueda expresar por otra equivalente en cuanto a la correspondencia entre la variable independiente y el valor de la función pero mas sencilla y, por lo tanto, de mas fácil manipulación. Lo anterior es lo que se conoce como ajuste de curvas, interpolación o cálculo de la ecuación de una curva. A continuación se muestra el método de ajuste de curvas a un polinomio como una Serie de Potencias o mediante procedimiento de interpolación como el de Newton y Lagrange. .SERIE DE POTENCIAS.

Prácticamente todas las funciones matemáticas se pueden expresar como un polinomio de grado n, es decir, mediante una expresión en serie de potencias.

Es más fácil encontrar el valor numérico de una función expandiéndola en una serie de potencia polinomial como la ecuación (1):

∑=

+++==n

i

ninii

nni xaxaxaaxaxf

0

2210 )1(...)(

y evaluando los coeficientes naa ....0 .

Las funciones logarítmicas, hiperbólicas y elípticas son casos puntuales.

- 1 -

Las series de potencias pueden usarse para ajustar un conjunto de datos tomando un número suficiente de términos. El número de términos está dado por el siguiente teorema:

Sí las enésimas diferencias divididas de una función tabulada son constantes cuando los valores de la variable independiente son tomadas en progresión aritmética, la función es un polinomio de grado n.

Ejemplo 1

44

33

2210)( xaxaxaxaaxf ++++= (2)

Tabla 1. Datos de la función

Punto 0 1 2 3 4 5x 1.0 1.1 1.2 1.3 1.4 1.5fx 5.000 5.785 6.763 7.971 9.451 11.25

Elabore una tabla de diferencias divididas determine los coeficientes del polinomio dado por la ecuación (2).

Las primeras diferencias divididas mediante los puntos (0), (1) y (1), (2), respectivamente, son:

[ ]

[ ] 7800.91.12.1

7852.57632.6,

8520.70.11.1

0000.57852.5,

21

10

=−−=

=−−=

xxf

xxf

La segunda diferencia dividida mediante los puntos (0), (1) y (2) es:

[ ] 6400.90.12.1

8520.77800.9,, 210 =

−−=xxxf

La Tabla 2 muestra los resultados correspondientes hasta la cuarta diferencia dividida.

Tabla 2. Diferencias divididas

x f(x) [ ] [ ] [ ] [ ]4321

iiii ffff ∆∆∆∆1 5.0000 7.852 9.640 6.200 2.000 1.1 5.7852 9.780 11.50 7.000 2.000 1.2 6.7632 12.08 13.60 7.800 1.3 7.9712 14.80 15.94 1.4 9.4512 17.98 1.5 11.250

Debe notarse que todas las diferencias divididas de cuarto orden tienen el mismo valor, independientemente de los argumentos que se usen para su cálculo, por lo tanto , la ecuación(2) se puede escribir en forma de series de potencias como un polinomio de cuarto orden.Para realizar los cálculos de diferencias divididas puede usarse el siguiente procedimiento codificado con MATLAB:

Procedimiento 1

x=[1.0 1.1 1.2 1.3 1.4 1.5];fx=[5.000 5.7852 6.7632 7.9712 9.4512 11.25];M=6; N= M-1;

for i=1:N T(i,1)= (fx(i+1)- fx(i))/(x(i+1)-x(i)); End

for j=2 :N for i=j : N T(i,j)= (T(i,j-1)- T(i-1,j-1))/(x(i+1)-x(i-j+1)); endend

T

- 2 -

Para encontrar los coeficientes 43210 ,,, ayaaaa del polinomio en series

de potencia de la ec(2), se escribe el siguiente procedimiento codificado con MATLAB:

Procedimiento 2

x=[1.0 1.1 1.2 1.3 1.4 1.5];fx =[5.00 5.7852 6.7632 7.9712 9.4512 11.25];

plot(x,fx,’o’)

a = polyfit (x, fx, 4);Y= polyval (a, x);

fprintf ( ‘ a0=%8.5f\n a1=%9.6f\n a2=%9.6f\n … a3=%9.6f\n a4=%9.6f\n’,a(5),a(4),a(3),a(2),a(1)) plot(x,fx,’o’,x,Y,’-‘)

Donde se obtiene que:

000.2

000.3

000.5

000.2

000.3

4

3

2

1

0

=−=

=−=

=

a

a

a

a

a

En la figura 1 se muestran los datos suministrado junto con el polinomio ajustado

Figura 1 Gráfica del polinomio ajustado.

FORMULA DE NEWTON EN DIFERENCIAS FINITAS HACIA ADELANTE.

La fórmula necesita una tabla de valores y0, y1, y2, .......yn para valores equidistantes x0, x1, x2, ..xn de la variable independiente x.Para usar la fórmula de Newton en diferencias finitas es de mucha ayuda construir una tabla de diferencias finitas.La tabla 3 es una tabla de diferencias finitas, para 3xy = Los valores numéricos están arriba y la nomenclatura está debajo.

Tabla 3 Diferencias finitas hacia adelante_______________________________________X 3xy = [ ] [ ] [ ] [ ]4321

iiii ffff

1.1 1.331 0.397 0.072 0.006 01.2 1.728 0.469 0.078 0.0061.3 2.197 0.547 0.081 01.4 1.744 0.631 0 0

- 3 -

1.5 3.375 0 0 0_______________________________________ x [ ] [ ] [ ]321

iii fffy

[ ] [ ] [ ]33102101000 ,,,,,,)( xxxxfxxxfxxfxfx

[ ] [ ] [ ][ ] [ ] [ ][ ] [ ][ ]

______________________________________

,)(

,,,)(

,,,,,,)(

,,,,,,)(

5444

4434333

54324323222

43213212111

xxfxfx

xxxfxxfxfx

xxxxfxxxfxxfxfx

xxxxfxxxfxxfxfx

La función tabulada debe ajustarse con un polinomio f(x) de n-ésimo grado, que se expresa por

hxxsxxxxhHaciendo

xxxxxxa

xxxxxxa

xxxxaxxaaxf

nn

/)(

)(..).)((

))()((

))(()()(

01201

110

2103

102010

−=−=−=−−−+−−−

+−−+−+=

−P

or derivación:[ ] [ ]

[ ] [ ]

[ ] )3(!

)1.(..)1(

!4

)3)(2)(1(

!3

)2)(1(

!2

)1()()(

43

210

ni

ii

ii

fn

nsss

fssss

fsss

fss

sfxfxfy

+−−

+−−−+−−+

+−++==

Siendo [ ] [ ] [ ] ,,, 321ii fff son la primera,

segunda y tercera diferencias finitas.

La fórmula es útil solo para valores puntuales, no para la ecuación de la curva total

Ejemplo 2

La velocidad de sedimentación de una suspensión, se relaciona con la concentración volumétrica del sedimento. Los datos y la curva para la sedimentación de una suspensión de precipitado de

carbonato de calcio se muestran en la figura 2. La graficación de la curva se deja como ejercicio para el lectorSe requiere:

1)Encontrar la ecuación de la curva que mejor se ajuste a los datos dados.

2)Calcular la velocidad de sedimentación para una concentración volumétrica de 2.5%.

.

Figura. 2 Datos de sedimentación.

Solución por Serie Potencias

Para encontrar el polinomio en serie de potencias, suponemos un polinomio de séptimo grado que se encuentra mediante el siguiente procedimiento codificado con MATLAB

Procedimiento 3

- 4 -

x= [ 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0];y= [0 3.2 4.8 4.25 3.23 2.87 2.75 2.70 2.65];plot(x,y,’o’)

Coef = polyfit(x,y,7);

X=1:0.1:8;

Y= polyval (Coef,X);

plot(x,y,’o’,X,Y)

fprintf (‘ a0=%9.6f\n a1=%9.6f\n a2=%9.6f\n …a3=%9.6f\n a4=%9.6f\n a5=%9.6f\n a6=%9.6f\n

a7=%9.6f\n’,a(8),a(7),a(6),a(5),a(4),a(3),a(2),a(1))

Donde los coeficientes del polinomio de séptimo grado son:

00040.027.3

0121.08871.3

1579.07060.1

023.10004.0

73

62

51

40

−=−===

−====

aa

aa

aa

aa

La ecuación de la curva es:

)4(

)(7

7

66

55

44

33

2210

xa

xaxaxaxaxaxaaxf

+

++++++=

La velocidad másica de concentración para una concentración volumétrica de 2.5%, se halla sustituyendo los coeficientes encontrados con el procedimiento 3 en la ecuación (4) para un valor de x =2.5.

Empleando los siguientes comandos de MATLAB:

Pol= [-0.0004 0.0121 -0.1579 1.023 –3.88710.0004];fx = polyval (Pol,2.5)

Obtenemos que: f(2.5) = 4.6783 hcmg 2/

Solución por la fórmula de Newton

Este problema se puede resolver utilizando la fórmula de Newton en diferencias finitas. Este método es válido solamente para calcular valores puntuales de la función y no para calcular la ecuación de la curva, por consiguiente, se calcula solamente el valor de la función para un valor de x = 2.5.

Se calculan las diferencias finitas que se resumen en la Tabla 4.

Tabla 4. Diferencias finitas

[ ] [ ] [ ] [ ] [ ]54321iiiii fffffyx

2.0 4.8 -0.55 -0.47 1.13 -1.55 1.803.0 4.25 –1.02 0.66 -0.42 0.25 -0.1504.0 3.23 –0.36 0.24 -0.17 0.105.0 2.87 –0.12 0.07 -0.076.0 2.75 –0.050 07.0 2.70 -0.05 08.0 2.65 0

h =1.0

5.00.1

0.25.2 =−=s

Aplicando la ecuación (3)

hcmg

f

2/7149.4

)55.1()2)(3)(4(

)35.0)(25.0)(15.0(5.0

)13.1()2)(3(

)25.0)(15.0(5.0

)47.0(2

)15.0(5.0)50.0)(5.0(8.4)5.2(

=

−−−−+

−−+

−−+−+=

Aunque la cuarta diferencia finita no es

- 5 -

constante, el resultado obtenido es satisfactorio. Es evidente a partir de éste ejemplo que tanto el polinomio en serie de potencias como la fórmula de Newton son bastante aproximadas al valor medido que es de 4.700.

Los cálculos anteriores se pueden realizar con el siguiente procedimiento codificado con MATLAB.

Procedimiento 4

x= [2.0 3.0 4.0 5.0 6.0 7.0 8.0];y= [4.8 4.25 3.23 2.87 2.75 2.70 2.65];N=7;

for i =1: N-1 f(i,1) = y(i+1) –y(i);end

for j=2: N-1 for i=j: N-1 f(i,j) = f(i,j-1) – f(i-1,j-1); endend

f

h= 1.0 ; xi = 2.5;s = (xi – x(1))/h ;yi = y(1) + s*f(1,1) + s*(s-1)/2*f(2,2) + s*(s-1)*(s-2)/(3*2)*f(3,3) + s*(s-1)*(s-2)*(s-3)/(4*3*2)*f(4,4) ;

fprintf(‘\n\n Resultado: 4º grado f(%4.2f) =... %6.2f \ n’, xi,yi )

FORMULA DE INTERPOLACION DE LAGRANGE

Muchas fórmulas de interpolación son aplicables solo cuando los valores de la variable independiente son dados en intervalos equidistantes. La fórmula de

Lagrange no tiene ésta limitación, pero solo utiliza datos que sean necesarios para aproximarse al valor correcto. Los datos donde los valores de x no son equidistantes, a menudo son resultados de observaciones experimentales o de análisis de datos históricos.Supóngase que se tiene una tabla de datos con cuatro pares de valores x y f(x)

i 0 1 2 3 . . n

n

n

fffffxf

xxxxxx

3210

3210

)(

Estos cuatro pares de datos es posible ajustarlos a una función cúbica. La fórmula de Lagrange para un polinomio de n-ésimo grado es

)5(

)()())((

)())((

)()())()((

)())()((

)()())((

)())((

)()())((

)())(()(

110

110

22321202

310

112101

20

002010

21

nnnnn

n

n

n

n

n

n

n

fxxxxxx

xxxxxx

fxxxxxxxx

xxxxxxxx

fxxxxxx

xxxxxx

fxxxxxx

xxxxxxxf

−

−

−−−−−−++

−−−−−−−−+

−−−−−−+

−−−−−−=

La fórmula de Lagrange se usa principalmente para :

(1) Calcular el valor de la variable independiente correspondiente a un valor dado de la función .

(2) Calcular cualquier valor de una función, cuando los valores dados de la variable independiente no son equidistantes.

Además de que la fórmula de Lagrange es tediosa, tiene una limitación muy seria,

- 6 -

cuando los valores no son tan cercanos unos a otros, los resultados tienden a ser indeseables. Sin embargo puede utilizarse cuando sea imposible utilizar otro método.Los siguientes ejemplos aplican la fórmula de Lagrange.

Ejemplo 3

Se desea estimar la densidad de una sustancia a una temperatura de 251º C a partir de los siguientes datos experimentales que se dan en la Tabla 5.

Tabla 5 Datos de Temperatura-Densidad

i 0 1 2

860902929,

37120594º,

3m

kg

CT

i

i

ρ

Como se dispone de tres datos, el orden de la fórmula de Lagrange es 2 y el cálculo de la densidad a 251 es dado por

3/5.890

)860()205371)(94371(

)205251)(94251(

)902()371205(94205

)371251)(94251(

)929()37194)(20594(

)371251)(205251()º251(

mkg

C

=−−−−+

−−−−+

−−−−=ρ

El siguiente procedimiento codifcado con MATLAB realiza los cálculos anteriores.

Procedimiento 5

X = [94 205 371];Y = [929 902 860];Xi= 251;Densidad =interp1(X,Y,Xi,’cubic’)

En la Tabla 6 se muestran las densidades en 3/mkg , de soluciones acuosas de ácido

sulfurico de diferentes concentraciones en % para un conjunto de temperaturas en ºC. Se desea calcular la densidad de una solución de ácido sulfúrico a una concentración del 40% y a una temperatura de 15 ºC.

Tabla 6 Tabulación de una función de dos variables ),( CTf=ρ

(%)

)(º

C

CT 10

30

60

100

5 1.0344 1.0281 1.0140 0.9888 20 1.1453 1.1335 1.1153 1.0885 40 1.3103 1.2953 1.2732 1.2446 70 1.6923 1.6014 1.5753 1.5417

Para una función polinómica de dos variables como éste caso, se puede aplicar la fórmula de Lagrange , tomando los datos de las densidades a una concentración del 40% y la temperatura como la variable independiente. El orden de la fórmula es de 1 y el cálculo de la densidad mediante la fórmula de Lagrange es:

3/3066.1

)12953()1030(

)1015(3103.1(

)3010(

)3015()º15(

mkg

C

=−−+

−−=ρ

El siguiente procedimiento codificado con MATLAB realiza los cálculos anteriores.

Procedimiento 6

x= [10 30];

y= [1.3103 1.2953];

xi = 15;

- 7 -

d = interp1(x,y,xi,’linear’)

FORMULA DE INTERPOLACION HACIA DELANTE DE DERIVADAS DE NEWTON

La fórmula de diferenciación de Newton para una estimación de f’(x) se obtiene

[ ] [ ] [ ] [ ] [ ]

−+−+−=′ 5

14321

5

1

4

1

3

1

2

11)( fffff

hxf iiii

(6)Derivaciones sucesivas se obtienen

[ ] [ ] [ ] [ ]

[ ] [ ] [ ]

[ ] [ ][ ] )9(21

)(

)8(4

7

2

31)(

)7(6

5

12

111)(

544

5433

54322

+−=

−+−=′′′

+−+−=′′

iiIV

iii

iiii

ffh

xf

fffh

xf

ffffh

xf

METODO DE DOUGLAS-AVAKIAN

Este método usa un polinomio de cuarto orden que se ajusta a siete puntos equidistantes por el método de mínimos cuadrados. El polinomio es

432 exdxcxbxay ++++=

Estos puntos son espaciados en intervalos iguales con las coordenadas escogidas, tal que, en x = 0 se encuentra el punto central de los siete. Los siete valores de x pueden escribirse como –3h, -2h, -h, 0, 2h y 3h. Por derivación,

)10(216

7

1512

397 3

0 h

yk

h

ky

dx

dy ∑∑ −=

Donde k representa el coeficiente de h en los valores de x, por ejemplo –3, -2 , -1, 0, 1, 2, 3.

Ejemplo 5

Una pasta de material cristalino se seca con aire, que se hace fluir por encima de ella . Para diseñar el sistema de secado, se obtuvieron los datos experimentales que se muestran en la figura 3. A partir de esto, calcule la velocidad de secado en 0.9h ,es decir, 9.0/ =dtdy , donde t es el tiempo en horas.

Figura 3 Curva de velocidad de secado.

Solución por la Fórmula de la derivada de Newton

Se divide parte de la curva en cinco subdivisiones comenzando en t=0.9 hora, como muestra la figura 3 y se elabora la Tabla de diferencias finitas ( Tabla 7)

- 8 -

Se elabora una tabla de diferencias finitas.x y [ ] [ ] [ ] [ ] [ ]54321

iiiii fffff 0.9 0.18335 -0.01995 0.0025 0.0003 0.00007 -.00021 1.0 0.1634 -0.01745 0.00280 0.00037 –0.00014 1.1 0.14595 –0.1465 0.00317 0.000231.2 0.1313 -0.001148 0.00340 0.11982 –0.08080.11174

+−−=

= 6

)0003.0(2

2

0025.001995.0

1.0

1

9.0tdx

dy

]120

)00021.0(6

24

)00007.0(6 −−− =

= osólidolbOHlb sec/2111.0 2−

Solución por el método de Douglas-Avakian.

Primero se preparó la Tabla 8, a partir del polinomio de cuarto orden ajustado los datos experimentales. Ecuación (11) con ayuda de Matlab Para hallar el polinomio ajustado de cuarto orden se utilizó MATLAB, obteniéndose el siguiente polinomio:

40.01958.0

1453.0119.00146.0 234

+−−+−=

x

xxxfx

(12)

Tabla 8 Datos de y = f(x)x f(x) k ky 3k y0.3 0.3313 -3 -0.9939 -8.94510.5 0.2798 -2 -0.5596 -2.23840.7 0.2291 -1 -0.2291 -0.2291

0.9 0.1833 0 0 01.1 0.1459 1 0.1459 0.1459 1.3 0.1198 2 0.2396 0.95841.5 0.1071 3 0.3219 2.8971 ∑ −= 0752.1

∑ −= 4112.7

La velocidad de secado se calcula con la ecuación (10), de la siguiente manera

osólidolbOlbH

dx

dy

t

sec/2106.0

)2.0)(216(

)4112.7(7

)2.0)(1512(

)0752.1)(397(

2

9.0

−=

−−−=

=

Comparando los resultados encontramos un valor de –0.2111 por el método de Newton y –0.2106 por el método de Douglas-Avakian. El valor medido es de –0.21. El método de Douglas-Avakian se basa en el método de mínimos cuadrados, por lo tanto, es un método inseguro.El siguiente procedimiento codificado con MATLAB realiza los cálculos anteriores donde se aplica ell método de Douglas-Avakian.

Procedimiento 7

function y = Douglas(y,k)x = [0.9 1.0 1.1 1.2 1.3 1.4] ;fx =[0.18335 0.1634 0.14595 0.1313 0.11982 0.11174] ;pol = polyfit (x, fx, 4);

xi = [0.3 0.5 0.7 0.9 1.1 1.3 1.5] ;

yi = polyval(pol,xi) ;k = [-3 -2 -1 0 1 2 3] ;y = yi ;

for i = 1 : 7 for j = 1 : 7 K(i,1) = k(i)*y(i); K(i,2) = k(i)^3*y(i);

- 9 -

endend

Ks= sum (K)

Derivada= 397*s(1)/(1512*0.2) . . . - 7*s(2)/(216*0.2) OTROS METODOS PARA AJUSTE DE CURVAS

Método de mínimos cuadrados. Este método se basa en la suposición, que la mejor curva representativa es aquella para la cual la suma de los cuadrados de los residuos (errores) es un mínimo. Los residuos son elevados al cuadrado para eliminar lo que concierne a su signo. Consultar el libro de Nieves-Domínguez página 1.362 Este método es mucho más complicado para polinomios de mayor grado y se usa para polinomios no mayor de segundo grado. Es menos seguro que la Fórmula interpolación de Newton y debe emplearse para correlacionar o encontrar el “mejor ajuste” de un conjunto de datos experimentales.

Fórmula de diferencia central de Stirling. Dos formas de la fórmula de Newton se usan para la interpolación cercana al comienzo y cercana al final de un conjunto de datos tabulados. La fórmula de Stirling es particularmente disponible para valores interpolados cercanos a la mitad de un conjunto de datos tabulados. Este método está explicado en el libro de Constantinides-Mostoufi, página 2176

Series de Taylor. Un método de expandir funciones en series de potencias es utilizando las series de Taylor. El último

término en la serie es el residuo o tamaño de error después de n términos y por lo tanto, la serie de Taylor tiene una ventaja sobre otros métodos, por que puede programarse en un computador, de tal manera que los términos se pueden agregar automáticamente hasta que el último término (término error) sea menor que el limite especificado. Una nota de precaución en el uso de todos los métodos de ajuste de curvas debe expresarse. La exactitud de la correlación entre los puntos de datos (xi,yi) se debe chequear.

BIBLIOGRAFIA

1. Nieves A y Domínguez F. Métodos numéricos aplicados a la ingeniería . 2ª Edición CECSA 2002.

2. Constantinides A y Mostoufi N Numerical methods for chemical engineers with MATLAB applications 1ª Edición Prentice-Hall 1999.

3. Gerald C.F y Wheatley P.O Analisis numérico con aplicaciones. 7ª Edición Pearson Educación 2000.

4. Nakamura S. Análisis numérico y visualización gráfica con MATLAB 1ª Edición Pearson Educación 1997.

- 10 -

CALCULO DE INTEGRALES POR INTEGRACION NUMERICA

El proceso de calcular el valor de una integral definida a partir de un conjunto de valores numéricos del integrando recibe el nombre de integración numérica. El integrando se representa por una fórmula de interpolación y la fórmula se integra entre los limites deseados.

Método de Simpson. Este método se puede resumir diciendo que se basa en la conexión de los puntos (xi,yi) por una series de parábolas.Las funciones de éste tipo son polinomios de segundo grado

2)( cxbxaxf ++= Hay un error inherente, por supuesto, si el polinomio es mayor de segundo grado. La fórmula final de la ecuación para la Regla 1/3 de Simpson es

]∫

+++

+++++=

−

−

b

a

nn

n

yyyy

yyyyh

y d x

)(2

)(4[3

242

1310

(13)La regla de Simpson sola es exacta para polinomios de primero y segundo grado. El

- 11 -

grado de la función es desconocida en muchas aplicaciones, por consiguiente , se debe calcular el error. El error se calcula por la siguiente ecuación:Error =

[ )(7)(490 11011 −+− +++−+− nnn yyyyyyh

])(8)(8 353242 −− ++++++ nn yyyyyy (14)Donde h = 6≥∆ nyxiMétodo trapezoidal compuesto. Consiste en dividir el intervalo[a , b] en n subintervalos y aproximar cada uno por un polinomio de primer grado, luego se aplica la fórmula trapezoidal a cada subintervalo y se obtiene el área de cada trapezoide, de tal modo que la suma de todas ellas da la aproximación al área bajo la curva de la función. La forma final de la ecuación para el método trapezoidal compuesto es:

[ ]∫ ++++++= −

b

a nn yyyyyyh

ydx )(22 13210

(15)Los siguiente dos ejemplo ilustran estos dos métodos. Una torre empacada absorbe un gas A de un gas de combustión. El gas de entrada a la torre contiene 10.5% molar de A y el gas de salida contiene 2.5% molar de A. Calcule el número de unidades de transferencia necesarias, OGN . Los datos se muestran en la tabla 6.

Tabla 6 Datos para el problema de unidades de transferencia. Datos Calculados de los datos y y* y – y*

*

1

yy −

0.015 )( 1−x 0.006342 0.008658 115.5 )( 1−y

0.025 )( 0x 0.014328 0.010672 93.7 )( 0y

0.035 )( 1x 0.02250 0.012500 80.0 )( 1y

0.045 . 0.031264 0.013736 72.8 .0.055 . 0.040141 0.014859 67.3 .0.065 . 0.049202 0.015798 63.3 0.075 0.058444 0.016556 60.40.085 . 0.067833 0.017167 58.25 .0.095 )( 7x 0.077425 0.017575 56.9

)( 7y

0.105 )( nx 0.087127 0.017873 55.95)( ny

0.115 )( 1+nx 0.096819 0.018181 55.0)( 1+ny

y* = Composición en equilibrio.

Primero resolvemos el problema aplicando el método 1/3 de Simpson. Suponiendo que la película gaseosa es la controlante, tenemos:

∫ ++=−

=)2(

)1( *3.6780(47.93[

3

01.0y

yOG yy

dyN

+ 60.4 +56.9) + 2 (72.8 + 63.3 + 58.25) + 55.95] = 5.3225 unidades de transf.Error =

] .000333.0)4.60

3.67(8)25.583.638.72(8)7.56

80(7)95.557.93(4555.115[90

01.0

transfdeunidades=+++++−+

+++−+−

El error es relativamente pequeño.

Por el método trapezoidal compuesto aplicamos la ecuación (15)

∫ ++==)2(

)1(8.7280(27.93[

2

01.0y

yOGN

+ 67.3 + 63.3 + 60.4 + 58.25 + 56.9) + 55.95 ] = 5.3377 unidades de transf.

Consideremos ahora una columna de destilación discontinua que contiene una mezcla de 50% molar de A en B, se destila hasta que la fracción molar de A en el calderin sea menor que 0.20.

- 12 -

Calcule la razón 0W

W Los datos se muestran

en la tabla 7. y se grafican en la figura 4.

Tabla 7 Datos para el problema de la columna de destilación discontinúa

WDWDWD xxxxxx

−− 1

0.549 0.129 )( 0x 0.420 2.38 )( 0y

0.691 0.191 )( 1x 0.500 2.00 )( 1y

0.793 0.253 )( 2x 0.540 1.85 )( 2y

0.806 0.314 . 0.492 1.83 .0.902 0.376 . 0.526 1.90 .0.928 0.438 )( 5x 0.490 2.04 )( 5y

0.950 0.50 )( nx 0.450 2.22 )( ny

Aplicando el método 1/3 de Simpson, tenemos

=−

= ∫ fx

xwD

w

xx

dxA

0

]4776.0739.0ln

739.022.2)04.285.1(2

)04.283.10.2(438.2[3

0618.0

00

=−=

=+++

+++

W

Wy

W

W

Fig 4 Gráfica de Xw vs 1/(XD- Xw)

Por el método trapezoidal compuesto, tenemos que

∫ ++=−

= wx

xwD

w

xx

dxA

0

0.2(238.2[2

0618.0

+ 1.85 + 1.83 + 1.90 + 2.04) + 2.22 ] = 0.7366

4787.0;7366.0ln00

=−=W

W

W

W

Se observa que los dos resultados son casi iguales debido a que el polinomio es de orden 3. El siguiente guión de MATLAB hace los cálculos de los dos problemas dados anteriormente.

x = input(‘Introduzca los valores de x = ’);y = input(‘Introduzca los valores de y =’);Area_1= trapz(x,y);Area_2= Simpson(x,y);fprintf (‘\ n Area_1(Método trapezoidal)=%9.4f’,Area_1)fprintf(‘\ n Area_2(Método 1/3 de Simpson)=%9.4f’,Area_2)

function A=Simpson ( x, y)puntos = length(x);if length(y) ~= puntos error(‘x y y no son de la misma longitud ‘) breakenddx = diff(x);if max(dx)-min(dx) > min(abs(x))/1000error ( ‘ x no son equidistantes’) breakendh= dx(1);if mod (puntos,2) == 0precaución (‘Agregue numeros de intervalos’) n= puntos – 1;else n= puntos;endy1 = y(2 : 2 : n – 1);y2 = y(3 : 2 : n –2 );A= (h/3)*(y(1) + 4*sum(y1) + 2* sum(y2) + y(n)) ;

- 13 -

if n ~= puntos A = A + (y(puntos) + y(n))* h/2;end.

BIBLIOGRAFIA

5. Nieves A y Domínguez F. Métodos numéricos aplicados a la ingeniería . 2ª Edición CECSA 2002.

6. Constantinides A y Mostoufi N Numerical methods for chemical

engineers with MATLAB applications 1ª Edición Prentice-Hall 1999.

7. Gerald C.F y Wheatley P.O Analisis numérico con aplicaciones. 7ª Edición Pearson Educación 2000.

8. Nakamura S. Análisis numérico y visualización gráfica con MATLAB 1ª Edición Pearson Educación 1997.

−

−

- 14 -