Vectorization And Matrix Multiplication By Scalars

I am new to python/numpy. I need to do the following calculation: for an array of discrete times t, calculate $e^{At}$ for a $2\times 2$ matrix $A$ What I did: def calculate(t_,x_

Solution 1:

A pure numpy calculation of t_ is (creates an array instead of a list):

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_ is np.ndarray. I made A a np.ndarray as well; yours is np.matrix, but your A_ will still be np.ndarray. np.matrix can only be 2d, where as A_ is 3d.

So t * A will be array elementwise multiplication, hence the broadcasting error, (10,) (10,2,2).

To do that elementwise multiplication right you need something like

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

---

But if you want matrix multiplication of the (10,) with (10,2,2), then einsum does it easily:

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

np.dot can't do it because its rule is 'last with 2nd to last'. tensordot can, but I'm more comfortable with einsum.

But that einsum expression makes it obvious (to me) that I can get the same thing from the elementwise *, by summing on the 1st axis:

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

Or (t[:,None,None]*A[None,:,:]).cumsum(axis=0) to get a 2x2 for each time.