提问人:rad 提问时间:3/3/2022 更新时间:3/3/2022 访问量:111
反正切函数的数值计算
Numerical evaluation of arctangent functions
问:
我在解释反正切函数的结果时遇到了困难。这种行为对于我遇到的所有实现都是一致的,因此我在这里将自己限制在 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 的每个倍数的“步骤”。
这看起来像是某种精度误差,但奇怪的是,对于我意想不到的角度。此外,这种行为也存在于普通功能中。atan
答:
请注意,您使用的是整数
因此,对于每对 (p,q),您将有给出相同角度的对。然后,对于 (p,q) 的倍数,是floor(sqrt(p**2 + q**2)/gcd(p,q)/r)
arctan(p,q)
gcd(p,q)
1
还要注意的是,对于 的倍数和奇数倍数,我们可以预测 的偶数倍数将多于奇数的多数项。这与我们在你的情节中看到的一致。p**2+q**2
1
pi/2
2
pi/4
pi/4
pi/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}')
在这里,您可以在图中看到几个与其他点具有相同角度的点。对于 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)
如果圆圈的大小增加,并且您使用足够宽的条柱进行直方图,则不会注意到变化,否则您将在直方图中看到此模式。
评论
[p, q]
r
sqrt(p^2+q^2)
p^2+q^2
pi/4
p
q
pi/4
评论
x = (rand(1, num)-0.5)*2*siz;
x = x+rand; y = y+rand