Erstaunlich großer Unterschied bei der Bewertung der trigonometrischen Identität mit NumPy

8

Laut Wolfram Alpha und dem Sage-Computeralgebrasystem gilt folgende Identität:

cos(arctan(l1l2d))=11+(l1l2)2d2

Als ich jedoch versuchte, es mit einem beliebigen Beispiel in NumPy zu verifizieren, bemerkte ich einen ziemlich großen Unterschied in den tatsächlichen Werten, die von beiden Seiten der Identität berechnet wurden. Ich habe den folgenden Code verwendet:

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

Die Auswertung der Ergebnisse für leftund rightder folgenden ergab:

    left = 0.99315060432287616
    right = 1.0

Es ist verlockend, dies einfach als numerischen Fehler abzuschreiben, aber da ich nur sehr wenig Erfahrung damit habe, wie groß numerische Fehler werden können, bin ich mir nicht so sicher. Ist das möglich oder fehlt mir etwas (offensichtliches)?

Daniel Eberts
quelle
2
Ihr rightist falsch eingegeben. es sollte sein right = 1/sqrt() Wenn ich die Formeln in meinen Ti-89 eingebe, bekomme ich eine Übereinstimmung mit 12 Stellen bei 0,99315 ...
Godric Seer
Ja, ich denke, er hätte das bemerkt, wenn der Quadratwurzelausdruck richtig ausgewertet worden wäre.
Michael Grant
Ja, als ich seinen Gesichtsausdruck betrat, bekam ich 1.007, was definitiv sichtbar gewesen wäre.
Godric Seer
1
Ja wirklich. Das ist peinlich. Danke, dass du es bemerkst.
Daniel Eberts
Habe diesen Fehler vor mir gemacht.
Michael Grant

Antworten:

10

Ich habe den Verdacht, dass Ihr Python-Ausdruck für die rechte Seite eine Ganzzahldivision und keine Gleitkommadivision durchführt. Infolgedessen ((l1-l2)**2)/(d**2)wird als Null ausgewertet, und der Term innerhalb der Quadratwurzel ist Eins.

Tatsächlich haben Sie das Gegenteil auch in Ihrem rechten Ausdruck vergessen, aber das ist nicht das erste 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
Michael Grant
quelle
6
Ich habe from __future__ import divisionmir angewöhnt, zu diesem Zweck alle meine Python-Skripte ganz oben zu verwenden.
Geoffrey Irving
@GeoffreyIrving, verwandle das in eine Antwort und ich wette, du bekommst einen ernsthaften Repräsentanten dafür. Guter Rat.
Michael Grant
7

Wie Michael C. Grant betonte, besteht das Problem darin, dass die Divisionsoperation eine ganzzahlige Division und keine Gleitkommadivision ist. In Versionen von Python vor 3 werden zwei Ganzzahlen durch /Runden auf eine Ganzzahl geteilt. Dieses Verhalten kann durch Hinzufügen geändert werden

from __future__ import division

oben in Ihrem Skript. Beachten Sie, dass diese Deklaration das Verhalten /nur für die aktuelle Datei ändert . Wenn Sie gelegentlich das ursprüngliche Verhalten wünschen, verwenden Sie //(doppelter Schrägstrich). Dies wird hier ausführlicher besprochen:

PEP 238 - Ändern des Abteilungsbetreibers

Beachten Sie, dass die meisten Sprachen, einschließlich C, eine ähnliche Semantik für die Division haben, die Semantik in Python jedoch verwirrender ist, da Ganzzahlen normalerweise nicht zu Floats heraufgestuft werden, wenn sie an Funktionen übergeben werden. Somit kann eine Funktion, die so geschrieben ist, als ob sie auf Floats arbeitet, tatsächlich auf ganzen Zahlen arbeiten.

Geoffrey Irving
quelle