Python 和matlab 关于DMD（动态模态分解）的实现和对比 21/06/08

参考资料：寇家庆，张伟伟. 动力学模态分解及其在流体动力学中的应用. 空气动力学报. 2018

 

通过numpy和matlab对于该论文中的典型算例1进行复现，并基于实数算例进行对比分析.

1 典型算例DMD方法代码和对比：

matlab:

X = [0.5,0.55-0.25i,0.48-0.55i,0.253-0.845i;1,1+0.1i,0.99+0.2i,0.97+0.299i;0.8,0.16-0.008i,0.0319-0.0032i,0.0064-0.001i];
Y = [0.55-0.25i,0.48-0.55i,0.253-0.845i,-0.1442-1.056i;1+0.1i,0.99+0.2i,0.97+0.299i,0.9401+0.396i;0.16-0.008i,0.0319-0.0032i,0.0064-0.001i,0.0013-0.0003i];
%X = [4,4,4,4; 3,9,27,81; 1,6,36,216];
%Y = [4,4,4,4; 9,27,81,243; 6,36,216,1476];

[u,s,v] = svd(X,'econ');
a1 = u' * Y *v * s^-1;
[V,D]= eig(a1);

其中，econ项将分解后的vt去除一列，从4*4矩阵变为4*3矩阵

python:

# python中特征值sigma是以列表形式出现，因此需要将其转化为对角阵

#numpy中linalg的SVD项目前似乎尚不支持econ项，因此需要手动进行去值。

# 两种工具进行SVD分解后的Vt项并不是互为转置，而是互为共轭转置。 因此，在某些情况下，U项的行列式结果会互为相反数，但不影响特征值的计算。

#最终进行eig特征值计算时，python中特征值项为V, 而matlab则为D

import numpy as np
from numpy import linalg as lg
# 条件矩阵
X = np.mat([[0.5, 0.55-0.25j, 0.48-0.55j, 0.253-0.845j],
           [1, 1+0.1j, 0.99+0.2j, 0.97+0.299j],
           [0.8, 0.16-0.008j, 0.319-0.0032j, 0.0064-0.001j]])
Y = np.mat([[0.55-0.25j,0.48-0.55j,0.253-0.845j,-0.1442-1.056j],
            [1+0.1j,0.99+0.2j,0.97+0.299j,0.9401+0.396j],
            [0.16-0.008j,0.0319-0.0032j,0.0064-0.001j,0.0013-0.0003j]])

u,sigma,vt = lg.svd(X)
#将sigma转为矩阵
sigma_ls = np.full((len(sigma),len(sigma)),0)
sigma_ls = sigma_ls.astype(np.float64)
for i in range(len(sigma_ls)) :
    for j in range(len(sigma_ls)):
        if i == j: 
            sigma_ls[i][j] = sigma[i]
sigma = sigma_ls 
#将vt去数据，相当于econ
vt_ls = np.full((X.shape[0],X.shape[1]),0)
vt_ls = vt_ls.astype(np.complex64)

for i in range(vt_ls.shape[0]):
     for j in range(vt_ls.shape[1]):
          vt_ls[i,j] = vt[i,j]
vt = vt_ls

#u^h
uh = u.T.conjugate()

#sigma的逆矩阵
sigma_1 = lg.inv(sigma)
#相似变换矩阵计算
A_1 = (uh * Y * vt.T.conjugate() * sigma_1)
v1,d1 = lg.eig(A_1)  #求解特征值，即为系统矩阵
print(np.around(v1,3))

结果对比：

matlab:

D =

   1.1000 - 0.5000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.1000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.2000 - 0.0100i

python:

[1.1  -0.5j   1.   +0.1j   0.167+0.005j]

完整复现算例设定的特征值

 

通过过程分析可知，对于复数矩阵，在SVD分解阶段两种工具差异即较大：

（需要注意：此处的两个vt项应互为共轭转置）

matlab:

u =
  -0.5215 + 0.0000i  -0.6396 + 0.0000i   0.5648 + 0.0000i
  -0.5171 - 0.6473i   0.3078 + 0.1550i  -0.1289 - 0.4222i
  -0.1582 - 0.1289i  -0.0873 + 0.6816i  -0.2449 + 0.6528i

