How to iteratively solve when equations are dependent on each other.

I have two equations, I know the varibles T_oQlast, K_dot_nom, and Sigma_ysrtM which are attached. The top equation is typically used to solve for T_oQlast, but in my case I need to solve for T_o_calc so I rearranged it as you can see. GAMMA is calculated based on T_o_calc.

What funtion can I used to solve for T_o_calc? I have thought about using an interative solution, select a starting GAMMA1 of 50, then calculate T_o_calc for that GAMMA1. Recacluate GAMMA1 for this new T_o_calc. But I can't figure out a way to loop the code until it converges on a solution.

T_o_calc = (((T_oQlast+273.15)*GAMMA1)/(GAMMA1-log(K_dot_nom/0.91)))-273.15;

GAMMA1 = (9.9 * exp((((T_o_calc+273.15)/190)^1.66)+(((Sigma_ysrtM)/722)^1.09)));

Heres the Vpasolve code I tied out. equLeft is the T_o_calc equation solved for GAMMA1, equRight is the equation for GAMMA1. This didn't work for some reason, it gave back no results.

syms T_o_calc

eqnLeft = -(273.15*log(K_dot_nom/0.91)-T_oQlast*log(K_dot_nom/0.91))/(T_o_calc-T_oQlast);

eqnRight = (9.9 * exp((((T_o_calc+273.15)/190)^1.66)+(((Sigma_ysrtM)/722)^1.09)));

T_o = vpasolve(eqnLeft == eqnRight, T_o_calc);

[Sam Chak]

Hi @John, Can you verify if the plots are correct before attempting to solve the problem?

load('Sigma_ysrtM.mat')
load('T_oQlast.mat')
load('K_dot_nom.mat')

syms T_o_calc
eqnLeft  = -(273.15*log(K_dot_nom) - T_oQ*log(K_dot_nom))/(T_o_calc + 546.3 - T_oQ);
eqnRight = (9.9*exp((((T_o_calc + 273.15)/190)^1.66) + (((Sigma_ysrtM)/722)^1.09)));

hold on
fplot(eqnLeft,  [-650, -100])
fplot(eqnRight, [-650, -100])
hold off

legend('eqnL', 'eqnR', 'location', 'northwest')
grid on

●

[John]

Sam, I made a typo in on of the equations it should now be fixed.

load('Sigma_ysrtM.mat')

load('T_oQlast.mat')

load('K_dot_nom.mat')

>> syms T_o_calc

eqnLeft = -(273.15*log(K_dot_nom) - T_oQ*log(K_dot_nom))/(T_o_calc - T_oQ);

eqnRight = (9.9*exp((((T_o_calc + 273.15)/190)^1.66) + (((Sigma_ysrtM)/722)^1.09)));

hold on

fplot(eqnLeft, [-650, -100])

fplot(eqnRight, [-650, -100])

hold off

legend('eqnL', 'eqnR', 'location', 'northwest')

grid on

vpasolve is returning this as the T_o_calc solution,

- 266.69623322637174915300510886183 + 300.47327450035187608909376093367i

where -144.5 should be the actual solution. I used solver in Excel to get the -145.5.

[Sam Chak]

Hi @John, Based on your corrected equations, there is an intersection around T_o_calc = -170. You can try out the approaches by @Star Strider and @Torsten.

load('Sigma_ysrtM.mat')
load('T_oQlast.mat')
load('K_dot_nom.mat')

syms T_o_calc
eqnLeft  = -(273.15*log(K_dot_nom) - T_oQ*log(K_dot_nom))/(T_o_calc - T_oQ);
eqnRight = (9.9*exp((((T_o_calc + 273.15)/190)^1.66) + (((Sigma_ysrtM)/722)^1.09)));

hold on
fplot(eqnLeft,  [-220, -140])
fplot(eqnRight, [-220, -140])
hold off

ylim([0, 100])
legend('eqnL', 'eqnR', 'location', 'northwest')
grid on

[John]

Third times the charm here, I broke things down to make sure each indivdual term was correct and found another error.

