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

Question

我正在尝试写一个Python程序，使用Tanh-sinh积分法来计算一个值：

虽然这个程序在每种情况下都能收敛到一个合理的值，没有错误，但它并没有收敛到正确的值（对于这个特定的积分，正确的值是π），我找不到问题出在哪里。

程序不是要求一个期望的精度，而是要求想要的函数评估次数，这样可以更方便地与简单的积分方法进行收敛比较。评估次数需要是一个奇数，因为使用的近似方法是：

有没有人能建议我可能做错了什么？

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

请查看这个回答在scicomp上的内容：

在使用tanh-sinh积分法时，有很多需要注意的地方。其中一个问题是，积分函数在区间边界附近的值需要非常精确地计算，距离要小于计算机的精度限制，比如在原始示例中提到的1.0 - 1.0e-20。当计算这个点时，它会被四舍五入为1.0，而在这个点上，函数f会出现奇异性，这样就会导致一些不可预知的情况。因此，你需要先对函数进行变换，使得奇异性位于0的位置。

以1 / sqrt(1 - x**2)为例，对于左右两个奇异性，可以用1 / numpy.sqrt(-x**2 + 2*x)来表示。使用tanh_sinh（这是我做的一个项目），你就可以得到：

import numpy
import tanh_sinh

# def f(x):
#    return 1 / numpy.sqrt(1 - x ** 2)

val, error_estimate = tanh_sinh.integrate_lr(
      lambda x: 1 / numpy.sqrt(-x**2 + 2*x),  # = 1 / sqrt(1 - (x-1)**2)
      lambda x: 1 / numpy.sqrt(-x**2 + 2*x),  # = 1 / sqrt(1 - (-(x-1))**2)
      2,  # length of the interval
      1.0e-10
      )
print(val, val - numpy.pi)

3.1415926533203944 -2.693987255497632e-10

Answer 2

我觉得问题的一部分可能跟范围和步长有关。我修改了代码，这样你可以单独输入范围和步长，并且重新写了一些数学计算。现在看起来能给出正确的答案。你可以试试输入5和0.1。

一个特别的问题是计算 math.cosh(0.5 * math.pi * math.sinh(k * h))，当 k * h 变大时，math.sinh(k * h) 会快速增长，这样计算它的 math.cosh 就会变得很困难。

def func(x):
#    return 1   # very simple test function
    # 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 max value for range \n")
    print "N =", N
h = input("Please the step size \n")
print "h =", h

k = -N
k_max = N
sum = 0
count = 0
print "k ", k , " " , k_max

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

    # Compute abscissa
    v = k
    u = math.pi * 0.5 * math.sinh(v)
    x_k = math.tanh(u)
    #print u
    # Compute weight 
    numerator = 0.5 * math.pi * math.cosh(v)
    csh = math.cosh(u)
    denominator = csh*csh
    w_k =  numerator / denominator
    print k, v, x_k , w_k
    sum += w_k * func(x_k)
    count += 1
    k += h      # note changed
res = sum * h
print "Integral =", sum * h

Answer 3

你需要明白，+1和-1是你要计算的函数的特殊点，当x接近+1或-1时，函数f(x)会变得非常大（趋向于正无穷）。所以，在这些边界点之外，你可以用你喜欢的积分方法来计算，但在这些点附近，你需要根据函数f(x)在这些点附近的表现来设计一个特别的积分方法。

下面是这个方法的大致步骤：

1. 先选一个很小的数，记作epsilon<<1。

2. 把要计算的积分I分成平滑部分和特殊部分：

   • I_smooth是指在区间[-1+epsilon, 1-epsilon]内的积分。
   • I_singular是指在区间[-1, -1+epsilon]和[1-epsilon, 1]内的积分。

3. 在区间[-1+epsilon, 1-epsilon]内，使用你常用的积分规则来计算I_smooth。

4. 在特殊点附近（比如x=1）进行局部展开：

   f(x) = 1/sqrt(1-x) * (a0 + a1*(1-x) + a2*(1-x)^2 + ...)
   
        = f0(x-1) + f1(x-1) + f2(x-1) + ..
   

   这实际上是对f(x)*sqrt(1-x)在x=1处进行的泰勒展开，前面乘以1/sqrt(1-x)。不过，你需要做一些数学运算来得到这个泰勒展开，除非你有Mathematica软件或者能找到相关的表格。

5. 每一项fn(x-1) = an*(1-x)^n/sqrt(1-x)都可以精确计算（这只是一个幂函数）。设Fn为fn在1-epsilon到1之间的积分。 你可以将I_singular = F0 + F1 + F2 + ...近似到你想要的精度。

6. 最后：

   I = I_smooth + I_singular  
   

注意：为了提高准确性，epsilon不应该选得太小，因为积分在这些特殊点附近会变得不稳定，反而应该增加泰勒展开的阶数。

Answer 4

使用多精度库 mpmath：

from mpmath import *

mp.dps = 100

h = mpf(2**-12);

def weights(k):
    num = mpf(0.5)*h*pi*cosh(k*h)
    den = cosh(mpf(0.5)*pi*sinh(k*h))**2
    return (num/den)

def abscissas(k):
    return tanh(mpf(0.5)*pi*sinh(k*h))

def f(x):
    return 1/sqrt(1 - mpf(x)**2)

N = 20000

result = 0
for k in range(-N, N+1):
    result = result + weights(k)*f(abscissas(k))

print result - pi

会得到错误信息

-3.751800610920472803216259350430460844457732874052618682441090144344372471319795201134275503228835472e-45

Answer 5

说实话，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），这个建议来自于这篇论文。

但要注意，这篇论文其实建议逐步调整步长——如果步长固定，随着kh的值增大，sinh或cosh会变得太大。可能更好的办法是根据那篇论文的方法来实现。

不过回到我们的问题，

2) 确保你设置足够的迭代次数，让积分真正收敛，也就是说，迭代次数要足够多，直到math.fabs(w_k * func(x_k)) < 1.0e-9。

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

例如，使用31111次迭代，计算出的π值是3.14159256208。

这是修改后的示例代码（注意我把sum替换成了thesum，因为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