Imaginary output from arccosine function when the input is close to -1

I am running a shell buckling simulation, where, I need to find the angle between two vectors. The first vector, m, is [0.000128276283345364, -0.000167881251734009, 7.36224543533651e-06], whereas, the second vector, n, is [-0.000128255846266013, 0.000167854504759308, -7.36107040874295e-06]. To find the angle between these vectors, I am using the following command:

m = [0.000128276283345364, -0.000167881251734009, 7.36224543533651e-06]
n = [-0.000128255846266013, 0.000167854504759308, -7.36107040874295e-06]
m_Mag = sqrt(m(1)^2 + m(2)^2 + m(3)^2);
n_Mag = sqrt(n(1)^2 + n(2)^2 + n(3)^2);
Numerator =  ( m(1) * n(1) + m(2) * n(2) + m(3) * n(3) );
Denominator = m_Mag * n_Mag;
Ratio = Numerator / Denominator
theta = acos(Ratio)

The angle should be 180 degrees or pi radians, but Matlab version 2022b is giving a complex number as an output. Is there a workaround this error?

[Chuguang Pan]

@Ali. The code is running normally in R2025a.

m = [0.000128276283345364, -0.000167881251734009, 7.36224543533651e-06];
n = [-0.000128255846266013, 0.000167854504759308, -7.36107040874295e-06];
m_Mag = sqrt(m(1)^2 + m(2)^2 + m(3)^2);
n_Mag = sqrt(n(1)^2 + n(2)^2 + n(3)^2);
Numerator =  ( m(1) * n(1) + m(2) * n(2) + m(3) * n(3) );
Denominator = m_Mag * n_Mag;
Ratio = Numerator / Denominator

Ratio = -1

theta = acos(Ratio)

theta = 3.1416

[Ali]

@Chuguang Pan Thanks for checking this in R2025a. I should clarify an important detail: in your test, the vectors m and n are defined explicitly with the numerical values I posted. In my actual workflow, these vectors are not hard-coded — they are generated inside a function as the result of several intermediate calculations in a shell buckling simulation.

If I hard-code the values exactly as written, I also get the correct angle, just like in your example. The issue appears only when the vectors come from the calculation chain inside the function.

[Stephen23]

@Ali: please upload both m and n in one .mat file, by clicking the paperclip button.

[Ali]

Please see the attached attached file, and thank you so much for your help.

[Chuguang Pan]

@Ali. I think the actual value of Ratio calculated by the inside function is less than -1, which causes the complex number of acos.

[Ali]

@Chuguang Pan I have uploaded the vectors m and n generated by the function in a previous comment. Can you please use them in your calculations? I think if you use these vectors, you will get Ratio = -1 and theta = 3.14159265358979 - 2.1073424255447e-08i.

[Chuguang Pan]

@Ali. I have used the vectors included in your NormalVectors.mat file, and calculated the theta as shown in the answer utilizing dot and vecnorm functions in R2024b.

Answer 1

[Chuguang Pan]

@Ali. After loading the m and n vectors from your attached NormalVectors.mat file. I find that the calculation result of theta is correct in R2024b. The codes are listed as follows, you can try this in your R2022b.

 load NormalVectors.mat m n 
Numerator = dot(m,n); % inner product of m and n
Denominator = vecnorm(m)*vecnorm(n); 
theta = acos(Numerator/Denominator)

theta = 3.1416

[Ali]

Unfortunately, I cannot use norm or vecnorm commands. To understand my reasoning, consider a vector

A = [1 + i 1 + i]

A = 
   1.0000 + 1.0000i   1.0000 + 1.0000i

If I use norm or vecnorm, the magnitude is

vecnorm(A)

ans = 2

norm(A)

ans = 2

However, I do not want this result. My whole formulation is based on the assumption that the magnitude of this A vector will be

sqrt( (1 + i)^2 + (1 + i)^2 )

ans = 1.4142 + 1.4142i

So, sorry but this solution does not work for my case. This is why I defined magnitude of the vector m as m_Mag = sqrt(m(1)^2 + m(2)^2 + m(3)^2), instead of m_Mag = norm(m) or m_Mag = vecnorm(m).

