利用matlab实现DMD动态模态分解(在一维信号或二维流场矢量中的应用）

（2020年9月更新1，增加了第二章对于DMD模态排序的方法介绍。另外增加了一个注释：%旧版本matlab可以用右边形式替换Time_DMD=(b*ones(1,N)).*Q;。 由于旧版本中的 $.\ast$不支持向量和矩阵之间的操作，所以特此注释。如果代码中还有类似的报错，请参照前面提到的方法进行自行修改即可。）

0 前言

DMD的全称为Dynamic Mode Decomposition，译为动态模态分解，是一种常见的数据降维方法，本文主要是尝试利用matlab对DMD方法进行实现。
本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

本文只代表个人的理解，如果有描述不准确或者错误的地方，欢迎在评论区指正。由于DMD可能涉及到了包括且不限于流体、机械、信号、控制、深度学习等各领域，所以请各位大佬们如果发现错误，或者有什么补充，欢迎在评论区指出，感谢。

期间涉及到一些POD分解与SVD分解的相关知识，建议先阅读下面这篇文章（自引）：
【利用matlab实现POD分解(在一维信号或二维流场矢量中的应用）】

文章利用到的参考文献如下：
[1]Modal Analysis of Fluid Flows: An Overview（Kunihiko Taira、Steven L. Brunton等人）
[2]Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems（J. Nathan Kutz，Steven L. Brunton等人）

0.1 特征根的计算与含义

matlab中特征值和特征向量都可以用eig()函数进行计算。

[V,D] = eig(A)

其中V为特征向量组成的矩阵，D为特征根所组成的矩阵。A、V和D三个矩阵之间的关系为：A* V = V* D。
其中V的每一列都为一个特征向量，特征值对应的D矩阵中对角线上的相应列。

从线性变换的角度来看，矩阵的特征根和特征向量代表了变换的大小和方向。

当特征值为实数时，大于1代表变换的增长，小于1代表变换的缩小，特征向量代表变换的方向。以矩阵A=[0.95,-0.35;-0.35,0.95]为例，其特征值为0.6和1.3，特征向量的方向一个为45度，一个-45°。矩阵A对于随机点阵的效果如下所示。可以看到在经过矩阵A相乘后，随机点阵发生了变形，变形的方向就是特征向量的方向，变形的大小就是特征根的大小。大于1的方向拉长，小于1的方向缩短。

当特征值为复数时，则代表旋转。以A=[0.866,0.5;-0.5,0.866]为例，特征值为 cos30° ± sin30° * i ，每次乘以矩阵A代表旋转30度（当然我这里为了作图好看，通过调整特征向量把图形变为正圆，一般情况是一个椭圆）。

1 DMD的基本思路

DMD分解的基本思想就是线性变换。比如定义一个信号Utx，x为一维信号上每一点的数值，t为随时间变化下，整个信号的变化。时间采样点为N个，空间采样点为m个。Utx的转置矩阵为Uxt。
$$U_{tx}=\left[\begin{array}{cccc}u\left(t_1,x_1\right)& u\left(t_1,x_2\right)& \cdots & u\left(t_1,x_m\right)\\ u\left(t_2,x_1\right)& u\left(t_2,x_2\right)& \cdots & u\left(t_2,x_m\right)\\ ⋮& ⋮& \ddots & ⋮\\ u\left(t_N,x_1\right)& u\left(t_N,x_2\right)& \cdots & u\left(t_N,x_m\right)\end{array}\right]$$
初始时刻t=1时，信号以列向量的形式可以记为Ux1：
$$U_{x1}=\left[\begin{array}{c}u\left(t_1,x_1\right)\\ u\left(t_1,x_2\right)\\ ⋮\\ u\left(t_1,x_m\right)\end{array}\right]$$
那么DMD方法认为，如果系统是一个线性系统，则我们可以找到一个矩阵A，使得：
$$\begin{gathered} U_{x2}=A\ast U_{x1} \\ {U}_{x3}=A\ast U_{x2}=A^2\ast U_{x1} \\ {U}_{x4}=A\ast U_{x3}=A^3\ast U_{x1} \end{gathered}$$
因此，我们只要知道初始的状态Ux1，以及系统的变换矩阵A，就可以知道系统之后任意时刻t的状态UxN。对于某些非线性系统，也可以找到一个近似的矩阵A，用来进行数据降维处理，然而效果肯定比线性系统的效果差。

因此，DMD分解的目的，就是把系统进行线性变换的分解。第k个时间步信号Uxk可以经过DMD分解变为
$$U_{x,k}=\sum _{i=1}^m{b}_i\cdot {\lambda }_i^k\cdot {\varphi }_i$$
其中b为初始值，$\lambda$为矩阵A的特征根，$\varphi$为对应的模态。
其中特征根可以利用衰减率$\mu$来替换：
$$\lambda =\exp (\mu \cdot \Delta t)$$
我们也可以把离散步k转换为时间t的形式：
$$U_{x,t}=\sum _{i=1}^m{b}_i\cdot e^{{\mu }_i\cdot t}\cdot {\varphi }_i$$

所以综上所述，我们只要知道矩阵A和初始时刻Ux1，就可以知道系统之后任意时刻的状态。同时根据矩阵A的特征值，还可以得到系统的增长率和变化频率。

我们设矩阵X为第1步至第N-1步组成的矩阵：
X = U x , t ( 1 ⋯ N − 1 ) = [ u ( x 1 , t 1 ) u ( x 1 , t 2 ) ⋯ u ( x 1 , t N − 1 ) u ( x 2 , t 1 ) u ( x 2 , t 2 ) ⋯ u ( x 2 , t N − 1 )