chemical equilibrium with Newton-Raphson Iterative Solver
xingyu
I am trying to solve the chemical equilibrium problem by using Newton-Raphson Iterative Solver. I download a custom solver and apply it to solve chemical equilibrium problem. I somehow can not get the answer,Please help what I'm doing wrong. the m-files are shown below:
function x=equilibrium(n,m,l,q,T,P,F,EGR)
% x13[CnHmOlNq+(n+m/4-l/2)/F*(O2+3.76N2+4.76EGR/(1-EGR)*(0.135CO2+0.124H2O+0.74N2))]
%=x1H+x2O+x3N+x4H2+x5OH+x6CO+x7NO+x8O2+x9H2O+x10CO2+x11N2
Tr=0.001*T;
k{1}=0.432168*log(Tr)-11.2464/Tr+2.67269-0.0745744*Tr+0.00242484*Tr*Tr; % 0.5H2=H
k{2}=0.310805*log(Tr)-12.9540/Tr+3.21779-0.0738336*Tr+0.00344645*Tr*Tr; % 0.5O2=O
k{3}=0.389716*log(Tr)-24.5828/Tr+3.14505-0.0963730*Tr+0.00585643*Tr*Tr; % 0.5N2=N
k{4}=-0.141784*log(Tr)-2.13308/Tr+0.853461+0.0355015*Tr-0.00310227*Tr*Tr;% 0.5O2+0.5H2=OH
k{5}=0.0150879*log(Tr)-4.70959/Tr+0.646096+0.00272805*Tr-0.00154444*Tr*Tr;% 0.5N2 + 0.5O2=NO
k{6}=-0.752364*log(Tr)+12.4210/Tr-2.60286+0.259556*Tr-0.0162687*Tr*Tr;% H2 + 0.5O2 =H2O
k{7}=-0.00415302*log(Tr)+14.8627/Tr-4.75746+0.124699*Tr-0.00900227*Tr*Tr;% CO + 0.5O2=CO2
for i=1:1:7,
k{i}=power(10,k{i});
end
r0=(n+m/4-l/2)/F;
R=r0*4.76*EGR/(1-EGR);
c1=k{1};c2=k{2};c3=k{3};c4=k{4};c5=k{5};c6=k{6};c7=k{7};c8=(P/101)^0.5;c9=r0;c10=R;c11=n;c12=m;c13=l;c14=q;
F = {'c1/c8*v1^0.5+2*v1+2*c6*v1*v3^0.5*c8+c4*v1^0.5*v3^0.5-(c12+c10*0.248)/(c11+c10*0.135)*(v2+c7*v2*v3^0.5*c8)';
'c2/c8*v3^0.5+c4*v1^0.5*v3^0.5+v2+c5*v3^0.5*v4^0.5+2*v3+c6*v1*v3^0.5*c8+2*c7*v2*v3^0.5*c8-(c13+c9*2+c10*0.394)/(c11+c10*0.135)*(v2+c7*v2*v3^0.5*c8)';
'c3*v4^0.5/c8+c5*v3^0.5*v4^0.5+2*v4-(c14+c9*2*3.76+c10*0.74*2)/(c11+c10*0.135)*(v2+c7*v2*v3^0.5*c8)'
'c1*v1^0.5/c8+c2*v3^0.5/c8+c3*v4^0.5/c8+v1+c4*v1^0.5*v3^0.5+v2+c5*v3^0.5*v4^0.5+v3+c6*v1*v3^0.5*c8+c7*v2*v3^0.5*c8+v4-1'};
dFdx={'0.5*c1/c8/v1^0.5+2+2*c6*v3^0.5*c8+0.5*c4*v3^0.5/v1^0.5','-(c12+c10*0.248)/(c11+c10*0.135)*(1+c7*v3^0.5*c8)','c6*v1/v3^0.5*c8+0.5*c4*v1^0.5/v3^0.5-(c12+c10*0.248)/(c11+c10*0.135)*0.5*c7*v2/v3^0.5*c8','0';
'0.5*c4/v1^0.5*v3^0.5+c6*v3^0.5*c8','1+2*c7*v3^0.5*c8-(c13+c9*2+c10*0.394)/(c11+c10*0.135)*(1+c7*v3^0.5*c8)','0.5*c2/c8/v3^0.5+0.5*c4*v1^0.5/v3^0.5+0.5*c5*v4^0.5/v3^0.5+2+0.5*c6*v1*c8/v3^0.5+c7*v2*c8/v3^0.5-(c13+c9*2+c10*0.394)/(c11+c10*0.135)*0.5*c7*v2/v3^0.5*c8','0.5*c5*v3^0.5/v4^0.5';
'0','-(c14+c9*2*3.76+c10*0.74*2)/(c11+c10*0.135)*(1+c7*v3^0.5*c8)','0.5*c5*v4^0.5/v3^0.5-(c14+c9*2*3.76+c10*0.74*2)/(c11+c10*0.135)*0.5*c7*v2/v3^0.5*c8)','0.5*c3/v4^0.5/c8+0.5*c5*v3^0.5/v4^0.5+2';
'0.5*c1/v1^0.5/c8+1+0.5*c4*v3^0.5/v1^0.5+c6*v3^0.5*c8','1+c7*v3^0.5*c8','0.5*c2/v3^0.5/c8+0.5*c4*v1^0.5/v3^0.5+0.5*c5*v4^0.5/v3^0.5+1+0.5*c6*v1/v3^0.5*c8+0.5*c7*v2/v3^0.5*c8','0.5*c3/v4^0.5/c8+1+0.5*c5*v3^0.5/v4^0.5+1'};
xi=[0.1;0.1;0.1;0.7];
tol=eps;
max_iter=1000;
x = nrsolve(F,dFdx,xi,tol,max_iter,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14);
function x = nrsolve(F,dFdx,xi,tol,max_iter,varargin)
%==========================================================================
% x = nrsolve( F, dFdx, xi, tol, max_iter, c1, c2.... )
%Newton-Raphson Iterative Solver nrsolve
%
%Author: John Fuller
% 703-404-7620
%
%An iterative solver for N equations with N unknown variables, up to N=8.
%Supports built-in Matlab functions defined within F and dFdx.
%--------------------------------------------------------------------------
%Inputs:
% F - Nx1 CELL array listing functions to be evaluated,
% to be formatted as follows:
% SYSTEM NRSOLVE
% Ex: F1 = x + y + 2*z = v1 + v2 + 2*v3 = 0
% F2 = (x - y)/2 = (v1 - v2)/2 = 0
% F3 = x + z = v1 + v3 = 0
%
% Translates to F = {'v1+v2+2*v3';'(v1-v2)/2';'v1+v3'}
% Variables must be named as v1, v2, v3... vN
% Constants must be named as c1, c2, c3... cM
%
% dFdx - NxN CELL array listing derivatives of functions to be evaluated,
% the same way F is defined. The {m,n} location of the array
% would define dF_m/dv_n. IE: The cell in row 2, column 1 would
% be derivative of F2 w.r.t. var 1. This is the string
% equivalent of the sensitivity matrix used in N-R solver. From
% the example above, row 2 column 1 would simply be '0.5' or '1/2'
%
% xi - Nx1 vector of initial guesses for unknown variables v1, v2, v3...
%
% tol - User-defined precision value at which to discontinue iteration.
% Defined as the minimum norm of difference vector dx at which to
% discontinue iteration and obtain solution.
%
% max_iter - maximum allowable iterations before termination.
%
% c1...cM - additional constant values that user would like to pass into
% nrsolve. See testnrsolve.m Example 5 for more information.
% These constants may be defined within F and dFdx in the same
% manner as v1, v2, etc.
%--------------------------------------------------------------------------
%Output:
% x - Nx1 vector of solved answers.
% x will be NaN if solver cannot resolve the equations in 1000
% iterations, and a warning will be generated.
%
%As with all iterative solvers, care needs to be taken when solving
%non-linear sets of equations. Particular attention must be dedicated to
%selecting a proper initial guess to avoid iteration failure.
%==========================================================================
%Calculations begin here
%Check inputs for errors.
if size(xi,2)~=1,
error('xi vector must be an Nx1 vector');
end
N=size(xi,1);
if N<1
error('xi vector must be an Nx1 vector where N>=1')
end
if tol<0
error('Precision tolerance may not be negative')
end
%Stop run if vector is larger than 8.
if N>8
error('Holy crap! More than 8 equations to solve??')
end
%Make sure cell arrays are properly formatted
if ~iscellstr(F) || size(F,1)~=N || size(F,2)~=1
error('Input F must be an Nx1 cell array')
end
if ~iscellstr(dFdx) || size(dFdx,1)~=N || size(dFdx,2)~=N
error('Input dFdx must be an NxN cell array')
end
%Assign variables in varargin to c1, c2, ... cM
M=nargin-5;
if M>0
for ii=1:M
if ~isnumeric(varargin{ii})
error('Invalid input for additional constant values.')
end
eval(['c',num2str(ii,'%i'),'=',num2str(varargin{ii}),';']);
end
end
%Append prefixes and suffixes to input function arrays
F_eval_prefix={'F1=';'F2=';'F3=';'F4=';'F5=';'F6=';'F7=';'F8='};
F_eval_suffix={';';';';';';';';';';';';';';';'};
F_eval_str_array=strcat(F_eval_prefix(1:N),F,F_eval_suffix(1:N));
dF_eval_prefix={'dF1dv1=','dF1dv2=','dF1dv3=','dF1dv4=','dF1dv5=','dF1dv6=','dF1dv7=','dF1dv8=';...
'dF2dv1=','dF2dv2=','dF2dv3=','dF2dv4=','dF2dv5=','dF2dv6=','dF2dv7=','dF2dv8=';...
'dF3dv1=','dF3dv2=','dF3dv3=','dF3dv4=','dF3dv5=','dF3dv6=','dF3dv7=','dF3dv8=';...
'dF4dv1=','dF4dv2=','dF4dv3=','dF4dv4=','dF4dv5=','dF4dv6=','dF4dv7=','dF4dv8=';...
'dF5dv1=','dF5dv2=','dF5dv3=','dF5dv4=','dF5dv5=','dF5dv6=','dF5dv7=','dF5dv8=';...
'dF6dv1=','dF6dv2=','dF6dv3=','dF6dv4=','dF6dv5=','dF6dv6=','dF6dv7=','dF6dv8=';...
'dF7dv1=','dF7dv2=','dF7dv3=','dF7dv4=','dF7dv5=','dF7dv6=','dF7dv7=','dF7dv8=';...
'dF8dv1=','dF8dv2=','dF8dv3=','dF8dv4=','dF8dv5=','dF8dv6=','dF8dv7=','dF8dv8='};
dF_eval_suffix={';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';';...
';',';',';',';',';',';',';',';'};
dF_eval_str_array=strcat(dF_eval_prefix(1:N,1:N),dFdx,dF_eval_suffix(1:N,1:N));
%Build strings for evaluation of functions and derivatives
F_num_eval_str='F_num=[F1;F2;F3;F4;F5;F6;F7;F8];';
F_num_eval_str=[F_num_eval_str(1:6+N*3),F_num_eval_str(end-1:end)];
dF_num_eval_str_full=['dF1dv1,dF1dv2,dF1dv3,dF1dv4,dF1dv5,dF1dv6,dF1dv7,dF1dv8,',...
'dF2dv1,dF2dv2,dF2dv3,dF2dv4,dF2dv5,dF2dv6,dF2dv7,dF2dv8,',...
'dF3dv1,dF3dv2,dF3dv3,dF3dv4,dF3dv5,dF3dv6,dF3dv7,dF3dv8,',...
'dF4dv1,dF4dv2,dF4dv3,dF4dv4,dF4dv5,dF4dv6,dF4dv7,dF4dv8,',...
'dF5dv1,dF5dv2,dF5dv3,dF5dv4,dF5dv5,dF5dv6,dF5dv7,dF5dv8,',...
'dF6dv1,dF6dv2,dF6dv3,dF6dv4,dF6dv5,dF6dv6,dF6dv7,dF6dv8,',...
'dF7dv1,dF7dv2,dF7dv3,dF7dv4,dF7dv5,dF7dv6,dF7dv7,dF7dv8,',...
'dF8dv1,dF8dv2,dF8dv3,dF8dv4,dF8dv5,dF8dv6,dF8dv7,dF8dv8,'];
F_eval_str=[];
dF_eval_str=[];
dF_num_eval_str='dFdx_num=[';
for ii=1:N
F_eval_str=[F_eval_str,F_eval_str_array{ii,1}]; %#ok<AGROW>
dF_eval_str=[dF_eval_str,dF_eval_str_array{ii,1:N}]; %#ok<AGROW>
dF_num_eval_str=[dF_num_eval_str,dF_num_eval_str_full(56*(ii-1)+1:56*(ii-1)+7*N-1),';']; %#ok<AGROW>
end
dF_num_eval_str=[dF_num_eval_str(1:end-1),'];'];
%By using eval on each string above, all F functions and derivatives will
%be calculated per iteration.
dx_norm=tol+1;
x=xi;
update_vars_str='v1=x(1);v2=x(2);v3=x(3);v4=x(4);v5=x(5);v6=x(6);v7=x(7);v8=x(8);';
update_vars_str=update_vars_str(1:N*8);
counter=0;
while dx_norm>tol
if counter>max_iter
fprint('Maximum iteration limit reached. x = NaN.\n');
break
end
%Update variables, solve functions, and derivatives at current x
eval([update_vars_str,F_eval_str,dF_eval_str,F_num_eval_str,dF_num_eval_str]);
dx=-dFdx_num\F_num;
x=x+dx;
dx_norm=norm(dx);
counter=counter+1;
end
x(abs(x)<eps)=0;
x(~isfinite(x))=NaN;
end
Nasser M. Abbasi
> Hi,
> I am trying to solve the chemical equilibrium problem by using Newton-Raphson Iterative Solver.
> I download a custom solver and apply it to solve chemical equilibrium problem. I somehow can not get the answer,
> Please help what I'm doing wrong. the m-files are shown below:
>
>
<SNIP> listing of m file.
Ok, lets analyze this a bit. You say "Please help what I'm doing wrong"
but you NEVER show what YOU did!
How is someone going to show what you did wrong, if you do not show what you did?
--Nasser
xingyu
Hi Nasser,
this custom solver works well for the simple problem such as :
-------------------------------------------------------------
%==========================================================================
%Solve an example problem of 3 equations, 3 unknowns and 1 solution:
% 2x + y + 3z = 1
% 2x + 6y + 8z = 3
% 6x + 8y + 18z = 5
clc
close all
fprintf('--------------------------------------------------------------\n')
fprintf('Example 1:\n')
fprintf('\nSolving linear system with 1 solution:\n')
fprintf('2x + y + 3z = 1\n')
fprintf('2x + 6y + 8z = 3\n')
fprintf('6x + 8y + 18z = 5\n')
fprintf('\n')
c1=2;
F = {'c1*v1+v2+3*v3-1';'2*v1+6*v2+8*v3-3';'6*v1+8*v2+18*v3-5'};
dFdx={'c1','1','3';'2','6','8';'6','8','18'};
xi=zeros(3,1);
tol=eps;
max_iter=1000;
fprintf('Initial guess = %s\n',num2str(xi','%7.2f'))
x = nrsolve(F,dFdx,xi,tol,max_iter,c1);
if ~isnan(x)
fprintf('Success!\n')
else
fprintf('No luck!\n')
end
fprintf('[x,y,z] = [%s]\n',num2str(x','%7.2f'))
fprintf('\n')
---------------------------------------------------------------------------------
the result is :
--------------------------------------------------------------
Example 1:
Solving linear system with 1 solution:
2x + y + 3z = 1
2x + 6y + 8z = 3
6x + 8y + 18z = 5
Initial guess = 0.00 0.00 0.00
Success!
[x,y,z] = [0.30 0.40 0.00]
However the error was shown when I try to execute this code:
n=1;m=4;l=0;q=0; T=2000; P=101; F=1; EGR=0;
>> x=equilibrium(n,m,l,q,T,P,F,EGR)
??? Error: Unbalanced or unexpected parenthesis or bracket.
Error in ==> nrsolve at 147
eval([update_vars_str,F_eval_str,dF_eval_str,F_num_eval_str,dF_num_eval_str]);
Error in ==> equilibrium at 32
x = nrsolve(F,dFdx,xi,tol,max_iter,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14);
dpb
...
> However the error was shown when I try to execute this code:
> n=1;m=4;l=0;q=0; T=2000; P=101; F=1; EGR=0;
>>> x=equilibrium(n,m,l,q,T,P,F,EGR)
> ??? Error: Unbalanced or unexpected parenthesis or bracket.
>
> Error in ==> nrsolve at 147
> eval([update_vars_str,F_eval_str,dF_eval_str,F_num_eval_str,dF_num_eval_str]);
>
...
That's a syntax error...look at what
[update_vars_str,F_eval_str,dF_eval_str,F_num_eval_str,dF_num_eval_str]
is. It will somewhere have an ill-formed string.
--
masum billah
can u get solution of this equilibrium equation. I need this equation urgently for my working purpose. If u have please give me.
my address is masu...@gmail.com
thanks
masum
"xingyu " <xuexi...@tamu.edu> wrote in message <imvo1l$fed$1...@fred.mathworks.com>...
vinitp...@gmail.com
I was working on a project titled Egr combustion.So I want to know from where u got this iterative solver.Its urgent.Reply ASAP