Python中的Adams-Bashforth方法

Question

我正在尝试用Python语言编写Adams-Bashforth方法( https://en.wikipedia.org/wiki/Linear_multistep_method#Two-step_Adams.E2.80.93Bashforth )。实际上，我的代码比Euler方法慢，并且给出的结果最差。你能帮我升级我的代码吗？

s1 = [1]
s2 = [-0.5, 1.5]
s3 = [5./12, -4./3, 23./12]
s4 = [-3./8, 37/24, -59./24, 55./24]
s5 = [251./720, -637./360, 109./30, -1387./360, 1901./720]


def adams_bashforth(f, ta, tb, xa, n, sn=s5):
    last_n = []
    h = (tb - ta) / float(n)
    t = ta
    x = xa
    # first n steps made by Euler method
    for i in range(len(sn)):
        last_n.append(h * f(t, x))
        x += last_n[-1]
        t += h
    # Adams-Bashforth method
    for i in range(n):
        x += h * sum([last_n[i] * sn[i] for i in range(len(sn))])
        last_n = last_n[1:]
        last_n.append(f(t, x))
        t += h
    return x

def f(t, x):
    return t * sqrt(x)

print adams_bashforth(f, 0, 1, 10, 1000, s5)

Answer 1

免责声明:我不熟悉Adams-Bashforth，所以我只是在重构您提供的Python代码。可以通过改变算法本身来获得更好的结果。

通过消除程序中的列表副本(last_n = last_n[1:])和附加，您可以获得不错的加速(在我的机器上，n=1M的速度为30-40%)。以下代码等效，并以循环方式使用SN (=len(sn))元素的列表：

def adams_bashforth(f, ta, tb, xa, n, sn=s5):
    h = (tb - ta) / float(n)
    t = ta
    x = xa
    SN = len(sn)
    # first n steps made by Euler method
    last_n = []
    for i in range(SN):
        y = h * f(t, x)
        last_n.append(y)
        x += y
        t += h
    # Adams-Bashforth method
    j = 0
    for i in range(n):
        ss = 0
        for s in sn:
            ss += last_n[j] * s
            j += 1
            if j >= SN: j = 0
        x += h * ss        
        last_n[j] = f(t, x)
        t += h
        j += 1
        if j >= SN: j = 0
    return x