# Matlab produces complex number where it cannot be produced

2 views (last 30 days)
SiSa on 23 Nov 2022
Commented: SiSa on 23 Nov 2022
I have a problem with following code: Why complex numbers are generated?
k=4/3;
P0=1.00*10^5;
Req=69.2*10^(-3);
b=0.0016;
R=0.75*10^(-3);
P_R=P0*((Req^3-b*Req^3)/(R^3-b*Req^3))^k;
P_R2=P0.*((Req^3-b*Req^3)./(R^3-b*Req^3))^k;
P_R3=P0*((Req.^3-b*Req.^3)/(R^3-b*Req.^3)).^k;
P_R4=P0.*((Req.^3-b.*Req.^3)/(R.^3-b.*Req.^3)).^k;

Torsten on 23 Nov 2022
Edited: Torsten on 23 Nov 2022
Because
(Req.^3-b.*Req.^3)/(R.^3-b.*Req.^3)
is negative.
If you use
k1 = 4;
k2 = 1/3;
P0=1.00*10^5;
Req=69.2*10^(-3);
b=0.0016;
R=0.75*10^(-3);
P_R=P0*(((Req^3-b*Req^3)/(R^3-b*Req^3))^k1)^k2
P_R = 5.3379e+08
P_R2=P0.*(((Req^3-b*Req^3)./(R^3-b*Req^3))^k1)^k2
P_R2 = 5.3379e+08
P_R3=P0*(((Req.^3-b*Req.^3)/(R^3-b*Req.^3)).^k1)^k2
P_R3 = 5.3379e+08
P_R4=P0.*(((Req.^3-b.*Req.^3)/(R.^3-b.*Req.^3)).^k1)^k2
P_R4 = 5.3379e+08
the problem will vanish.
Torsten on 23 Nov 2022
Edited: Torsten on 23 Nov 2022
In school, I learned that x^y is only defined for x>0. : (-2)^3=-8 and (-8)^(1/3)=-2
ok, I should have said x^y for y not an integer.
(-8)^(1/3)
ans = 1.0000 + 1.7321i
1 + sqrt(3)*1i
ans = 1.0000 + 1.7321i
SiSa on 23 Nov 2022
Thank you!
I changed my code to evalute correct values:
if ((Req^3-b*Req^3)/(R.^3-b.*Req^3))<0
dP_Rdt=P0*k1*k2*((b*Req^3-Req^3)./(R.^3-b.*Req^3)).^(k2).*(3*R.^2.*v.*(Req^3-b*Req^3))./((R.^3-b.*Req^3).^2);
else
dP_Rdt=-P0*k1*k2*((Req^3-b*Req^3)./(R.^3-b.*Req^3))...
.^(k1*k2-1).*(3*R.^2.*v.*(Req^3-b*Req^3))./((R.^3-b.*Req^3).^2);
end