Dynamic Mode Decomposition Example



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def generate_data():
    s1 = np.sin(np.linspace(0, 50, 1000))
    s2 = np.cos(np.linspace(0, 20, 1000))
    x = np.tile(s1, (9,1))
    x = (x.T*np.array([1,2,3,4,5,4,3,2,1])).T
    y = np.tile(s2, (9,1))
    y = (y.T*np.array([1,2,3,4,5,4,3,2,1])).T
    x = np.concatenate((x, np.zeros((16,1000))), axis=0)
    y = np.concatenate((np.zeros((16,1000)), y), axis=0)
    out = x + y
    out = out + np.random.randn(*out.shape)*0.4
    return out

Make some data



In [2]:

    
data = generate_data()
data.shape
plt.plot(data.T);

PCA



In [3]:

    
X = data[:,:-1];
Y = data[:, 1:];
[U, S, Vh] = np.linalg.svd(X);
plt.plot(S,'*-');



In [4]:

    
for i in range(3):
    plt.subplot(3,1,i+1);
    plt.plot(np.dot(X.T,U[:,i]));

Using DMD with r=2 (=nk)



In [5]:

    
nk = 2
A_tilda = U[:,:nk].T.dot(Y).dot(Vh[:nk,:].T)/S[:nk]



In [6]:

    
A_evals, W = np.linalg.eig(A_tilda)



In [7]:

    
phi = Y.dot(Vh[:nk, :].T).dot(W)/S[:nk]

Solving for z



In [8]:

    
z = np.linalg.pinv(phi).dot(X[:,0])

Reconstruction



In [9]:

    
Xhat0 = phi.dot(z)



In [14]:

    
Xhat4 = phi.dot(np.diag(A_evals**4)).dot(z)



In [29]:

    
np.corrcoef(np.real(Xhat4), X[:,5]) # Not sure if this should be 5th or 6th sample









    Out[29]:





array([[ 1.        ,  0.96520607],
       [ 0.96520607,  1.        ]])

Spectrum



In [10]:

    
dt = 0.001 # sampling rate of 1KHz
f = np.imag(np.log(A_evals))/(2*np.pi)/dt
P = np.diag(phi.T.dot(phi))



In [30]:

    
plt.stem(f, P);