Einführung

Ich konnte mir Fourier-Transformation oder Wavelet-Transformation nur als Methode zur Untersuchung von Zeitänderungen vorstellen, aber ich habe gelernt, dass es eine gute Methode namens Dynamic Mode Decomposition (DMD) gibt, mit der Moden sowohl in zeitlicher als auch in räumlicher Richtung extrahiert werden können. Ich werde es hier aufzeichnen, um mein Verständnis zu vertiefen.

Der Inhalt entspricht fast dem der Website, auf die ich verwiesen habe, und ich werde eine leicht modifizierte Version der Google-Übersetzung schreiben. Der Code für das Zeichnen von Grafiken befindet sich nicht in der Referenzquelle, daher habe ich ihn hinzugefügt. Sie sollten in der Lage sein, alles zu tun, von der Ausführung von DMD bis zum Zeichnen des Diagramms.

Referenz

Der Einfachheit halber werden wir die DMD des 3D-Vektorfeldes weglassen und nur eine einfache 1D-Skalarfunktion betrachten. Dynamic Mode Decomposition in Python

Ich wusste nicht, was SVD ist, also bezog ich mich darauf. Über die Beziehung zwischen PCA und SVD

Dynamische Moduszerlegung

Dynamic Mode Decomposition (DMD) ist eine relativ junge mathematische Innovation, die dynamische Systeme unter anderem in Bezug auf kohärente Strukturen lösen oder approximieren kann, die im Laufe der Zeit wachsen, zerfallen und / oder vibrieren. Die kohärente Struktur wird als DMD-Modus bezeichnet. Für jeden DMD-Modus ist eine entsprechende Zeitdynamik für einen einzelnen eindeutigen Wert definiert.

Mit anderen Worten, DMD transformiert ein dynamisches System in eine Überlagerung von Moden, deren Dynamik von Eigenwerten dominiert wird.

Überraschenderweise ist das mathematische Verfahren zur Unterscheidung von DMD-Moden von Eigenwerten rein linear, aber das System selbst kann nicht linear sein. Obwohl hier nicht behandelt, hat die Behauptung, dass nichtlineare Systeme als eine Menge von Moden- und Eigenwertpaaren beschrieben werden können, gute Gründe. Weitere Informationen finden Sie im Artikel zur Integration des Koopman-Operators und von DMD [^ 1] [^ 2] [^ 3].

DMD ist nicht nur ein nützliches Diagnosewerkzeug zur Analyse des internen Verhaltens eines Systems, sondern kann auch zur Vorhersage des zukünftigen Zustands des Systems verwendet werden. Alles was Sie brauchen ist der Modus und die eindeutigen Werte. Mit geringem Aufwand können Sie Modi und Eigenwerte kombinieren, um jederzeit eine Funktion zu generieren, die sich dem Systemzustand annähert.

Offizielle Definition

Betrachten Sie n Datensätze. ${(x_0,y_0),(x_1,y_1),\dots (x_n,y_n)\}$

Wobei $ x_i $ und $ y_i $ Spaltenvektoren der Größe $ m $ sind. Definieren Sie zwei $ m \ times n $ Matrizen. $X=[x_0\ x_1\ \dots\ x_n],\quad Y=[y_0\ y_1\ \dots\ y_n]$

Wenn Sie den Operator $ A $ wie folgt definieren $A=YX^\dagger$

Wobei $ X ^ \ dagger $ eine pseudoinverse Matrix von $ X $ ist und die dynamische Modenzerlegung von $ (X, Y) $ durch die Eigenwertzerlegung von $ A $ gegeben ist. Das heißt, der DMD-Modus und die Eigenwerte sind die Eigenvektoren und Eigenwerte von $ A $.

Die Definition von Tu et al. [^ 2] oben ist als exakte DMD bekannt. Dies ist derzeit die am häufigsten verwendete Definition und kann auf jeden Datensatz angewendet werden, der bestimmte Anforderungen erfüllt. In diesem Artikel interessiert mich am meisten, ob $ A $ die Formel $ y_i = Ax_i $ (wahrscheinlich ungefähr) für alle $ i $ erfüllt. Oder genauer: $Y=AX$

