DMD in Python one dimension

Introduction

Dynamic Mode Decomposition (DMD) was announced in 2008 [^ 1] [^ 2]. DMD is a decomposition method that combines principal component analysis and Fourier transform, and is a method that can separate data into dimensional and temporal features.

The commentary articles that I referred to are below.

• http://www.pyrunner.com/weblog/2016/07/25/dmd-python/ -<https://iqujack-lequina.hatenablog.com/entry/2018/05/20/ Memorandum of Understanding on Dynamic Mode Decomposition>
• https://qiita.com/Miyabi1456/items/702f62c9bcd9c063d664
• https://qiita.com/harmegiddo/items/84552c32f4b75c26878a

The eigenvalues and eigenvectors of DMD are complex numbers, and the complex components of the eigenvalues represent the decomposition in the time direction. However, there is a problem that it cannot be separated well if it is a one-dimensional standing wave. The explanation is in the following article.

-<https://iqujack-lequina.hatenablog.com/entry/2018/12/02/ Standing wave DMD>

I don't think that DMD will be used for 1D data, but since the Python code has not been released, I will describe it for reference. Here are some examples that cannot be decomposed well and some that can. The base paper is Tu 2014 [^ 3], and the code is DMD Book. It is based on the Matlab code of.

environment

• Windows 10 home

• Anaconda(Python 3.7.6)

• Numpy(1.18.1)

Disassembly target

Find the period from the numerical data of the one-dimensional sine wave. The period is $ 1j $ for complex numbers. $ y = \sin(x) $

DMD example that doesn't work

In the example that doesn't work, the eigenvalues are real and the complex numbers don't come out because the data is only one-dimensional. Therefore, the periodic component cannot be decomposed.

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

t = np.arange(0, 10.01, 0.01)
x = np.sin(t)
dt = t[2] - t[1]

#Data matrix
X1 = x[:-1]
X2 = x[1:]

#SVD
U,S,Vh = LA.svd(X1[np.newaxis,:],False)
V = Vh.T

#Atilde with A reduced by left singular vector
Atilde = np.dot(np.dot(np.dot(U.T, X2[np.newaxis,:]), V), LA.inv(np.diag(S)))

#Find the eigenvalues and eigenvectors of Atilde
Lam, W = LA.eig(Atilde)

print("eigenvalue:",Lam)#Lam:1.It becomes 00025541 and is a real number

#Find the eigenvector of A from the eigenvector of Atilde
Phi = np.dot(np.dot(np.dot(X2[np.newaxis,:], V), LA.inv(np.diag(S))), W)

#Discrete to continuous exp(**)of**Seeking
Omega = np.log(Lam)/dt

#Restore the original function from Omega and Phi.
b = np.dot(LA.pinv(Phi * np.exp(Omega * t)).T, x)
x_dmd = b * Phi * np.exp(Omega * t)

plt.plot(t, x_dmd[0,:])
plt.show()

DMD example that works

It is introduced in Explanatory article, but the data is shifted in the time direction. Make another dimension and make it two-dimensional data. By doing this, the complex component appears in the eigenvalue, and $ \ sin (x) $ can be found in the numerical data. That's it.

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

t = np.arange(0, 10.01, 0.01)
x = np.sin(t)
dt = t[2] - t[1]

#Data matrix
#The main difference is here. It is two-dimensional.
X_aug1 = np.array([x[:-2], x[1:-1]])
X_aug2 = np.array([x[1:-1], x[2:]])

#SVD
U,S,Vh = LA.svd(X_aug1,False)
V = Vh.conj().T

#Atilde with A reduced by left singular vector
#U takes the complex conjugate and makes it a transposed matrix.
Atilde = np.dot(np.dot(np.dot(U.conj().T, X_aug2), V), LA.inv(np.diag(S)))

#Find the eigenvalues and eigenvectors of Atilde
Lam, W = LA.eig(Atilde)

print("eigenvalue:", Lam)#Lam:[0.9995+0.00999983j 0.9995-0.00999983j]Complex number

#Find the eigenvector of A from the eigenvector of Atilde
Phi = np.dot(np.dot(np.dot(X_aug2, V), LA.inv(np.diag(S))), W)

#Discrete to continuous exp(**)of**Seeking
Omega = np.diag(np.log(Lam)/dt)

print("Omega:", Omega)#Just in case, it will be 1j

#Restore the original function from Omega and Phi.
#What you are doing is the same as in 1D, but the writing style is different.
b = np.dot(LA.pinv(Phi), X_aug1[:,0])
x_dmd = np.zeros([2,len(t)], dtype='complex')
for i, _t in enumerate(t):
    x_dmd[:,i] = np.dot(np.dot(Phi, np.exp(Omega * _t)), b)

plt.plot(t, x_dmd[0,:].real)
plt.show()

References