如何优化python中的时间步进算法？

我有一段代码，它为所谓Lorenz95 model（由Ed Lorenz于1995年发明）增加了一个时间步长。它通常作为40-变量模型实现，并显示混沌行为。我已经编写了时间步的算法如下：如何优化python中的时间步进算法？

class Lorenz: 
    '''Lorenz-95 equation''' 

    global F, dt, SIZE 
    F = 8 
    dt = 0.01 
    SIZE = 40 

    def __init__(self): 
     self.x = [random.random() for i in range(SIZE)] 

    def euler(self): 
     '''Euler time stepping''' 
     newvals = [0]*SIZE 
     for i in range(SIZE-1): 
      newvals[i] = self.x[i] + dt * (self.x[i-1] * (self.x[i+1] - self.x[i-2]) - self.x[i] + F) 
     newvals[SIZE-1] = self.x[SIZE-1] + dt * (self.x[SIZE-2] * (self.x[0] - self.x[SIZE-3]) - self.x[SIZE-1] + F) 
     self.x = newvals

此功能欧拉不慢，但不幸的是，我的代码需要做出一个非常大量的调用它。有没有一种方法可以编写时间步进以使其运行更快？

非常感谢。

来源

2013-05-08 Neill Bowler

到底是在做什么'global'？ – 2013-05-08 15:44:30

你可以展开'for'循环并且明确它所做的一系列计算 - 但是没有得到解决它们每个都必须完成的事实。可能要考虑写一个C扩展来完成实际的数字处理。 – martineau 2013-05-08 16:14:02

正如我们都知道'global'的使用是懒惰编程的标志。 ;-)我想确保类“Lorenz”的所有实例的数据数组大小相同。因此，迫使SIZE成为这个课程的全球化是确保这一点的最简单方法。 – 2013-05-09 08:20:45

至少有两种可能的优化：以更智能的方式工作（算法改进）并加快工作速度。

在算法侧，您使用的Euler method，这是一阶方法（因此全局误差正比于步长大小），并且具有一个短小稳定区域。也就是说，效率不高。
另一方面，如果您使用标准的CPython实现，这种代码将会非常缓慢。为了解决这个问题，你可以试着在PyPy下运行它。它的即时编译器可以使数字代码运行速度提高100倍。您也可以编写自定义的C或Cython扩展。

但有一个更好的办法。求解常微分方程组是非常普遍的，因此scipy是Python中的核心科学库之一，它包装了快速的，经过测试的Fortran库来解决它们。通过使用scipy，您可以同时获得算法改进（因为集成商将拥有更高的订单）和快速实施。

解决洛伦茨95模型的一组的扰动的初始条件是这样的：

import numpy as np 


def lorenz95(x, t): 
    return np.roll(x, 1) * (np.roll(x, -1) - np.roll(x, 2)) - x + F 

if __name__ == '__main__': 
    import matplotlib.pyplot as plt 
    from scipy.integrate import odeint 
    SIZE = 40 
    F = 8 
    t = np.linspace(0, 10, 1001) 
    x0 = np.random.random(SIZE) 
    for perturbation in 0.1 * np.random.randn(5): 
     x0i = x0.copy() 
     x0i[0] += perturbation 
     x = odeint(lorenz95, x0i, t) 
     plt.plot(t, x[:, 0]) 
    plt.show()

和输出（设定np.random.seed(7)，你可以是不同的）是很好的混乱。初始条件下的小扰动（仅在他的坐标之一中）产生了非常不同的解决方案： Lorenz-95 dynamical system

但是，它真的比欧拉时间步进快吗？对于dt = 0.01它看起来几乎快了三倍，但解决方案除了在最开始时不匹配。 Euler vs odeint

如果dt减小，由欧拉方法提供的解决方案变得越来越类似于odeint解决方案，但它需要更长的时间。请注意较小的dt，后来的欧拉解决方案是如何松动odeint解决方案的。最精确的欧拉解决方案花费了600倍的时间来计算t = 6时的解决方案，而不是t = 10时的解决方案。请参阅完整脚本here。 Euler vs odeint

最后，这个系统是如此不稳定，我甚至不认为odeint解决方案在所有绘图时间都是准确的。

来源

2013-05-08 17:23:40 jorgeca

感谢@jorgeca的建议。你的代码做了一些重要的改变： 1.使用'np.roll'意味着整个数组可以在一次扫描中更新，所以我现在有一个新的欧拉时间步进方案：'def euler2（self ）：self.x = self.x + dt *（np.roll（self.x，1）*（np.roll（self.x，-1） - np.roll（self.x，2）） - self .x + F）'。这比原始实施速度快大约3倍。 2.使用'odeint'，而不是粗糙的欧拉时间步。这可能更稳定，但速度比较慢（比欧拉例程慢2倍）。我想我将不得不使用编译代码。 – 2013-05-09 15:02:51

@NeillBowler你的欧拉例程更快的结果是虚假的，除非在一开始。我会在明天展开我的回答，告诉你，但总结是“在解决方案中实现同样的错误，欧拉时间步长将比odeint多出几个数量级的时间”。 – jorgeca 2013-05-09 19:52:58

@NeillBowler我已经扩展了答案。 – jorgeca 2013-05-11 21:42:19

如何优化python中的时间步进算法？

回答

相关问题