Offensichtlich ist $ X $ eine Menge von Eingabevektoren und $ Y $ ist eine Menge von entsprechenden Ausgabevektoren. Diese spezielle Interpretation von DMD ist sehr leistungsfähig, da sie eine bequeme Möglichkeit bietet, dynamische Stämme mit unbekannten maßgeblichen Gleichungen zu analysieren (und vorherzusagen). Wir werden später über dynamische Systeme sprechen.

Es gibt mehrere Sätze, die dieser Definition von DMD folgen [^ 2]. Einer der nützlicheren Sätze ist, dass $ Yv = 0 ist, wenn $ X $ und $ Y $ linear konsistent sind (mit anderen Worten, wenn $ Xv = 0 $ für den Vektor $ v $ ist). Nur wenn $ auch erfüllt ist, ist $ Y = AX $ vollständig erfüllt. Das Testen der linearen Konsistenz ist relativ einfach, wie wir später sehen werden. Kurz gesagt, lineare Konsistenz ist keine zwingende Voraussetzung für die Verwendung von DMD. Selbst wenn die DMD-Zerlegung von $ A $ die Formel $ Y = AX $ nicht vollständig erfüllt, ist sie die Methode mit dem kleinsten Quadrat und minimiert den Fehler in der $ L ^ 2 $ -Norm.

DMD-Algorithmus

Auf den ersten Blick scheint die Eigenwertzerlegung von $ A = YX ^ \ dagger $ kein großes Problem zu sein. Wenn die Größen $ X $ und $ Y $ stimmen, funktioniert es ein paar Mal, die Pinv- und Eig-Methoden von Numpy oder MATLAB aufzurufen. Das Problem entsteht, wenn $ A $ wirklich groß ist. Da $ A $ $ m \ mal m $ ist, wird die Eigenwertzerlegung umständlich, wenn $ m $ (die Anzahl der Signale in jeder Zeitabtastung) sehr groß ist.

Glücklicherweise können Sie mit Hilfe des genauen DMD-Algorithmus das Problem in kleinere Teile zerlegen.

Berechnen Sie die SVD (Singular Value Decomposition) von $ X $ und führen Sie bei Bedarf gleichzeitig eine Kürzung auf niedriger Ebene durch: $X=U\Sigma V^*$
Projizieren Sie die Matrix $ A $ auf $ U $, um $ \ tilde A $ zu berechnen: $\tilde A=U^* AU=U^*YV\Sigma^{-1}$
Berechnen Sie den Eigenwert $ \ lambda_i $ von $ \ tilde A $ und den Eigenvektor $ w_i : $\tilde AW=W\Lambda$$
Rekonstruieren Sie die Eigenwertzerlegung von $ A $ aus $ W $ und $ \ Lambda $. Der eindeutige Wert von $ A $ entspricht dem eindeutigen Wert von $ \ tilde A $. Der Eigenvektor von $ A $ ist in den Spalten $ \ Phi $ angegeben. $A\Phi=\Phi\Lambda,\quad \Phi=YV\Sigma^{-1}W$

Eine detailliertere Erklärung der Ableitung des Algorithmus finden Sie in Lit. [^ 1] [^ 2]. Es kann auch interessant sein festzustellen, dass $ \ Phi = UW $ eine alternative Ableitung von $ \ Phi $ ist, die als projizierter DMD-Modus bezeichnet wird. Dieser Artikel verwendet nur den exakten DMD-Modus.

Schauen wir uns Schritt für Schritt den Algorithmus in Python an. Beginnen Sie mit der Installation und dem Import aller benötigten Pakete.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import dot, multiply, diag, power
from numpy import pi, exp, sin, cos, cosh, tanh, real, imag
from numpy.linalg import inv, eig, pinv
from scipy.linalg import svd, svdvals
from scipy.integrate import odeint, ode, complex_ode
from warnings import warn

