• Author / Uploaded
• mdmorello

• Categories
• Análisis numérico
• Integral
• Conceptos matemáticos
• Matemática aplicada
• Algoritmos

Citation preview

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

-1-

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

f ( xi ) = ∑ an x n = a0 + a1 xi + a2 xi2 + ...an xin

(1)

i =0

y evaluando los coeficientes a0 .. ..an . Las funciones logarítmicas, hiperbólicas y elípticas son casos puntuales.

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 f ( x) = a0 + a1 x + a2 x 2 + a3 x 3 + a4 x 4

(2)

Tabla 1. Datos de la función Punto 0 1 2 3 4 5 x 1.0 1.1 1.2 1.3 1.4 1.5 fx 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: 5.7852 − 5.0000 = 7.8520 1.1 −1.0 6.7632 − 5.7852 = = 9.7800 1.2 −1.1

f [ x0 , x1 ] = f [ x1 , x2 ]

x f(x) 1 5.0000 1.1 5.7852 1.2 6.7632 1.3 7.9712 1.4 9.4512 1.5 11.250

[1]

[ 2]

[ 3]

[ 4]

∆ fi ∆fi ∆fi ∆ fi 7.852 9.640 6.200 2.000 9.780 11.50 7.000 2.000 12.08 13.60 7.800 14.80 15.94 17.98

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

La segunda diferencia dividida mediante los puntos (0), (1) y (2) es: f [ x0 , x1 , x2 ] =

Tabla 2. Diferencias divididas

9.7800 − 7.8520 = 9.6400 1.2 −1.0

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

-2-

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)); end end T

Para

encontrar

a0 , a1 , a2 , a3 y

los

coeficientes

a4 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 11.25];

9.4512

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

Figura 1 Gráfica del polinomio ajustado.

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: a0 = 3.000

a1 = −2.000 a2 = 5.000 a3 = −3.000 a4 = 2.000

En la figura 1 se muestran los datos suministrado junto con el 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 y = x 3 Los valores numéricos están arriba y la nomenclatura está debajo. Tabla 3 Diferencias finitas hacia adelante _______________________________________ [ 1] f i [ 2 ] f i [ 3] fi [ 4] X y = x 3 fi 1.1 1.331 0.397 0.072 0.006 0 1.2 1.728 0.469 0.078 0.006 1.3 2.197 0.547 0.081 0 1.4 1.744 0.631 0 0

-3-

1.5 3.375 0 0 0 carbonato de calcio se muestran en la figura _______________________________________ 2. La graficación de la curva se deja como [ 1] [ 2] [ 3] ejercicio para el lector fi fi fi x y x0 f ( x0 ) f [ x0 , x1 ] f [ x0 , x1 , x2 ] f [ x0 , x1 , x3 , x3 ] Se requiere: x1 f ( x1 ) f [ x1 , x2 ] x2 f ( x2 ) f [ x2 , x3 ] x3 f ( x3 ) f [ x3 , x4 ]

x4 f ( x4 ) f [ x4 , x5 ]

f [ x1 , x2 , x3 ]

f [ x2 , x3 , x4 ]

f [ x3 , x4 , x4 ]

f [ x1 , x2 , x3 , x4 ]

1)Encontrar la ecuación de la curva que

f [ x2 , x3 , x4 , x5 ] mejor se ajuste a los datos dados.

2)Calcular la velocidad de sedimentación para una concentración volumétrica de __________ __________ __________ ________ 2.5%. La función tabulada debe ajustarse con un polinomio f(x) de n-ésimo grado, que se expresa por

.

f ( x) = a0 + a1 ( x − x0 ) + a2 ( x − x0 )( x − x1 ) + a3 ( x − x0 )( x − x1 )( x − x2 ) +

P

an ( x − x0 )( x − x1 ). . . ( x − xn −1 ) Haciendo

h = x1 − x0 = x2 − x1

s = ( x − x0 ) / h

or derivación: s ( s −1) [ 2 ] fi + 2! s ( s −1)( s − 2)( s −3) [ 4 ] + fi + 4!

y = f ( x ) = f ( x0 ) + sf i [1] + s ( s −1)( s − 2) [ 3] fi 3! s ( s −1) ...( s − n +1) [ n ] fi n!

+

[1]

[2]

(3) [ 3]

f i , f , f i , son la primera, Siendo segunda y tercera diferencias finitas.

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

Figura. 2 Datos de sedimentación. Solución por Serie Potencias

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

-4-

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

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

Pol= [-0.0004 0.0121 3.8871 0.0004]; fx = polyval (Pol,2.5)

-0.1579 1.023 –

Obtenemos que: f(2.5) = 4.6783 g / cm 2 h

Coef = polyfit(x,y,7);

Solución por la fórmula de Newton

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: a 0 = 0.0004

a 4 = 1.023

a1 = 1.7060

a 5 = − 0.1579

a 2 = 3.8871

a 6 = 0.0121

a 3 = − 3.27

a 7 = − 0.00040

La ecuación de la curva es: f ( x) = a0 + a1 x + a2 x 2 + a3 x 3 + a4 x 4 + a5 x 5 + a6 x 6 + a7 x 7

