scipy.integrate.solve_ivp矢量化 [英] scipy.integrate.solve_ivp vectorized

问题描述

尝试将向量化选项用于solve_ivp，奇怪的是，它抛出一个错误，y0必须为一维. MWE:

Trying to use the vectorized option for solve_ivp and strangely it throws an error that y0 must be 1 dimensional. 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，)；它可以具有形状(n，).然后乐趣必须返回形状为(n，)的array_like.或者，它可以具有形状(n，k)；那么fun必须返回一个形状为(n，k)的array_like，即每一列对应于y中的单个列.这两个选项之间的选择由矢量化参数确定(请参见下文).向量化的实现可以通过有限差分(对于刚性求解器而言)可以更快地逼近雅可比行列式.

The calling signature is fun(t, y). Here t is a scalar, and there are two options for the ndarray y: It can either have shape (n,); then fun must return array_like with shape (n,). Alternatively it can have shape (n, k); then fun must return an array_like with shape (n, k), i.e. each column corresponds to a single column in y. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows a faster approximation of the Jacobian by finite differences (required for stiff solvers).

错误:

ValueError:y0必须为一维.

Python 3.6.8

Python 3.6.8

scipy.版本 '1.2.1'

scipy.version '1.2.1'

推荐答案

此处vectorize的含义有些混乱.这并不意味着y0可以是2d，而是传递给函数的y可以是2d.换句话说，如果求解器愿意，可以一次在多个点上对func进行评估.多少点取决于求解器，而不是您.

The meaning of vectorize here is a bit confusing. It doesn't mean that y0 can be 2d, but rather that y as passed to your function can be 2d. In other words that func may be evaluated at multiple points at once, if the solver so desires. How many points is up to the solver, not you.

更改f以在每次调用时显示形状为y:

Change the f to show the shape a y at each call:

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的调用，从而加快执行速度.在这种情况下，没关系.

With False, the y passed to f is (2,), 1d; with True it is (2,1). I'm guessing it could be (2,2) or even (2,3) if the solver method so desires. That could speed up the execution, with fewer calls to f. In this case, it doesn't matter.

quadrature具有类似的vec_func布尔参数:

quadrature has a similar vec_func boolean parameter:

具有使用scipy进行矢量输入

相关的错误/问题讨论:

A related bug/issue discussion:

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

这篇关于scipy.integrate.solve_ivp矢量化的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！