#Das ist mein Lieblings
plt.rcParams["font.family"] = "Times New Roman"      #Stellen Sie die gesamte Schriftart ein
plt.rcParams["xtick.direction"] = "in"               #Nach innen die x-Achsen-Skalierungslinie
plt.rcParams["ytick.direction"] = "in"               #Y-Achsen-Skalierungslinie nach innen
plt.rcParams["xtick.minor.visible"] = True           #Hinzufügen einer x-Achsen-Hilfsskala
plt.rcParams["ytick.minor.visible"] = True           #Hinzufügen einer Hilfsskala der y-Achse
plt.rcParams["xtick.major.width"] = 1.5              #Linienbreite der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.width"] = 1.5              #Linienbreite der Hauptskalenlinie der y-Achse
plt.rcParams["xtick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der x-Achse
plt.rcParams["ytick.minor.width"] = 1.0              #Linienbreite der Hilfsskalenlinie der y-Achse
plt.rcParams["xtick.major.size"] = 10                #Länge der x-Achsen-Hauptskalenlinie
plt.rcParams["ytick.major.size"] = 10                #Länge der Hauptmaßstablinie der y-Achse
plt.rcParams["xtick.minor.size"] = 5                 #Länge der Hilfslinie der x-Achse
plt.rcParams["ytick.minor.size"] = 5                 #Länge der Hilfsskala der y-Achse
plt.rcParams["font.size"] = 14                       #Schriftgröße
plt.rcParams["axes.linewidth"] = 1.5                 #Gehäusedicke

Lassen Sie uns einige Spieldaten generieren. Beachten Sie, dass Sie in der Praxis die maßgebliche Gleichung der Daten nicht unbedingt kennen. Hier erstellen wir einige Gleichungen, um einen Datensatz zu erstellen. Vergessen Sie nach der Generierung der Daten, dass sie vorhanden sind.

#Definieren Sie Zeit und Raum
x = np.linspace(-10, 10, 100)
t = np.linspace(0, 6*pi, 80)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)

#Erstellen Sie drei Zeit- und Raummuster
f1 = multiply(20-0.2*power(Xm, 2), exp((2.3j)*Tm))
f2 = multiply(Xm, exp(0.6j*Tm))
f3 = multiply(5*multiply(1/cosh(Xm/2), tanh(Xm/2)), 2*exp((0.1+2.8j)*Tm))

#Kombinieren Sie Signale, um eine Datenmatrix zu erstellen
D = (f1 + f2 + f3).T

#Erstellen Sie eine DMD-E / A-Matrix
X = D[:,:-1]
Y = D[:,1:]

fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111, projection="3d")
surf = ax.plot_surface(Xm, Tm, np.real(D.T), cmap="gray", linewidth=0)
ax.set_xlabel("x")
ax.set_ylabel("t")
fig.colorbar(surf)

Berechnen Sie nun die SVD für $ X $. Die erste Variable von primärem Interesse ist der Singularwert von $ X $, $ \ Sigma $. Wenn Sie eine $ X $ SVD erhalten, können Sie "Hochenergie" -Modi extrahieren und die Abmessungen Ihres Systems mit der richtigen orthogonalen Zerlegung (Proper Orthogonal Decomposition, POD) reduzieren. Ein Blick auf die Singularität bestimmt die Anzahl der abzuschneidenden Modi.

#Eingabematrix SVD
U2,Sig2,Vh2 = svd(X, False)

fig, ax = plt.subplots()
ax.scatter(range(len(Sig2)), Sig2, label="singular values")
ax.legend()

In Anbetracht der obigen Singularität können wir schließen, dass die Daten drei wichtige Modi haben. Schneiden Sie daher die SVD so ab, dass nur diese Modi enthalten sind. Dann baue $ \ tilde A $ und finde seine Eigenwertzerlegung.

#Rang 3 abgeschnitten
r = 3
U = U2[:,:r]
Sig = diag(Sig2)[:r,:r]
V = Vh2.conj().T[:,:r]

# A~Bauen
Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig))
mu,W = eig(Atil)

def circle(r=1):
    x,y = [],[]
    for _x in np.linspace(-180,180,360):
        x.append(np.sin(np.radians(_x)))
        y.append(np.cos(np.radians(_x)))
    return x,y

c_x,c_y = circle(r=1)

