Problem using fsolve() for a Chemical equilibrium

Hey all,

I'm trying to solve a chemical equilibrium problem using fsolve(), but I'm having problems making it work. Hope someone here can help me.

I'm interested in obtaining the molar fraction (y) of the spicies (CH4,H2O...) in the equilibrium. At the begining I only have 2 mols of CH4 and 3 mols of H2O, that then react and leave a misture of the five spicies.

For the equations I have:

function Y = fun1(x)

%CH4, H2O --> CH4, H2O, CO, CO2, H2

%"i" [CH4, H2O, CO, CO2, H2]
%"k" [C, H, O]

A = [2,14,3]; %Initial mols of "k"
a = [1, 0, 4; 0, 1, 2; 1, 1, 0; 1, 2, 0; 0, 0, 2]; %Number of "k" atoms in spicie "i"
DelG = [19720, -192420, -200240, -395790, 0]; %Delta G0 @ 1000K (J/mol)
R = 8.314; %J/mol.K
T = 1000; %K

n = x(1:length(a)); % number of mol of specie "i"
Lamb = x(length(a)+1:length(x)); %Lagrange operators

for k=1:length(A) %Atomic balance
    for i=1:length(a)
        t(i) = n(i)*a(i,k);
        AB(k) = sum(t)-A(k);
    end
end

for i=1:length(a) %Minimization of "Del G"
    for k=1:length(A)
        y(i) = n(i)/sum(n);
        Eq(i) = (DelG(i)/(R*T)) + log(y(i)) + (sum(Lamb(k)*a(i,k))/(R*T));
    end
end

Y = [Eq, AB];

end

For the fsolve() part I have:

n0 = [2, 3, 0, 0, 0]; %initial guess for the number of mols for each spicie
Lamb0 = [0, 0, 0]; %initial guess for the Lagrange operators
x0 = [n0, Lamb0];
[x,fval] = fsolve('fun1', x0, []);
n = x(1:5);

for i=1:length(n)
    y(i) = n(i)/sum(n);
end

As a result I keep getting the following message:

Error using trustnleqn (line 28)
Objective function is returning undefined values at initial point. FSOLVE cannot continue.

Error in fsolve (line 368)
    [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...

Error in Questao3 (line 9)
[x,fval] = fsolve('funcaoQ3_2', x0, []);

I've already checked the equations and even tryed to write them 1by1, but the same thing occurs.

Hope someone can help, and thanks in advance.

Answer 1

[Walter Roberson]

        y(i) = n(i)/sum(n);
        Eq(i) = (DelG(i)/(R*T)) + log(y(i)) + (sum(Lamb(k)*a(i,k))/(R*T));

Your first two n values are non-zero but the other 3 are 0. so n(i) will be 0 for i=3 so y(3)=0 and log(y(3)) will be -inf

[João Lazari]

Thaks Walter,

So many hours studing to miss the log(0).

I've just tried other sets of initial values and now it 'works'.

No solution found.

fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the default value of the function tolerance.

<stopping criteria details>


y =

  -1.0802 + 0.5981i   2.8978 - 0.0633i   0.1940 + 1.2523i   2.4587 - 0.3812i  -4.4703 - 1.4059i

[Walter Roberson]

You have 8 equations in 6 variables.

[ (2*X8)/4157 + log(X1/(X1 + X2 + X3 + X4 + X5)) + 9860/4157,
  X8/4157 + log(X2/(X1 + X2 + X3 + X4 + X5)) - 96210/4157,
  log(X3/(X1 + X2 + X3 + X4 + X5)) - 6779233742668001/281474976710656,
  log(X4/(X1 + X2 + X3 + X4 + X5)) - 6699842496530583/140737488355328,
  X8/4157 + log(X5/(X1 + X2 + X3 + X4 + X5)),
  X1 + X3 + X4 - 2,
  X2 + X3 + 2*X4 - 14, 4*X1 + 2*X2 + 2*X5 - 3]

This has no X6 and no X7.

In the tests I am doing at the moment to find the best (real-valued) fit I can find to all of the equations, I am finding that X8 is becoming large, and X1 and X5 are becoming quite small. log() of quite small can be fairly negative, and that is balancing X8 becoming large in the first and fifth equations. For example

[X1, X2, X3, X4, X5, X8] =
[5.15508711453230437e-21, 0.660422184314845673, 1.75619892254994348, 4.86725538886841491, 2.27046406365568618e-09, 96940.8426038390608]

[João Lazari]

Hi Walter,

I've tried with the code were all the equations are writen 1by1 as follows, and obtained the result I expected

function Y = fun2(x)

%"i" [CH4, H2O, CO, CO2, H2]
%"k" [C, H, O]

Ak = [2,14,3]; 
a = [1, 0, 4; 0, 1, 2; 1, 1, 0; 1, 2, 0; 0, 0, 2];
DelG = [19720, -192420, -200240, -395790, 0]; %Del G0 @ 1000k (J/mol)
R = 8.314; %J/mol.K
T = 1000; %K

y(1)=x(1)/sum(x(1:5));
y(2)=x(2)/sum(x(1:5));
y(3)=x(3)/sum(x(1:5));
y(4)=x(4)/sum(x(1:5));
y(5)=x(5)/sum(x(1:5));


Y(1) = (DelG(1)/(R*T)) + log(y(1)) + (x(6) + 4*x(8))/(R*T);
Y(2) = (DelG(2)/(R*T)) + log(y(2)) + (x(7) + 2*x(8))/(R*T);
Y(3) = (DelG(3)/(R*T)) + log(y(3)) + (x(6) + x(7))/(R*T);
Y(4) = (DelG(4)/(R*T)) + log(y(4)) + (x(6) + 2*x(7))/(R*T);
Y(5) = (DelG(5)/(R*T)) + log(y(5)) + (2*x(8))/(R*T);
Y(6) = x(1) + x(3) + x(4) - Ak(1);
Y(7) = 4*x(1) +  2*x(2) + 2*x(5) - Ak(2);
Y(8) = x(2) + x(3) + 2*x(4) - Ak(3);

end

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>


y =

    0.0196    0.0980    0.1743    0.0371    0.6711

Now I'm trying to find what I missed in the previous code (I already finded some, but it isn't quite right)

    function Y = fun1(x)

%Ordem dos "i" [CH4, H2O, CO, CO2, H2]
%Ordem dos "k" [C, H, O]

A = [2,14,3]; 
a = [1, 0, 4; 0, 1, 2; 1, 1, 0; 1, 2, 0; 0, 0, 2];
DelG = [19720, -192420, -200240, -395790, 0]; %Delta G0 @ 1000K (J/mol)
R = 8.314; %J/mol.K
T = 1000; %K

n = x(1:length(a)); 
Lamb = x(length(a)+1:length(x)); 

for i=1:length(a)    
    for k=1:length(A)
        y(i) = n(i)/sum(n);
        z(i,k) = Lamb(k)*a(i,k);
        Eq(i) = (DelG(i)/(R*T)) + log(y(i)) + (sum(z(i,:))/(R*T));
    end
end

for k=1:length(A)
    for i=1:length(a)
        t(i,k) = (n(i)*a(i,k));
        BA(k) = sum(t(:,k)) - A(k);
    end
end


Y = [Eq, BA];

end

[Walter Roberson]

I recommend against using length(a) . a is a 2D array, and length(a) is defined as:

if isempty(a)
    length is 0
else
    length is max(size(a))
end

a is 5 x 3 so length() for it should be 5, but it is not clear that is that is the value you are expecting.