Solve complex chamical equilibrium using Newton-Raphson

2015 6 19
1 个回答
更新时间：2015 7 16
9 次查看（30 天）

1 个投票

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

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2015-6-19

0 个投票

help fsolve
Best wishes
Torsten.

2 个评论
显示 无 隐藏 无

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
Torsten 2015-7-16
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Parallel Computing Toolbox 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by