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

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

总有琐碎的autojit解决方案：

def forward_euler(initial_x, initial_y, num_steps, h):

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

     x[0, 0] = initial_x[0] # initial condition 1st solution
     y[0, 0] = initial_y[0]

     x[0, 1] = initial_x[1]  # initial condition 2nd solution
     y[0, 1] = initial_y[1]

     for step in xrange(int(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

时间安排：

from numba import autojit
jit_forward_euler = autojit(forward_euler)

%timeit forward_euler([10,0], [5,0.0000000001], 1E4, 0.01)
1 loops, best of 3: 385 ms per loop

%timeit jit_forward_euler([10,0], [5,0.0000000001], 1E4, 0.01)
100 loops, best of 3: 3.51 ms per loop

Answer 2

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

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

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

您可以使用可变的默认参数来实现这一点，该参数在定义函数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_tmp 视为函数 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))