我尝试了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积分进行比较,并发现对于相同的精度,这种自适应方法需要使用更多的点来计算积分。
这是否仅仅是由于其固有的局限性?使用此自适应算法的优势是什么?有没有让它更快更准确的方法?