how to solving equation with trials

hi all, i'm a student and i've a assingment and i'm kind a new with matlab.

i have an equation with 2 variables, its like r=(x,T)

i want x and y values to make r=0 and then i will plot a graph to T to x. in more specificly its a rate equation for ammonia synthesis and i'm trying to plot a conversion vs temperature graph.

i try solve with solve(r,T). with this command , i think my equation is kind a complicated and matlab gives an warning; "Warning: Explicit solution could not be found. "

so i'm trying to solve with trials, i want T changes in a to b, and x changes in c to d and which values gives r=0 (10^-3 will be enough to accept 0) display to me this x and T values or even plot a graph.

i know i should use 2 'for loop' like

for T=a:1:b and for x=c:0.01:d

but i dont know how to arrange to r=0 (or 0<=r<=10^-3) and display x,y for r=0.

Answer 1

[Star Strider]

The fminsearch function might work, but assuming T is temperature and you want to solve for different values of x at defined values of T that make r=0, I suggest fzero. You did not post your function so I cannot run code and test it. (I do not know what r returns as output, or how it works.)

[murat]

like i said i'm new with matlab and i really don't know how to add fzero to this equation.

-----------------code

function rNH3=fml(FN2,FH2,FNH3,Ftot)

syms x T 

% x : conversion

% P : total pressure

P=150; % atm

R=1.9872041; % kcal/kmol*K

% fugacity coefficents from literature;

% H2 from Cooper, Shaw and Wones

% N2 and NH3 from Cooper and Newton

fcH2=exp(exp(-3.8402*T^(0.125)+0.541)*P-exp(-0.1263*T^(0.5)-15.98)*P^(2)+300*exp(-0.011901*T-5.941)*(exp(-P/300)-1));

fcN2=0.93431737+0.3101804*10^(-3)*T+0.295896*10^(-3)*P-0.270727*10^(-6)*T^(2)+0.47752507*10^(-6)*P^(2);

fcNH3=0.1438996+0.2028538*10^(-2)*T-0.4487672*10^(-3)*P-0.1142945*10^(-5)*T^(2)+0.2761216*10^(-6)*P^(2);

% The equation of Gillespie and Beattie for equilibrium constant

% logK=(-2.691122*log(T)-50.51925*T+1.848863*10^-7*T.^2+2001.6./T+2.689)

K=10^((-2.691122*log10(T))-(5.519265*10^(-5)*T)+(1.848863*10^(-7)*T^(2))+(2001.6/T)+2.689);

% f0i : pure component fugacity

% f0i=fci*P 

f0N2=P*fcN2;

f0H2=P*fcH2;

f0NH3=P*fcNH3;

% yi : mol fractions 

% N2 + 3 H2 <---> 2 NH3

% from stoichiometry of reaction, mol fractions can be found

yN2=FN2*(1-x)/(Ftot-2*FN2*x);

yH2=(FH2-3*FN2*x)/(Ftot-2*FN2*x);

yNH3=(FNH3+2*FN2*x)/(Ftot-2*FN2*x);

% fi : fugacity

% fi= yi*f0i

fN2=f0N2*yN2;

fH2=f0H2*yH2;

fNH3=f0NH3*yNH3;

% rate formation of ammonia; kmol/h*m^3catbed

rNH3=1.7698*10^(15)*exp(-40765/(R*T))*((K^(2)*fN2*fH2^(1.5)/fNH3)-(fNH3/fH2^(1.5)));
-------------------------------

rNH3 is r in my main question.

T=540:1:800 and x=0:0.0001:1

can you teach me how pls.

btw u can call function fml(1,3,0.13,4.35)

[Star Strider]

In your function line:

function rNH3=fml(FN2,FH2,FNH3,Ftot) syms x T

why ‘syms x T’? Why syms and why in the function line? The x and T you mentioned in your original question do not appear at all in your fml argument list. The only input arguments listed are: (FN2,FH2,FNH3,Ftot).

Also, is rNH3 a scalar (I hope) or a vector?

EDIT — After your last edit (to get the syms statement out of the function line), I still do not understand why x and T fail to appear in the argument list. Declaring them as Symbolic variables (invoking the Symbolic Math Toolbox) is not going to add anything to your code.

I believe what you need is:

function rNH3=fml(x,T,FN2,FH2,FNH3,Ftot)

and eliminate the syms statement.

Also, please use the [{}Code] button to format your code.

[murat]

it isnt in the function line when i c/p it wrote that like and didnt see that mistake.

syms because for define ? it isnt working without it and from documents i read/watch always using syms. and i wrote a code for another reaction equation like that it was more simple equation and with solve command i could get T=f(x) function. so that's kind a why i wrote like this.

this is my basic code i want add 2 'for loop' for T=540:1:800 and x=0:0.0001:1 to rNH3=0 and yes its scalar. but i don't know where/how to add rNH3=0 and so i didnt add loops too.

Edit: but if i use function rNH3=fml(x,T,FN2,FH2,FNH3,Ftot) i should give x and T everytime. but i want plot x and T graph and i need a function like T=f(x) or trial and error values for rNH3=0. i couldn't get a T=f(x) because of warning in my main question, i'm trying trial and error. i think it should be like

rNH3=fml(FN2,FH2,FNH3,Ftot) 

for T=540:1:800

for x=0:0.0001:1

.

.

.

but how could i add rNH3=0 ?

Edit2: i think i understand what u mean, syms is look like unnecessary when using 'for loop'.

[Star Strider]

Knowing that, and adding x and T to the argument list as I noted in my last comment, and assuming FN2,FH2,FNH3,Ftot are defined, I suggest:

FN2 = 3; FH2= 5; FN3 = 7; Ftot = 13;
Tv = 540:1:800;
TrNH3 = [];
for k1 = 1:length(Tv)
    T = Tv(k1);
    f = @(x) fml(x, T, FN2, FH2, FNH3, Ftot);
    rNH3 = fzero(f, 1);
    TrNH3 = [TrNH3;  T  rNH3];
end

The TrHN3 output matrix is a (261x2) array with T as the first column and rNH3 as the second. I tested this loop with a different but similar test function, and it works. Be sure to substitute the correct values for FN2, FH2, FNH3, Ftot for the ones I posted here. The Tv vector has the temperatures. The loop chooses one to pass to fml for each iteration of the loop.

EDIT — You do not add rNH3=0. In the loop, you supply a value of T to fml at each step of the loop, then you let fzero vary x at each value of T until it finds the x that makes rNH3=0 (or as close to zero as the fzero function can approximate it).

I also suggest you vectorize your statements in your function to allow for element-by-element operations. This means substituting (.*) for (*), and similar substitutions resulting in ./ and .^ in place of / and ^.

[murat]

thank you very much.

[Star Strider]

My pleasure!