Astoundingly large difference when evaulating trigonometric identity with NumPy

Question

According to Wolfram Alpha and the Sage computer algebra system, the following identity holds: $$ \cos\left(\arctan\left(\frac{l_1-l_2}{d}\right)\right) = \frac{1}{\sqrt{1 + \frac{(l_1-l_2)^2}{d^2}}} $$

However, when I tried to verify it with an arbitrary example in NumPy, I noticed a rather large difference in the actual values computed by both sides of the idenitity. I have used the following code:

    l1 = 10; l2 = 8; d = 17
    from numpy import arctan2, cos, sin, sqrt

    alpha = arctan2((l1-l2),d)
    left = cos(alpha)

    right = sqrt(1 + ((l1-l2)**2)/(d**2))

Evaluating the results for left and right yielded the following:

    left = 0.99315060432287616
    right = 1.0

It is tempting to write this off as simply being a numerical error, but since I have very little experience in how large numerical errors can get, I am not so sure. Is this possible or am I missing something (obvious)?

your right is incorrectly entered. it should be right = 1/sqrt() When I enter the formulas into my Ti-89 I get a match out to 12 digits at 0.99315... — Godric Seer, Apr 26 '13 at 15:39
Yeah, I think he would have caught that if the square root expression had been evaluated correctly. — Michael Grant, Apr 26 '13 at 15:40
yeah, when I entered his expression I got 1.007, which definitely would have been visible. — Godric Seer, Apr 26 '13 at 15:42

score 10 · Accepted Answer · answered Apr 26 '13 at 15:37

10

I have a suspicion that your Python expression for the right-hand side is doing integer division, not floating-point division. As a result, ((l1-l2)**2)/(d**2) is being evaluated as zero, and the term inside the square root is one.

In fact, you forgot the reciprocal in your right-hand expression as well, but that's not the first problem...

In MATLAB:

>> l1 = 10; l2 = 8; d = 17;
>> alpha = atan2((l1-l2),d)
alpha =
    0.1171
>> left = cos(alpha)
left =
    0.9932
>> right = 1/sqrt(1 + ((l1-l2)^2)/(d^2))
right =
    0.9932

answered Apr 26 '13 at 15:37

Michael Grant

2,062
11
24

6

I've gotten into the habit of using from __future__ import division at the top of all my python scripts for this purpose. – Geoffrey Irving Apr 26 '13 at 16:28
@GeoffreyIrving, turn that into an answer and I'll bet you get some serious rep for that. Good advice. – Michael Grant Apr 27 '13 at 15:52

score 7 · Answer 2 · answered Apr 27 '13 at 19:32

As Michael C. Grant pointed out, the problem is that the division operation is integer division, not floating point division. In versions of Python before 3, dividing two integers using / rounds down to an integer. This behavior can be changed by adding

from __future__ import division

at the top of your script. Note that this declaration changes the behavior of / only for the current file. If you want the original behavior occasionally, use // (double slash). This is discussed in more detail here:

PEP 238 -- Changing the Division Operator

Note that while most languages including C have similar semantics for division, the semantics are more confusing in Python because integers are not usually promoted to floats when passed to functions. Thus, a function written as if it was operating on floats may in fact be operating on integers.

Astoundingly large difference when evaulating trigonometric identity with NumPy

2 Answers2