Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
574 views
in Technique[技术] by (71.8m points)

python - scipy.integrate.solve_ivp vectorized

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'

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

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


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

1.4m articles

1.4m replys

5 comments

57.0k users

...