收敛到错误值的tanh-sinh正交数值积分方法 [英] Tanh-sinh quadrature numerical integration method converging to wrong value

问题描述

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

I'm trying to write a Python program to use Tanh-sinh quadrature to compute the value of:

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

but although the program converges to a sensible value with no errors in every case, it's not converging to the correct value (which is pi for this particular integral) and I can't find the problem.

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

Instead of asking for a desired level of accuracy, the program asks for the number of function evaluations wanted, to make comparisons of convergence with simpler integration methods easier. The number of evaluations needs to be an odd number as the approximation used is

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

Can anyone suggest what I might have done wrong?

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 具有数值积分功能

For what it's worth, Scipy has numerical integration functions

例如

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)

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

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

where the second number in the tuple is the absolute error.

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

That said, I was able to get your program to work with some tuning (although this is just an initial attempt):

1) 按照 这篇论文

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

But note - the paper actually suggests altering the step size iteratively - with a fixed step size you will reach a point where the sinh or cosh grow too large for large enough values of kh. It would probably be better to attempt an implementation based on that paper's approach.

但是坚持手头的问题，

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

2) Make sure you set enough iterations for the integration to really converge , i.e. enough iterations that math.fabs(w_k * func(x_k)) < 1.0e-9

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

With these tunings, I was able to get the integration to converge to the correct value of pi to 4 significant figures using > 30000 iterations.

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

For instance with 31111 iterations, the value of pi computed was 3.14159256208

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

Example code with modifications (note I replaced sum with thesum, sum is the name of a Python built-in function):

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

这篇关于收敛到错误值的tanh-sinh正交数值积分方法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！