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

问题描述

尝试对solve_ivp 使用矢量化选项，奇怪的是它抛出了y0 必须是一维的错误.MWE:

from scipy.integrate import solve_ivp将 numpy 导入为 np导入数学定义 f(t, y):θ = math.pi/4火腿 = np.array([[1,0],[1,np.exp(-1j*theta*t)]])return-1j * np.dot(ham,y)定义主():y0 = np.eye(2,dtype= np.complex128)t0 = 0tmax = 10**(-6)sol=solve_ivp( lambda t,y :f(t,y),(t0,tmax),y0,method='RK45',vectorized=True)打印(sol.y)如果 __name__ == '__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'

解决方案

vectorize 这里的意思有点混乱.这并不意味着 y0 可以是 2d，而是传递给您的函数的 y 可以是 2d.换句话说，func 可以一次在多个点进行计算，如果求解器愿意的话.多少分取决于求解者，而不是你.

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

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

示例调用:

在[47]中:integre.solve_ivp(f,(t0,tmax),[1j,0],method='RK45',vectorized=False)(2,)(2,)(2,)(2,)(2,)(2,)(2,)(2,)出[47]:消息:求解器成功到达积分区间的末端."非传染性疾病:8尼耶夫:0零:0溶胶:无状态:0成功:正确t: 数组([0.e+00, 1.e-06])t_events:无y: 数组([[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:

在[48]中:integre.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)出[48]:消息:求解器成功到达积分区间的末端."非传染性疾病:8尼耶夫:0零:0溶胶:无状态:0成功:正确t: 数组([0.e+00, 1.e-06])t_events:无y: 数组([[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 布尔参数:

标量值函数的数值求积与使用scipy的向量输入

相关的错误/问题讨论:

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

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()

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).

Error :

ValueError: y0 must be 1-dimensional.

Python 3.6.8

scipy.version '1.2.1'

解决方案

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.

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)

A sample call:

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]])

Same call, but with 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]])

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 has a similar vec_func boolean parameter:

Numerical Quadrature of scalar valued function with vector input using scipy

A related bug/issue discussion:

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

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