Matlab produces complex number where it cannot be produced

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;

Answer 1

[Torsten]

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.

[SiSa]

Thank you for your fast response. The problem was vanished and I am happy :) but I want to mention that mathematically there is no difference between phrases 1 and 2:

phrase 1: P_R=P0*((Req^3-b*Req^3)/(R^3-b*Req^3))^k;

phrase 2: P_R=P0*(((Req^3-b*Req^3)/(R^3-b*Req^3))^k1)^k2 with k1=4 , k2=1/3

It means that matlab calculates a signle equation in two different ways that has no mathematical justification.

[Torsten]

In school, I learned that x^y is only defined for x>0.

At university, I learned that x^y is computed as exp(y*log(x)) which gives a complex result for x<0.

So (x^a)^b = x^(a*b) does not necessarily hold for x<0.

What would you expect as result if k = pi in your equations from above ?

[SiSa]

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);

[Steven Lord]

In theory there is no difference between theory and practice. In practice there is.

[Torsten]

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]

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