麦克劳林级数利用欧拉的思想,使用适当的多项式来近似函数的值。多项式显然与这样的函数不同cos(x)
因为它们在某个时刻都趋向于无穷大,而cos
没有。 100 阶多项式在零的每一侧最多可以近似函数的 50 个周期。由于 50 * 2pi can't近似cos(1000)
.
为了更接近合理的解决方案,多项式的阶数必须至少为x / pi
。您可以尝试计算 300+ 阶的多项式,但由于浮点数的有限精度和阶乘的巨大性,您很可能会遇到一些主要的数值问题。
相反,使用周期性cos(x)
并将以下内容添加为函数的第一行:
x %= 2.0 * math.pi
您还需要限制多项式的阶数,以避免因阶乘太大而无法放入浮点数的问题。此外,您可以而且应该通过增加先前的结果来计算阶乘,而不是在每次迭代时从头开始。这是一个具体的例子:
import math
def cos(x, i=30):
x %= 2 * math.pi
c = 2
n = 0
f = 2
for i in range(i):
if i % 2 == 0:
n += x**c / f
else:
n -= x**c / f
c += 2
f *= c * (c - 1)
return 1 - n
>>> print(cos(5), math.cos(5))
0.28366218546322663 0.28366218546322625
>>> print(cos(1000), math.cos(1000))
0.5623790762906707 0.5623790762907029
>>> print(cos(1000, i=86))
...
OverflowError: int too large to convert to float
您可以通过注意到增量乘积来进一步摆脱数值瓶颈x**2 / (c * (c - 1))
。对于更大的范围来说,这将保持良好的界限i
比你可以用直接阶乘支持:
import math
def cos(x, i=30):
x %= 2 * math.pi
n = 0
dn = x**2 / 2
for c in range(2, 2 * i + 2, 2):
n += dn
dn *= -x**2 / ((c + 1) * (c + 2))
return 1 - n
>>> print(cos(5), math.cos(5))
0.28366218546322675 0.28366218546322625
>>> print(cos(1000), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, i=86), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, i=1000), math.cos(1000))
0.5623790762906709 0.5623790762907029
请注意,经过某个点,无论执行多少次循环,结果都不会改变。这是因为现在dn
正如欧拉的意图,收敛到零。
您可以使用此信息进一步改进您的循环。由于浮点数的精度有限(具体来说,尾数为 53 位),因此您可以在以下情况下停止迭代:|dn / n| < 2**-53
:
import math
def cos(x, conv=2**-53):
x %= 2 * math.pi
c = 2
n = 1.0
dn = -x**2 / 2.0
while abs(n / dn) > conv:
n += dn
c += 2
dn *= -x**2 / (c * (c - 1))
return n
>>> print(cos2(5), math.cos(5))
0.28366218546322675 0.28366218546322625
>>> print(cos(1000), math.cos(1000))
0.5623790762906709 0.5623790762907029
>>> print(cos(1000, 1e-6), math.cos(1000))
0.5623792855306163 0.5623790762907029
>>> print(cos2(1000, 1e-100), math.cos(1000))
0.5623790762906709 0.5623790762907029
参数conv
不仅仅是上的界限|dn/n|
。由于以下项交换符号,因此它也是结果整体精度的上限。