scipy.integrate.solve_ivp 矢量化

Question

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

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”

Answer 1

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

更改以在每次调用时f显示形状 a ：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，则y传递给的f是 (2,), 1d；True 则为 (2,1)。我猜如果求解器方法需要的话，它可能是 (2,2) 甚至 (2,3)。这可以加快执行速度，同时减少对f. 在这种情况下，这并不重要。

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

使用 scipy 对向量输入的标量值函数进行数值求积

相关错误/问题讨论：

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