fig, ax = plt.subplots(figsize=(8,8))
ax.plot(c_x, c_y, c="k", linestyle="dashed")
ax.scatter(np.real(mu[0]), np.imag(mu[0]), label="1st")
ax.scatter(np.real(mu[1]), np.imag(mu[1]), label="2nd")
ax.scatter(np.real(mu[2]), np.imag(mu[2]), label="3rd")
ax.set_aspect("equal")
ax.set_xlabel(r"$\it{Re}\,\mu$")
ax.set_ylabel(r"$\it{Im}\,\mu$")
ax.legend()

Jeder eindeutige Wert von $ \ Lambda $ gibt Auskunft über das dynamische Verhalten des entsprechenden DMD-Modus. Die genaue Interpretation hängt von der Art der Beziehung zwischen $ X $ und $ Y $ ab. Bei Differenzgleichungen können viele Schlussfolgerungen gezogen werden. Wenn der Eigenwert einen Imaginärteil ungleich Null hat, tritt im entsprechenden DMD-Modus eine Vibration auf. Wenn der Eigenwert innerhalb des Einheitskreises liegt, wird der Modus gedämpft. Wenn der Eigenwert außen liegt, wächst der Modus. Wenn die Eigenwerte genau zum Einheitskreis passen, wächst oder fällt der Modus weder.

Erstellen Sie als Nächstes den genauen DMD-Modus.

#DMD-Modus erstellen
Phi = dot(dot(dot(Y, V), inv(Sig)), W)

fig, ax = plt.subplots()
ax.plot(x, np.real(Phi[:,0]), label="1st mode")
ax.plot(x, np.real(Phi[:,1]), label="2nd mode")
ax.plot(x, np.real(Phi[:,2]), label="3rd mode")
ax.set_xlabel("x")
ax.legend()

Die Spalte $ \ Phi $ ist der oben dargestellte DMD-Modus. Sie sind eine konsistente Struktur, die innerhalb des Systems entsprechend der unterschiedlichen Zeitdynamik wächst / abfällt / vibriert. Vergleichen Sie die Kurve im obigen Diagramm mit der rotierenden, sich entwickelnden Form im ursprünglichen 3D-Oberflächendiagramm. Sie werden Ähnlichkeiten feststellen.

Dies ist das technische Ende des DMD-Algorithmus. Mit einer Eigenwertzerlegung von $ A $ und einem grundlegenden Verständnis der Natur des Systems $ Y = AX $ können wir eine Matrix $ \ Psi $ erstellen, die der zeitlichen Entwicklung des Systems entspricht. Um den folgenden Code vollständig zu verstehen, lesen Sie die Differenzgleichungsfunktion $ x (t) $ im nächsten Abschnitt.

#Berechnung der Zeitentwicklung
b = dot(pinv(Phi), X[:,0])
Psi = np.zeros([r, len(t)], dtype='complex')
for i,_t in enumerate(t):
    Psi[:,i] = multiply(power(mu, _t/dt), b)

fig = plt.figure(figsize=(10,10))
ax1,ax2,ax3 = fig.add_subplot(311), fig.add_subplot(312), fig.add_subplot(313)
plt.subplots_adjust(wspace=0.5, hspace=0.5)
ax1.set_title("1st mode"), ax2.set_title("2nd mode"), ax3.set_title("3rd mode")

ax1.plot(t, np.real(Psi[0]), label=r"$\it{Re}\,\Psi$")
ax1.plot(t, np.imag(Psi[0]), label=r"$\it{Re}\,\Psi$")
ax2.plot(t, np.real(Psi[1])), ax2.plot(t, np.imag(Psi[1]))
ax3.plot(t, np.real(Psi[2])), ax3.plot(t, np.imag(Psi[2]))
ax3.set_xlabel("t")

fig.legend()

Die obigen drei Diagramme zeigen die Zeitdynamik der drei DMD-Modi. Beachten Sie, dass alle drei vibrieren. Außerdem scheint der zweite Modus exponentiell zu wachsen. Dies wird durch das Eigenwertdiagramm bestätigt.

Um eine Annäherung an die ursprüngliche Datenmatrix zu erstellen, multiplizieren Sie einfach $ \ Phi $ mit $ \ Psi $. In diesem speziellen Fall entspricht die Annäherung genau dem Original.

#DMD-Rekonstruktionsberechnung
D2 = dot(Phi, Psi)
np.allclose(D, D2) # True

Dynamisches System