(4)

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:

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 resumen en la Tabla 4.

finitas que se

Tabla 4. Diferencias finitas x 2.0 3.0 4.0 5.0 6.0 7.0 8.0

y f i [ 1] fi [ 2] 4.8 -0.55 4.25 –1.02 3.23 –0.36 2.87 –0.12 2.75 –0.050 2.70 -0.05 2.65 0

f i [ 3] -0.47 0.66 0.24 0.07 0 0

fi [ 4 ] fi [ 5] 1.13 -1.55 1.80 -0.42 0.25 -0.150 -0.17 0.10 -0.07

h =1.0 s=

2.5 − 2.0 = 0 .5 1 .0

Aplicando la ecuación (3)

0.5(0.5 −1) (−0.47 ) 2 0.5(0.5 −1)( 0.5 −2) + (1.13 ) (3)( 2) 0.5(0.5 −1)( 0.5 −2)( 0.5 −3) + (−1.55 ) ( 4)( 3)( 2)

f (2.5) = 4.8 + (0.5)( −0.50 ) +

= 4.7149 g / cm 2 h

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.

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

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;

0

2

3

..

x0

x1

x2

x3

 xn

f ( x)

f0

f1

f2

f3

 fn

n

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

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

f ( x) =

( x − x1 )( x − x2 )( x − xn ) ( f0 ) ( x0 − x1 )( x0 − x2 )( x0 − xn )

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

+

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

( x − x0 )( x − x2 )( x − xn ) ( f1 ) ( x1 − x0 )( x1 − x2 )( x1 − xn )

( x − x0 )( x − x1 )( x − x3 )( x − xn ) ( f2 ) ( x2 − x0 )( x2 − x1 )( x2 − x3 )( x2 − xn ) ( x − x0 )( x − x1 )( x − xn −1 ) ( fn ) ( xn − x0 )( xn − x1 )( xn − xn −1 ) (5)

La fórmula de Lagrange principalmente para : (1)

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

1

x

DE INTERPOLACION DE

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

-6-

(2)

se

usa

Calcular el valor de la variable independiente correspondiente a un valor dado de la función . 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,

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

Ti , º C

ρi ,

kg m3

1

2

94

205

371

929

902

860

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 ( 251 − 205 )( 251 − 371 ) (929 ) (94 − 205 )( 94 − 371 ) (251 − 94 )( 251 − 371 ) + (902 ) 205 − 94 (205 − 371 ) (251 − 94 )( 251 − 205 ) + (860 ) (371 − 94 )( 371 − 205 )

ρ( 251 º C ) =

En la Tabla 6 se muestran las densidades en kg / m3 , 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 ρ = f (T , C ) T (º C ) C (%)

5 20 40 70

10 1.0344 1.1453 1.3103 1.6923

30 1.0281 1.1335 1.2953 1.6014

60 1.0140 1.1153 1.2732 1.5753

100 0.9888 1.0885 1.2446 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: ρ(15 º C ) =

(15 −30 ) (15 −10 ) (1.3103 + (12953 ) (10 −30 ) (30 −10 )

=1.3066 kg / m3

= 890 .5 kg / m3

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

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

Procedimiento 5 x= [10 30];

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

y= [1.3103 1.2953]; xi = 15;

-7-

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

d = interp1(x,y,xi,’linear’) FORMULA DE INTERPOLACION HACIA DELANTE DE DERIVADAS DE NEWTON

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

1  [1] 1 [ 2 ] 1 [ 3] 1 [ 4 ] 1 [ 5] f i − f i + f i − f i + f1 h  2 3 4 5

(6) Derivaciones sucesivas se obtienen

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, dy / dt = 0.9 , donde t es el tiempo en horas.

5 [ 5]  [ 2]  [ 3] 11 [ 4 ]  f i − f i + 12 f i − 6 f i +     1  3 7  f ′′′( x) = 3  f i [ 3] − f i [ 4 ] + f i [ 5 ] −  h  2 4  1 f IV ( x) = 4 f i [ 4 ] − 2 f i [ 5] +  h f ′′( x) =

1 h2

[

(7 ) (8)

]

(9)

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 y = a + bx + cx 2 + dx 3 + ex 4

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.

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)

Por derivación, 397 ∑ ky 7 ∑ k 3 y  dy  −   = 1512 h 216 h  dx  0

Figura 3 Curva de velocidad de secado.

(10)

-8-

Se elabora una tabla de diferencias finitas. [ 1] f i [ 2 ] f i [ 3] f i [ 4 ] f i [ 5] x y fi

0.9 1.1 1.3 1.5

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.00023 1.2 0.1313 -0.001148 0.00340 0.11982 –0.0808 0.11174

0.1833 0.1459 0.1198 0.1071

0 1 2 3

0 0.1459 0.2396 0.3219

∑= −7.4112

0 0.1459 0.9584 2.8971

∑= −1.0752

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

