• Author / Uploaded
• Jimenena2

Citation preview

METODO GRAFICO RIGUROSO DE DESTILACION BINARIA Parra Milián, R. Universidad Nacional Pedro Ruiz Gallo Facultad de Ingeniería Química e Industrias Alimentarias

RESUMEN El presente trabajo se desarrolló, con la finalidad de aplicar el modelo UNIFAC en la destilación binaria de metanol-agua a 760 mmHg utilizando el método grafico riguroso de Ponchon y Savarit. Las ecuaciones fueron transcritas en un lenguaje de programación conocida MATLAB, siendo ejecutadas al mismo tiempo. Se determino que tanto el modelo UNIFAC así como el método grafico pueden elaborarse sin restricción alguna. PALABRAS CLAVES Equilibrio líquido vapor. Ley de Raoult modificada, Método grafico Ponchon-Savarit ABSTRACT This research was made with the aim to apply UNIFAC model into methanol-water binary distillation at 760 mmHg using Ponchon and Savarit rigorous graphic method. The equations were written in a programming language called MATLAB, being executed simultaneously. It was determined that both the UNIFAC model and Ponchon and Savarit rigorous graphic method may be made without restriction. KEY WORDS Liquid vapor equilibrium, Modified Law of Raoult, Ponchon-Savarit graphic method. INTRODUCCION El método de destilación grafica simplificada de McCabe-Thiele, asume que los flujos molares tanto de líquido descendente como vapor ascendente permanecen constantes en cualquier sección de la columna. De esta manera realiza solo el intercambio de materia en cada etapa de la columna, despreciando el intercambio de calor entre aquellos. El método de destilación grafica riguroso de Ponchon y Savarit toma en cuenta la variación de los flujos molares de líquido descendente como vapor ascendente. La única consideración que toma en cuenta es que la columna opere adiabáticamente. REVISION DE LITERATURA EVL a presiones bajas y moderadas. Cuando se opera en sistema de presiones bajas o moderadas, para el cálculo tanto del punto de burbuja como de rocío se realizan con la ley de Raoult: Ley de Raoult

P

P

Desafortunadamente esta ley nos introduce ciertos errores en el cálculo del equilibrio liquido-vapor, por lo cual se introduce un coeficiente , que es el coeficiente de actividad, dando de esta manera la ley de Raoult modificada Ley de Raoult modificada

P

P

Se han desarrollado modelos que pueden predecir los coeficientes de actividad con precisión (Van Laar, Margules, Wilson, NRTL, UNIQUAC y UNIFAC). Este último se basa en la interacción que tienen los grupos funcionales, un conjunto de sumatorias que su ejecución a mano es laborioso. Al no obtener laos parámetros binarios de los demás modelos se opto por la aplicación del modelo UNIFAC. La expresión matemática se presenta en el texto de Prausnitz, siendo complejo en su transcripción, se despejo término a término y se llego a la forma simplificada que presenta Van Ness. Método Grafico de Ponchon y Savarit QD

El método grafico de Ponchon y Savarit, es un método grafico riguroso que combina tanto los balances de materia como los balances entálpicos, y que podrá ser aplicable a cualquier sistema binario a separar, sin que se tenga restricción alguna.

D xD VM+N

LD M+N

n+1

Zona de Rectificación Balance de materia:

V

L

Vn

D

n

Balance de materia del componente más volátil: Balance de energia: Donde

Vy

L

DM

Q

VH

L

VH

L

Ln-1

x

Vn-1

Dx

h

n-1

Q

Dh

Ln

Dh

a+1 M+N

F xF

h

DM

a

a-1

Despejando tenemos: ! "#

V m+1

La-1

Vm

Lm+1

m+1

"#

m

Zona de Agotamiento Balance de materia: L V R Balance de materia del componente más volátil: Rx&

Balance de energia: Donde

L% x%

V% y%

1

L% h% Q& V% H% Rh& RM& Rh& ' Q & L% h% V% H% RM&

Despejando tenemos: (

(

!)"#

)

)

)"#

+

Calculo de MD:

M

*

Calculo de MR:

M&

H. ' *

Finalmente:

m-1

Q

VR

L1

QR

)

)

+

1- H ' * - h /

/

- 0x . ' x & 1

D0M ' h 1 , Q&

, donde

R0h& ' M& 1

R= LD/D

R xR

