why does similar code generate complex numbers?

Question

IB00 2020-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573925-why-does-similar-code-generate-complex-numbers

评论： Walter Roberson 2020-8-2

In the following code, the array 'f1' contains real numbers while the array 'f2' contains complex numbers. From what I understand and am expecting, both arrays should contain identical data, consisting only of real numbers. But this is not what I get. Instead, the array 'f2', generated using exponential notation, contains complex numbers and the question is why is this happening?

x = [-27 -8 -1 1 8 27];
f1   = nthroot(x, 3);
f2   = x .^ (1/3);

The array 'f1' contains the expected results:

-3 -2 -1 1 2 3

For the array 'f2', here are the numbers generated:

1.5000 + 2.5981i 1.0000 + 1.7321i 0.5000 + 0.8660i 1.0000 + 0.0000i 2.0000 + 0.0000i 3.0000 + 0.0000i

The array 'f2' unexpectedly contains complex numbers. The last 3 numbers give the correct/expected result but the first 3 numbers are not the correct/expected result. That is, the first 3 numbers in 'f2' are NOT the cubed root of the first 3 numbers in 'x'. Why is this happening?

In trying to debug/understand the problem, when I manually enter the negative numbers of 'x' in the command window then I get the expected results:

-27 ^ (1/3) ---> ans = -3

-8 ^ (1/3) ---> ans = -2

-1 ^ (1/3) ---> ans = -1

However, when I use array syntax in the command window and select the same negative numbers from 'x' then I get unexpected and incorrect results:

x(1) ^ (1/3) ---> ans = 1.5000 + 2.5981i

x(2) ^ (1/3) ---> ans = 1.0000 + 1.7321i

x(3) ^ (1/3) ---> ans = 0.5000 + 0.8660i

I have no experience working with complex numbers. The question I am asking is why do I get incorrect results (i.e. in complex number format) when using exponential notation in my script?

Thanks to all who respond

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

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2020-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573925-why-does-similar-code-generate-complex-numbers#answer_473962

When x is negative, x^n is defined by MATLAB as exp(n*log(x)) . log(x) when x is negative, is log(-x) + pi * 1i . n times that is n*log(-x) + n*pi*1i . exp() of that is going to be complex except when n is an integer (in which case n*pi becomes an exact multiple of pi)

Thus for x negative, f2 = x .^ (1/3); is going to be exp(log(-x)/3 + pi/3*1i) which is going to be complex valued.

nthroot() is defined in terms of real roots, with it being an error for the value to be negative and the power to not be an odd integer.

3 个评论
显示 1更早的评论隐藏 1更早的评论

IB00 2020-8-2

编辑：IB00 2020-8-2

Operator precedence! Good catch. My mistake. Thanks Walter.

I just saw your first post. You wrote, "When x is negative, x^n is defined by MATLAB as exp(n*log(x))".

1) Does Matlab define x^n the same way when x is positive? I assume yes but, given I'm in a bit over my head, I am asking just to confirm.

2) Would you be able to refer me to the section in the Matlab documentation where I can read more about how "x^n is defined when x is negative"?

Thanks again Walter

Walter Roberson 2020-8-2

If I recall correctly, MATLAB actually uses log() and exp() when the base is negative, even for integer exponents . At some point I saw someone demonstrate this based on precision, a case where log() / exp() exactly matched the observed result but repeated multiplication did not exactly match.. but that would have been a number of years ago now.

MATLAB does not definitely use log() and exp() when the base is positive. I believe that it does so when the exponent is not an integer, but for the case where the exponent is a hard-coded integer, it might not. People have done some timing tests, which I have mostly forgotten the details of now. I seem to recall that X.^2 has timing consistent with being implemented as X*X . My memory seems to think that X.^4 had timing consistent with X*X*X*X, and not consistent with temp = X*X; temp*temp. My memory is telling me that the precision results were consistent with repeated multiplication, not with log() / exp() for positive integer expontnets. Someone submitted a FEX contribution to do faster integer exponents by decomposing the power into binary.

I do not know what the algorithm is for nthroot; I have seen cases where nthroot(x, N) did not exactly equal x.^(1/N), and my memory is telling me that x.^(1/N) in that case was exactly equal to the log/exp solution.

请先登录，再进行评论。

Answer 2

John D'Errico 2020-8-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573925-why-does-similar-code-generate-complex-numbers#answer_473935

编辑：John D'Errico 2020-8-1

What, for example, are the THREE cube roots of a negative number? We can find them all as solutions to the equation

syms x
xsol = solve(x^3 == -4,'maxdegree',3)
xsol =
                        -2^(2/3)
  2^(2/3)*((3^(1/2)*1i)/2 + 1/2)
 -2^(2/3)*((3^(1/2)*1i)/2 - 1/2)
vpa(xsol)
ans =
                          -1.5874010519681994748
 0.79370052598409973738 + 1.3747296369986026264i
 0.79370052598409973738 - 1.3747296369986026264i

So -4 has 3 cube roots. We can plot them in the complex plane.

plot(double(xsol),'o')
xline(0);
yline(0);

Which one has the smallest phase angle in the complex plane?

angle(xsol)
ans =
    pi
  pi/3
 -pi/3

So the second root it found is the one that lives at an angle of pi/3 in the complex plane, or 60 degrees for the degree fanatics. This is the root with the smallest complex phase angle, and is the one chosen by power.

NTHROOT explicitly looks for the real root.

nthroot(-4,3)
ans =
      -1.5874

It tells you that.

 help nthroot
 nthroot Real n-th root of real numbers.
 
    nthroot(X, N) returns the real Nth root of the elements of X.
    Both X and N must be real, and if X is negative, N must be an odd integer.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

IB00 2020-8-2

Thanks for helping out. I may be wrong but I believe you have not understood the question I was asking. I have just edited my original post in order to specfically explain what the problem is. If you have a chance, can you review my edit. Thanks.

请先登录，再进行评论。