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
%             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=','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

On 3/30/2011 10:08 AM, xingyu wrote:

<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

"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;

??? 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);

On 3/30/2011 12:32 PM, xingyu wrote:
...

...

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.

--

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.
my address is masum0...@gmail.com
thanks
masum

...

read more »