Inyección directa de vapor: Si uno de los componentes de la mezcla binaria a separar fuese agua, se es factible la aplicación directa de vapor en sustitución de la caldera. Para la determinación de vapor de agua se despejan las ecuaciones mostradas, métodos diferentes realizan la determinación gráficamente. Balance de materia F V D R B. de materia del componente más volátil Fx. Dx Rx& Balance de energia FH. Vλ DM Rh& 0 / λ1 ( 0!( λ1 / λ1 ( 0!( λ1

D

Despejando:

F *0

Calculo del punto (xc,Mc) /

3

/

3

!( λ (

Igualando y despejando tenemos:

x5 M5

Donde:

-

3

4

3

λ 4 7 89/ ; : 8:/ 7 89/ 6 ; : 8:/

λ 6

< 8λ 6 ( ; :(

λ

*

!( λ (

- x5

RESULTADOS Se elaboró el archivo Meth_Ponchon_Savarit.m con la capacidad de leer tanto datos calculados como experimentales para el cálculo de etapas por el método de Ponchon-Savarit. Se utilizo con fines teóricos datos calculados mediante la ecuación de Raoult modificada, calculando una nube de 10000 puntos utilizando el modelo UNIFAC para el cálculo de los coeficientes de actividad, siendo el caso más sencillo en estudio el sistema binario metanol-agua a 760 mmHg. Caso en estudio: A la presión de 760 mmHg se separa una mezcla binaria de metanol-agua ingresando a su punto de burbuja, con 36.0% mol de metanol con una alimentación de 100 Kmol/h. Se obtiene un producto de destilado de composición molar de 91.5% y de fondos de 2.87%. El destilado se va a condensar totalmente y el reflujo se va a regresar en su punto de burbuja. Se va a utilizar una relación de reflujo externa de 1.17. Se estudia la separación en dos etapas: Con condensador total y reboiler, Con condensador total e inyección de vapor de agua. En esta sección solo se presentan un número limitado de puntos que se muestran en la tabla n°01. Tabla n°01: Datos Entálpicos y de Equilibrio liquido-vapor calculados para el sistema Metanol-Agua a 760 mmHg Temperatura (K) 373.1521 360.8420 355.0072

Fracción mol de metanol en el liquido 0.0000 0.1000 0.2000

Fracción mol de metanol en el vapor 0.0000 0.4230 0.5795

hl (Mcal/kgmol) 1.8050 1.5967 1.5023

Hg (Mcal/kgmol) 11.2292 11.1612 10.9095

351.3263 348.5849 346.3233 344.3358 342.5203 340.8215 339.2081 337.6608

0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000

0.6692 0.7336 0.7864 0.8333 0.8770 0.9188 0.9597 1.0000

1.4454 1.4046 1.3718 1.3432 1.3172 1.2928 1.2694 1.2466

10.7105 10.5494 10.4110 10.2867 10.1717 10.0632 9.9596 9.8599

Fuente: Autor Diagrama Entalpico para el sistema binario Metanol-Agua a 760 mmHg 12

Entalpia especifica molar H,h ( Mcal/Kmol )

Hg vs y 10

8

6

4

2 hl vs x 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fraccion molar de metanol x,y

Grafica n°01: Diagrama Entalpia vs composición

Calculo del número de Etapas por el Método grafico de Ponchon-Savarit Con condensador total y reboiler: Metodo grafico Ponchon-Savarit para el sistema binario Metanol-Agua a 760 mmHg 22 20

(x ,M ) D

D

18

Entalpia especifica molar H,h ( Mcal/Kmol )

16 14 12

7

6

5

4

3

2

10

1

8 6 4 2 0 -2 -4 -6 -8 (x ,M )

-10 -12

R

0

R

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fraccion molar de metanol x,y

Grafica n°02: Método de destilación grafica Ponchon-Savarit con condensador total

De la grafica se ha obtenido 7 platos teóricos, de los cuales 1 corresponde al hervidor. La entrada de la alimentación se realizara en el plato n°4. Obteniendo: MD= 20.0983 Mcal/Kmol, MR= -9.7303 Mcal/Kmol Destilado= 37.3801 Kmol/h, Fondos= 62.6199 Kmol/h. Calculando los flujos de calor del condensador y del reboiler: QD= 703.9552 Mcal/h, QR = 717.1834 Mcal/h Con condensador total e inyección de vapor de agua. Metodo grafico Ponchon-Savarit para el sistema binario Metanol-Agua a 760 mmHg 22 20

(x ,M ) D

D

Entalpia especifica molar H,h ( Mcal/Kmol )

18 16 14 12

7

6

5

4

3

2

10

1

8 6 4 2 0 -2 -4 -6 -8

(x ,M ) C

-10

0

0.1

C

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fraccion molar de metanol x,y

Grafica n°03: Método grafico Ponchon-Savarit con condensador total e inyección de vapor

De la grafica se ha obtenido 7 platos teóricos, con la inyección de vapor, siendo la entrada de la alimentación en el plato n°4. Obteniendo: MD= 20.0983 Mcal/Kmol, MC= -8.7464 Mcal/Kmol Calor latente del vapor de agua: 12.00 Mcal/Kmol Calculando los flujos de calor del condensador: QD= 663.7287 Mcal/h y el flujo de vapor: 65.9640 Kmol/h, Destilado= 35.2441 Kmol/h, Fondos= 130.7199 Kmol/h CONCLUSIONES Se comprobó que la ejecución del método grafico riguroso de Ponchon-Savarit así como el cálculo de los coeficientes de actividad por el modelo UNIFAC puede realizarse en el lenguaje de programación MATLAB. BIBLIOGRAFIA Henley S. “Operaciones de separación por etapas de equilibrio en ingeniería química”. Ed. Reverte. 2000 Martínez P. y Rus E. “Operaciones de separación en ingeniería química-Métodos de cálculo”. Ed. Prentice Hall. 2004

Perry R.H., and Green D.W., “Perry's Chemical Engineer's Handbook”. McGrawHill, 8th edition, 2008. Poling B.E., Prausnitz J.M. and O’Connell J.P., “The properties of Gases and Liquids”. McGrawHill, 5th edition, 2001. Smith J.M., Van Ness H.C., and Abbott M.M., “Introduction to Chemical Engineering Thermodynamics”. McGrawHill, 7th edition, 2005. Yaws, Carl. “Chemical Properties Handbook”. McGraw-Hill. 1999.

Archivo Meth_Ponchon_Savarit.m disp(' Pära condensador total y una sola alimentacion') data=load('datos_UNIFAC_Ideal.txt'); T=data(:,1);x=data(:,2);y=data(:,3);h=data(:,4);H=data(:,5); xd=0.915; %Composicion de destilado 0.9 xb=0.0287; %Composicion de fondos 0.0008 zf=0.36; %Composicion de alimentacion 0.4507 q= 1; F=100;% Kmol/h D=F*(zf-xb)/(xd-xb);W=F-D; figure(1); hold on; grid off % une las graficas. plot(x,h,'-',y,H,'-'); axis on; box on; set(gca,'YMinortick','on','XMinortick','on') hd=interp1(x,h,xd); Hd=interp1(x,H,xd); yb=interp1(x,y,xb);Hb=interp1(y,H,yb); yf=interp1(x,y,zf); Hf=interp1(y,H,yf); hf=interp1(x,h,zf); hb=interp1(x,h,xb); MDm=(Hf-hf)/(yf-zf)*(xd-zf)+hf;Rm=(MDm-Hf)/(Hf-hf); R=1.17; MD=((R+1)*Hd-R*hd); MR=hf-((MD-hf)/(xd-zf))*(zf-xb); set(line([zf zf] ,[0 hf]),'Color',[.1 .1 0],'Linestyle',':') %para zf set(line([xd xd] ,[0 MD]),'Color',[.1 .1 0],'Linestyle',':') %para xd set(line([xb xb] ,[Hb MR]),'Color',[.1 .1 .1],'Linestyle',':') %para xb set(line([zf xd] ,[hf MD]),'Color',[1 0 0],'Linestyle','-') set(line([zf xb] ,[hf MR]),'Color',[1 0 0],'Linestyle','-') set(line([0 1] ,[0 0]),'Color',[0 0 0],'Linestyle','-') %Zona de enriquecimiento i=1; xp(1)=xd; yp(1)=xd; count1=0; tic while (xp(i)>zf) xp(i+1)=interp1(y,x,yp(i)); %interpolacion H1(i)=interp1(y,H,yp(i)); h1(i+1)=interp1(x,h,xp(i+1)); set(line([xp(i+1) yp(i)] ,[h1(i+1) H1(i)]),'Color',[0 .5 .5],'LineWidth',1.5) if (xp(i+1)>zf) set(line([xp(i+1) xd] ,[h1(i+1) MD]),'Color',[1 .1 1],'Linestyle','--','LineWidth',1) end %Iterar para encontrar y% ypp(i+1)=.5*yp(i);er(i+1)=1; hh=0; while abs(er(i+1))>1e-12 He(i+1)=interp1(y,H,ypp(i+1)); mm(i+1)=(xp(i+1)-xd)/(h1(i+1)-MD); yp1(i+1)=(He(i+1)+xp(i+1)/mm(i+1)-h1(i+1))*mm(i+1); er(i+1)=ypp(i+1)-yp1(i+1); ypp(i+1)=yp1(i+1); hh=hh+1; end ii=num2str(i); text(yp(i),H1(i),ii,'FontSize',10,'VerticalAlignment','bottom'); yp(i+1)=ypp(i+1); i=i+1; count1=count1+1; end % Zona de agotamiento count2=count1; while (xp(i)>xb) %Iterar para encontrar y% ypp(i+1)=.5*yp(i);er(i+1)=1; hh=0; while abs(er(i+1))>1e-12 He(i+1)=interp1(y,H,ypp(i+1)); nn(i)=(xp(i)-xb)/(h1(i)-MR); yp1(i+1)=(He(i+1)+xp(i)/nn(i)-h1(i))*nn(i);

er(i+1)=ypp(i+1)-yp1(i+1); ypp(i+1)=yp1(i+1); hh=hh+1; end yp(i+1)=ypp(i+1); xp(i+1)=interp1(y,x,yp(i+1)); %interpolacion H1(i+1)=interp1(y,H,yp(i+1)); h1(i+1)=interp1(x,h,xp(i+1)); set(line([xp(i+1) yp(i+1)] ,[h1(i+1) H1(i+1)]),'Color',[0 .5 .5],'LineWidth',1.5) set(line([yp(i+1) xb] ,[H1(i+1) MR]),'Color',[1 .1 1],'Linestyle','--','LineWidth',1) ii=num2str(i); text(yp(i+1),H1(i+1),ii,'FontSize',10,'VerticalAlignment','bottom'); i=i+1; count2=count2+1; end toc [M]=get(gca,'Ylim'); M(1)=M(1)-2;M(2)=M(2)-3; M=[M(1) M(2)];set(gca,'Ylim',M); [N]=get(gca,'Ytick'); N(1)=M(1); N(end)=M(2); N=[N(1):2:N(end)];set(gca,'Ytick',N); fprintf('Numero de platos en la columna: %3d\n',(count2)-1) title('Metodo grafico Ponchon-Savarit para el sistema binario Metanol-Agua a 760 mmHg','fontsize',12) xlabel('Fraccion molar de metanol x,y','fontsize',12); ylabel('Entalpia especifica molar H,h ( Mcal/Kmol )','fontsize',12) text(xd,MD,' (x_D,M_D)','FontSize',9,'VerticalAlignment','cap'); text(xb,MR,' (x_R,M_R)','FontSize',9,'VerticalAlignment','middle') Qd=D*(MD-hd);QR=W*(hb-MR); fprintf('MD= %7.4f Mcal/Kmol, MR= %7.4f Mcal/Kmol, QD= %7.4f Mcal/h, QR= %7.4f Mcal/h \n',MD,MR,Qd,QR) fprintf('F= %7.4f Kmol/h, D= %7.4f Kmol/h, W= %7.4f Kmol/h \n',F,D,W)

Archivo Meth_Ponchon_Savarit_2.m clear all; clc; clf disp(' Pära condensador total e inyeccion de vapor de agua') data=load('datos_UNIFAC_Ideal.txt'); T=data(:,1);x=data(:,2);y=data(:,3);h=data(:,4);H=data(:,5); xd=0.915; %Composicion de destilado 0.9 xb=0.0287; %Composicion de fondos 0.0008 zf=0.36; %Composicion de alimentacion 0.4507 q= 1; F=100;% Kmol/h Lambda=12; % Calor latente de vaporizacion del vapor de agua Mcal/Kmol figure(1); hold on; grid off % une las graficas. plot(x,h,'-',y,H,'-');axis on; box on set(gca,'YMinortick','on','XMinortick','on') hd=interp1(x,h,xd); Hd=interp1(x,H,xd); yb=interp1(x,y,xb);Hb=interp1(y,H,yb); yf=interp1(x,y,zf); Hf=interp1(y,H,yf); hf=interp1(x,h,zf);hb=interp1(x,h,xb); MDm=(Hf-hf)/(yf-zf)*(xd-zf)+hf; Rm=(MDm-Hf)/(Hf-hf); R=1.17; %Razon de reflujo R=L/D calcular % R=Rm*2; %Relacion de reflujo R=L/D calcular o ingresar directo MD=((R+1)*Hd-R*hd); MR=hf-((MD-hf)/(xd-zf))*(zf-xb); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% m1=(MD-hf)/(xd-zf); m2=(hb-Lambda)/xb; xc=(MD-Lambda-m1*xd)/(m2-m1); Mc=Lambda+m2*xc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(line([0 xc] ,[Lambda Mc]),'Color',[.1 .1 0],'Linestyle',':') set(line([zf zf] ,[0 hf]),'Color',[.1 .1 0],'Linestyle',':') %para zf set(line([xd xd] ,[0 MD]),'Color',[.1 .1 0],'Linestyle',':') %para xd set(line([zf xd] ,[hf MD]),'Color',[1 0 0],'Linestyle','-') set(line([zf xc] ,[hf Mc]),'Color',[1 0 0],'Linestyle','-') set(line([0 1] ,[0 0]),'Color',[0 0 0],'Linestyle','-') %Zona de enriquecimiento

i=1; xp(1)=xd; yp(1)=xd; count=0; tic while (xp(i)>zf) xp(i+1)=interp1(y,x,yp(i)); %interpolacion H1(i)=interp1(y,H,yp(i)); h1(i+1)=interp1(x,h,xp(i+1)); set(line([xp(i+1) yp(i)] ,[h1(i+1) H1(i)]),'Color',[0 .5 .5],'LineWidth',1.5) if (xp(i+1)>zf) set(line([xp(i+1) xd] ,[h1(i+1) MD]),'Color',[1 .1 1],'Linestyle','--','LineWidth',1) end %Iterar para encontrar y% ypp(i+1)=.5*yp(i);er(i+1)=1; hh=0; while abs(er(i+1))>1e-12 He(i+1)=interp1(y,H,ypp(i+1)); mm(i+1)=(xp(i+1)-xd)/(h1(i+1)-MD); yp1(i+1)=(He(i+1)+xp(i+1)/mm(i+1)-h1(i+1))*mm(i+1); er(i+1)=ypp(i+1)-yp1(i+1); ypp(i+1)=yp1(i+1); hh=hh+1; end ii=num2str(i); text(yp(i),H1(i),ii,'FontSize',10,'VerticalAlignment','bottom'); yp(i+1)=ypp(i+1); i=i+1; count=count+1; end % Zona de agotamiento while (xp(i)>xb) %Iterar para encontrar y% ypp(i+1)=.5*yp(i);er(i+1)=1; hh=0; while abs(er(i+1))>1e-12 He(i+1)=interp1(y,H,ypp(i+1)); nn(i)=(xp(i)-xc)/(h1(i)-Mc); yp1(i+1)=(He(i+1)+xp(i)/nn(i)-h1(i))*nn(i); er(i+1)=ypp(i+1)-yp1(i+1); ypp(i+1)=yp1(i+1); hh=hh+1; end yp(i+1)=ypp(i+1); xp(i+1)=interp1(y,x,yp(i+1)); %interpolacion H1(i+1)=interp1(y,H,yp(i+1)); h1(i+1)=interp1(x,h,xp(i+1)); set(line([xp(i+1) yp(i+1)] ,[h1(i+1) H1(i+1)]),'Color',[0 .5 .5],'LineWidth',1.5) set(line([yp(i+1) xc] ,[H1(i+1) Mc]),'Color',[1 .1 1],'Linestyle','--','LineWidth',1) ii=num2str(i); text(yp(i+1),H1(i+1),ii,'FontSize',10,'VerticalAlignment','bottom'); i=i+1; count=count+1; end toc [M]=get(gca,'Ylim'); M(1)=M(1);M(2)=M(2)-3; M=[M(1) M(2)];set(gca,'Ylim',M); [N]=get(gca,'Ytick'); N(1)=M(1); N(end)=M(2); N=[N(1):2:N(end)];set(gca,'Ytick',N); fprintf('Numero de platos en la columna: %3d\n',count) title('Metodo grafico Ponchon-Savarit para el sistema binario Metanol-Agua a 760 mmHg','fontsize',12) xlabel('Fraccion molar de metanol x,y','fontsize',12); ylabel('Entalpia especifica molar H,h ( Mcal/Kmol )','fontsize',12) text(xd,MD,' (x_D,M_D)','FontSize',9,'VerticalAlignment','cap'); text(xc,Mc,' (x_C,M_C) ','FontSize',9,'VerticalAlignment','middle') D=F*((hf-Lambda)*xb-(hb-Lambda)*zf)/((MD-Lambda)*xb-(hb-Lambda)*xd); Qd=D*(MD-hd); W=(F*zf-D*xd)/xb; Vv=D+W-F; fprintf('MD= %7.4f Mcal/Kmol, MC= %7.4f Mcal/Kmol, QD= %7.4f Mcal/h \n',MD,Mc,Qd) fprintf('F= %7.4f Kmol/h, D= %7.4f Kmol/h, W= %7.4f Kmol/h, V= %7.4f Kmol/h \n',F,D,W,Vv)

El Archivo gamma_UNIFAC.m se utilizó para el cálculo de Temperatura de burbuja. El algoritmo de ejecución se encuentra en el texto de Van Ness. Archivo gamma_UNIFAC.m vka=Vkdata('metanol'); vkb=Vkdata('agua'); vka=vka'; vkb=vkb'; data=load('Grupos_Subgrupos.txt'); Grupo=1; Subgrupo=2; R_k=3; Q_k=4; G_data=data(:,Grupo);G_m=G_data(:); Sg_data=data(:,Subgrupo);S_k=Sg_data(:); R_data=data(:,R_k);R=R_data(:); Q_data=data(:,Q_k);Q=Q_data(:); a=load('a_m_k.txt'); vki=[vka vkb]; vt=vka'+vkb'; K=find(vt~=0); nn1=length(vki(1,:)); T= ; %Temperatura en K xx1=0:0.0001:1; % Componente mas volátil. xx1=xx1'; xx2=1-xx1; xx=[xx1 xx2]; nn2=length(xx1); T=T*ones(nn2,1); sumae1=zeros(nn2,1);sumae2=zeros(nn2,1); sumad1=zeros(nn2,1); sumaf=zeros(nn2,nn1); z=10; tic for i=1:nn1 % Calculo de ri y qi r(i)=0; q(i)=r(i); for k=K r(i)=vki(k,i)*R(k)+r(i); q(i)=vki(k,i)*Q(k)+q(i); end end for j=1:nn1 % Calculo de eki for k=K e(k,j)=vki(k,j)*Q(k)/q(j); end end e(K,:); for KK=1:nn2; for k=K % Calculo de thetak for h=1:nn1 p(k,h,KK)=xx(KK,h)*q(h)*e(k,h); d(KK,h)=xx(KK,h)*q(h); end theta1(KK,k)=sum(p(k,:,KK))/sum(d(KK,:)); end theta1(KK,K); a(G_m(K),G_m(K)); for i=K % Calculo de Sk S(KK,i)=0; for j=K S(KK,i)=theta1(KK,j)*exp(-a(G_m(j),G_m(i))/T(KK))+S(KK,i); end end S(KK,K); for i=1:nn1 % Calculo de beta beta(i,K,KK)=0; for k=K for k2=K beta(i,k,KK)=e(k2,i)*exp(-a(G_m(k2),G_m(k))/T(KK))+beta(i,k,KK); end end end beta(:,K,KK); for w=1:nn1

sumae1(KK)=r(w)*xx(KK,w)+sumae1(KK); sumae2(KK)=q(w)*xx(KK,w)+sumae2(KK); end for j=1:nn1 JJ(KK,j)=r(j)/sumae1(KK); LL(KK,j)=q(j)/sumae2(KK); end for kk=1:nn1 % Calculo de lngammaC lngammaC(KK,kk)=1-JJ(KK,kk)+log(JJ(KK,kk))-(z/2)*q(kk)*(1JJ(KK,kk)/LL(KK,kk)+log(JJ(KK,kk)/LL(KK,kk))); end lngammaC(KK,:); for ff=1:nn1 % Calculo de lngammaR for k=K sumaf(KK,ff)=theta1(KK,k)*beta(ff,k,KK)/S(KK,k)- e(k,ff)*log(beta(ff,k,KK)/S(KK,k))+sumaf(KK,ff); end lngammaR(KK,ff)=q(ff)*(1-sumaf(KK,ff)); end lngammaR(KK,:); for jj=1:nn1 gamma(KK,jj)=exp(lngammaC(KK,jj)+lngammaR(KK,jj)); end fprintf('%3d %7.16f %7.14f\n',KK,gamma(KK,1),gamma(KK,2)); gamma(KK,:); end toc