在 Python 3.7 中评估积分；奇怪的行为答案

【问题标题】：Evaluating integrals in Python 3.7; strange behaviour在 Python 3.7 中评估积分；奇怪的行为
【发布时间】：2018-12-09 16:06:51
【问题描述】：

我编写了一些代码来近似积分：使用 Python 3.7，但是发生了一些奇怪的行为，这给了我错误的结果。我通过以下方式得出公式：

让然后和对于 n = 1,2,3...

这在我的代码中实现：

import numpy as np

I = 1
for n in range(1,21):
    I = 2*(np.log(2))**n - n*I

这应该导致 I = 0.0000419426270488826，但我的代码给了我 50.40429353428721。我一直试图弄清楚使用打印语句发生了什么：

print("Iteration: ",n)
print("first half: ",2*(np.log(2))**n)
print("second half: ", n*I)
print("New I: ",I, "\n")

您可以看到等式的后半部分在第 17 次迭代中变为负数，但我不明白为什么，因为 I 和 n 都应该是正数。我的猜测是，这就是问题开始出现的地方。有谁知道为什么结果不正确，我的假设是否正确？我正在使用 Mac OS X el Capitan 10.11.6

【问题讨论】：

看起来n*I 在n >= 17 之后大于2*(np.log(2))**n。这将使I 成为负数。
不应该这样。在计算第 17 次迭代中的新 I 时，应使用第 16 次迭代中的 I。第 16 次迭代中的 I 为 0.0006564050553949539。这乘以 17 应该是一个正数。
来自您的代码I = 2*(np.log(2))**n - n*I。即使n*I是一个正整数，如果它大于2*(np.log(2))**n，也会产生一个负数。
我知道一个正数减去一个更大的正数会给我一个负数。但这并不能解释为什么在 python 中一个正数乘以一个正数会给我一个负数。
嗯...您在哪里看到正数乘以正数等于负数？也许这个电子表格会有所帮助：i.stack.imgur.com/y3Xha.png

标签： python-3.x for-loop math integral

【解决方案1】：

我认为这是因为溢出错误。如果您使用Decimal 来提高浮点精度，它将解决问题。代码如下：

import numpy as np
from decimal import *
getcontext()

n = 20
print('n=',n)
def In(n):
    if n==0:
        return 1
    else:
        return 2*Decimal(2).ln()**n - n*In(n-1)

print('I{}={}'.format(n,float(In(n))))

# check with an existing function of Scipy
import scipy.integrate as integrate
result = integrate.quad(lambda x: np.log(x)**n, 1, 2)
print('I{}={}'.format(n,result))

这是输出：

n= 20
I20=4.1942535120159466e-05
I20=(4.1942627048882645e-05, 7.29235597161084e-18)

顺便说一句，公式需要通过将ln(x)替换为ln(2)来更正。

更新：问题的名称可能是loss of significance

【讨论】：

这似乎可以解决问题，谢谢！我有一个小问题：为什么这称为溢出错误？当数字太大而无法以当前精度表示时，是否会溢出？我认为问题可能是失去了意义，因为两个小数字被以有限的精度减去。越来越多的迭代会导致 Python 中的值越来越多的差异。
是的，我认为你是对的，我已经更新了答案。显然，当值越来越小时，减去两个小数可以删除有效数字。如果这样重复很多次，就会导致结果出现更多的差异。