矢量化前向欧拉方法的微分方程组

Question

我正在数值求解一阶微分方程组的x(t).该系统是:

dx/dt = y
dy/dt = -x - a*y(x ^ 2 + y ^ 2 -1)

我已经实现了Forward Euler方法来解决这个问题,如下所示:

def forward_euler():
    h = 0.01
    num_steps = 10000

    x = np.zeros([num_steps + 1, 2]) # steps, number of solutions
    y = np.zeros([num_steps + 1, 2])
    a = 1.

    x[0, 0] = 10. # initial condition 1st solution
    y[0, 0] = 5.

    x[0, 1] = 0.  # initial condition 2nd solution
    y[0, 1] = 0.0000000001

    for step in xrange(num_steps):
        x[step + 1] = x[step] + h * y[step]
        y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1))

    return x, y

现在我想进一步向量化代码并将x和y保持在同一个数组中,我已经提出了以下解决方案:

def forward_euler_vector():
    num_steps = 10000
    h = 0.01

    x = np.zeros([num_steps + 1, 2, 2]) # steps, variables, number of solutions
    a = 1.

    x[0, 0, 0] = 10. # initial conditions 1st solution
    x[0, 1, 0] = 5.  

    x[0, 0, 1] = 0.  # initial conditions 2nd solution
    x[0, 1, 1] = 0.0000000001

    def f(x): 
        return np.array([x[1],
                         -x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)])

    for step in xrange(num_steps):
        x[step + 1] = x[step] + h * f(x[step])

    return x

问题:forward_euler_vector()有效,但是这是最好的矢量化方法吗？我问,因为矢量化版本在我的笔记本电脑上运行速度大约20毫秒:

In [27]: %timeit forward_euler()
1 loops, best of 3: 301 ms per loop

In [65]: %timeit forward_euler_vector()
1 loops, best of 3: 320 ms per loop

Answer 1

@Ophion 评论很好地解释了发生的事情。array()对inside的调用f(x)引入了一些开销，从而消除了在表达式 中使用矩阵乘法的好处h * f(x[step])。

正如他所说，您可能有兴趣查看scipy.integrate一组不错的数值积分器。

为了解决代码矢量化的问题，您希望避免每次调用时都重新创建数组f。您希望初始化该数组一次，并在每次调用时返回修改后的数组。这类似于staticC/C++ 中的变量。

您可以使用可变的默认参数来实现此目的，该默认参数在定义 function 时解释一次f(x)，并且具有本地作用域。由于它必须是可变的，因此您将其封装在单个元素的列表中：

 def f(x,static_tmp=[empty((2,2))]): 
    static_tmp[0][0]=x[1]
    static_tmp[0][1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)
    return static_tmp[0]

通过对代码的修改，数组创建的开销消失了，并且在我的机器上我获得了一个小小的改进：

%timeit forward_euler()        #258ms
%timeit forward_euler_vector() #248ms

这意味着用 numpy 优化矩阵乘法的增益相当小，至少在手头的问题上是这样。

f您可能还想立即删除该函数，在 for 循环中执行其操作，从而消除调用开销。然而，默认参数的这个技巧也可以应用于scipy更通用的时间积分器，其中您必须提供一个函数f。

编辑：正如 Jaime 所指出的，另一种方法是将其视为static_tmpfunction 的属性f，并在声明该函数之后但在调用它之前创建它：

 def f(x): 
    f.static_tmp[0]=x[1]
    f.static_tmp[1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)
    return f.static_tmp
 f.static_tmp=empty((2,2))