In diesem Artikel werden nur zwei Interpretationen des Ausdrucks $ Y = AX $ betrachtet. Die erste Interpretation ist, wo $ A $ die Differenzgleichung definiert $x_{i+1}=Ax_i$

In diesem Fall erhöht der Operator $ A $ den dynamischen Systemzustand $ x_i $ um einen Zeitschritt. Angenommen, Sie haben eine Zeitreihe $ D . $D=[x_0\ x_1\ \dots\ x_{n+1}]$$

Wobei $ x_i $ ein Spaltenvektor ist, der den Dimensionsstatus $ m $ des Systems zum Zeitschritt $ i $ definiert. Dann können Sie $ X $ und $ Y $ wie folgt definieren: $X=[x_0\ x_1\ \dots\ x_{n}],\quad Y=[x_1\ x_2\ \dots\ x_{n+1}]$

Auf diese Weise entspricht jedes Paar von $ X $ - und $ Y $ -Spaltenvektoren einer einzelnen Iteration der Differenzgleichung, typischerweise: $Y=AX$

Verwenden Sie DMD, um die intrinsische Zerlegung von $ A \ Phi = \ Phi \ Lambda $ zu ermitteln. Im DMD-Modus und mit Eigenwerten können Sie $ Y = AX $ einfach in eine Funktion konvertieren, die durch die diskrete Zeititeration $ k $ des Zeitschritts $ \ Delta t $ definiert ist. $x_k=\Phi\Lambda^k\Phi^\dagger x_0$

Die entsprechende Funktion für die kontinuierliche Zeit $ t $ ist $x(t)=\Phi\Lambda^{t/\Delta t}\Phi^\dagger x(0)$

Was wirklich erstaunlich ist, ist, dass ich eine explizite Funktion in der Zeit definiert habe, indem ich nur die Daten verwendet habe. Dies ist ein gutes Beispiel für eine gleichungslose Modellierung.

In der zweiten Interpretation von $ Y = AX $, die in diesem Artikel betrachtet wird, definiert $ A $ das Differentialgleichungssystem. $\dot x=Ax$

In diesem Fall berechnet der Operator $ A $ die lineare Ableitung des Vektors $ x_i $ in Bezug auf die Zeit. Die Matrizen $ X $ und $ Y $ bestehen aus $ n $ Abtastwerten von Vektorfeldern. Die $ i $ -te Spalte von $ X $ ist der Positionsvektor $ x_i $. Die $ i $ -te Spalte von $ Y $ ist der Geschwindigkeitsvektor $ \ dot x_i $.

Nach der Berechnung der DMD ist die Funktion der Zeit der vorherigen sehr ähnlich (dh der Differenzgleichung). Wenn $ x (0) $ eine Anfangsbedingung ist und $ t $ eine kontinuierliche Zeit ist $x(t)=\Phi\text{exp}(\Lambda t)\Phi^\dagger x(0)$

Hilfsmethode

Kombinieren Sie den DMD-Code der Einfachheit halber zu einer Methode, definieren Sie mehrere Hilfsmethoden, um die lineare Konsistenz zu überprüfen, und bestätigen Sie die Lösung.

def nullspace(A, atol=1e-13, rtol=0):
    # from http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
    A = np.atleast_2d(A)
    u, s, vh = svd(A)
    tol = max(atol, rtol * s[0])
    nnz = (s >= tol).sum()
    ns = vh[nnz:].conj().T
    return ns

def check_linear_consistency(X, Y, show_warning=True):
    # tests linear consistency of two matrices (i.e., whenever Xc=0, then Yc=0)
    A = dot(Y, nullspace(X))
    total = A.shape[1]
    z = np.zeros([total, 1])
    fails = 0
    for i in range(total):
        if not np.allclose(z, A[:,i]):
            fails += 1
    if fails > 0 and show_warning:
        warn('linear consistency check failed {} out of {}'.format(fails, total))
    return fails, total

def dmd(X, Y, truncate=None):
    U2,Sig2,Vh2 = svd(X, False) # SVD of input matrix
    r = len(Sig2) if truncate is None else truncate # rank truncation
    U = U2[:,:r]
    Sig = diag(Sig2)[:r,:r]
    V = Vh2.conj().T[:,:r]
    Atil = dot(dot(dot(U.conj().T, Y), V), inv(Sig)) # build A tilde
    mu,W = eig(Atil)
    Phi = dot(dot(dot(Y, V), inv(Sig)), W) # build DMD modes
    return mu, Phi

