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

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)

Answer 1

[Walter Roberson]

    [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.

[roham farhadi]

thank you for your answer , but this equations extracted from thermodynamics law and I don't understand why there is no physical solution

[Walter Roberson]

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]

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

Answer 2

[Konrad]

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

[roham farhadi]

Hi,

thank you for your answer , but when i tried to use abs function in that line , i received this ereor :

Error in sym/subsref (line 881)

R_tilde = builtin('subsref',L_tilde,Idx);

Error in Untitled77 (line 29)

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

[Konrad]

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]

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.