收敛到错误值的 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

值得一提的是，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) 按照this paper 的建议，将步长 h 设置为 0.0002（大约 1/2^12）

但请注意 - 该论文实际上建议迭代地改变步长 - 使用固定的步长，您将达到一个点，即 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

Answer 2

使用多精度库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 3

我认为部分问题可能是由于范围和步长。我修改了代码 所以你可以分别输入范围和步长并重写一些数学。它似乎给出了正确的答案。尝试例如 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 4

你必须意识到 +1 和 -1 是你的被积函数的奇异点，f(x)-->+infinity 和 x-->+1,-1。因此，您可以使用您最喜欢的求积公式远离边界点，但您必须根据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 5

在 scicomp 上查看 this answer：

在 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