反正切函数的数值计算

Numerical evaluation of arctangent functions

提问人:rad 提问时间:3/3/2022 更新时间:3/3/2022 访问量:111

问:

我在解释反正切函数的结果时遇到了困难。这种行为对于我遇到的所有实现都是一致的,因此我在这里将自己限制在 NumPy 和 MATLAB 上。

这个想法是有随机放置的点的圆圈。目标是表示它们在极坐标系中的位置,并且由于它们是均匀分布的,因此我希望 θ 角(使用函数计算)也沿区间 -π 随机分布......π.atan2

以下是 MATLAB 的代码:

stp = 2*pi/2^8;
siz = 100;
num = 100000000;

x = randi([-siz, siz], [1, num]);
y = randi([-siz, siz], [1, num]);
m = (x.^2+y.^2) < siz^2;

[t, ~] = cart2pol(x(m), y(m));
figure()
histogram(t, -pi:stp:pi);

这里是 Python 和 NumPy:

import numpy as np
import matplotlib.pyplot as pl

siz = 100
num = 100000000

rng = np.random.default_rng()
x = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
y = rng.integers(low=-siz, high=siz, size=num, endpoint=True)
m = (x**2+y**2) < siz**2

t = np.arctan2(y[m], x[m]);
pl.hist(t, range=[-np.pi, np.pi], bins=2**8)
pl.show()

在这两种情况下,我都得到了这样的结果,人们可以很容易地看到 π/4 的每个倍数的“步骤”。

Atan2 results

这看起来像是某种精度误差,但奇怪的是,对于我意想不到的角度。此外,这种行为也存在于普通功能中。atan

numpy matlab precision atan2 atan

评论

0赞 Cris Luengo 3/3/2022
正如 Bob 在下面所说,这是一个离散化问题。您可以使用 创建非离散样本,也可以将域的来源设置在随机位置:。x = (rand(1, num)-0.5)*2*siz;x = x+rand; y = y+rand

答:

3赞 Bob 3/3/2022 #1

请注意,您使用的是整数

因此,对于每对 (p,q),您将有给出相同角度的对。然后,对于 (p,q) 的倍数,是floor(sqrt(p**2 + q**2)/gcd(p,q)/r)arctan(p,q)gcd(p,q)1

还要注意的是,对于 的倍数和奇数倍数,我们可以预测 的偶数倍数将多于奇数的多数项。这与我们在你的情节中看到的一致。p**2+q**21pi/22pi/4pi/4pi/4

让我们用整数坐标绘制位于半径为 10 的圆内的点。

import numpy as np
import matplotlib.pyplot as plt
from collections import Counter
def gcd(a,b):
    if a == 0 or b == 0:
        return max(a,b)
    while b != 0:
        a,b = b, a%b
    return a;
R = 10
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
t = np.linspace(0, np.pi / 2)
plt.figure(figsize=(6, 6))
plt.plot(x, y, 'o')
plt.plot(R * np.cos(t), R * np.sin(t))

lines = Counter((xi / gcd(xi,yi), 
                 yi / gcd(xi,yi)) for xi, yi in zip(x,y))
plt.axis('off')
for (x,y),f in lines.items():
    if f != 1:
        r = np.sqrt(x**2 + y**2)
        plt.plot([0, R*x/r], [0, R*y/r], alpha=0.25)
        plt.text(R*1.03*x/r, R*1.03*y/r, f'{int(y)}/{int(x)}: {f}')

enter image description here

在这里,您可以在图中看到几个与其他点具有相同角度的点。对于 45 度有 7 个点,对于 90 的倍数有 10 个点。许多点都有独特的角度。 基本上,你有很多角度,只有很少的点,还有一些角度击中了很多点。

但总体而言,这些点相对于角度分布几乎均匀。在这里,我绘制了几乎是一条直线的累积频率(如果分布是单向分布的,则会是什么样子),并且 bin 频率形成一些三角形分形模式。

R = 20
x,y = np.indices((R+1, R+1))
m = (x**2 + y**2) <= R**2
x,y = x[m], y[m]
plt.figure(figsize=(6,6))
plt.subplot(211)
plt.plot(np.sort(np.arctan2(x,y))*180/np.pi, np.arange(len(x)), '.', markersize=1)

plt.subplot(212)
plt.plot(np.arctan2(x,y)*180/np.pi, np.gcd(x,y), '.', markersize=4)

enter image description here

如果圆圈的大小增加,并且您使用足够宽的条柱进行直方图,则不会注意到变化,否则您将在直方图中看到此模式。

enter image description here

评论

0赞 rad 3/3/2022
看来你是对的,谢谢!至少如果我尝试浮动,它会给我预期的结果。但是,我很难理解你的推理。你会这么善良并扩展它吗?该货币对的定义间隔是多少(从第二段开始,它看起来必须位于正方形边缘 [-1..1;-1..1])?不等于 ?对于奇数和偶数倍,如何给出不同的结果,因为 ang 的平方在所有情况下都必须相等 (for)?谢谢。[p, q]rsqrt(p^2+q^2)p^2+q^2pi/4pqpi/4