Векторизация и умножение матриц на скаляры

Я новичок в python/numpy. Мне нужно сделать следующий расчет: для массива дискретных времен t вычислить $e^{At}$ для матрицы $2\times 2$ $A$

Что я сделал:

def calculate(t_,x_0,v_0,omega_0,c):

    # define A
    a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
    A =np.matrix([[a_11,a_12], [a_21, a_22]]) 
    print A
    # use vectorization 
    temps = np.array(t_)
    A_ = np.array([A  for k in  range (1,n+1,1)])
    temps*A_
    x_=scipy.linalg.expm(temps*A)
    v_=A*scipy.linalg.expm(temps*A)
    return x_,v_

n=10
omega_0=1
c=1
x_0=1
v_0=1
t_ = [float(5*k*np.pi/n)  for k in  range (1,n+1,1)]
x_, v_ = calculate(t_,x_0,v_0,omega_0,c)

Однако я получаю эту ошибку при умножении A_ (массив, содержащий n раз A ) и temps (содержащий время, для которого я хочу вычислить exp(At) :

ValueError: операнды не могли передаваться вместе с формами (10,) (10,2,2)

Насколько я понимаю векторизацию, каждый элемент в A_ будет умножаться на элемент с тем же индексом из temps; но я думаю, что я не понимаю это правильно. Любая помощь/комментарии очень ценятся


python numpy array-broadcasting
person PerelMan    schedule 05.02.2018    source источник
comment
В стороне, как я только что заметил это; есть ли какая-то особая причина использовать np.matrix? Создатели Numpy в основном хотят, чтобы его не существовало, поэтому, как правило, лучше использовать многомерный np.array, за исключением очень избранного количества операций.   -  person roganjosh    schedule 06.02.2018
comment
@Azevedo Я добавил экспоненту. Я не указал это в вопросе, потому что хотел изолировать проблему   -  person PerelMan    schedule 06.02.2018
comment
@roganjosh, цель этой строки - сделать векторизацию. Я думал, что мне нужно предоставить 2 массива numpy с одинаковым количеством элементов, чтобы иметь возможность выполнять векторизованное умножение. Не было особой причины использовать np.matrix, я просто использовал его, потому что это было первое решение, которое я нашел, поскольку я новичок в numpy.   -  person PerelMan    schedule 06.02.2018

Ответы (2)


arrow_upward
3
arrow_downward

Чистый numpy расчет t_ (создает массив вместо списка):

In [254]: t = 5*np.arange(1,n+1)*np.pi/n
In [255]: t
Out[255]: 
array([ 1.57079633,  3.14159265,  4.71238898,  6.28318531,  7.85398163,
        9.42477796, 10.99557429, 12.56637061, 14.13716694, 15.70796327])

In [256]: a_11,a_12, a_21, a_22=0,1,-omega_0^2,-c
In [257]: a_11
Out[257]: 0
In [258]: A = np.array([[a_11,a_12], [a_21, a_22]]) 
In [259]: A
Out[259]: 
array([[ 0,  1],
       [-3, -1]])
In [260]: t.shape
Out[260]: (10,)
In [261]: A.shape
Out[261]: (2, 2)
In [262]: A_ = np.array([A  for k in  range (1,n+1,1)])
In [263]: A_.shape
Out[263]: (10, 2, 2)

A_ это np.ndarray. Я тоже сделал A np.ndarray; у вас np.matrix, но у вас A_ все равно будет np.ndarray. np.matrix может быть только 2d, тогда как A_ — 3d.

Таким образом, t * A будет поэлементным умножением массива, отсюда и ошибка вещания (10,) (10,2,2).

Чтобы правильно выполнить это поэлементное умножение, вам нужно что-то вроде

In [264]: result = t[:,None,None]*A[None,:,:]
In [265]: result.shape
Out[265]: (10, 2, 2)

Но если вам нужно матричное умножение (10,) на (10,2,2), то einsum сделает это легко:

In [266]: result1 = np.einsum('i,ijk', t, A_)
In [267]: result1
Out[267]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

np.dot не может этого сделать, потому что его правило "последний со 2-м до последнего". tensordot может, но мне удобнее einsum.

Но это выражение einsum делает очевидным (для меня), что я могу получить то же самое из поэлементного *, просуммировав по 1-й оси:

In [268]: (t[:,None,None]*A[None,:,:]).sum(axis=0)
Out[268]: 
array([[   0.        ,   86.39379797],
       [-259.18139392,  -86.39379797]])

Или (t[:,None,None]*A[None,:,:]).cumsum(axis=0), чтобы каждый раз получать 2x2.

person hpaulj    schedule 06.02.2018
comment
Вы можете использовать точку с осью вращения. Кроме того, небольшая придирка, но я уверен, что класс называется np.ndarray, а не np.array. +1 - person Mad Physicist; 06.02.2018

arrow_upward
0
arrow_downward

Это то, что я бы сделал.

import numpy as np
from scipy.linalg import expm
A = np.array([[1, 2], [3, 4]])
for t in np.linspace(0, 5*np.pi, 20):
    print(expm(t*A))

Здесь нет попытки векторизации. Функция expm применяется к одной матрице за раз и, безусловно, занимает большую часть времени вычислений. Не нужно беспокоиться о стоимости умножения A на скаляр.

person Community    schedule 05.02.2018