[Chuguang Pan]

@Ali. I find that the Ratio is not equal to -1, it is less than -1. The reason is the floating error induced by sqrt function. The vector magnitudes calculated by vecnorm and sqrt are not equal.

[Chuguang Pan]

@Ali. I have tried another calculation method which first normalizing m and n vectors, I find this calculation process is correct.

load NormalVectors.mat m n
mMag = sqrt(sum(m.^2));
nMag = sqrt(sum(n.^2));
theta1 = acos(dot(m/mMag,n/nMag)); % first normalize and then make inner product

theta2 = acos(dot(m,n)/(mMag*nMag)); % original method
isequal(theta1,theta2) % theta1 and theta2 are not equal, theta1 is correct!

ans = logical
   0

[Ali]

@Chuguang Pan Your proposed solution based on first normalizing the vectors and then evaluating the angle between the vectors works perfectly. Thank you so much for your help.

[Paul]

The proposed solution is not guaranteed to work.

rng(101);
v1 = randn(1,3);
v1 = v1/norm(v1); % normalized before taking the dot product.
v2 = -v1;

By construction, the dot product of v2 and v1 should be exactly -1.

format long e
d = dot(v2,v1)

d = 
    -1.000000000000000e+00

Even though the dot product looks like -1, the display showing zeros past the decimal point indicates the numerical result is not what it seems. In fact

d < -1

ans = logical
   1

and

acos(d)

ans = 
      3.141592653589793e+00 - 2.107342425544702e-08i

We can see that d differs from -1 in the least significant bit

format hex
[d,-1]

ans = 1×2

bff0000000000001 bff0000000000000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The example above uses norm despite the objection in this comment, which I don't understand.

[Chuguang Pan]

There is an exsiting answer about how to calculate the angle between two vectors. This linked answer maybe a better solution.

https://ww2.mathworks.cn/matlabcentral/answers/101590-how-can-i-determine-the-angle-between-two-vectors-in-matlab

[Paul]

The solution on that linked page is also robust against small rounding errors in dot(v2,v1) because atan2 is not limited to arguments with absolute value less than 1.

rng(101);
v1 = randn(1,3);
v1 = v1/norm(v1); % normalized before taking the dot product.
v2 = -v1;
dot(v2,v1) < -1

ans = logical
   1

atan2(norm(cross(v2,v1)),dot(v2,v1))

ans = 3.1416

Answer 2

[David Goodmanson]

Hi Ali,

The acos and (m dot n) approach is not a good way to go about this at all. It's inaccurate when n is almost equal to -m as you have. It seems like a good idea to normalize m and n, so assume that has been done. Then

m.n  = cos(alpha)           alpha = acos(m.n)

and alpha is very close to -pi. Let

alpha = -pi + theta 

where theta is small, and is the amount by which alpha falls short of -pi.

Define q = -n so that m and q are amost identical. Then a much better method, which treats small differences in linear fashion, shows that

d = m-q
theta = |d|

In your case

theta = 9.7767e-09.             % which is tiny
alpha = -pi + 9.7767e-09

Details:

First of all, consider

m.n  = cos(alpha)           alpha = acos(m.n),

Since m.n is almost exactly -1, alpha is almost exactly - pi. But if you take a look at acos, acos(arg) is a rapidly varying function of its argument when that argument is close to -1 or 1. And at -1 or 1, acosh has infinite slope. So a small change in arg leads to a resulting angle that can easily be off.

Now it's easy to show that

m.(-n) = cos(-pi + alpha)

Let

-pi + alpha = theta

where theta is small. Also so define a new vector q = -n, so q and m are almost identical. Then

m.q  = cos(theta)         theta = acos(m.q)        % angle between m and q

This does not help yet, because here the argument of acos is almost 1, and acos has the same problem as before. But now you have the small quantity

d = m-q

available as the third side of a triangle whose other two sides have length 1 and angle theta between them. So

|d| = 2*sin(theta/2)             theta = 2*asin(|d|/2) 

asin for small argument is basically linear and has no issues like acos does. For very small |d|,

theta = |d|