The old Google Groups will be going away soon, but your browser is incompatible with the new version.
chemical equilibrium with Newton-Raphson Iterative Solver
 There are currently too many topics in this group that display first. To make this topic appear first, remove this option from another topic. There was an error processing your request. Please try again. Standard view   View as tree
 5 messages

From:
To:
Cc:
Followup To:
Subject:
 Validation: For verification purposes please type the characters you see in the picture below or the numbers you hear by clicking the accessibility icon.

More options Mar 30 2011, 1:08 pm
Newsgroups: comp.soft-sys.matlab
From: "xingyu " <xuexingy...@tamu.edu>
Date: Wed, 30 Mar 2011 17:08:05 +0000 (UTC)
Local: Wed, Mar 30 2011 1:08 pm
Subject: chemical equilibrium with Newton-Raphson Iterative Solver
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:

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=r 0;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*c 7*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.1 35)*(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.39 4)/(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*c 8)','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,c 14);

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
%             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=','dF2d v8=';...
'dF3dv1=','dF3dv2=','dF3dv3=','dF3dv4=','dF3dv5=','dF3dv6=','dF3dv7=','dF3d v8=';...
'dF4dv1=','dF4dv2=','dF4dv3=','dF4dv4=','dF4dv5=','dF4dv6=','dF4dv7=','dF4d v8=';...
'dF5dv1=','dF5dv2=','dF5dv3=','dF5dv4=','dF5dv5=','dF5dv6=','dF5dv7=','dF5d v8=';...
'dF6dv1=','dF6dv2=','dF6dv3=','dF6dv4=','dF6dv5=','dF6dv6=','dF6dv7=','dF6d v8=';...
'dF7dv1=','dF7dv2=','dF7dv3=','dF7dv4=','dF7dv5=','dF7dv6=','dF7dv7=','dF7d v8=';...
'dF8dv1=','dF8dv2=','dF8dv3=','dF8dv4=','dF8dv5=','dF8dv6=','dF8dv7=','dF8d v8='};
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,dF1 dv8,',...
'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

To post a message you must first join this group.
You do not have the permission required to post.
More options Mar 30 2011, 1:15 pm
Newsgroups: comp.soft-sys.matlab
From: "Nasser M. Abbasi" <n...@12000.org>
Date: Wed, 30 Mar 2011 10:15:17 -0700
Local: Wed, Mar 30 2011 1:15 pm
Subject: Re: chemical equilibrium with Newton-Raphson Iterative Solver
On 3/30/2011 10:08 AM, xingyu wrote:

> 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,

<SNIP> listing of m file.

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

To post a message you must first join this group.
You do not have the permission required to post.
More options Mar 30 2011, 1:32 pm
Newsgroups: comp.soft-sys.matlab
From: "xingyu " <xuexingy...@tamu.edu>
Date: Wed, 30 Mar 2011 17:32:20 +0000 (UTC)
Local: Wed, Mar 30 2011 1:32 pm
Subject: Re: chemical equilibrium with Newton-Raphson Iterative Solver
"Nasser M. Abbasi" <n...@12000.org> wrote in message <imvof3\$pm...@speranza.aioe.org>...

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,c 14);

To post a message you must first join this group.
You do not have the permission required to post.
More options Mar 30 2011, 1:49 pm
Newsgroups: comp.soft-sys.matlab
From: dpb <n...@non.net>
Date: Wed, 30 Mar 2011 12:49:54 -0500
Local: Wed, Mar 30 2011 1:49 pm
Subject: Re: chemical equilibrium with Newton-Raphson Iterative Solver
On 3/30/2011 12:32 PM, xingyu wrote:
...

> 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.

--

To post a message you must first join this group.
You do not have the permission required to post.
More options May 11 2012, 2:32 am
Newsgroups: comp.soft-sys.matlab
From: "masum billah" <masum0...@gmail.com>
Date: Fri, 11 May 2012 06:32:30 +0000 (UTC)
Local: Fri, May 11 2012 2:32 am
Subject: Re: chemical equilibrium with Newton-Raphson Iterative Solver
Dear sir,
can u get solution of this equilibrium equation. I need this equation urgently for my working purpose. If u have please give me.
thanks
masum

...