Tanh-sinh 求积数值积分方法收敛到错误值

我正在尝试编写一个 Python 程序来使用 Tanh-sinh 求积来计算以下值：

但是，尽管程序收敛到一个合理的值，在每种情况下都没有错误，但它没有收敛到正确的值（对于这个特定的积分来说是 pi ），我找不到问题。

该程序不要求所需的准确度，而是要求所需的函数评估数量，以便更容易地比较与更简单的积分方法的收敛性。评估次数必须是奇数，因为使用的近似值为

谁能建议我可能做错了什么？

import math

def func(x):
    # Function to be integrated, with singular points set = 0
    if x == 1 or x == -1 :
        return 0
    else:
        return 1 / math.sqrt(1 - x ** 2)

# Input number of evaluations
N = input("Please enter number of evaluations \n")
if N % 2 == 0:
    print "The number of evaluations must be odd"
else:
    print "N =", N  

# Set step size
h = 2.0 / (N - 1)
print "h =", h

# k ranges from -(N-1)/2 to +(N-1)/2
k = -1 * ((N - 1) / 2.0)
k_max  = ((N - 1) / 2.0)
sum = 0

# Loop across integration interval
while k < k_max + 1:

    # Compute abscissa
    x_k = math.tanh(math.pi * 0.5 * math.sinh(k * h))

    # Compute weight
    numerator = 0.5 * h * math.pi * math.cosh(k * h)
    denominator = math.pow(math.cosh(0.5 * math.pi * math.sinh(k * h)),2)
    w_k =  numerator / denominator

    sum += w_k * func(x_k)

    k += 1

print "Integral =", sum

就其价值而言，Scipy 具有数值积分函数

例如，

from scipy import integrate
check = integrate.quad(lambda x: 1 / math.sqrt(1 - x ** 2), -1, 1)
print 'Scipy quad integral = ', check

给出结果

Scipy四积分 = (3.141592653589591, 6.200897573194197e-10)

其中元组中的第二个数字是绝对误差。

也就是说，我能够通过一些调整让您的程序正常工作（尽管这只是初步尝试）：

1) 按照建议将步长 h 设置为 0.0002（大约 1/2^12）这张纸 http://www.davidhbailey.com/dhbpapers/dhb-tanh-sinh.pdf

但请注意 - 论文实际上建议迭代地改变步长 - 使用固定的步长，您将达到一个点，即 sinh 或 cosh 对于足够大的值来说变得太大kh。尝试基于该论文的方法进行实现可能会更好。

但坚持手头的问题，

2) 确保设置足够的迭代以使积分真正收敛，即足够的迭代 math.fabs(w_k * func(x_k))

通过这些调整，我能够使用 > 30000 次迭代将积分收敛到 pi 的正确值（4 位有效数字）。

例如，迭代次数为 31111 次，计算出的 pi 值为 3.14159256208

修改后的示例代码（注意我用 thesum 替换了 sum，sum 是 Python 内置函数的名称）：

import math

def func(x):
    # Function to be integrated, with singular points set = 0
    if x == 1 or x == -1 :
        return 0
    else:
        return 1 / math.sqrt(1 - x ** 2)

# Input number of evaluations
N = input("Please enter number of evaluations \n")
if N % 2 == 0:
    print "The number of evaluations must be odd"
else:
    print "N =", N  

# Set step size
#h = 2.0 / (N - 1)
h=0.0002 #(1/2^12)
print "h =", h

# k ranges from -(N-1)/2 to +(N-1)/2
k = -1 * ((N - 1) / 2.0)
k_max  = ((N - 1) / 2.0)
thesum = 0

# Loop across integration interval
actual_iter =0
while k < k_max + 1:

    # Compute abscissa
    x_k = math.tanh(math.pi * 0.5 * math.sinh(k * h))

    # Compute weight
    numerator = 0.5 * h * math.pi * math.cosh(k * h)
    dcosh  = math.cosh(0.5 * math.pi * math.sinh(k * h))
    denominator = dcosh*dcosh
    #denominator = math.pow(math.cosh(0.5 * math.pi * math.sinh(k * h)),2)
    w_k =  numerator / denominator

    thesum += w_k * func(x_k)
    myepsilon = math.fabs(w_k * func(x_k))
    if actual_iter%2000 ==0 and actual_iter > k_max/2:
        print "Iteration = %d , myepsilon = %g"%(actual_iter,myepsilon)


    k += 1
    actual_iter += 1

print 'Actual iterations = ',actual_iter
print "Integral =", thesum