Newton Raphson multivariate xy equilibrium

Hello, I am developing a code to find the equilibrium between xy and the behavior of the temperature and each of the components. The 2 variables I have to find are y component and T the temperature. I think the code is good but obviously have some issues since it isn’t returning the right answer. If you could help me I will be really grateful.

clc
clear
close all
%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);

%Inicializacion;
for i = 1: n
    iter=0;
    itermax=1000;
    error=1000;
    tol=1e-6;
    % Inicialización  
    m0= [0.01 ; 400]';
    % Recorrido del método NR multivariado
    while error>tol && iter<itermax
        f0= raoult(m0,x(i)) ;
        j0=jacobiano(m0,x(i)) ;
        delta=-(j0\f0);
        xnew=m0+delta;
        % Criterio del ciclo
        error=norm (delta) ;
        % Reemplazo nuevo valor
        m0=xnew;
    end
    yresp(i)= xnew(1); 
    Tresp(i) = xnew(2);
end

figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'

figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno' 

function j = jacobiano(m,x1)

n = length(x1);
h = 1e-9;
for i=1:n
    dh = zeros(n,1);
    dh(i) = h;
    derv = (raoult(m+dh,x1))-raoult(m,x1)/(h);
    j(:,i) = derv;
end
end

CONSTANTES 
function resp = Antoine(T) 

%Benceno 
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno 
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;

resp (:,1) = 10^(A(1)-B(1)/(T+C(1)));

resp (:,2) = 10^(A(2)-B(2)/(T+C(2)));

end 

function gamma= wilson (T,x)
x1 = x; 
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884; 
R = 8.314; 
z = 2;   

beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);

delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;

lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T)); 
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));

gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
end 

function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
gm = wilson(x1,T);
resp(:,1) = x1.*gm(1).*Psat(1)- y1*P;
resp(:,2) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end

Answer 1

[Torsten]

function main
%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);
m0= [0.01 ; 400]';
for i=1:n
    fun = @(m0)raoult(m0,x(i)) ;
    sol = fsolve(fun,m0);
    yresp(i)= sol(1); 
    Tresp(i)= sol(2);  
    m0 = sol;
end
figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'
figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno' 
end
function resp = Antoine(T) 
%Benceno 
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno 
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;
resp (:,1) = 10^(A(1)-B(1)/(T+C(1)));
resp (:,2) = 10^(A(2)-B(2)/(T+C(2)));
end 
function gamma= wilson (T,x)
x1 = x; 
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884; 
R = 8.314; 
z = 2;   
beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);
delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;
lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T)); 
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));
gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
end 
function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
%gm = wilson(x1,T)
gm(1:2) = 1.0;
resp(:,1) = x1.*gm(1).*Psat(1)- y1*P;
resp(:,2) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end

I think you should check the calculation of the activity coefficients.

[Juliana Quintana]

Hello, thank you so much for answering so fast. I forgot to mention that I can’t use fsolve or similar functions. Also thanks for the advise I will check the calculation of gamma. If you could also help me with the think I just mention it will be amazing. Thank you again. Have a nice day.

[Torsten]

I think you should use "fsolve" until the physical problems are solved.

Then you can switch to your own Newton code.

[Juliana Quintana]

I just fixed the physical problems and I think it’s giving the proper answer. Now I’ll change to the newton code. Also do you have any recommendation por the jacobian function . Thank you

[Torsten]

%Vector de x
n=100;
x=linspace(0,1,n);
%Vector para guardas las variables
yresp = zeros(n,1);
Tresp = zeros(n,1);
%Inicializacion;
for i = 1: n
    iter=0;
    itermax=1000;
    error=1000;
    tol=1e-6;
    % Inicialización  
    m0= [0.01 ; 400];
    % Recorrido del método NR multivariado
    while error>tol && iter<itermax
        f0= raoult(m0,x(i)) ;
        j0=jacobiano(m0,x(i)) ;
        delta=-(j0\f0);
        xnew=m0+delta;
        % Criterio del ciclo
        error=norm (delta) ;
        % Reemplazo nuevo valor
        m0=xnew;
    end
    yresp(i) = xnew(1); 
    Tresp(i) = xnew(2);
end
figure (1)
plot(x,Tresp, 'r',yresp,Tresp,'b')
legend ('x' ,'y')
xlabel 'x,y, Benceno'
ylabel 'T (K)'
figure (2)
plot(x,yresp)
xlabel 'x'
ylabel 'y'
legend 'Benceno' 
end
function j = jacobiano(m,x1)
y = raoult(m,x1);
h = 1e-8;
for i=1:numel(m)
    m(i) = m(i) + h;
    yh = raoult(m,x1);
    j(:,i) = (yh-y)/h;
    m(i) = m(i) - h;
end  
end
CONSTANTES 
function resp = Antoine(T) 
%Benceno 
A(1)= 6.01905; B(1) = 1204.637; C(1)=53.081;
%Tolueno 
A(2)=6.08426; B(2)= 1347.2; C(2)=53.363;
resp(1,:) = 10^(A(1)-B(1)/(T+C(1)));
resp(2,:) = 10^(A(2)-B(2)/(T+C(2)));
end 
function gamma= wilson (T,x)
x1 = x; 
x2 = 1-x1;
%Benceno(1)
%Tolueno (2)
v251 = 90.4; v252 = 104.9 ;
vb1=96; vb2=118.2;
delta251=18.8; delta252=18.7;
tb1=80.09; tb2=110.622;
eps12=0.0851; eps21 = -0.0884; 
R = 8.314; 
z = 2;   
beta1 =(vb1-v251)/(tb1-25);
beta2= (vb2-v252)/(tb2-25);
v1= v251+beta1*((T-273.15)-25);
v2= v252+beta2*((T-273.15)-25);
delta1= (v251/v1)*delta251;
delta2= (v252/v2)*delta252;
lambda11 = -(2/z)*v1*(delta1)^2 ;
lambda12 = -(1-eps12)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda21 = -(1-eps21)*(2/z)*((v1*v2)^0.5)*(delta2*delta1);
lambda22 = -(2/z)*v2*(delta2)^2 ;
A21= (v1/v2)*exp(-(lambda21-lambda22)/(R*T)); 
A12=(v2/v1)*exp(-(lambda12-lambda11)/(R*T));
gamma(1) = exp(-log(x1+A12*x2)+(x2*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
gamma(2) = exp(-log(x2+A21*x1)-(x1*((A12/(x1+(A12*x2)))-(A21/((A21*x1+x2)))))); 
end 
function resp = raoult(m,x1)
P = 101.325 ;
y1 = m(1);
T = m(2);
Psat = Antoine(T);
%gm = wilson(x1,T);
gm(1:2) = 1.0;
resp(1,:) = x1.*gm(1).*Psat(1)- y1*P;
resp(2,:) = (1-x1).*gm(2).*Psat(2)- (1-y1)*P;
end

I used your code for the version with activity coefficients equal to 1.

You will have to replace "wilson" with your corrected version.

[Juliana Quintana]

I just made al the changes and it works perfecty. Thank you for your help and time. Amazing work :)