def check_dmd_result(X, Y, mu, Phi, show_warning=True):
    b = np.allclose(Y, dot(dot(dot(Phi, diag(mu)), pinv(Phi)), X))
    if not b and show_warning:
        warn('dmd result does not satisfy Y=AX')

Einschränkungen

DMD weist einige bekannte Einschränkungen auf. Erstens kann die Unveränderlichkeit von Translation und Rotation nicht besonders gut gehandhabt werden. Zweitens kann eine vorübergehende Zeitoperation insgesamt fehlschlagen. Das folgende Beispiel veranschaulicht diese Probleme.

1. Translationale Invarianz

Der folgende Datensatz ist sehr einfach. Es besteht aus einem einzelnen Modus (Gauß), der sich im Verlauf des Systems entlang der räumlichen Domäne verschiebt. Sie könnten denken, dass DMD dies sauber handhabt, aber das Gegenteil passiert. SVD erhält viele Werte anstelle einer einzelnen, genau definierten Singularität.

#Definieren Sie Zeit und Raum
x = np.linspace(-10, 10, 50)
t = np.linspace(0, 10, 100)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)

#Erstellen Sie Daten in einem einzigen Modus, der sich sowohl räumlich als auch zeitlich bewegt
D = exp(-power((Xm-Tm+5)/2, 2))
D = D.T

#Eingabe- / Ausgabematrix extrahieren
X = D[:,:-1]
Y = D[:,1:]
check_linear_consistency(X, Y)

U2,Sig2,Vh2 = svd(X, False)

fig = plt.figure(figsize=(10,5))
ax1,ax2 = fig.add_subplot(121), fig.add_subplot(122)
ax1.set_xlabel("x"), ax1.set_ylabel("t")
ax1.imshow(D.T, aspect=0.5)
ax2.scatter(range(len(Sig2)), Sig2)

Das Diagramm links zeigt die zeitliche Variation des Systems. Das Diagramm rechts zeigt die Singularwerte. Wir fanden heraus, dass fast 10 DMD-Modi erforderlich waren, um das System genau zu approximieren. Betrachten Sie die folgende Darstellung. Vergleichen Sie die wahre Dynamik mit einer unterschiedlichen Anzahl überlappender Modi.

def build_dmd_modes(t, X, Y, r):
    """
DMD-Rekonstruktion
    """
    mu, Phi = dmd(X, Y, truncate=r)
    b = dot(pinv(Phi), X[:,0])
    Psi = np.zeros([r, len(t)], dtype='complex')
    for i,_t in enumerate(t):
        Psi[:,i] = multiply(power(mu, _t/dt), b)
    return dot(Phi, Psi)

fig = plt.figure(figsize=(10,10))
axes = []
for i in range(9):
    axes.append(fig.add_subplot(331+i))
plt.subplots_adjust(wspace=0.5, hspace=0.5)

for i,ax in enumerate(axes):
    ax.set_title("{} modes".format(i+1))
    ax.imshow(np.real(build_dmd_modes(t, X, Y, r=i+1).T), aspect=0.5)

2. Vorübergehendes Zeitverhalten

In diesem letzten Beispiel wird ein Datensatz untersucht, der temporäre Zeitdynamik enthält. Insbesondere wird angezeigt, wann Gauß vorhanden ist und wann es in den Daten nicht vorhanden ist. Leider kann DMD diese Daten nicht genau zerlegen.

#Definieren Sie Zeit und Raum
x = np.linspace(-10, 10, 50)
t = np.linspace(0, 10, 100)
dt = t[2] - t[1]
Xm,Tm = np.meshgrid(x, t)

#Erstellen Sie Daten in einem einzigen temporären Modus
D = exp(-power((Xm)/4, 2)) * exp(-power((Tm-5)/2, 2))
D = D.astype('complex').T

#Eingabe- / Ausgabematrix extrahieren
X = D[:,:-1]
Y = D[:,1:]
check_linear_consistency(X, Y)

