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

Question

我正在尝试编写一个 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

Answer 1

就其值(value)而言，Scipy 具有数值积分函数

例如，

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

给出结果

Scipy quad integral = (3.141592653589591, 6.200897573194197e-10)

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

也就是说，我能够让您的程序进行一些调整(尽管这只是初步尝试):

1) 按照 this paper 的建议将步长 h 设置为 0.0002(大约 1/2^12)

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

但坚持手头的问题，

2) 确保为集成设置足够的迭代次数以真正收敛，即 math.fabs(w_k * func(x_k)) < 1.0e-9 的足够迭代次数

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

例如有 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