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

5 visualizzazioni (ultimi 30 giorni)
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.
  9 Commenti
Walter Roberson
Walter Roberson il 18 Apr 2021
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.
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.)

Accedi per commentare.

Risposta accettata

Walter Roberson
Walter Roberson il 17 Apr 2021
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.
  1 Commento
Gabriela Flores Miranda
Gabriela Flores Miranda il 17 Apr 2021
Modificato: Gabriela Flores Miranda il 18 Apr 2021
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!

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Loops and Conditional Statements in Help Center e File Exchange

Prodotti


Release

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by