could someone please tell me why my code returns complex answers while it is about room temperature ?

3 次查看（过去 30 天）
roham farhadi
roham farhadi 2021-8-12
编辑： Walter Roberson 2021-8-12
clc , clear
% inputs :
ho = 34 ;
k1 = 25;
cp1 = 10 ;
c1 = 12682.7 ;
ro = 3410.9 ;
q = 0 ;
desired_time=25; %sec
dt = 0.5 ;
n = desired_time/dt ;
wall_thickness = 219; %mm
number_of_nodes = 5 ;
dx = ( wall_thickness / (number_of_nodes - 1) )*10^-3 ;
outside_temperature = 10 ;
To = sym(outside_temperature) ;
initial_room_temperature = 25 ;
initial_wall_temperature = 20 ;
% equations :
T_p_1_n = zeros(1,number_of_nodes+1)+ initial_wall_temperature ;
T_p_1_n(1,number_of_nodes+1) = initial_room_temperature ;
T_p_1 = sym(T_p_1_n);
T_p_a = sym(zeros(n,number_of_nodes+1)) ;
T_p = sym('T_p%d' ,[1,number_of_nodes+1],'real');
eqn = zeros(1,number_of_nodes+1) ;
eqn = sym(eqn) ;
for o=1:n
eqn(1) = ho * (To - T_p(1)) - (k1 / dx)*(T_p(1) - T_p(2)) - ((ro * cp1 * (dx/2))/dt) * (T_p(1) -T_p_1(1)) ; % outside of the wall node
for i = 2:number_of_nodes-1 %inside wall nodes
eqn(i) = (k1/dx) * (T_p(i-1) - T_p(i)) - (k1/dx)*(T_p(i) - T_p(i+1)) - ((ro * cp1 * dx )/dt) * (T_p(i)-T_p_1(i)) ;
end
eqn(number_of_nodes) = (k1/dx) / ( T_p(number_of_nodes-1) - T_p(number_of_nodes)) - (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes) - T_p(number_of_nodes+1)) - ((ro * cp1 * (dx/2))/dt)*(T_p(number_of_nodes)-T_p_1(number_of_nodes)) ; % inside of the wall node
eqn(number_of_nodes+1) = q - (28.9805) * (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes+1)-T_p(number_of_nodes)) - (c1/dt) * (T_p(number_of_nodes+1)-T_p_1(number_of_nodes+1)) ; % room air equation
[x1 , x2 , x3 , x4 , x5 , x6] = vpasolve( eqn,symvar(eqn) )
T_p_1 = [x1 , x2 , x3 , x4 , x5 , x6]
end
double(T_p_1)

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2021-8-12
[x1 , x2 , x3 , x4 , x5 , x6] = vpasolve( eqn,symvar(eqn) )
if you breakpoint at that line, and run to there, then you can do
syms(symvar(eqn));
sol1_4 = solve(eqn(1:4), [T_p1, T_p2, T_p3, T_p4]) %easy
neqn = simplify(subs(eqn(5:end), sol1_4))
It is then possible to solve neqn(1) for T_p5 -- or at least Maple can do it. The result is T_p5 expressed in terms of the roots of a polynomial of degree 56 with non-zero coefficients at degree 56, 50, 31, 25, and 0.
Now you could hope for a moment to factor into a polynomial of degree 31 times a polynomial of degree 25, but that would get you entries with degree 56, 31, 25, and 0... but no entries of degree 50. And if you try to include (polynomial of degree 25)^2 then you will start to get additional powers; conjugate roots can get you degree 50 but not the 25 needed to add to 31 to get 56. Anyhow, you probably cannot factor.
Once you have the root() form of the polynomial of degree 56, you can substitute it into neqn(2) to get the equation that needs to hold for T_p6 . You are not going to be able to solve that symbolically. If you try to solve it numerically... you have problems. If you plot the imaginary part, it looks like there might be only one place the imaginary part is 0, near 20.749... but at that location, the real part is over 100000.
It appears that all of the roots for T_p6 are imaginary.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2021-8-12
Ran in:
clc , clear
% inputs :
ho = 34 ;
k1 = 25;
cp1 = 10 ;
c1 = 12682.7 ;
ro = 3410.9 ;
q = 0 ;
desired_time=25; %sec
dt = 0.5 ;
n = desired_time/dt ;
wall_thickness = 219; %mm
number_of_nodes = 5 ;
dx = ( wall_thickness / (number_of_nodes - 1) )*10^-3 ;
outside_temperature = 10 ;
To = sym(outside_temperature) ;
initial_room_temperature = 25 ;
initial_wall_temperature = 20 ;
% equations :
T_p_1_n = zeros(1,number_of_nodes+1)+ initial_wall_temperature ;
T_p_1_n(1,number_of_nodes+1) = initial_room_temperature ;
T_p_1 = sym(T_p_1_n);
T_p_a = sym(zeros(n,number_of_nodes+1)) ;
T_p = sym('T_p%d' ,[1,number_of_nodes+1],'real');
eqn = zeros(1,number_of_nodes+1) ;
eqn = sym(eqn) ;
for o=1:1 %n
eqn(1) = ho * (To - T_p(1)) - (k1 / dx)*(T_p(1) - T_p(2)) - ((ro * cp1 * (dx/2))/dt) * (T_p(1) -T_p_1(1)) ; % outside of the wall node
for i = 2:number_of_nodes-1 %inside wall nodes
eqn(i) = (k1/dx) * (T_p(i-1) - T_p(i)) - (k1/dx)*(T_p(i) - T_p(i+1)) - ((ro * cp1 * dx )/dt) * (T_p(i)-T_p_1(i)) ;
end
eqn(number_of_nodes) = (k1/dx) / ( T_p(number_of_nodes-1) - T_p(number_of_nodes)) - (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes) - T_p(number_of_nodes+1)) - ((ro * cp1 * (dx/2))/dt)*(T_p(number_of_nodes)-T_p_1(number_of_nodes)) ; % inside of the wall node
eqn(number_of_nodes+1) = q - (28.9805) * (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes+1)-T_p(number_of_nodes)) - (c1/dt) * (T_p(number_of_nodes+1)-T_p_1(number_of_nodes+1)) ; % room air equation
%{
[x1 , x2 , x3 , x4 , x5 , x6] = vpasolve( eqn,symvar(eqn) )
T_p_1 = [x1 , x2 , x3 , x4 , x5 , x6]
%}
eqn(1:4).'
syms(symvar(eqn));
sol1_4 = solve(eqn(1:4), [T_p1, T_p2, T_p3, T_p4]) %easy
subs([T_p1; T_p2; T_p3; T_p4], sol1_4)
neqn = simplify(subs(eqn(5:end), sol1_4))
obj(T_p5, T_p6) = simplify(sum(neqn.^2))
fsurf(obj, [-1000 1000 -100 100])
xlabel('T_p5'); ylabel('T_p6')
vpa(obj(0,0))
vpa(obj(-200,0))
end
ans = 
sol1_4 = struct with fields:
T_p1: [1×1 sym] T_p2: [1×1 sym] T_p3: [1×1 sym] T_p4: [1×1 sym]
ans = 
neqn = 
obj(T_p5, T_p6) = 
ans = 
403524066470.01354894148238674629
ans = 
The sol1_4 part shows that you can solve the first 4 equations for the first 4 variables in an entirely linear way. There is no ambiguity there about possibly having chosen the wrong roots of a polynomial or exponential or whatever. When you look at the output of eqn(1:4) it is obvious that that portion is linear and should be solved with linear results.
That reduces you down to two nonlinear equations in two unknowns, as shown in neqn.
Each of the resulting equation is of the form of an expression that is implicitly required to be 0 for there to be a solution for the equations.
When you have two expressions, X and Y that are each individually required to be 0, then it follows that X^2 must be 0 and Y^2 must be 0, and that X^2+Y^2 must be 0 at the point of solution. But where X and Y might individually be positive or negative, X^2 and Y^2 cannot be negative and X^2+Y^2 cannot be negative (well, not unless there are complex values involved). So to see if there is a possible solution for the pair of equations, you can examine the surface X^2+Y^2 and if it ever reaches 0 then at those points there is a potential solution -- and if X^2+Y^2 never reaches 0 then it follows that the two equations cannot be solved together (even though there might be solutions for each of them.)
We then proceed to fsurf() the X^2+Y^2 and we can see that it has a minima near T_p6 == 0 and it might be tilting along T_p5 direction. But notice that it vanishes; examining neqn we see that T_p5 must be greater than T_p6 or else we get complex values.
If you look closely at the first neqn you see we have a division by an expression. Do we have a discontinuity?
critical5 = solve(-219/children(neqn(1), 2)/100000)
vpa(critical5)
ans = 
19.999842239001277067602483203082
Yes, we do. What effect does that have on the values?
fsurf(obj, [19 20 0 20]); xlabel('T_p5'); ylabel('T_p6')
Effectively we have two different systems, before and after the critical value, and after the critical value the results get larger. What about before the critical point?
fsurf(obj, double([-20, critical5*0.99, -20, critical5])); xlabel('T_p5'); ylabel('T_p6')
... Just not good enough.
Walter Roberson
Walter Roberson 2021-8-12
编辑：Walter Roberson 2021-8-12
Now, these plots do depend upon MATLAB having solved equations correctly... namely the linear equations in the first 4 parts. But even if it got the exact values of the linear solutions wrong, it would have to have gotten the form of the solutions wrong for it to have ruined the possibility of a solution for the remaining two equations. And my work with the equations in Maple shows that MATLAB did not get the form of the solutions wrong for the first four equations.
The surface plots do not depend upon MATLAB being able to solve the remaining two equations: only upon it being able to evaluate the equations at points and graph the results correctly. MATLAB is not without its flaws in graphing equations that require high precision because they have terms that nearly balance, but this is not one of those situations.
At this point, we are faced with one of the following possibilities:
  • the laws of thermodynamics might be wrong
  • fundamental Mathematics might be wrong
  • you might have made a mistake in deriving the equations
  • you might have made a mistake in transcribing the equations into MATLAB
  • MATLAB and Maple might both have major major mistakes in calculating the form of solutions to simple linear equations
  • MATLAB and Maple might both have significant problems in evaluating the objective function obj() that is calculated here

请先登录，再进行评论。

更多回答（1 个）

Konrad
Konrad 2021-8-12
Hi,
One quick guess:
in your equation (eqn(number_of_nodes) = ...) you are raising to the power of .24, which means you are taking the 4.166..th root. If the term inside the parentheses
(T_p(number_of_nodes) - T_p(number_of_nodes+1))^.24
is negative, you get a complex result.
Best, Konrad
  3 个评论
显示 1更早的评论
Konrad
Konrad 2021-8-12
You didn't explain what your code is supposed to do, so it's difficult to say how to change it to make it work. I don't know if simply using abs() is the right way...
However, the problem is that now vpasolve() does not find a solution anymore, so T_p_1 is empty and trying to index T_p_1 in the next iteration errors.
roham farhadi
roham farhadi 2021-8-12
sorry for bad explanations ,
this code use explicit finite-difference method to find transient room air and walls temperature so the basic is in each time step a system of equations must be solved to find T_ P , and then T_P_1 (which is initial conditions ) must be replaced with T_p for the next time step , and this continued to reach the desired time (like 25 sec ) , this is what the main loop do.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by