问题描述
我使用 Taylors series 来计算一个数的 cos,对于小数,函数返回准确的结果,例如 cos(5)
给出 0.28366218546322663
。但是对于较大的数字,它会返回不准确的结果,例如 cos(1000)
给出 1.2194074101485173e+225
def factorial(n):
c = n
for i in range(n-1,-1):
c *= i
return c
def cos(x,i=100):
c = 2
n = 0
for i in range(i):
if i % 2 == 0:
n += ((x**c) / factorial(c))
else:
n -= ((x**c) / factorial(c))
c += 2
return 1 - n
我尝试使用 round(cos(1000),8)
把它仍然返回一个用科学记数法写的数字 1.2194074101485173e+225
和 e+ 部分。 math.cos(1000)
给出 0.5623790762907029
,我怎样才能四舍五入我的数字,使它们与 math.cos 方法相同?
解决方法
麦克劳林级数使用欧拉的思想通过适当的多项式来近似函数的值。多项式显然与 cos(x)
之类的函数不同,因为它们都在某个时刻趋于无穷大,而 cos
则不然。一个 100 阶多项式最多可以在零的每一侧近似函数的 50 个周期。由于 50 * 2pi 不能近似 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),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),1e-6),math.cos(1000))
0.5623792855306163 0.5623790762907029
>>> print(cos2(1000,1e-100),math.cos(1000))
0.5623790762906709 0.5623790762907029
参数 conv
不仅仅是 |dn/n|
的边界。由于以下项切换符号,因此它也是结果整体精度的上限。
返回的数字只是一个数字;在您打印它之前,它没有符号感。如果您希望控制值的打印,假设您像这样打印
print(cos(1000))
然后我们可以使用format strings来控制输出
print("{:f}".format(cos(1000)))
如果您使用的是 Python 3.6 或更新版本,我们甚至可以将它直接 interpolate 到字符串文字中。
print(f"{cos(1000):f}")
您可以阅读以上链接以查看有关格式迷你语言的更多详细信息(两个功能之间的语言相同)。例如,如果您想打印特定数量的小数位,您也可以提出要求。我们可以精确打印三位小数如下
print("{:.3f}".format(cos(1000)))
print(f"{cos(1000):.3f}")
但是,正如 Mad Physicist 所指出的,您的代码也存在更多数学问题,因此我强烈建议您也阅读他的回答。