Solving Simultaneous Equations need HELP ASAP

I have the following three equations:

Q_room = h_inside*A_roof*(T_house - T_rin) + E*o*A_roof*((T_house+273)^4 - (T_rin+273)^4)

Q_roof = (k_roof*A_roof*(T_rout - T_infinity))/l_roof

Q_outside = h_outside*A_roof*(T_rout - T_infinity) + E*o*A_roof*((T_rout+273)^4 - (T_sky+273)^4)

*NOTE: All these equations are equal to each other.

I need to solve for T_rin, T_rout. I have tried using the solve command and have only gotten lines of equations with no actual answer. I'm using the student version of matlab.

Answer 1

[Walter Roberson]

No go. The solution involves the roots of a 16th order polynomial. There is no possible closed form algebraic solution for this equation.

I put in some random numbers in to the equations and solved for the roots; I got 4 values for T_rin that were completely real and which were very very close to the house temperature I used; there were 4 other roots that had effectively negligible imaginary components that were also very very close. The other 8 roots of the polynomial reversed the imaginary and real portions.

I had no good idea as to what temperature to use for T_infinity so I tried a variety of temperatures; with the coefficients I used for the other values, the difference due to T_infinity changes was very small, and decreased as T_infinity increased.

[Aaron]

I think you can avoid the 16th order poly if you set Q_roof = Q_room and if you set Q_roof = Q_outside. This will give you 2 equations with 2 unknowns. I'm just not sure how to do this in matlab....

[Walter Roberson]

No, that 16th order polynomial was *with* the values set to be the same.

To set the values to be the same,

solve([Q_room-Q_roof, Q_room-Q_outside], T_rin, T_rout)

This will not give you two equations in two unknown: this will give you two quartic equations in 12 unknowns. Four possible solutions for the first quartic times four possible solutions for the second quartic leads you to 16 possible solutions.

[Aaron]

Is there anyway to give matlab a starting point to and then have it plug in values until the values of the all the Q's are all equal? The other variables are given so I have the numerical value for them. Its just T_rin and T_rout that I need to calculate. As long as I can get the Q's pretty close to equal it'll be good enough.

[Walter Roberson]

A polynomial of even order that has strictly real roots, must have an even number of real roots. The mathematics can't tell the difference between them.

If you are fortunate, your coefficients might be such that you end up with exactly two real roots: if you do, then one will be negative kelvin degrees, which you can then eliminate. But I _suspect_ you will instead end up with four real roots, two of which express negative kelvin degrees, but the other two being mathematically undecidable between.

If you post the numeric values of each of the constants, I will run it through.

[Aaron]

First off thank you for your help.

Here is my code for the 1st part to set everything up. You can just copy and paste it in. The starting value of V is 10.

%Entering in givins

t_roof = 15E-2; %Thickness of the roof [m]

k_roof = 2; %Thermal conductivity of the roof [W/m*K]

w_roof = 15; %Width of the roof [m]

l_roof = 20; %Length of the roof [m]

h_inside = 5; %Heat transfer co. of the inner surface [W/m*K]

T_infinity = 10; %Temperature of the ambient air [C]

T_sky = 100; %Temperature of the sky [K]

T_house = 20; %Temperature of the interior of the house [C]

E = 0.9; %Emissivity of both surfaces of the roof

o = 5.67E-8; %Stefan-Boltmann constant [W/m^2*K^4]

T_rin = sym('T_rin'); %Temperature of the inner surface of the roof

T_rout = sym('T_rout'); %Temperature of the out surface of the roof

%Inputting the parameters of the air

p = 1.246; %Density of air [kg/m^3]

cp = 1006; %Specific heat of air at constant pressure [J/kg*K]

k_air = 0.02439; %Thermal conductivity of the air [W/m*k]

u = 1.778E-5; %Dynamic viscosity of the air [kg/m*s]

%Asking the user to input the velocity of the wind

V = input('the velocity of the wind in m/s ')

%Beginning of Calculations

A_roof = w_roof*l_roof; %Calculate the area of the roof

L_c = l_roof; %Charateristic length of the roof

Re = (p*V*L_c)/u; %Calculate the Reynolds number

Pr = (u*cp)/k_air; %Calculate the Prantl number

Nu= (0.037*Re^0.8 - 871)*Pr^(1/3);

h_outside = (Nu*k_air)/L_c;

[Walter Roberson]

It'll take me a few minutes to convert that to the form I can use in the other package, but I'm working on it.

[Walter Roberson]

The following _might_ be directly usable within a MuPad notebook. I converted all of the floating point constants to fractions to prevent premature evaluation at low precision.

#Entering in givins

t_roof := 15/10^2; #Thickness of the roof [m]

k_roof := 2; #Thermal conductivity of the roof [W/m*K]

w_roof := 15; #Width of the roof [m]

l_roof := 20; #Length of the roof [m]

h_inside := 5; #Heat transfer co. of the inner surface [W/m*K]

T_infinity := 10; #Temperature of the ambient air [C]

T_sky := 100; #Temperature of the sky [K]

T_house := 20; #Temperature of the interior of the house [C]

E := 9/10; #Emissivity of both surfaces of the roof

o := 567/10^10; #Stefan-Boltmann constant [W/m^2*K^4]

#T_rin = sym('T_rin'); #Temperature of the inner surface of the roof

#T_rout = sym('T_rout'); #Temperature of the out surface of the roof