s =
    2.3981         0         0
         0    0.9074         0
         0         0    0.2808

v =
  -0.3771 - 0.3129i  -0.0902 + 0.7717i  -0.1512 + 0.3563i
  -0.3724 - 0.3119i  -0.0528 + 0.0801i   0.3386 - 0.5897i
  -0.3738 - 0.3456i   0.0262 - 0.2628i   0.1749 - 0.2190i
  -0.3453 - 0.3815i   0.2004 - 0.5266i  -0.3940 + 0.3925i

python:

u:  [[-0.5127+0.j     -0.6781+0.j      0.5266+0.j    ]
 [-0.5132-0.6352j  0.3336+0.1338j -0.0701-0.4462j]
 [-0.2036-0.1699j -0.0813+0.6359j -0.3028+0.6534j]]

s:  [2.4329 0.871  0.2764]

vt:  [[-0.3833+0.3169j -0.3658+0.3045j -0.3887+0.3547j -0.3365+0.3688j]
 [-0.0809-0.7377j -0.0506-0.0367j  0.0041+0.1201j  0.2192+0.6188j]
 [-0.1775-0.2771j  0.4387+0.7432j -0.0164-0.2512j -0.256 -0.134j ]
 [ 0.2767+0.1308j -0.0164+0.1565j -0.6949-0.4029j  0.4722+0.1156j]]

但均可通过U,s,vt对X进行还原

2.实数算例DMD代码和结果对比

选取一个符合特征要求的X和Y矩阵，见代码

matlab:

%X = [0.5,0.55-0.25i,0.48-0.55i,0.253-0.845i;1,1+0.1i,0.99+0.2i,0.97+0.299i;0.8,0.16-0.008i,0.0319-0.0032i,0.0064-0.001i];
%Y = [0.55-0.25i,0.48-0.55i,0.253-0.845i,-0.1442-1.056i;1+0.1i,0.99+0.2i,0.97+0.299i,0.9401+0.396i;0.16-0.008i,0.0319-0.0032i,0.0064-0.001i,0.0013-0.0003i];
X = [4,4,4,4; 3,9,27,81; 1,6,36,216];
Y = [4,4,4,4; 9,27,81,243; 6,36,216,1476];

[u,s,v] = svd(X,'econ');
a1 = u' * Y *v * s^-1;
[V,D]= eig(a1);

python:

import numpy as np
from numpy import linalg as lg

#假设 A = [[1,0,0],
#          [0,3,0],
#          [0,0 6]]
#初始条件为[4,3,1]T

X = np.mat([[4,4,4,4],
            [3,9,27,81],
           [1,6,36,216]])

Y = np.mat([[4,4,4,4],
            [9,27,81,243],
            [6,36,216,1476]])

u,sigma,vt = lg.svd(X)

sigma_ls = np.full((len(sigma),len(sigma)),0)
sigma_ls = sigma_ls.astype(np.float64)
for i in range(len(sigma_ls)) :
    for j in range(len(sigma_ls)):
        if i == j: 
            sigma_ls[i][j] = sigma[i]
sigma = sigma_ls 

vt_ls = np.full((X.shape[0],X.shape[1]),0)
vt_ls = vt_ls.astype(np.complex64)

for i in range(vt_ls.shape[0]):
     for j in range(vt_ls.shape[1]):
          vt_ls[i,j] = vt[i,j]
vt = vt_ls

#u^h
uh = u.T.conjugate() 

#sigma的逆矩阵
sigma_1 = lg.inv(sigma)

A_1 = (uh * Y * vt.T.conjugate() * sigma_1)
v1,d1 = lg.eig(A_1)
print(np.around(A_1,4))
print(np.around(v1,3))

结果对比：

matlab:

D =
    7.7886         0         0
         0    3.0000         0
         0         0    1.0000

python:

[7.789+0.j 3.   +0.j 1.   +0.j]

结果一致，可能是样本数过少，估计值偏差稍大。

结论：

对于应用型算例来说，二者具备相同精度和功能。但python编写时，应考虑是否符合主流算法说明，避免踩坑。