scipy.integrate.solve_ivp векторизация - PullRequest
Новая
       10

scipy.integrate.solve_ivp векторизация

0 голосов
/ 04 марта 2019

Попытка использовать векторизованную опцию для solve_ivp и странным образом выдает ошибку, что y0 должен быть одномерным.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()

Вызывающая подпись забавная (t, y).Здесь t - скаляр, и есть два варианта ndarray y: он может иметь форму (n,);тогда fun должен вернуть array_like с shape (n,).В качестве альтернативы он может иметь форму (n, k);тогда fun должен возвращать array_like с shape (n, k), то есть каждый столбец соответствует одному столбцу в y.Выбор между двумя вариантами определяется векторизованным аргументом (см. Ниже).Векторизованная реализация позволяет быстрее приближать якобиан с помощью конечных разностей (требуется для жестких решателей).

Ошибка:

ValueError: y0 должно быть одномерным.

Python 3.6.8

scipy. версия '1.2.1'

1 Ответ

0 голосов
/ 05 марта 2019

Значение vectorize здесь немного сбивает с толку.Это не значит, что y0 может быть 2d, но скорее всего y, передаваемый вашей функции, может быть 2d.Другими словами, что func может быть оценено в нескольких точках одновременно, если решатель этого пожелает.Сколько очков зависит от решателя, а не от вас.

Измените f, чтобы при каждом вызове отображать форму 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

...