如何高精度积分正弦余弦乘积函数？Scipy积分精度问题求助

我尝试对正弦（sin）和余弦（cos）类函数的乘积进行积分。当结果量级较大（如1e-2_{1e-4）时，结果匹配度良好，但量级极小（如1e-11}1e-32）的结果存在较大相对误差，有时甚至符号相反。

使用的代码如下：

import numpy as np
import scipy.integrate as spi

# Define the integrand function
def f0(x,n):
    if n == 1:
        return (1-x)
    elif n == 2:
        return x
    else:
        return np.sin((n - 2) * np.pi * x)


def f00(x, m, n):
    return f0(x,m)*f0(x,n)

# Define the dimensions of the tensor
NMHT_0 = 12

# Create the tensor
I00 = np.zeros((NMHT_0, NMHT_0))

# Perform the integration
for m in range(NMHT_0):
    for n in range(NMHT_0):
        result1, error1 = spi.quad(f00, 0, 1, args=(m, n))
        
        
        I00[m, n] = result1      

我最初使用scipy.integrate.quad_vec函数进行积分，但所有小量级结果均不匹配；替换为scipy.integrate.quad函数后，仅部分小量级结果精度改善，多数仍不匹配。请问有什么方法可以提升积分精度，或是需要尝试其他库？

1. 利用正交函数解析解彻底避免数值误差 #

你的积分函数属于正交函数系（线性函数+正弦函数）的乘积，大部分情况可以直接用解析公式计算，完全消除数值误差：

• 当m=1且n=1：$\int_0^1 (1-x)^2 dx = 1/3$
• 当m=1且n=2：$\int_0^1 (1-x)x dx = 1/6$
• 当m=2且n=2：$\int_0^1 x^2 dx = 1/3$
• 当m≥3且n≥3：若m≠n，$\int_0^1 \sin((m-2)\pi x)\sin((n-2)\pi x)dx = 0$（正交性）；若m=n，结果为$1/2$
• 当m=1且n≥3：$\int_0^1 (1-x)\sin((n-2)\pi x)dx = 1/[(n-2)\pi]^2$
• 当m=2且n≥3：$\int_0^1 x\sin((n-2)\pi x)dx = 1/[(n-2)\pi]^2$

直接用这些公式填充I00，精度能达到机器精度，示例代码：

import numpy as np

NMHT_0 = 12
I00 = np.zeros((NMHT_0, NMHT_0))

# 填充解析解
for m in range(NMHT_0):
    for n in range(NMHT_0):
        if m == 1 and n == 1:
            I00[m,n] = 1/3
        elif m ==1 and n ==2:
            I00[m,n] = 1/6
        elif m ==2 and n ==2:
            I00[m,n] =1/3
        elif m >=3 and n >=3:
            if m ==n:
                I00[m,n] = 1/2
            else:
                I00[m,n] =0
        elif m ==1 and n >=3:
            k = n-2
            I00[m,n] = 1/(k*np.pi)**2
        elif m ==2 and n >=3:
            k =n-2
            I00[m,n] =1/(k*np.pi)**2
        # 对称情况直接复用结果
        elif n ==1 and m >=3:
            k =m-2
            I00[m,n] =1/(k*np.pi)**2
        elif n ==2 and m >=3:
            k =m-2
            I00[m,n] =1/(k*np.pi)**2

2. 调整scipy.quad的精度参数 #

如果一定要用数值积分，可以给quad传入更高精度的控制参数，强制提升计算精度：

• 设置epsabs和epsrel为更小的阈值（比如1e-14），要求积分结果的绝对误差和相对误差都低于该值
• 增加limit参数，允许quad使用更多的采样点迭代计算

修改后的积分代码：

result1, error1 = spi.quad(f00, 0, 1, args=(m, n), epsabs=1e-14, epsrel=1e-14, limit=100)

3. 尝试其他数值积分方法 #

高斯-勒让德求积法（fixed_quad） #

对于光滑函数（你的积分函数都是光滑的），高斯求积法的精度远高于普通自适应积分，指定足够多的节点数即可：

result1, _ = spi.fixed_quad(f00, 0, 1, args=(m, n), n=50)

高密度采样的梯形法/辛普森法 #

用numpy生成足够密的采样点，再用梯形法或辛普森法积分，虽然是低阶方法，但采样点足够多时精度也能满足要求：

x = np.linspace(0, 1, 100000)  # 生成10万个采样点
y = f00(x, m, n)
result1 = np.trapz(y, x)  # 梯形法积分
# 或者用辛普森法：from scipy.integrate import simpson; result1 = simpson(y, x)

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