#DMD-Zersetzung
r = 1
mu,Phi = dmd(X, Y, r)
check_dmd_result(X, Y, mu, Phi)

#Zeitentwicklung
b = dot(pinv(Phi), X[:,0])
Psi = np.zeros([r, len(t)], dtype='complex')
for i,_t in enumerate(t):
    Psi[:,i] = multiply(power(mu, _t/dt), b)

fig = plt.figure(figsize=(10,10))
ax1,ax2 = fig.add_subplot(221),fig.add_subplot(222)
ax3,ax4 = fig.add_subplot(223),fig.add_subplot(224)

ax1.imshow(np.real(D), aspect=2.0), ax1.set_title("True")
ax2.imshow(np.real(dot(Phi, Psi).T), aspect=0.5), ax2.set_title("1-mode approx.")
ax3.plot(x, np.real(Phi)), ax3.set_xlabel("x"), ax3.set_title("mode1")
ax4.plot(t, np.real(Psi[0])), ax3.set_xlabel("t"), ax4.set_title("time evol. of mode 1")

Das DMD identifiziert den Modus korrekt, kann jedoch das Verhalten der Zeit nicht vollständig identifizieren. Dies kann verstanden werden, wenn man bedenkt, dass das Zeitverhalten der DMD-Zeitreihe von den Eigenwerten abhängt. Der Eigenwert kann nur eine Kombination aus exponentiellem Wachstum (Realteil des Eigenwerts) und Schwingung (Imaginärteil) charakterisieren.

Das Interessante an diesem System ist, dass die ideale Zerlegung aus einer Einzelmodus-Überlagerung (wie gezeigt) mit unterschiedlichen Eigenwerten bestehen kann. Stellen Sie sich einen einzelnen Modus vor, multipliziert mit vielen orthogonalen Sinus- und Cosinus-Linearkombinationen (Fourier-Reihen), die sich der tatsächlichen Zeitdynamik annähern. Leider kann eine einzelne SVD-basierte DMD-Anwendung nicht denselben DMD-Modus mehrmals mit unterschiedlichen eindeutigen Werten generieren.

Darüber hinaus ist zu beachten, dass selbst wenn das Zeitverhalten als große Anzahl von Eigenwerten korrekt extrahiert werden kann, die Vorhersagefunktion der Lösung ohne ein vollständiges Verständnis des Übergangsverhaltens selbst unzuverlässig ist. Temporäres Verhalten ist von Natur aus nicht dauerhaft.

DMD mit mehreren Auflösungen (mrDMD) versucht, das Problem des temporären Zeitbetriebs durch rekursives Anwenden von DMD zu verringern.

Fazit

Trotz seiner Einschränkungen ist DMD ein sehr leistungsfähiges Werkzeug zur Analyse und Vorhersage dynamischer Systeme. Alle Datenwissenschaftler mit unterschiedlichem Hintergrund sollten ein gutes Verständnis für DMD und dessen Anwendung haben. Der Zweck dieses Artikels ist es, die Theorie hinter DMD und praktische Python-Codebeispiele bereitzustellen, die mit realen Daten verwendet werden können. Wir haben die formale Definition von DMD überprüft, Schritt für Schritt durch den Algorithmus gegangen und einige einfache Anwendungsfälle ausprobiert, einschließlich fehlgeschlagener. Wir hoffen, dass Sie dadurch ein klareres Verständnis dafür erhalten, wie DMD für Forschungs- oder Ingenieurprojekte gilt.

DMD hat viele Erweiterungen. Zukünftige Arbeiten können Beiträge zu einigen dieser Erweiterungen enthalten, z. B. DMD mit mehreren Auflösungen (mrDMD) und DMD mit geringer Auflösung (sDMD).

Zusammenfassung

Wenn Sie die zeitliche Entwicklung eines zweidimensionalen Arrays, z. B. eines Hochgeschwindigkeitskamerabilds, DMD möchten, reduzieren Sie das zweidimensionale Array in ein eindimensionales Array, und der obige Code funktioniert. Ich werde ein Beispiel für DMD der Simulation des Kalman-Wirbels mit CFD hinzufügen.

[PYTHON] Einführung in die dynamische Moduszerlegung