我已经尝试过由 Newman (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方法有什么优势？有什么方法可以使其更快，更准确？
最佳答案
您的代码正在完善，以满足每个子间隔中的容错能力。它还使用了低阶积分规则。这两方面的改进可以显着减少功能评估的次数。
与其单独考虑每个子间隔中的误差，不如使用更高级的代码来计算所有子间隔中的总误差并进行细化，直到总误差低于所需的阈值。根据子间隔对总误差的贡献，选择子间隔进行细化，首先对较大的误差进行细化。通常，优先级队列用于快速选择子间隔进行细化。
高阶集成规则可以精确地集成更复杂的功能。例如，您的代码基于Simpson规则，该规则对于度为3的多项式都是精确的。更高级的代码可能会使用对度更高的多项式(例如10-15)精确的规则。
从实践的角度来看，最简单的方法是使用实现上述想法的固定例程，例如scipy.integrate.quad。除非您对要集成的内容有特别的了解，否则您不可能做得更好。
Romberg集成需要在等距的点进行评估。如果您可以随时评估函数，则对于“平滑”(类似于多项式)函数，其他方法通常更准确。而且，如果您的函数在任何地方都不平滑，那么自适应代码会做得更好，因为它可以专注于消除非平滑区域中的错误。