Solving Non-Linear Equations in Matlab by Iteration

I am doing my matlab project and I am stuck in this part. This is a part of my code:

dhO2out= r*(O2a*(taf-tref)+(O2b/2)*(taf^2 - tref)+O2c*log(taf-tref)+O2d*(taf-tref)*log(taf-tref)-(taf-tref)+(O2e/3)*(taf^3 - tref));

dhN2out= r*(N2a*(taf-tref)+(N2b/2)*(taf^2 - tref)+N2c*log(taf-tref)+N2d*(taf-tref)*log(taf-tref)-(taf-tref)+(N2e/3)*(taf^3 - tref));

dhCO2 = r*(CO2a*(taf-tref)+(CO2b/2)*(taf^2 - tref)+CO2c*log(taf-tref)+CO2d*(taf-tref)*log(taf-tref)-(taf-tref)+(CO2e/3)*(taf^3 - tref));

dhH2O = r*(H2Oa*(taf-tref)+(H2Ob/2)*(taf^2 - tref)+H2Oc*log(taf-tref)+H2Od*(taf-tref)*log(taf-tref)-(taf-tref)+(H2Oe/3)*(taf^3 - tref));

syms taf

hp = dhO2out+dhN2out+dhCO2+dhH2O == sum (v10);

s = vpasolve (hp,taf,'Random',true);

disp ('The adiabatic flame temperature is (K): ')

disp (s)

I am using the function vpasolve, but it's giving me unspecified answer. I asked my professor and he told me that I can use iteration (like Newton raphson method) to find taf. But how do I do this in matlab? - Pls help, THANK YOU.

Answer 1

[Torsten]

First step:

Give numerical values to all variables except "taf"

Second step:

Plot dhO2out+dhN2out+dhCO2+dhH2O - sum (v10) as a function of taf to approximately locate the zero:

dhO2out= @(taf) r*(O2a*(taf-tref)+(O2b/2)*(taf.^2 - tref)+O2c*log(taf-tref)+O2d*(taf-tref).*log(taf-tref)-(taf-tref)+(O2e/3)*(taf.^3 - tref))
dhN2out= @(taf) r*(N2a*(taf-tref)+(N2b/2)*(taf.^2 - tref)+N2c*log(taf-tref)+N2d*(taf-tref).*log(taf-tref)-(taf-tref)+(N2e/3)*(taf.^3 - tref));
dhCO2 = @(taf) r*(CO2a*(taf-tref)+(CO2b/2)*(taf.^2 - tref)+CO2c*log(taf-tref)+CO2d*(taf-tref).*log(taf-tref)-(taf-tref)+(CO2e/3)*(taf.^3 - tref)); 
dhH2O = @(taf) r*(H2Oa*(taf-tref)+(H2Ob/2)*(taf.^2 - tref)+H2Oc*log(taf-tref)+H2Od*(taf-tref).*log(taf-tref)-(taf-tref)+(H2Oe/3)*(taf.^3 - tref));
hp = @(taf) dhO2out(taf)+dhN2out(taf)+dhCO2(taf)+dhH2O(taf) - sum (v10);
taf = tref+1:tref+500;
plot(taf,hp(taf))

Third step:

Locate taf from the graph, set it as taf0 and use fzero to solve:

taf0 = ...;
taf = fzero(hp,taf0)

[Veronica Gabriel]

Thank you for the response!

My code has already generated the graph, but I am not sure on how I can locate the taf from the graph?

[Torsten]

Can't you see a point where the graph of

dhO2out(taf)+dhN2out(taf)+dhCO2(taf)+dhH2O(taf) - sum (v10)

crosses the x-axis ? This is taf.

[Veronica Gabriel]

Hi, thank you! I already found it - but it is now giving me an error when I run it.

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -84 is -3111.2844+26.17844i.)
Check function or try again with a different starting value.

taf =

   NaN

[Torsten]

Where does the plot of

dhO2out(taf)+dhN2out(taf)+dhCO2(taf)+dhH2O(taf) - sum (v10) 

cross the x-axis when you do the plot as I wrote you should ?