#Inputting the parameters of the air

p := 1246/1000; #Density of air [kg/m^3]

cp := 1006; #Specific heat of air at constant pressure [J/kg*K]

k_air := 2439/100000; #Thermal conductivity of the air [W/m*k]

u := 1778/10^8; #Dynamic viscosity of the air [kg/m*s]

#Asking the user to input the velocity of the wind

#V = input('the velocity of the wind in m/s ')

V := 10;

#Beginning of Calculations

A_roof := w_roof*l_roof; #Calculate the area of the roof

L_c := l_roof; #Charateristic length of the roof

Ren := (p*V*L_c)/u; #Calculate the Reynolds number

Pr := (u*cp)/k_air; #Calculate the Prantl number

Nu := (37/1000*Ren^(4/5) - 871)*Pr^(1/3);

h_outside := (Nu*k_air)/L_c;

Q_room := h_inside*A_roof*(T_house - T_rin) + E*o*A_roof*((T_house+273)^4 - (T_rin+273)^4);

Q_roof := (k_roof*A_roof*(T_rout - T_infinity))/l_roof;

Q_outside := h_outside*A_roof*(T_rout - T_infinity) + E*o*A_roof*((T_rout+273)^4 - (T_sky+273)^4);

solve([Q_room-Q_roof,Q_room-Q_outside],[T_rin, T_rout]);

As I post this, my system is still generating the solution(s)

[Aaron]

When I used this solution it gave me the following error:

??? Undefined function or variable 'T_rin'.

Error in ==> project2 at 33

Q_room = h_inside*A_roof*(T_house - T_rin) +

E*o*A_roof*((T_house+273)^4 - (T_rin+273)^4);

[Walter Roberson]

You did not paste it in to a Mupad Notebook!

http://www.mathworks.com/help/toolbox/symbolic/mupad.html

The text I posted would be pasted *exactly* as-is, := and all.

[Aaron]

Sorry, I've never used this tool box before. I have copied and pasted the exact text into a new mupad. How do I run the program?

[Walter Roberson]

My system is still computing the exact algebraic answers.

The real-valued numeric answers, to 20 decimal places, are:

{T_rin = 19.74391918, T_rout = 35.93498828}

{T_rin = 30.60180844, T_rout = -1094.691378}

{T_rin = -820.1754549, T_rout = 35.93498828}

{T_rin = -824.1010342, T_rout = -1094.691378}

You probably are only interested in the first of these.

I have the general form here expressed in terms of V, but it is a bit messy to copy and explain.

The more exact solution (in terms of fractions) is taking a long time to calculate. I need to go and do other things now and will likely not be able to work on this problem further for a number of hours.

[Walter Roberson]

Sorry, I don't know how you execute in mupad. Possibly just by pressing return.

To get answers in a timely fashion, you may wish to start with

solve(evalf([Q_room-Q_roof,Q_room-Q_outside])

Or perhaps it would be vpa() instead of evalf()

[Aaron]

how did you calculate those values then? did you call the mupad from matlab? I'm a little confused.

[Walter Roberson]

I didn't calculate in MuPad, I calculated in Maple, which has nearly identical syntax and calls for this particular purpose.

The calculation with the exact values filled up 2 Gb of RAM before my system suspended it due to lack of memory, so I do not currently have available anything of higher precision than the previous answers.

[Walter Roberson]

Here is the more accurate solution; evaluate to however many decimal places you need, using whatever method you feel appropriate to maintain internal precision:

T_rin = roots([439472356653575249361, ...

1919615253862816689208848, ...

3930412232284117171155116280, ...

5179585966471049050051618140720, ...

4989652385189759765375690944201176, ...

3701216626436076128923822822819761120, ...

2167738513269420675763523358690452634816, ...

1014175222962471235825403347115223039097920, ...

375521286131149755649059452793869503124493456, ...

107765264413221130466648910711313765537500086400, ...

22729447005300698735785330717673666230367109659520, ...

2981106516667655067061115576212724030895091084534272, ...

112895710200000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) - 37762200000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) * 1780000000^(4/5) * 127^(1/5) + 155495669455565081740318293689101535474615698750607616, ...

123282115538400000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) - 41236322400000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) * 1780000000^(4/5) * 127^(1/5) - 36147832241668783475876325899847545451769037650384097280, ...

1526511069296699278844302786120960097093079818771165798400 + 50484026312974800000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) - 16886274022800000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) * 1780000000^(4/5) * 127^(1/5), ...

-6773301872149600000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) * 1780000000^(4/5) * 127^(1/5) + 20249792788961413600000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) - 25590823497477413523929210576045879829337048228280000512000, ...

154289193741448180667268985246955382748223801505444700160000 + 142556479583312000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3) * 1780000000^(4/5) * 127^(1/5) - 426193786542357392000000000000000000000000000000000 * 447167^(1/3) * 609750^(2/3)])

The observant will note that it is impossible to get a numerically meaningful answer from this without evaluating well beyond 59 decimal places.

T_rout is

-(5103/10000000000)*T_rin^4 - (1393119/2500000000)*T_rin^3 - (1140964461/5000000000)*T_rin^2 - (228827765951/2500000000)*T_rin + 242054864161/125000000

Answer 2

[Sean de Wolski]

once you get the lines of equations from solve use

subs

and

double

to get meaningful values.

[Aaron]

could you be a little more specific on how you would use those two functions I'm having a little trouble with them.