1  0.0025 2(0.0003 )  dy  (397 )( −1.0752 ) 7( −7.4112 ) = − 0.01995 − +  dy     = −   0.1  2 6  dx t =0.9 (1512 )( 0.2) ( 216 )( 0.2)  dx t =0.9 = −0.2106 lbH 2O / lb sólido sec o 6(0.00007 ) 6( −0.00021 ) ]= − − 24 120 = − 0.2111 lb H 2O / lb sólido sec o

Solución por el método de DouglasAvakian. 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: fx = −0.0146 x 4 +0.119 x 3 −0.1453 x 2 −0.1958 x +0.40

(12)

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 DouglasAvakian. 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] ;

Tabla 8 Datos de y = f(x) x f(x) k ky 0.3 0.3313 -3 -0.9939 0.5 0.2798 -2 -0.5596 0.7 0.2291 -1 -0.2291

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

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

-8.9451 -2.2384 -0.2291

-9-

end end K s= 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 362 .1 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 ConstantinidesMostoufi, página 176 2 Series de Taylor. Un método de expandir funciones en series de potencias es utilizando las series de Taylor. El último

- 10 -

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

CALCULO DE INTEGRALES INTEGRACION NUMERICA

POR

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 f ( x ) = a + bx + cx 2

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 (13) La regla de Simpson sola es exacta para polinomios de primero y segundo grado. El b

h

∫ y d x= 3 [ y + 4( y + y +  + y a

0

1

3

2 ( y2 + y 4 +  yn − 2 ) + yn

- 11 -

n−1

]

)+

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 = h [ y−1 + yn +1 − 4( y0 + yn ) + 7( y1 + yn −1 ) 90 8( y2 + y4 +  + yn −2 ) + 8( y3 + y5 +  yn −3 ) ]

−

(14) Donde h = ∆xi y n ≥ 6 Mé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: h

∫ ydx = 2 [ y b

a

0

+ 2( y1 + y2 + y3 +  + yn −1 ) + yn

]

(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, N OG . 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 y − y* 0.015 ( x−1 )

0.025 ( x0 )

0.014328

0.010672

93.7

0.02250

0.012500

80.0

( y0 )

0.035 ( x1 ) ( y1 )

0.045 . 0.055 . 0.065 . 0.075 0.085 . 0.095 ( x7 )

0.031264 0.013736 0.040141 0.014859 0.049202 0.015798 0.058444 0.016556 0.067833 0.017167 0.077425 0.017575

72.8 . 67.3 . 63.3 60.4 58.25 . 56.9

( y7 )

0.105 ( xn ) 0.087127

0.017873

55.95

0.115 ( xn +1 ) 0.096819 0.018181

55.0

( yn ) ( yn +1 )

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: N OG = ∫

y ( 2)

y (1)

dy 0.01 = [93 .7 + 4(80 + 67 .3 * y −y 3

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

0.01 [115 .5 + 55 − 4(93 .7 + 55 .95 ) + 7(80 + 90 + 56 .7) −8(72 .8 + 63 .3 + 58 .25 ) +8(67 .3 + + 60 .4) ] = 0.000333 unidades de transf .

El error es relativamente pequeño. Por el método trapezoidal compuesto aplicamos la ecuación (15) N OG = ∫

y ( 2)

y (1)

=

0.01 [93 .7 + 2(80 + 72 .8 2

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

0.008658

115.5

( y−1 )

- 12 -

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.

Calcule la razón

W

W0

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 xD

0.549 0.691 0.793 0.806 0.902 0.928 0.950

xW

xD − xW

1 xD − xW

0.129 ( x0 ) 0.420 0.191 ( x1 ) 0.500 0.253 ( x2 ) 0.540 0.314 . 0.492 0.376 . 0.526 0.438 ( x5 ) 0.490 0.50 ( xn ) 0.450

2.38 ( y0 ) 2.00 ( y1 ) 1.85 ( y2 ) 1.83 . 1.90 . 2.04 ( y5 ) 2.22 ( yn )

Aplicando el método 1/3 de Simpson, tenemos xf dx w A=∫ = x0 x − x D w 0.0618 [ 2.38 + 4 ( 2.0 +1.83 + 2.04 ) 3 + 2 (1.85 + 2.04 ) + 2.22 ] = 0.739 W W ln = − 0.739 y = 0.4776 W0 W0

Fig 4 Gráfica de Xw vs 1/(XD- Xw) Por el método trapezoidal compuesto, tenemos que

- 13 -

A=∫

dx w 0.0618 = [2.38 + 2 (2.0 + x D − xw 2 + 1.85 + 1.83 + 1.90 + 2.04) + 2.22 ]

xw

x0

= 0.7366 W ln = −0.7366 W0

;

W = 0.4787 W0

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 ‘) break end dx = diff(x); if max(dx)-min(dx) > min(abs(x))/1000 error ( ‘ x no son equidistantes’) break end h= dx(1); if mod (puntos,2) == 0 precaución (‘Agregue numeros de intervalos’) n= puntos – 1; else n= puntos; end y1 = y(2 : 2 : n – 1); y2 = y(3 : 2 : n –2 ); A= (h/3)*(y(1) + 4*sum(y1) + 2* sum(y2) + y(n)) ;

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

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.

−

−

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

- 14 -