Matlab produces complex number where it cannot be produced

Question

SiSa 2022-11-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1860743-matlab-produces-complex-number-where-it-cannot-be-produced

评论： SiSa 2022-11-23

采纳的回答： Torsten

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2022-11-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1860743-matlab-produces-complex-number-where-it-cannot-be-produced#answer_1109623

编辑：Torsten 2022-11-23

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

SiSa 2022-11-23

编辑：SiSa 2022-11-23

在 MATLAB Online 中打开

Thank you Torsten:

In school, I learned that x^y is only defined for x>0. : (-2)^3=-8 and (-8)^(1/3)=-2

I do not know how Matlab calculates my eqautions but what I know is that complex numbers are produced when we calculate (even) roots of a negative real number. For odd numbers there is no reason to have a complex number : (-8)^(1/3)=-2.

example with negative x: ((-2)^2)^2= ((-2)*(-2))*((-2)*(-2))=(-2)*(-2)*(-2)*(-2)= (-2)^4=(-2)^(2*2)

I dont know what should I expect when k=pi. :) Do you have experience with Keller Miksis equation for bubble dynamics? I am trying to solve a second order nonlinear ODE with varying parameters using RK method.

and now I am facing the same problem and I do not know how should I formulate my eqaution in a format, that matlab can evaluate it correctly!

Thank you again for your support.

K1=0;
L1=713018093.860669
k1=4;
k2 = 1/3;
k=k1*k2;
P0=1.00*10^5;
Req=69.2*10^(-3);
b=0.0016;
R=0.75*10^(-3);  
v=0
P_R=P0*(((Req^3-b*Req^3)/(R^3-b*Req^3))^k1)^k2
dP_Rdt2=-P0*k*((Req^3-b*Req^3)./((R+K1./2).^3-b.*Req^3))...
  .^(k-1).*(3*(R+K1./2).^2.*(v+L1./2).*(Req^3-b*Req^3))...
 ./(((R+K1./2).^3-b.*Req^3).^2);

Torsten 2022-11-23

编辑：Torsten 2022-11-23

在 MATLAB Online 中打开

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), e.g., is already complex-valued:

(-8)^(1/3)
ans = 1.0000 + 1.7321i
1 + sqrt(3)*1i
ans = 1.0000 + 1.7321i

SiSa 2022-11-23

在 MATLAB Online 中打开

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 