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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With