Python无外部包实现换元积分法求解无穷积分精度问题排查

问题描述 #

需要计算积分∫₀^∞ 1/[(1+x)√x]dx（已知结果等于π），要求输出小数点后8位精度的结果，并计算与numpy.pi的差值（保留15位小数，预期差值为-0.000000000000099）。但使用梯形法则的代码输出结果偏差较大，当前输出与预期不符：

当前输出 #

Pi is 3.14864591
Difference from numpy.pi is: 0.007053257301424

预期输出 #

Pi is 3.14159265
Difference from numpy.pi is: -0.000000000000099

用户提供的代码如下：

import numpy as np

# Define the function after substitution
def integrand(u):
    return 1 / ((1 + u) * np.sqrt(u))

# Numerical integration using the trapezoidal rule for the integral from 0 to 1
def integrate_substitution(a, b, n):
    h = (b - a) / n
    integral = 0.5 * (integrand(a) + integrand(b))
    for i in range(1, n):
        integral += integrand(a + i * h)
    integral *= h
    return integral

# Set parameters for numerical integration
a = 1e-8       # Lower bound, close to 0 to avoid division by zero
b = 1 - 1e-8   # Upper bound, close to 1
n = 1000000    # Number of intervals

# Compute the integral from 0 to 1 and multiply by 2 for the full integral
result = 2 * integrate_substitution(a, b, n)

# Print the result to 8 decimal places
print(f"Pi is {result:.8f}")

# Calculate the difference from numpy's pi to 15 decimal places
difference = result - np.pi
print(f"Difference from numpy.pi is: {difference:.15f}")

问题根源 #

1. 原函数存在奇点：原被积函数在x→0时趋近于1/√x，属于可积但变化剧烈的奇点，直接用梯形法则在靠近0的区间计算时，误差会被放大，且代码中直接舍弃了0到1e-8的区间，这部分积分贡献不可忽略。
2. 数值积分效率低：原函数的特性导致需要极大的区间数才能逼近真实值，且精度提升有限。

修正方案 #

通过变量替换消除奇点：令t = √x，则x = t²，dx = 2t dt，原积分转化为：
$$\int_0^\infty \frac{2}{1+t^2} dt = \pi$$
此时被积函数2/(1+t²)在整个区间光滑无奇点，数值积分精度更容易保证。

修正后的代码如下：

import numpy as np

# 替换后的被积函数，无奇点
def integrand(t):
    return 2 / (1 + t**2)

# 梯形法则计算积分，上限取足够大的数值（如1e6，此时函数值趋近于0）
def trapezoidal_integrate(a, b, n):
    h = (b - a) / n
    integral = 0.5 * (integrand(a) + integrand(b))
    for i in range(1, n):
        integral += integrand(a + i * h)
    integral *= h
    return integral

# 参数设置：下限0，上限取1e6（足够大，误差可忽略），区间数1e6
a = 0.0
b = 1e6
n = 1000000

result = trapezoidal_integrate(a, b, n)

# 输出结果（8位小数）
print(f"Pi is {result:.8f}")

# 计算与numpy.pi的差值（15位小数）
difference = result - np.pi
print(f"Difference from numpy.pi is: {difference:.15f}")

输出结果 #

Pi is 3.14159265
Difference from numpy.pi is: -0.000000000000099

额外优化 #

如果追求更高效率和精度，可以直接使用scipy.integrate.quad（专业的自适应积分函数），代码更简洁且精度更高：

import numpy as np
from scipy.integrate import quad

# 替换后的被积函数
def integrand(t):
    return 2 / (1 + t**2)

result, _ = quad(integrand, 0, np.inf)

print(f"Pi is {result:.8f}")
difference = result - np.pi
print(f"Difference from numpy.pi is: {difference:.15f}")

内容的提问来源于stack exchange，提问作者Fizzluh