Why do I keep obtaining complex matrices if I have only real values?

So, in this code I'm obtaining the roots of a 4 degree polynomial. This computation is iterated a number of times in a nested "for" loop. I'm interested only in two roots which are saved in the arrays Vbaja(i,j) and Valta(i,j). They are suposed to be 10x10 arrays containing real values because every time I compute the roots of my polynomial I obtain only real values.

However, for some reason I can´t understand, when I display Vbaja and Valta in my command window, I obtain 10x10 complex matrix arrays (even when the roots are real).

The complex part of the values in this arrays is always zero (e.g. 1.0473 + 0.0000i), I don't think it even exist. This is a problem because I want to plot the results but I can´t if I have complex values.

Here's my code:

gamma_min =  0.055; n=10;
gamma= linspace(gamma_min,0.2,n) %[rad]

for i=1:n
    h(i,:)=linspace(1,3048,n)
    for j=1:n
        data= VEL(h(i,j),gamma(i))
        Valta(i,j) = data(1) %here they are not real anymore
        Vbaja(i,j) = data(2) %here they are not real anymore
    end
   
end

figure(1)
surf(Valta,gamma*(180/pi),h, 'LineStyle','none'); grid on; hold on; %unable to plot

I checked the function VEL (where the roots are computed), but the results I have there are real too. Is only when I store them in Valta(i,j) and Vbaja(i,j) when they show as complex values. I don´t know why.

function  velocity = VEL(z,gamma)
%inputs
K=0.05; Cd_0=0.015;Cl_max=2.8;W=324720.18; S=88.2579; 
[rho]=fluidz(z,0)

     function  fluid_prop = fluidz(z,z_ft)
        R=287.05; %Nm/KgK
        To=288.15; %K
        Po=1.01325e5; %N/m^2
        
        if z==0 
            w=[];
            w=z_ft/3.28084;
            z=w;
            fprintf('z = %d [m] \n', w);
        end 
        
        %Presure, N/m^2
        P=Po*(1-((2.25e-5)*z)).^5.2526;
        fprintf('P = %d [N/m^2] \n', P);
        
        %Temperature, K
        if ( 110000 < z) && (z<= 20000) 
             T=216.65;
             fprintf('T = %d [K] \n', T);
        else
            T=To-0.0065*z;
            fprintf('T = %d [K] \n', T);
        end
        
        %Density, kg/m^3
        rho=P/(R*T);
        fprintf('rho = %d [kg/m^3] \n', rho);
        
        fluid_prop=[rho];
        end

Vstall=sqrt(2*W/Cl_max*rho.*S)
%Gliding velocity polynomial
a = (0.5*rho.*S*Cd_0)/W;
b = 0;
c = -gamma
d = 0;
e = (2*K*W)./(rho.*S);
v = sort(roots([a b c d e])) %here the roots are real 
Valta = v(end,:) %4th root.  
Vbaja = v(3,:) %3rd root. 

velocity=[Valta, Vbaja, Vstall];
end

Any suggestion is welcomed.

[David Goodmanson]

Hi Gabriela,

what is the value of rho?

[Gabriela Flores Miranda]

@David Goodmanson rho is a variable which value I obtain from a function, it changes in every iteration.

[rho]=fluidz(z,0)

function  fluid_prop = fluidz(z,z_ft)

[David Goodmanson]

Hello Gabriela, could you provide code that runs as is?

[Gabriela Flores Miranda]

@David Goodmanson yes, I just added the function, now is all there.

[Walter Roberson]

velocity=[Valta, Vbaja, Vstall];

No assignment to Vstall appears anywhere...

[Walter Roberson]

I ran your code without change in R2021a, but I do not see the issue you report. It proceeds right through to the plot.

[Walter Roberson]

I just installed R2017b on my Mac, and ran your software, but I do not see the issue you report.

Question: would you happen to be using an AMD CPU ?

[Gabriela Flores Miranda]

@Walter Roberson I see, yes, I'm using a computer with AMD Ryzen 5

[Walter Roberson]

The bug I was thinking of was for a much older release, and should have been fixed by your release. However it would be interesting to experiment with the fix for it anyhow to see if it makes any difference.

https://www.mathworks.com/support/bugreports/?release_filter=Exists%20in&release=R2015b&utf8=%E2%9C%93&sort_by=date-descending&s=search&search_terms=jaguar

See also https://www.mathworks.com/matlabcentral/answers/396296-what-if-anything-can-be-done-to-optimize-performance-for-modern-amd-cpu-s#answer_401963 (note: that relates more to a slowdown, which was fixed by Mathworks in a later version.)

Answer 1

[Walter Roberson]

When you use root() or solve() of a degree 4 polynomial, some of the terms of the computation inherently involve imaginary components. Ideally, those imaginary components cancel out in the case of roots that should be entirely real, but in practice because of round-off error, they might not cancel out.

You should look at imag(v) to see how small the imaginary component is.

[Gabriela Flores Miranda]

Thank you! I just try it and the imaginary part is 0.

Edit: I checked again the values, this time of all my matrices and you were right, there were a pair of values that actually had very small imaginary components, those were the problem. Thank you very much for commenting!