如何在Python中实现自适应海恩法？

Question

我正在尝试实现一个叫做海恩法的函数。不过，我想让它具有自适应的特点。一般来说，比如说矩形法，如果你想用自适应的方式，就需要比较从a到b的区域面积，看看这个区域的面积是否和从a到a和b的中间点，再加上从这个中间点到b的面积相等。如果它们的差别在一个可以接受的范围内，就返回这个面积的总和；如果不相等，就要在这两个区间上继续进行计算。

这是我目前写的代码：

import math    

def k(x,y):
    return math.sin(x+y)

''' h is the step size '''
def Heun(f,x,y,h,end):
    while x < end:
        f0=f(x,y)
        z=y+h*f0
        x+=h
        f1=f(x,z)
        y+=h*0.5*(f0+f1)
    return y 

def AdaptDiff(diffF,f,x,y,h,end,tol):
    if abs(1.-x/end) < tol:
        return y
    y1=diffF(f,x,y,1,x+h)
    y_=diffF(f,x,y,1,x+h/2.)
    y2=diffF(f,x+h/2.,y_,1,x+h)
    if abs(1.-y1/y2) > tol:
        return AdaptDiff(diffF,f,x+h/2.,y2,h/2.,end,tol)
    return AdaptDiff(diffF,f,x+h,y1,h*2.,end,tol)

print AdaptDiff(Heun,k,0,1,1/2.,1,1e-10)

但是我遇到了一个错误，提示说最大递归深度超出了限制。有没有人知道我哪里出错了？我应该得到的结果是1.73513792333。

谢谢。

Answer 1

Python在递归编程方面不是特别好，因为它默认的递归深度有限，而且内存也有限。Python为了保留调用栈的信息，不支持尾调用优化。

你可以通过sys模块来查看你的递归限制，我的限制是1000：

>>> import sys
>>> sys.getrecursionlimit()
1000

你也可以设置自己的递归限制，但这样做不太建议，因为可能会导致内存耗尽，从而让你的系统崩溃：

sys.setrecursionlimit(limit)

你应该考虑把你的函数改成循环的方式来实现。