如何改善这种自适应梯形法则? [英] how to improve this adaptive trapezoidal rule?

问题描述

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

I have attempted an exercise from the computational physics written by Newman and written the following code for an adaptive trapezoidal rule. When the error estimate of each slide is larger than the permitted value, it divides that portion into two halves. I am just wondering what else I can do to make the algorithm more efficient.

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

Besides, I used a few simple formula to compare this to Romberg integration, and found that for the same accuracy, this adaptive method uses many more point to calculate the integral.

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

Is it just because of its inherent limitations? Any advantages of using this adaptive algorithm instead of the Romberg method? any ways to make it faster and more accurate?

推荐答案

您的代码正在完善，以满足每个子间隔中的容错能力.它还使用了低阶积分规则.这两个方面的改进可以显着减少功能评估的次数.

Your code is refining to meet an error tolerance in each individual subinterval. It's also using a low-order integration rule. Improvements in both of these can significantly reduce the number of function evaluations.

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

Rather than considering the error in each subinterval separately, more advanced codes compute the total error over all the subintervals and refine until the total error is below the desired threshold. Subintervals are chosen for refinement according to their contribution to the total error, with larger errors being refined first. Typically a priority queue is used to quickly chose the subinterval for refinement.

高阶集成规则可以准确地集成更复杂的功能.例如，您的代码基于Simpson规则，该规则对于度为3的多项式都是精确的.更高级的代码可能会使用对度为高的多项式(例如10-15)精确的规则.

Higher-order integration rules can integrate more complicated functions exactly. For example, your code is based on Simpson's rule, which is exact for polynomials of degree up to 3. A more advanced code will probably use a rule that's exact for polynomials of much higher degree (say 10-15).

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

From a practical point of view, the simplest thing is to use a canned routine that implements the above ideas, e.g., scipy.integrate.quad. Unless you have particular knowledge of what you want to integrate, you're unlikely to do better.

Romberg集成需要在等距的点进行评估.如果您可以随时评估函数，则对于平滑"(类似于多项式)函数，其他方法通常更准确.而且，如果您的函数在任何地方都不平滑，那么自适应代码会做得更好，因为它可以专注于消除非平滑区域中的错误.

Romberg integration requires evaluation at equally-spaced points. If you can evaluate the function at any point, then other methods are generally more accurate for "smooth" (polynomial-like) functions. And if your function is not smooth everywhere, then an adaptive code will do much better because it can focus on beating down the error in the non-smooth regions.

这篇关于如何改善这种自适应梯形法则?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！