Solve complex chamical equilibrium using Newton-Raphson

I have 9 equation, such as this figure below,. I have to solve this equation to get value n1, n2, n3, n4, n5, n6. Can you help me? Danke. I have MATLAB R2013a

Answer 1

[Torsten]

help fsolve

Best wishes

Torsten.

[rifa'atul mahmudhah]

thankyou Torsten.

I just make the program like this: I have two m.file,... cobafun.m for input the parameter, and gas_func.m for write the solver.

.............there is my cobafun.m...............

function [fx] = cobafun(x)

nH2 = x(1);

nC = x(2);

nCO2 = x(3);

nH2O = x(4);

nCO = x(5);

nCH4 = x(6);

lamda1 = x(7);

lamda2 = x(8);

lamda3 = x(9);

R = 0.008314; %kJ/(mol.K)

%that will be variated

a = 0.5; %molar ratio

T = 773; %Temperature

%total of mol

ntotal = nH2 + nC + nCO2 + nH2O + nCO + nCH4;

%balance atom equation

fx (7,1) = 2*nH2 + 4*nCH4 + 2*nH2O - 1.67 + 2*a; %hydrogen (H) balance

fx (8,1) = nCO2 + nCO + nC + nCH4 - 1; %carbon (C) balance

fx (9,1) = 2*nH2 + nCO + nH2O - 0.75 - a; %oksygen (O) balance

%Gibbs equation

fx (1,1) = log(nH2/ntotal)*R*T + (2*lamda1);

fx (2,1) = log(nC/ntotal)*R*T + (lamda2);

fx (3,1) = -395.4928817 + log(nCO2/ntotal)*R*T + (lamda2 + 2*lamda3);

fx (4,1) = -205.0850656 + log(nH2O/ntotal)*R*T + (2*lamda1 + lamda3);

fx (5,1) = -180.3948228 + log(nCO/ntotal)*R*T + (lamda2 + lamda3);

fx (6,1) = -4.922967331 + log(nCH4/ntotal)*R*T + (lamda2 + 4*lamda1);

......................there is my gas_func.m.............

function gasf_func

clear, clc;

x0 = [0.5 0.5 0.5 0.1 0.3 0.1 1 1 1]; %initial guess (still confuse how to take the value)

disp('variable values at the initial estimate');

f0=cobafun(x0);

disp(' variabel value function value');

for i=1:size(x0,2);

    disp([' x' int2str(i) '                 ' num2str(x0(i)) '             
' num2str(f0(i))])

end

options=optimset('Display','final','MaxIter',10000,'MaxFunEvals',100000 );

xsolv=fsolve(@cobafun,x0,options);

disp('variable at the solution');

fsolv=cobafun(real(xsolv));

disp(' variable value function value');

for i=1:size(x0,2);

    disp([' x' int2str(i) '                 ' num2str(real(xsolv(i))) '             
' num2str(fsolv(i))])

end

and I get the output like this:...

variable values at the initial estimate

 variabel              value              function value

 x1                 0.5             -6.9093

 x2                 0.5             -7.9093

 x3                 0.5             -401.4022

 x4                 0.1             -221.3378

 x5                 0.3             -190.5871

 x6                 0.1             -19.1757

 x7                 1             0.93

 x8                 1             0.4

 x9                 1             0.15

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.

criteria details

variable at the solution

 variable               value              function value

 x1                 -0.50911             141.9083+20.19014i

 x2                 -0.51078             140.6942+20.19014i

 x3                 0.93946             -245.9691

 x4                 0.64581             -59.1916

 x5                 0.45233             -38.0249

 x6                 -1.0177             145.0134+20.19014i

 x7                 1.2029             -4.4675

 x8                 1.1708             -1.1367

 x9                 2.4567             -1.1701

>>

I get wrong result, I think. because I have negative result,.. ald also, maybe I do a mistake to take the initial point.

can you help me, what must I do?

[Torsten]

You won't get negative values for your variables if you solve in variables squared, e.g. set nH2 = x(1)^2 instead of nH2 = x(1) (same for the other molar masses).

You won't get complex values for your variables if you take exp of your equations. Thus instead of solving ln(a)=b, solve a=exp(b).

Best wishes

Torsten.