scipy.integrate.solve_ivp 矢量化

Posted 2023-03-12

【中文标题】scipy.integrate.solve_ivp 矢量化

【英文标题】：scipy.integrate.solve_ivp vectorized

【发布时间】：2019-07-26 05:26:33

【问题描述】：

尝试对solve_ivp 使用矢量化选项，但奇怪的是它会抛出一个错误，即y0 必须是一维的。 MWE：

from scipy.integrate import solve_ivp
import numpy as np
import math

def f(t, y):
    theta = math.pi/4
    ham = np.array([[1,0],[1,np.exp(-1j*theta*t)]])
    return-1j * np.dot(ham,y)


def main():

    y0 = np.eye(2,dtype= np.complex128)
    t0 = 0
    tmax = 10**(-6)
    sol=solve_ivp( lambda t,y :f(t,y),(t0,tmax),y0,method='RK45',vectorized=True)
    print(sol.y)

if __name__ == '__main__':
    main()

调用签名是 fun(t, y)。这里 t 是一个标量，对于 ndarray y 有两个选项：它可以具有形状 (n,)；那么 fun 必须返回具有形状 (n,) 的 array_like。或者，它可以具有形状 (n, k)；那么 fun 必须返回一个形状为 (n, k) 的 array_like，即每一列对应于 y 中的单个列。两个选项之间的选择由矢量化参数确定（见下文）。矢量化实现允许通过有限差分更快地逼近雅可比行列式（刚性求解器需要）。

错误：

ValueError：y0 必须是一维的。

Python 3.6.8

scipy.版本 '1.2.1'

【参考方案1】：

这里vectorize的含义有点混乱。这并不意味着y0 可以是二维的，而是传递给您的函数的y 可以是二维的。换句话说，func 可以同时在多个点进行评估，如果求解器愿意的话。多少分取决于求解器，而不是您。

更改f 以在每次调用时显示y 的形状：

def f(t, y):
    print(y.shape)
    theta = math.pi/4
    ham = np.array([[1,0],[1,np.exp(-1j*theta*t)]])
    return-1j * np.dot(ham,y)

一个示例调用：

In [47]: integrate.solve_ivp(f,(t0,tmax),[1j,0],method='RK45',vectorized=False) 
(2,)
(2,)
(2,)
(2,)
(2,)
(2,)
(2,)
(2,)
Out[47]: 
  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 8
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.e+00, 1.e-06])
 t_events: None
        y: array([[0.e+00+1.e+00j, 1.e-06+1.e+00j],
       [0.e+00+0.e+00j, 1.e-06-1.e-12j]])

同样的电话，但使用vectorize=True：

In [48]: integrate.solve_ivp(f,(t0,tmax),[1j,0],method='RK45',vectorized=True)  
(2, 1)
(2, 1)
(2, 1)
(2, 1)
(2, 1)
(2, 1)
(2, 1)
(2, 1)
Out[48]: 
  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 8
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([0.e+00, 1.e-06])
 t_events: None
        y: array([[0.e+00+1.e+00j, 1.e-06+1.e+00j],
       [0.e+00+0.e+00j, 1.e-06-1.e-12j]])

如果为 False，则传递给 f 的 y 是 (2,), 1d； True 为 (2,1)。如果求解器方法如此需要，我猜它可能是 (2,2) 甚至 (2,3)。这可以加快执行速度，减少对f 的调用。在这种情况下，没关系。

quadrature 有一个类似的vec_func 布尔参数：

Numerical Quadrature of scalar valued function with vector input using scipy

相关的错误/问题讨论：

https://github.com/scipy/scipy/issues/8922

【讨论】：

非常感谢。这有助于消除混乱。因此，如果我想对大量一维数组重复该过程，是循环遍历我的最佳选择吗？如果您能指出一些比遍历每个数组更快的方法，那就太好了。我有 (243) 个大小为 243 的一维数组。因此以这种方式消耗大量时间。再次感谢您！

目前未得到答复，链接返回此处：What does the letter k mean in the documentation of solve_ivp function of Scipy?

@uhoh，正如我所展示的，如果vectorized，k 至少为 1，并且可能更大。除了确保您的函数返回类似尺寸的结果之外，您可以或必须做的事情并不多。如果您遇到大于 1 的情况（例如在 jacobian 评估中），请发布答案供大家学习。

@hpaulj 谢谢，我期待y0 可以是二维的。当solve_ivp 第一次被宣布时，我很兴奋，因为我读过一些东西，关于维度不仅仅是一维的，但现在我开始使用它，我发现情况并非如此。

@uhoh，如果每个初始值数组的最优 ODE 步进不同，那么solve_ivp 无法以 2d 方式执行计算。每个都必须单独解决。或者你可以让y0n*k 长。对于某些问题，这样可能会更快，对于其他k 大小的迭代 n y0 可能是。您希望 2d 案例如何工作？