如何改善这种自适应梯形法则？

Question

我已经尝试过由纽曼（Newman）编写的计算物理学进行练习，并为自适应梯形规则编写了以下代码。当每张幻灯片的误差估计值大于允许值时，它将把该部分分成两半。我只是想知道我还能做些什么来使算法更有效。

xm=[]
def trap_adapt(f,a,b,epsilon=1.0e-8):    
    def step(x1,x2,f1,f2):
        xm = (x1+x2)/2.0
        fm = f(xm)
        h1 = x2-x1
        h2 = h1/2.0
        I1 = (f1+f2)*h1/2.0
        I2 = (f1+2*fm+f2)*h2/2.0
        error = abs((I2-I1)/3.0) # leading term in the error expression
        if error <= h2*delta:
            points.append(xm) # add the points to the list to check if it is really using more points for more rapid-varying regions
            return h2/3*(f1 + 4*fm + f2)
        else:
            return step(x1,xm,f1,fm)+step(xm,x2,fm,f2) 
    delta = epsilon/(b-a)
    fa, fb = f(a), f(b)  
    return step(a,b,fa,fb)

此外，我使用了一些简单的公式将其与Romberg积分进行比较，发现对于相同的精度，该自适应方法使用更多的点来计算积分。

仅仅是因为其固有的局限性吗？使用这种自适应算法代替Romberg方法有什么优势？有什么方法可以使其更快，更准确？

Answer 1

您的代码正在完善，以满足每个单独子区间的容错能力。它还使用低阶积分规则。这两者的改进可以显着减少功能评估的数量。

更高级的代码不是单独考虑每个子区间中的误差，而是计算所有子区间的总误差并进行细化，直到总误差低于所需阈值。根据子区间对总误差的贡献来选择细化，首先细化较大的误差。通常，优先级队列用于快速选择要细化的子区间。

高阶积分规则可以精确地积分更复杂的函数。例如，您的代码基于辛普森规则，该规则对于次数最多为 3 的多项式来说是精确的。更高级的代码可能会使用对于次数更高的多项式（例如 10-15）精确的规则。

从实用的角度来看，最简单的事情是使用实现上述想法的固定例程，例如 scipy.integrate.quad。除非您对想要集成的内容有特定的了解，否则您不太可能做得更好。

Romberg 积分需要在等距点进行评估。如果您可以在任意点评估函数，那么对于“平滑”（类似多项式）函数，其他方法通常更准确。如果您的函数并非处处平滑，那么自适应代码会做得更好，因为它可以专注于消除非平滑区域中的错误。