syms T_o_calc

Sigma_ysrtM = sigma_ysrt * 6.895;

K_dot_nom_Metric = table2array(inputs(1,10)); %attached is this value

lnK = log(K_dot_nom_Metric);

Term1 = 273.15 * lnK;

Term2 = T_oQ * lnK;

Term3 = ((Sigma_ysrtM)/722)^1.09;

eqnLeft = -(Term1+Term2)/(T_o_calc-T_oQ);

eqnRight = 9.9*exp((((T_o_calc+273.15)/190)^1.66) + (Term3));

T_o = vpasolve(eqnLeft == eqnRight, T_o_calc);

Answer is -145.3 if you plot it, which matches my excel solver, but vpasolve doesn't work still.

Answer 1

[Sam Chak]

Hi @John

You can specify the range so that 'vpasolve' searches the solution in that region.

load('Sigma_ysrtM.mat')
load('T_oQlast.mat')
load('K_dot_nom_Metric')

syms T_o_calc

lnK      = log(K_dot_nom_Metric);
Term1    = 273.15 * lnK;
Term2    = T_oQ * lnK;
Term3    = ((Sigma_ysrtM)/722)^1.09;
eqnLeft  = -(Term1+Term2)/(T_o_calc-T_oQ);
eqnRight = 9.9*exp((((T_o_calc+273.15)/190)^1.66) + (Term3));

hold on
fplot(eqnLeft,  [-200, -130])
fplot(eqnRight, [-200, -130])
hold off
ylim([0, 100])
legend('eqnL', 'eqnR', 'location', 'northwest')
grid on

●

sol = vpasolve(eqnLeft == eqnRight, T_o_calc, [-150 140])

sol = 

Answer 2

[Torsten]

Sigma_ysrtM = load("Sigma_ysrtM.mat").Sigma_ysrtM;
T_oQlast = load("T_oQlast.mat").T_oQ;
K_dot_nom = load("K_dot_nom.mat").K_dot_nom;
T_o_calc = fsolve(@(T_o_calc)fun(T_o_calc,Sigma_ysrtM,T_oQlast,K_dot_nom),300)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

T_o_calc = -123.4789

function res = fun(T_o_calc,Sigma_ysrtM,T_oQlast,K_dot_nom)
  GAMMA1 = (9.9 * exp((((T_o_calc+273.15)/190)^1.66)+(((Sigma_ysrtM)/722)^1.09)));
  res = T_o_calc - ((((T_oQlast+273.15)*GAMMA1)/(GAMMA1-log(K_dot_nom/0.91)))-273.15);
end

[John]

T_o_calc = fsolve(@(T_o_calc)fun(T_o_calc,Sigma_ysrtM,T_oQlast,K_dot_nom_Metric),300)

function res = fun(T_o_calc,Sigma_ysrtM,T_oQlast,K_dot_nom_Metric)

GAMMA1 = (9.9 * exp((((T_o_calc+273.15)/190)^1.66)+(((Sigma_ysrtM)/722)^1.09)));

res = T_o_calc - ((((T_oQlast+273.15)*GAMMA1)/(GAMMA1-log(K_dot_nom_Metric)))-273.15);

end

I get -123.4790, but it should be aroun -145.3

Answer 3

[John]

syms T_o_calc

Sigma_ysrtM = sigma_ysrt * 6.895;

K_dot_nom_Metric = table2array(inputs(1,10)); %Load from attached

lnK = log(K_dot_nom_Metric);

Term1 = 273.15 * lnK;

Term2 = T_oQ * lnK;

Term3 = ((Sigma_ysrtM)/722)^1.09;

eqnLeft = -(Term1+Term2)/(T_o_calc-T_oQ);

eqnRight = 9.9*exp((((T_o_calc+273.15)/190)^1.66) + (Term3));

T_o = vpasolve(eqnLeft == eqnRight, T_o_calc,[-500,500]); %ADDED bounds

Returns correct value and NOT an imaginary number which is the issue I had with the origianl code.