Multiple Equations in Matlab Solver

Hi everyone,

I'm trying to calculate heat change through a solar panel. In order to do this I need 4 equations that I then use an energy balance to solve.

The 4 equations are:

%Heat Flux between glass & sky
Qsky = Eg .* StefanBoltzmann .*(TEa^4 - TEsky^4);
%Heat Flux loss due to forced wind convection
Qwind = Hwind*(TEg - TEa);
%Heat Flux to PV due to radiation
Qrpv = (StefanBoltzmann*(TEpv^4 - TEg^4))/(1/Epv + 1/(Eg-1));
%Conductive tranfer to PV layer
Qpv = Keva*(TEpv - TEg);

The unknown here is TEpv, the code I'm currently using to find this is:

syms Qsky Qwind Qrpv Qpv
%Those 4 equations in here
eqn = Qsky + Qwind == Qrpv + Qpv + J*Ag 
S = solve(eqn,TEpv)

When I use this however no value is sent to S.

The variables Eg, StefanBoltzmann, etc. are defined earlier in the program, they are constants that will be used for more steps along the way, do I have to define them in the syms line? If anone knows how to fix this it would be greatly appreciated.

Many Thanks,

Ieuan

Answer 1

[John D'Errico]

No. You have ONE equation. It just happens that you defined it in 4 pieces, then you added them all up.

If you defined them as constants earlier, then you do NOT want to define those other variables in the syms call. If you did that, it would overwrite the values for those constants. Similarly, putting this line in your code,

syms Qsky Qwind Qrpv Qpv

again overwrites the variables you just created. That just undid all the work you did in the previous lines. So DELETE that use of syms. If you don't believe me, consider this example:

X = 3

X = 3

So at this point, X is defined. It has the value 3. MATLAB tells us that.

syms X
X

X = 

X

Here X no longer has the value 3. It is just a simple symbolic parametre, with comepletely unknown value.

Anyway, when you use that call to solve, what does it show? (Nothing is sent.) MATLAB does not use mail in any form. Why do you think nothing happens?

If all of the other parameters are known, then you have what is merely a 4th degree polynomial in the unknown, so it should be solvable.

[Ieuan Price-davies]

Thank you, have removed the sys and am not getting a value, however it is in terms of z, do you know what this means as a real result?

%Ambient Temp
TEa = 27;
%Sky Temp
TEsky = 0.0552 * TEa^1.5;
%Glass Temp
TEg = 35;
%Convective Heat Transfer coeffciecnt due to wind flow over PVT
Hwind = 5;
%Global Solar Irradiance on PVT
J = 600;

%Heat Flux between glass & sky
Qsky = Eg * StefanBoltzmann*(TEa^4 - TEsky^4);
%Heat Flux loss due to forced wind convection
Qwind = Hwind*(TEg - TEa);
%Heat Flux to PV due to radiation
Qrpv = (StefanBoltzmann*(TEpv^4 - TEg^4))/(1/Epv + 1/(Eg-1));
%Conductive tranfer to PV layer
Qpv = Keva*(TEpv - TEg);

%Energy Balance
eqn = Qsky + Qwind == Qrpv + Qpv + J*Ag;
S = solve(eqn,TEpv);

When I run this now S is given as:

S =
root(z^4 - (122781528554610775556096*z)/2754084132809577 + 54673524965266606675210225/19278588929667039, z, 1)
root(z^4 - (122781528554610775556096*z)/2754084132809577 + 54673524965266606675210225/19278588929667039, z, 2)
root(z^4 - (122781528554610775556096*z)/2754084132809577 + 54673524965266606675210225/19278588929667039, z, 3)
root(z^4 - (122781528554610775556096*z)/2754084132809577 + 54673524965266606675210225/19278588929667039, z, 4)

Do you know how I can turn this into an actual value?

[John D'Errico]

Ok. Now you are getting a result, and have hit a common point of confusion for people using solve. There are several ways to force MATLAB to produce a result there.

TEa = 27;
%Sky Temp
TEsky = 0.0552 * TEa^1.5;
%Glass Temp
TEg = 35;
%Convective Heat Transfer coeffciecnt due to wind flow over PVT
Hwind = 5;
%Global Solar Irradiance on PVT
J = 600;
%Heat Flux between glass & sky
Qsky = Eg * StefanBoltzmann*(TEa^4 - TEsky^4);
%Heat Flux loss due to forced wind convection
Qwind = Hwind*(TEg - TEa);
%Heat Flux to PV due to radiation
Qrpv = (StefanBoltzmann*(TEpv^4 - TEg^4))/(1/Epv + 1/(Eg-1));
%Conductive tranfer to PV layer
Qpv = Keva*(TEpv - TEg);
%Energy Balance
eqn = Qsky + Qwind == Qrpv + Qpv + J*Ag;
S = solve(eqn,TEpv);

However, I cannot reproduce your result here, in an attempt to show you the alternatives. Please learn to go SLOWLY, as the code you show cannot be used. That makes it easier for us to help you.

This next line fails, because in the folllowing line, Eg is undefined. What is Eg? Perhaps you meant TEg.

Qsky = Eg * StefanBoltzmann*(TEa^4 - TEsky^4);

And then we see that StefanBoltzmann is undefined. If I look further, I see Keva make its first appearance, never having been defined.

But if I assume you were to fix all of those problems, then you could have used solve, as:

S = solve(eqn,TEpv,'maxdegree',4);

That tells MATLAB you understand the result is just a 4th degree polynomial. Since I cannot use your code as an example, consider this simple example:

syms x
P = 3*x^4 + 2*x^3 + 5
P =
3*x^4 + 2*x^3 + 5

S = solve(P)
S =
root(z^4 + (2*z^3)/3 + 5/3, z, 1)
root(z^4 + (2*z^3)/3 + 5/3, z, 2)
root(z^4 + (2*z^3)/3 + 5/3, z, 3)
root(z^4 + (2*z^3)/3 + 5/3, z, 4)

As you see, it produces the same result with roots in it. I could have done this instead, however:

S = solve(P,'maxdegree',4)

S =

- (((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)) - (- 20*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + 2*((79^(1/2)*10i)/27 + 10/27)^(1/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - 9*((79^(1/2)*10i)/27 + 10/27)^(2/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + (2*5^(1/2)*6^(1/2)*(1 + 79^(1/2)*1i)^(1/2))/9)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/4)) - 1/6

- (((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)) + (- 20*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + 2*((79^(1/2)*10i)/27 + 10/27)^(1/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - 9*((79^(1/2)*10i)/27 + 10/27)^(2/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + (2*5^(1/2)*6^(1/2)*(1 + 79^(1/2)*1i)^(1/2))/9)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/4)) - 1/6

(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)) - (- 20*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + 2*((79^(1/2)*10i)/27 + 10/27)^(1/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - 9*((79^(1/2)*10i)/27 + 10/27)^(2/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - (2*5^(1/2)*6^(1/2)*(1 + 79^(1/2)*1i)^(1/2))/9)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/4)) - 1/6

(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)) + (- 20*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) + 2*((79^(1/2)*10i)/27 + 10/27)^(1/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - 9*((79^(1/2)*10i)/27 + 10/27)^(2/3)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/2) - (2*5^(1/2)*6^(1/2)*(1 + 79^(1/2)*1i)^(1/2))/9)^(1/2)/(6*((79^(1/2)*10i)/27 + 10/27)^(1/6)*(((79^(1/2)*10i)/27 + 10/27)^(1/3) + 9*((79^(1/2)*10i)/27 + 10/27)^(2/3) + 20)^(1/4)) - 1/6

And the result is a very nasty looking mess of radicals and possibly complex numbers. That is easily reparied using vpa. Double would also have worked.

vpa(S)

ans =

- 0.99536423828141728061460845652417 + 0.76956680335850447298817463829837i

- 0.99536423828141728061460845652417 - 0.76956680335850447298817463829837i

0.66203090494808394728127512319084 + 0.78395127329590816022141746812957i

0.66203090494808394728127512319084 - 0.78395127329590816022141746812957i

double would also have worked.

format long g
double(S)

ans =

-0.995364238281417 + 0.769566803358504i

-0.995364238281417 - 0.769566803358504i

0.662030904948084 + 0.783951273295908i

0.662030904948084 - 0.783951273295908i

[Ieuan Price Davies]

Thanks, I have defined thos values seperately as there are many constants in use in this model. I appreciate your help with the solver though - really helped cheers!