如何改进该自适应梯形规则？

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积分进行比较，并发现对于相同的精度，这种自适应方法需要使用更多的点来计算积分。

这是否仅仅是由于其固有的局限性？使用此自适应算法的优势是什么？有没有让它更快更准确的方法？

- Physicist

---

不清楚该算法的具体功能，但我想使用numpy或cython加速可能会大大提高性能。你怎么调用这个函数？ - Padraic Cunningham

使用中点求积法，而不是梯形求积法。与一开始可能会预期的相反，从数值上讲，中点求积法是更好的选择。 - Nico Schlömer

Answer 1

您的代码正在精细调整，以满足每个子区间中的误差容限。它还使用低阶积分规则。在这两方面的改进可以显著减少函数计算的数量。

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

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

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

Romberg积分需要在等间距点处进行评估。如果您可以在任意点评估函数，则其他方法通常对于“平滑”（类似于多项式）函数更准确。如果您的函数并非在所有地方都是平滑的，则自适应代码将做得更好，因为它可以专注于降低非平滑区域中的误差。

- bg2b