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.0000i0.0000 + 0.0000i   1.0000 + 0.1000i   0.0000 + 0.0000i0.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.6528is =2.3981         0         00    0.9074         00         0    0.2808v =-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]TX = 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         00    3.0000         00         0    1.0000

python:

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

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

结论：

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

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

利用matlab实现DMD动态模态分解(在一维信号或二维流场矢量中的应用）
利用matlab实现DMD动态模态分解(在一维信号或二维流场矢量中的应用) 0 前言 0.1 特征根的计算与含义 1 DMD的基本思路 2 一维DMD算法 3 二维DMD算法 4 总结 (2020年9 ...
动态模态分解 DMD | 做高维时间序列数据短时预测
通过这篇博客您将收获: 熟悉动态模态分解(DMD)的关键原理和重要的数学推导: 掌握利用 DMD 做多元时间序列预测任务的技术: 相关的有价值的资料分享,用于补充学习和拓展. CSDN 叶庭云:htt ...
lmd matlab 信号处理程序,LMD经验模态分解matlab程序.doc
LMD经验模态分解matlab程序 LMD经验模态分解matlab程序--原味的曾经也用滑动平均写过LMD,其实滑动平均的EMD才是原汁原味的居于均值分解. 分享给有需要的人,程序写的不好,只是希望 ...
【Python】这篇文章能让你明白经验模态分解（EMD）——EMD在python中的实现方法
暂时打断一下滤波专题,插播一条EMD在python中实现方法的文章. 本篇是Mr.看海:这篇文章能让你明白经验模态分解(EMD)--EMD在MATLAB中的实现方法的姊妹篇,也就是要在python中实 ...
动态模式分解（DMD）
本文地址:https://goodgoodstudy.blog.csdn.net/article/details/111873026 文章目录动态模式分解正式描述 DMD算法动态系统参考文献 ...
经验模态分解与Python调用实例
经验模态分解与Python调用实例经验模态分解 Python代码实现经验模态分解如果需要对一个信号进行降噪的话,我们通常会尝试将一个时域信号变换到不同的域空间,以期将信号中的噪声和有效信号在该域 ...
使用EMD【经验模态分解】对一维波形信号进行滤波去噪以及Python实现代码[emd eemd ceemdan]
使用EMD[经验模态分解]对一维波形信号进行滤波去噪以及Python实现代码 EMD[ Emprical Mode Decomposition]经验模态分解方法此处不再过多用赘述, 该信号处理方法可以 ...
matlab中使用VMD(变分模态分解)
最近我们被客户要求撰写关于VMD(变分模态分解)的研究报告,包括一些图形和统计输出. 拨号音信号的变模分解创建一个以4 kHz采样的信号,类似于拨打数字电话的所有键.将信号另存为MATLAB®时间数 ...
CEEMDAN：完全噪声辅助聚合经验模态分解（matlab）——学习笔记3
CEEMDAN:完全噪声辅助聚合经验模态分解--学习笔记3 从EMD到CEEMDAN 1.EMD EMD算法将基于原始信号的局部特征时间尺度,将原始信号分解为特征模态函数,即将其分解为从高频到低频的一 ...

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

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

最新文章

热门文章